R help - Apr 2010 - Problem with Loops

Hello,

I am trying to write a script with the end goal of graphing power (y) as a
result of sample size (x) at a variety of effects sizes.  I am new to loops,
and I think my problem is there.  Here's the script, which is modified from
the script found at the bootom of
http://www.statmethods.net/stats/power.html.  ANy help would be much
appreciated!mtesche at ualberta.ca

library(pwr)
library(lattice)

# range of sample sizes
r <- seq(5,300,5)
nr <- length(r)

# effect sizes
p <- seq(.4,1.2,.2)
np <- length(p)

# obtain power
power <- array(numeric(nr*np), dim=c(nr,np))
for (i in 1:np){
  for (j in 1:nr){
    result <- pwr.t.test(n = r[j], d = p[i],
    sig.level = .05, power = NULL,
    alternative = "two.sided")
    power[j,i] <- ceiling(result$power)
  }
}

# set up graph
xrange <- range(r)
yrange <- round(range(power))
colors <- rainbow(length(p))
plot(xrange, yrange, type="n",
  xlab="Sample Size (n)",
  ylab="Power" )

# add power curves
for (i in 1:np){
  lines(r, power[,i], type="l", lwd=2, col=colors[i])
}

# add annotation (grid lines, title, legend)
abline(v=0, h=seq(0,yrange[2],50), lty=2, col="grey89")
abline(h=0, v=seq(xrange[1],xrange[2],.02), lty=2,
   col="grey89")
title("Power Estimation for Student's T-Test\Sample size
  Sig=0.05 (Two-tailed)")
legend("topright", title="Power", as.character(p),
   fill=colors)

-- 
View this message in context:
http://n4.nabble.com/Problem-with-Loops-tp1794288p1794288.html
Sent from the R help mailing list archive at Nabble.com.

Hello,

mtesche wrote:
> Hello,
> 
> I am trying to write a script with the end goal of graphing power (y) as a
> result of sample size (x) at a variety of effects sizes.  I am new to
loops,
> and I think my problem is there.  Here's the script, which is modified
from
> the script found at the bootom of
> http://www.statmethods.net/stats/power.html.  ANy help would be much
> appreciated!mtesche at ualberta.ca
> 
> library(pwr)
> library(lattice)
> 
> # range of sample sizes
> r <- seq(5,300,5)
> nr <- length(r)
> 
> # effect sizes
> p <- seq(.4,1.2,.2)
> np <- length(p)
> 
> # obtain power
> power <- array(numeric(nr*np), dim=c(nr,np))
> for (i in 1:np){
>   for (j in 1:nr){
>     result <- pwr.t.test(n = r[j], d = p[i],
>     sig.level = .05, power = NULL,
>     alternative = "two.sided")
>     power[j,i] <- ceiling(result$power)
>   }
> }
First problem.  Why are you calling ceiling, which will return 1 for any value 
of a power?  What were you thinking that line was doing?