[R] (no subject)

[Kalia Schuster]

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]]