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? 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 would appreciate the help, Marco MPM Boks, MD PhD, Department of Psychiatry, B01.206 University Medical Centre Utrecht, PO box 85500, 3508 GA Utrecht The Netherlands. Tel: +31 30 2506370 Fax: +31 30 2505509 Email: m.p.m.boks@umcutrecht.nl <mailto:m.p.m.boks@umcutrecht.nl> [[alternative HTML version deleted]]
You need to debug your function. If you put some 'cat(p, ind)' statements: for (ind in 1:500) {cat(p, ind, '\n') res[ind,]<-powerb(x[c(1:(500-ind),(500+ind):1000)],4,0.05,0.20,0.1,m)} Here is the place where it seemed to 'stall' on me: 0.3 250 0.3 251 0.3 252 0.3 253 0.3 254 0.3 255 0.3 256>It looks like at this point there was probably a loop in your function. You will have to learn how to use debug to trace through what it is doing at this point. It had produced the graphs for the other values of 'p'. So there must be a 'bug' there. 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? > > > 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 would appreciate the help, > > Marco > > > MPM Boks, MD PhD, > > Department of Psychiatry, B01.206 > > University Medical Centre Utrecht, > > PO box 85500, 3508 GA Utrecht > > The Netherlands. > > Tel: +31 30 2506370 > > Fax: +31 30 2505509 > > Email: m.p.m.boks@umcutrecht.nl <mailto:m.p.m.boks@umcutrecht.nl> > > > > [[alternative HTML version deleted]] > > ______________________________________________ > 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 >-- Jim Holtman Cincinnati, OH +1 513 646 9390 (Cell) +1 513 247 0281 (Home) What is the problem you are trying to solve? [[alternative HTML version deleted]]
On 7/3/06, Boks, M.P.M. <M.P.M.Boks at 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)) }