Greetings List -
About a week ago I posted a message requesting assistance in cleaning up
some R-code for automating power curves at varying sampling effort, and
just wanted to share the informative responses I received. Thanks to all
of the responders and the list for contributing to information sharing.
There were two equally successful approaches, I'll copy the first, which
involves loops, for its simplicity in understanding.
rm(list=ls())
x<-c(4,5,5,6,8,7,9,9,7,7,8,6,5,4,6,8) # x is my vector of
samples, in this case the total number of coral species in 16 quadrat
tosses
delta<-mean(x)*0.2 # delta is my desired level
of change to detect in the power analysis
quad<-2:16 # quad is simply the number
of quadrats sampled, starting with 2 and ending with the total
maxquad<-16
MEAN<-STDEV<-p<-rep(0,maxquad-1) # MEAN is your "mean"
vector,STDEV is the "stdev" vector, and p is your "power" vector.
for (i in 2:maxquad){
MEAN[i-1]<-mean(sample(x=x,size=i,replace=FALSE))
STDEV[i-1]<-sd(sample(x=x,size=i,replace=FALSE))
p[i-1]<-power.t.test(i,sd=STDEV[i-1],delta)$power
}
## the plot:
par(mfrow=c(3,1))
plot(quad,MEAN)
plot(quad,STDEV)
abline(h=(mean(x)*0.3), lty=3)
plot(quad,p,ylab="POWER")
The second approach, which I was rightfully informed is computationally
much faster in R, identifies a matrix, then fills it. However, for this
exercise, computation time is not really a factor.
x<-c(35, 27, 32, 38, 36); delta<-mean(x)*0.3; quad<-1:5; maxquad<-5
rSamples=apply(X=matrix(nrow=maxquad,ncol=maxquad,data=x),MARGIN=2,FUN="sample",x,size=maxquad,replace=F)
rSamples[lower.tri(rSamples,diag=F)]=NA
rMeans=colMeans(rSamples,na.rm=T)
rSD=sd(rSamples,na.rm=T)
Cheers,
Peter
--
Peter Houk, PhD
Chief Biologist
Pacific Marine Resources Institute
www.pacmares.com
www.micronesianfishing.com
