From: Simon Knapp <sleepingwell_at_gmail.com>

Date: Mon 03 Jul 2006 - 23:10:55 EST

}

return(c(sx,n1))

}

repeat{

}

return(c(sx,n1))

}

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Mon Jul 03 23:19:08 2006

Date: Mon 03 Jul 2006 - 23:10:55 EST

On 7/3/06, Boks, M.P.M. <M.P.M.Boks@umcutrecht.nl> wrote:

*>
**>
**> Hi,
**>
*

> I am new in R and stumbled on a problem my (more experienced) friends

*> can not help with with. Why isnt this code working?
**>
**> The function is working, also with the loop and the graph appears,
**>
**> only when I build another loop around it (for different values of p) ,
**> R stays in a loop?
**>
**> Can't it take more then 2 loops in one program?
*

You can have as many as you like (as far as I know)

*>
**>
*

> powerb<-function(x,sp2,a,b,b1,m)

*> { sx<-(sum(x^2)-(sum(x)^2)/length(x))/length(x)
**> n0<-ceiling((((qnorm(1-(a/2))+qnorm(1-b))/b1)^2)*(((m+1)/m)*sp2/sx))
**> repeat
**> {
**> n1<-ceiling((((qt(1-(a/2),n0-4)+qt(1-b,n0-4))/b1)^2)*(((m+1)/m)*sp2/sx))
**> if(n0==n1) break
**> n0<-n1
**> }
**> return(c(sx,n1))
**> }
**>
**> x<-rnorm(1000,0,1)
**> x<-x[order(x)]
**>
**> res<-matrix(0,1000,2)
**>
**>
**> #use the function and plot for different values of ind and p
**> for ( p in c(0.05,0.10,0.15,0.20,0.25,0.30,0.40,0.50))
**> { risk<-p*(2-p)
**> nonrisk<-(1-p)^2
**> m<-nonrisk/risk
**>
**> for (ind in 1:500)
**> {res[ind,]<-powerb(x[c(1:(500-ind),(500+ind):1000)],4,0.05,0.20,0.1,m)}
**>
**> plot(res[,1],res[,2],type="p",ylab="n",xlab="var(x)",main="b=0.1,power=0
**> .80,alpha=0.05,dominant met p=0.25")}
**>
*

I modified your function as follows:

powerb<-function(x,sp2,a,b,b1,m){ sx<-(sum(x^2)-(sum(x)^2)/length(x))/length(x) n0<-ceiling((((qnorm(1-(a/2))+qnorm(1-b))/b1)^2)*(((m+1)/m)*sp2/sx)) repeat{ n1<-ceiling((((qt(1-(a/2),n0-4)+qt(1-b,n0-4))/b1)^2)*(((m+1)/m)*sp2/sx)) if(n0==n1) break else cat("p = ", p, ", ind = ", ind, ", n0 = ", n0, ", n1 = ", n1, "\n", sep='') n0<-n1

}

return(c(sx,n1))

}

and got, for example, the following output:

....

p = 0.2, ind = 278, n0 = 2288, n1 = 2289 p = 0.2, ind = 278, n0 = 2289, n1 = 2288 p = 0.2, ind = 278, n0 = 2288, n1 = 2289 p = 0.2, ind = 278, n0 = 2289, n1 = 2288 p = 0.2, ind = 278, n0 = 2288, n1 = 2289 p = 0.2, ind = 278, n0 = 2289, n1 = 2288 p = 0.2, ind = 278, n0 = 2288, n1 = 2289 p = 0.2, ind = 278, n0 = 2289, n1 = 2288...

so your in an infinite loop. Maybe you can check the difference between n0 and n1 instead, or perhaps better, do something like:

powerb<-function(x,sp2,a,b,b1,m){ sx<-(sum(x^2)-(sum(x)^2)/length(x))/length(x) n0<-ceiling((((qnorm(1-(a/2))+qnorm(1-b))/b1)^2)*(((m+1)/m)*sp2/sx))n0.1 <- n1.1 <- Inf

repeat{

n1<-ceiling((((qt(1-(a/2),n0-4)+qt(1-b,n0-4))/b1)^2)*(((m+1)/m)*sp2/sx)) if(n0==n1) break if(n0.1==n1 && n1.1==n0) break n0.1 <- n0 n1.1 <- n1 n0<-n1

}

return(c(sx,n1))

}

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Mon Jul 03 23:19:08 2006

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.1.8, at Tue 04 Jul 2006 - 02:14:48 EST.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*