Greetings,
I'm slowly becoming familiar with R and have written some (long winded)
code to automate the simple generation of power curves to help monitoring
programs define a desired level of sampling that matches their populations
in question. My question for the list is whether someone might be able to
assist in a quick review of the code, and to identify how/where to place a
second loop to reduce. Thanks much for any insight, my goal is to make a
refined version of this simple code freely available for monitoring
programs dealing with sample design issues.
A simple example of this code:
##First part is to define 5 samples of % cover that were estimated from
quadrat surveys
x<-c(35, 27, 32, 38, 36); delta<-mean(x)*0.3; quad<-2:5; maxquad<-5
##Second part is to generate the mean values associated with varying
sampling effort (n=2 to 5 quadrats). Here, I wrote (long winded) code to
take every possible combination of each sampling effort, and provide the
overall mean.
rdraw2<-mean(sample(x=x, size=2, replace=FALSE)); for (t in 1:maxquad)
rdraw2[t+1] <-{rdraw2[t]<-mean(sample(x=x, size=2, replace=FALSE))};
rdraw3<-mean(sample(x=x, size=3, replace=FALSE)); for (t in 1:maxquad)
rdraw3[t+1] <-{rdraw3[t]<-mean(sample(x=x, size=3, replace=FALSE))};
rdraw4<-mean(sample(x=x, size=4, replace=FALSE)); for (t in 1:maxquad)
rdraw4[t+1] <-{rdraw4[t]<-mean(sample(x=x, size=4, replace=FALSE))};
rdraw5<-mean(sample(x=x, size=5, replace=FALSE)); for (t in 1:maxquad)
rdraw5[t+1] <-{rdraw5[t]<-mean(sample(x=x, size=5, replace=FALSE))}
##Third part is to generate the standard deviation values associated with
varying sampling effort (n=2 to 5 quadrats). Here, I wrote (long winded)
code to take every possible combination of each sampling effort, and
provide the overall standard deviation.
rdrawsd2<-sd(sample(x=x, size=2, replace=FALSE)); for (t in 1:maxquad)
rdrawsd2[t+1] <-{rdrawsd2[t]<-sd(sample(x=x, size=2, replace=FALSE))};
rdrawsd3<-sd(sample(x=x, size=3, replace=FALSE)); for (t in 1:maxquad)
rdrawsd3[t+1] <-{rdrawsd3[t]<-sd(sample(x=x, size=3, replace=FALSE))};
rdrawsd4<-sd(sample(x=x, size=4, replace=FALSE)); for (t in 1:maxquad)
rdrawsd4[t+1] <-{rdrawsd4[t]<-sd(sample(x=x, size=4, replace=FALSE))};
rdrawsd5<-sd(sample(x=x, size=5, replace=FALSE)); for (t in 1:maxquad)
rdrawsd5[t+1] <-{rdrawsd5[t]<-sd(sample(x=x, size=5, replace=FALSE))}
##Fourth part is to produce a single vector for standard deviations and
means that were produced above for input into power analysis
stdev<-c(mean(rdrawsd2), mean(rdrawsd3), mean(rdrawsd4), mean(rdrawsd5));
mean<- c(mean(rdraw2), mean(rdraw3), mean(rdraw4), mean(rdraw5));
##Fifth part is to run the power analysis for each level of sampling, then
make a vector for the resultant power at each level of effort.
p2<-power.t.test(n=2, sd=mean(rdrawsd2), delta=delta);
p3<-power.t.test(n=3, sd=mean(rdrawsd3), delta=delta);
p4<-power.t.test(n=4, sd=mean(rdrawsd4), delta=delta);
p5<-power.t.test(n=5, sd=mean(rdrawsd5), delta=delta);
power<-c(p2$power, p3$power, p4$power, p5$power);
##Finally, the last part is to produce three simple graphs that help to
understand whether sampling is sufficient, to isolate upon steady means,
standard deviation, and resultant power. Clearly power will alway continue
to increase with sample size, however, I'd like to estimate when standard
deviations saturate out...
plot(quad, power); par(mfrow=c(3,1)); plot(quad,mean); {plot(quad,stdev);
abline(h=(mean(x)*0.3), lty=3)}; plot(quad,power)
Assistance and discussion is greatly appreciated,
Peter
--
Peter Houk, PhD
Chief Biologist
Pacific Marine Resources Institute
www.pacmares.com
www.micronesianfishing.com
