Hello,
I'm new to R. I'm trying to run a power analysis in a for loop to find an
ideal sample size. The situation is I am doing counts of fish at 24
different sites and rating those sites as 1 or 2 for good or poor habitat.
So the counts are continuous but the habitat rating isn't. When I try to
run a negative binomial general linearized model with given values for mu:
exp1.302 and exp-0.32 and a given theta: exp(0.75), I get an error message
before I can run the whole loop and plot the power of the sample size
against sample. I used a template my teacher showed in class, but since he
used a poisson dist., I'm stuck on the glm.nb.
Please help! This is an assignment, so any suggestions would be wonderful
and I don't expect anything more. But please remember, I am a beginner and
only have a crude knowledge of R at this point.
Much appreciated.
#Power analysis
rm(list=ls())
graphics.off()
set.seed(1)
alpha=0.05
m=200
Ns = seq(2,40,2)
nNs = length(Ns)
NAs=rep(NA,nNs)
results= data.frame(n=Ns, power=NAs)
for(j in 1:nNs) {
n=Ns[j]
reject_count = 0
for(i in 1:m ) {
mu_poor=exp(-0.32)
mu_suit=exp(1.302)
mu=c(rep(mu_poor,n/2),rep(mu_suit,n/2))
theta=exp(0.175)
count=rnbinom(n=n,size=theta,mu=mu)
hab_poor=1
hab_suit=2
habitat=sample(hab_poor:hab_suit, size=n, replace=T)
model = glm.nb(count~habitat,link="log")
summary(model)
pval = summary(model)$coeff[2,4]
pval
reject = pval < alpha
if (reject) reject_count = reject_count+1
}
reject_rate = reject_count/m
power = reject_rate
typeIIerror=1-reject_rate
results[j,]$power=power
}
plot (x=Ns, y=results$power, xlab="Sample size", ylab="Power")
lines(Ns,rep(0.95,nNs))
par(lty=2)
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
