I'm trying to convert an S-Plus program to R. Since I'm a SAS
programmer I'm not facile is either S-Plus or R, so I need some help. All I
did was convert the underscores in S-Plus to the assignment operator <-.
Here are the first few lines of the S-Plus file:
sshc _ function(rc, nc, d, method, alpha=0.05, power=0.8,
tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5)
{
### for method 1
if (method==1) {
ne1 _ ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01)
return(ne=ne1)
}
My translation looks like this:
sshc<-function(rc, nc=500, d=.5, method=3, alpha=0.05, power=0.8,
tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5)
{
### for method 1
if (method==1) {
ne1<-ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01)
return(ne=ne1)
}
The program runs without throwing errors, but I'm not getting any ourput in
the console. This is where it should be, right? I think I have this set up
correctly. I'm using method=3 which only requires nc and d to be
specified. Any ideas why I'm not seeing output?
Here is the entire output:
> ## sshc.ssc: sample size calculation for historical control studies
> ## J. Jack Lee (jjlee@mdanderson.org) and Chi-hong Tseng
> ## Department of Biostatistics, Univ. of Texas M.D. Anderson Cancer Center
> ##
> ## 3/1/99
> ## updated 6/7/00: add loess
> ##------------------------------------------------------------------
> ######## Required Input:
> #
> # rc number of response in historical control group
> # nc sample size in historical control
> # d target improvement = Pe - Pc
> # method 1=method based on the randomized design
> # 2=Makuch & Simon method (Makuch RW, Simon RM. Sample size
considerations
> # for non-randomized comparative studies. J of Chron Dis 1980;
3:175-181.
> # 3=uniform power method
> ######## optional Input:
> #
> # alpha size of the test
> # power desired power of the test
> # tol convergence criterion for methods 1 & 2 in terms of sample
size
> # tol1 convergence criterion for method 3 at any given obs Rc in terms of
difference
> # of expected power from target
> # tol2 overall convergence criterion for method 3 as the max absolute
deviation
> # of expected power from target for all Rc
> # cc range of multiplicative constant applied to the initial values ne
> # l.span smoothing constant for loess
> #
> # Note: rc is required for methods 1 and 2 but not 3
> # method 3 return the sample size need for rc=0 to (1-d)*nc
> #
> ######## Output
> # for methdos 1 & 2: return the sample size needed for the experimental
group (1 number)
> # for given rc, nc, d, alpha, and power
> # for method 3: return the profile of sample size needed for given nc,
d, alpha, and power
> # vector $ne contains the sample size corresponding to
rc=0, 1, 2, ... nc*(1-d)
> # vector $Ep contains the expected power corresponding
to
> # the true pc = (0, 1, 2, ..., nc*(1-d)) / nc
> #
> #------------------------------------------------------------------
> sshc<-function(rc, nc=500, d=.5, method=3, alpha=0.05, power=0.8,
+ tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5)
+ {
+ ### for method 1
+ if (method==1) {
+ ne1<-ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01)
+ return(ne=ne1)
+ }
+ ### for method 2
+ if (method==2) {
+ ne<-nc
+ ne1<-nc+50
+ while(abs(ne-ne1)>tol & ne1<100000){
+ ne<-ne1
+ pe<-d+rc/nc
+ ne1<-nef(rc,nc,pe*ne,ne,alpha,power)
+ ## if(is.na(ne1))
print(paste('rc=',rc,',nc=',nc,',pe=',pe,',ne=',ne))
+ }
+ if (ne1>100000) return(NA)
+ else return(ne=ne1)
+ }
+ ### for method 3
+ if (method==3) {
+ if (tol1 > tol2/10) tol1<-tol2/10
+ ncstar<-(1-d)*nc
+ pc<-(0:ncstar)/nc
+ ne<-rep(NA,ncstar + 1)
+ for (i in (0:ncstar))
+ { ne[i+1]<-ss.rand(i,nc,d,alpha=.05,power=.8,tol=.01)
+ }
+ plot(pc,ne,type='l',ylim=c(0,max(ne)*1.5))
+ ans<-c.searchd(nc, d, ne, alpha, power, cc, tol1)
+ ### check overall absolute deviance
+ old.abs.dev<-sum(abs(ans$Ep-power))
+ ##bad<-0
+ print(round(ans$Ep,4))
+ print(round(ans$ne,2))
+ lines(pc,ans$ne,lty=1,col=8)
+ old.ne<-ans$ne
+ ##while(max(abs(ans$Ep-power))>tol2 & bad==0){ #### unnecessary ##
+ while(max(abs(ans$Ep-power))>tol2){
+ ans<-c.searchd(nc, d, ans$ne, alpha, power, cc, tol1)
+ abs.dev<-sum(abs(ans$Ep-power))
+ print(paste(" old.abs.dev=",old.abs.dev))
+ print(paste(" abs.dev=",abs.dev))
+ ##if (abs.dev > old.abs.dev) { bad<-1}
+ old.abs.dev<-abs.dev
+ print(round(ans$Ep,4))
+ print(round(ans$ne,2))
+ lines(pc,old.ne,lty=1,col=1)
+ lines(pc,ans$ne,lty=1,col=8)
+ ### add convex
+ ans$ne<-convex(pc,ans$ne)$wy
+ ### add loess
+ ###old.ne<-ans$ne
+ loess.ne<-loess(ans$ne ~ pc, span=l.span)
+ lines(pc,loess.ne$fit,lty=1,col=4)
+ old.ne<-loess.ne$fit
+ ###readline()
+ }
+ return(ne=ans$ne, Ep=ans$Ep)
+ }
+ }>
> ## needed for method 1
> nef2<-function(rc,nc,re,ne,alpha,power){
+ za<-qnorm(1-alpha)
+ zb<-qnorm(power)
+ xe<-asin(sqrt((re+0.375)/(ne+0.75)))
+ xc<-asin(sqrt((rc+0.375)/(nc+0.75)))
+ ans<- 1/(4*(xc-xe)^2/(za+zb)^2-1/(nc+0.5)) - 0.5
+ return(ans)
+ }> ## needed for method 2
> nef<-function(rc,nc,re,ne,alpha,power){
+ za<-qnorm(1-alpha)
+ zb<-qnorm(power)
+ xe<-asin(sqrt((re+0.375)/(ne+0.75)))
+ xc<-asin(sqrt((rc+0.375)/(nc+0.75)))
+ ans<-(za*sqrt(1+(ne+0.5)/(nc+0.5))+zb)^2/(2*(xe-xc))^2-0.5
+ return(ans)
+ }> ## needed for method 3
> c.searchd<-function(nc, d, ne, alpha=0.05, power=0.8,
cc=c(0.1,2),tol1=0.0001){
+ #---------------------------
+ # nc sample size of control group
+ # d the differece to detect between control and experiment
+ # ne vector of starting sample size of experiment group
+ # corresonding to rc of 0 to nc*(1-d)
+ # alpha size of test
+ # power target power
+ # cc pre-screen vector of constant c, the range should cover the
+ # the value of cc that has expected power
+ # tol1 the allowance between the expceted power and target power
+ #---------------------------
+ pc<-(0:((1-d)*nc))/nc
+ ncl<-length(pc)
+ ne.old<-ne
+ ne.old1<-ne.old
+ ### sweeping forward
+ for(i in 1:ncl){
+ cmin<-cc[1]
+ cmax<-cc[2]
+ ### fixed cci<-cmax bug
+ cci <-1
+ lhood<-dbinom((i:ncl)-1,nc,pc[i])
+ ne[i:ncl]<-(1+(cci-1)*(lhood/lhood[1])) * ne.old1[i:ncl]
+ Ep0 <-Epower(nc, d, ne, pc, alpha)
+ while(abs(Ep0[i]-power)>tol1){
+ if(Ep0[i]<power) cmin<-cci
+ else cmax<-cci
+ cci<-(cmax+cmin)/2
+ ne[i:ncl]<-(1+(cci-1)*(lhood/lhood[1])) * ne.old1[i:ncl]
+ Ep0<-Epower(nc, d, ne, pc, alpha)
+ }
+ ne.old1<-ne
+ }
+ ne1<-ne
+ ### sweeping backward -- ncl:i
+ ne.old2<-ne.old
+ ne <-ne.old
+ for(i in ncl:1){
+ cmin<-cc[1]
+ cmax<-cc[2]
+ ### fixed cci<-cmax bug
+ cci <-1
+ lhood<-dbinom((ncl:i)-1,nc,pc[i])
+ lenl <-length(lhood)
+ ne[ncl:i]<-(1+(cci-1)*(lhood/lhood[lenl]))*ne.old2[ncl:i]
+ Ep0 <-Epower(nc, d, cci*ne, pc, alpha)
+ while(abs(Ep0[i]-power)>tol1){
+ if(Ep0[i]<power) cmin<-cci
+ else cmax<-cci
+ cci<-(cmax+cmin)/2
+ ne[ncl:i]<-(1+(cci-1)*(lhood/lhood[lenl]))*ne.old2[ncl:i]
+ Ep0<-Epower(nc, d, ne, pc, alpha)
+ }
+ ne.old2<-ne
+ }
+ ne2<-ne
+ ne<-(ne1+ne2)/2
+ #cat(ccc*ne)
+ Ep1<-Epower(nc, d, ne, pc, alpha)
+ return(ne=ne, Ep=Ep1)
+ }> ###
> vertex<-function(x,y)
+ { n<-length(x)
+ vx<-x[1]
+ vy<-y[1]
+ vp<-1
+ up<-T
+ for (i in (2:n))
+ { if (up)
+ { if (y[i-1] > y[i])
+ {vx<-c(vx,x[i-1])
+ vy<-c(vy,y[i-1])
+ vp<-c(vp,i-1)
+ up<-F
+ }
+ }
+ else
+ { if (y[i-1] < y[i]) up<-T
+ }
+ }
+ vx<-c(vx,x[n])
+ vy<-c(vy,y[n])
+ vp<-c(vp,n)
+ return(vx=vx,vy=vy,vp=vp)
+ }> ###
> convex<-function(x,y)
+ {
+ n<-length(x)
+ ans<-vertex(x,y)
+ len<-length(ans$vx)
+ while (len>3)
+ {
+ #cat("x=",x,"\n")
+ #cat("y=",y,"\n")
+ newx<-x[1:(ans$vp[2]-1)]
+ newy<-y[1:(ans$vp[2]-1)]
+ for (i in (2:(len-1)))
+ {
+ newx<-c(newx,x[ans$vp[i]])
+ newy<-c(newy,y[ans$vp[i]])
+ }
+ newx<-c(newx,x[(ans$vp[len-1]+1):n])
+ newy<-c(newy,y[(ans$vp[len-1]+1):n])
+ y<-approx(newx,newy,xout=x)$y
+ #cat("new y=",y,"\n")
+ ans<-vertex(x,y)
+ len<-length(ans$vx)
+ #cat("vx=",ans$vx,"\n")
+ #cat("vy=",ans$vy,"\n")
+ }
+ return(wx=x,wy=y)}> ###
> Epower<-function(nc, d, ne, pc = (0:((1 - d) * nc))/nc, alpha = 0.05)
+ {
+ #-------------------------------------
+ # nc sample size in historical control
+ # d the increase of response rate between historical and experiment
+ # ne sample size of corresonding rc of 0 to nc*(1-d)
+ # pc the response rate of control group, where we compute the
+ # expected power
+ # alpha the size of test
+ #-------------------------------------
+ kk <- length(pc)
+ rc <- 0:(nc * (1 - d))
+ pp <- rep(NA, kk)
+ ppp <- rep(NA, kk)
+ for(i in 1:(kk)) {
+ pe <- pc[i] + d
+ lhood <- dbinom(rc, nc, pc[i])
+ pp <- power1.f(rc, nc, ne, pe, alpha)
+ ppp[i] <- sum(pp * lhood)/sum(lhood)
+ }
+ return(ppp)
+ }>
> # adapted from the old biss2
> ss.rand<-function(rc,nc,d,alpha=.05,power=.8,tol=.01)
+ {
+ ne<-nc
+ ne1<-nc+50
+ while(abs(ne-ne1)>tol & ne1<100000){
+ ne<-ne1
+ pe<-d+rc/nc
+ ne1<-nef2(rc,nc,pe*ne,ne,alpha,power)
+
+ ## if(is.na(ne1))
print(paste('rc=',rc,',nc=',nc,',pe=',pe,',ne=',ne))
+ }
+ if (ne1>100000) return(NA)
+ else return(ne1)
+ }> ###
> power1.f<-function(rc,nc,ne,pie,alpha=0.05){
+ #-------------------------------------
+ # rcnumber of response in historical control
+ # ncsample size in historical control
+ # ne sample size in experitment group
+ # pietrue response rate for experiment group
+ # alphasize of the test
+ #-------------------------------------
+
+ za<-qnorm(1-alpha)
+ re<-ne*pie
+ xe<-asin(sqrt((re+0.375)/(ne+0.75)))
+ xc<-asin(sqrt((rc+0.375)/(nc+0.75)))
+ ans<-za*sqrt(1+(ne+0.5)/(nc+0.5))-(xe-xc)/sqrt(1/(4*(ne+0.5)))
+ return(1-pnorm(ans))
+ }
[[alternative HTML version deleted]]
Hi Scott, I am not familiar with S-Plus (though many aspects are quite similar to R). I will say that your function looks approximately correct. I am not familiar with the ss.rand function. I searched, and found some things that I suspect are similar in the packages MBESS, but without knowing more about it from S-Plus, it is tough to make a testable example. Do you have access to S-Plus? Can you provide more information about this function, what it does, what is like, etc.? There are some active members of this list who are quite familiar with S-Plus so one of them may be more insightful. Cheers, Josh On Tue, Oct 4, 2011 at 6:53 PM, Scott Raynaud <scott.raynaud at yahoo.com> wrote:> I'm trying to convert an S-Plus program to R.? Since I'm a SAS programmer I'm not facile is either S-Plus or R, so I need some help.? All I did was convert the underscores in S-Plus to the assignment operator <-.? Here are the first few lines of the S-Plus file: > > sshc _ function(rc, nc, d, method, alpha=0.05, power=0.8, > ???????????? tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5) > { > ### for method 1 > if (method==1) { > ne1 _ ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01) > return(ne=ne1) > ?????????????? } > > > My?translation looks like this: > > sshc<-function(rc, nc=500, d=.5, method=3, alpha=0.05, power=0.8, > ????????????? tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5) > { > ### for method 1 > if (method==1) { > ?ne1<-ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01) > ?return(ne=ne1) > ?????????????? } > > The program runs without throwing errors, but I'm not getting any ourput in the console.? This is where it should be, right?? I think I have this set up correctly.? I'm using method=3 which only requires nc and d to be specified.? Any ideas why I'm not seeing output? > > Here is the entire output: > >> ## sshc.ssc: sample size calculation for historical control studies >> ## J. Jack Lee (jjlee at mdanderson.org) and Chi-hong Tseng >> ## Department of Biostatistics, Univ. of Texas M.D. Anderson Cancer Center >> ## >> ## 3/1/99 >> ## updated 6/7/00: add loess >> ##------------------------------------------------------------------ >> ######## Required Input: >> # >> # rc???? number of response in historical control group >> # nc???? sample size in historical control >> # d????? target improvement = Pe - Pc >> # method 1=method based on the randomized design >> #??????? 2=Makuch & Simon method (Makuch RW, Simon RM. Sample size considerations >> #????????? for non-randomized comparative studies. J of Chron Dis 1980; 3:175-181. >> #??????? 3=uniform power method >> ######## optional Input: >> # >> # alpha? size of the test >> # power? desired power of the test >> # tol??? convergence criterion for methods 1 & 2 in terms of sample size >> # tol1?? convergence criterion for method 3 at any given obs Rc in terms of difference >> #????????? of expected power from target >> # tol2?? overall convergence criterion for method 3 as the max absolute deviation >> #????????? of expected power from target for all Rc >> # cc???? range of multiplicative constant applied to the initial values ne >> # l.span smoothing constant for loess >> # >> # Note:? rc is required for methods 1 and 2 but not 3 >> #??????? method 3 return the sample size need for rc=0 to (1-d)*nc >> # >> ######## Output >> # for methdos 1 & 2: return the sample size needed for the experimental group (1 number) >> #??????????????????? for given rc, nc, d, alpha, and power >> # for method 3:????? return the profile of sample size needed for given nc, d, alpha, and power >> #??????????????????? vector $ne contains the sample size corresponding to rc=0, 1, 2, ... nc*(1-d) >> #??????????????????? vector $Ep contains the expected power corresponding to >> #????????????????????? the true pc = (0, 1, 2, ..., nc*(1-d)) / nc >> # >> #------------------------------------------------------------------ >> sshc<-function(rc, nc=500, d=.5, method=3, alpha=0.05, power=0.8, > +????????????? tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5) > + { > + ### for method 1 > + if (method==1) { > + ne1<-ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01) > + return(ne=ne1) > +??????????????? } > + ### for method 2 > + if (method==2) { > + ne<-nc > + ne1<-nc+50 > + while(abs(ne-ne1)>tol & ne1<100000){ > + ne<-ne1 > + pe<-d+rc/nc > + ne1<-nef(rc,nc,pe*ne,ne,alpha,power) > + ## if(is.na(ne1)) print(paste('rc=',rc,',nc=',nc,',pe=',pe,',ne=',ne)) > + } > + if (ne1>100000) return(NA) > + else return(ne=ne1) > + } > + ### for method 3 > + if (method==3) { > + if (tol1 > tol2/10) tol1<-tol2/10 > + ncstar<-(1-d)*nc > + pc<-(0:ncstar)/nc > + ne<-rep(NA,ncstar + 1) > + for (i in (0:ncstar)) > + { ne[i+1]<-ss.rand(i,nc,d,alpha=.05,power=.8,tol=.01) > + } > + plot(pc,ne,type='l',ylim=c(0,max(ne)*1.5)) > + ans<-c.searchd(nc, d, ne, alpha, power, cc, tol1) > + ### check overall absolute deviance > + old.abs.dev<-sum(abs(ans$Ep-power)) > + ##bad<-0 > + print(round(ans$Ep,4)) > + print(round(ans$ne,2)) > + lines(pc,ans$ne,lty=1,col=8) > + old.ne<-ans$ne > + ##while(max(abs(ans$Ep-power))>tol2 & bad==0){? #### unnecessary ## > + while(max(abs(ans$Ep-power))>tol2){ > + ans<-c.searchd(nc, d, ans$ne, alpha, power, cc, tol1) > + abs.dev<-sum(abs(ans$Ep-power)) > + print(paste(" old.abs.dev=",old.abs.dev)) > + print(paste("???? abs.dev=",abs.dev)) > + ##if (abs.dev > old.abs.dev) { bad<-1} > + old.abs.dev<-abs.dev > + print(round(ans$Ep,4)) > + print(round(ans$ne,2)) > + lines(pc,old.ne,lty=1,col=1) > + lines(pc,ans$ne,lty=1,col=8) > + ### add convex > + ans$ne<-convex(pc,ans$ne)$wy > + ### add loess > + ###old.ne<-ans$ne > + loess.ne<-loess(ans$ne ~ pc, span=l.span) > + lines(pc,loess.ne$fit,lty=1,col=4) > + old.ne<-loess.ne$fit > + ###readline() > + } > + return(ne=ans$ne, Ep=ans$Ep) > +??????????????? } > + } >> >> ## needed for method 1 >> nef2<-function(rc,nc,re,ne,alpha,power){ > + za<-qnorm(1-alpha) > + zb<-qnorm(power) > + xe<-asin(sqrt((re+0.375)/(ne+0.75))) > + xc<-asin(sqrt((rc+0.375)/(nc+0.75))) > + ans<- 1/(4*(xc-xe)^2/(za+zb)^2-1/(nc+0.5)) - 0.5 > + return(ans) > + } >> ## needed for method 2 >> nef<-function(rc,nc,re,ne,alpha,power){ > + za<-qnorm(1-alpha) > + zb<-qnorm(power) > + xe<-asin(sqrt((re+0.375)/(ne+0.75))) > + xc<-asin(sqrt((rc+0.375)/(nc+0.75))) > + ans<-(za*sqrt(1+(ne+0.5)/(nc+0.5))+zb)^2/(2*(xe-xc))^2-0.5 > + return(ans) > + } >> ## needed for method 3 >> c.searchd<-function(nc, d, ne, alpha=0.05, power=0.8, cc=c(0.1,2),tol1=0.0001){ > + #--------------------------- > + # nc???? sample size of control group > + # d????? the differece to detect between control and experiment > + # ne???? vector of starting sample size of experiment group > + #??? corresonding to rc of 0 to nc*(1-d) > + # alpha? size of test > + # power? target power > + # cc? pre-screen vector of constant c, the range should cover the > + #??? the value of cc that has expected power > + # tol1?? the allowance between the expceted power and target power > + #--------------------------- > + pc<-(0:((1-d)*nc))/nc > + ncl<-length(pc) > + ne.old<-ne > + ne.old1<-ne.old > + ### sweeping forward > + for(i in 1:ncl){ > + cmin<-cc[1] > + cmax<-cc[2] > + ### fixed cci<-cmax bug > + cci <-1 > + lhood<-dbinom((i:ncl)-1,nc,pc[i]) > + ne[i:ncl]<-(1+(cci-1)*(lhood/lhood[1])) * ne.old1[i:ncl] > + Ep0 <-Epower(nc, d, ne, pc, alpha) > + while(abs(Ep0[i]-power)>tol1){ > + if(Ep0[i]<power) cmin<-cci > + else cmax<-cci > + cci<-(cmax+cmin)/2 > + ne[i:ncl]<-(1+(cci-1)*(lhood/lhood[1])) * ne.old1[i:ncl] > + Ep0<-Epower(nc, d, ne, pc, alpha) > + } > +? ne.old1<-ne > + } > + ne1<-ne > + ### sweeping backward -- ncl:i > + ne.old2<-ne.old > + ne???? <-ne.old > + for(i in ncl:1){ > + cmin<-cc[1] > + cmax<-cc[2] > + ### fixed cci<-cmax bug > + cci <-1 > + lhood<-dbinom((ncl:i)-1,nc,pc[i]) > + lenl <-length(lhood) > + ne[ncl:i]<-(1+(cci-1)*(lhood/lhood[lenl]))*ne.old2[ncl:i] > + Ep0 <-Epower(nc, d, cci*ne, pc, alpha) > + while(abs(Ep0[i]-power)>tol1){ > + if(Ep0[i]<power) cmin<-cci > + else cmax<-cci > + cci<-(cmax+cmin)/2 > + ne[ncl:i]<-(1+(cci-1)*(lhood/lhood[lenl]))*ne.old2[ncl:i] > + Ep0<-Epower(nc, d, ne, pc, alpha) > + } > +? ne.old2<-ne > + } > + ne2<-ne > + ne<-(ne1+ne2)/2 > + #cat(ccc*ne) > + Ep1<-Epower(nc, d, ne, pc, alpha) > + return(ne=ne, Ep=Ep1) > + } >> ### >> vertex<-function(x,y) > + { n<-length(x) > + vx<-x[1] > + vy<-y[1] > + vp<-1 > + up<-T > + for (i in (2:n)) > + { if (up) > + { if (y[i-1] > y[i]) > + {vx<-c(vx,x[i-1]) > +? vy<-c(vy,y[i-1]) > +? vp<-c(vp,i-1) > +? up<-F > + } > + } > + else > + { if (y[i-1] < y[i]) up<-T > + } > + } > + vx<-c(vx,x[n]) > + vy<-c(vy,y[n]) > + vp<-c(vp,n) > + return(vx=vx,vy=vy,vp=vp) > + } >> ### >> convex<-function(x,y) > + { > + n<-length(x) > + ans<-vertex(x,y) > + len<-length(ans$vx) > + while (len>3) > + { > + #cat("x=",x,"\n") > + #cat("y=",y,"\n") > + newx<-x[1:(ans$vp[2]-1)] > + newy<-y[1:(ans$vp[2]-1)] > + for (i in (2:(len-1))) > + { > +? newx<-c(newx,x[ans$vp[i]]) > + newy<-c(newy,y[ans$vp[i]]) > + } > + newx<-c(newx,x[(ans$vp[len-1]+1):n]) > + newy<-c(newy,y[(ans$vp[len-1]+1):n]) > + y<-approx(newx,newy,xout=x)$y > + #cat("new y=",y,"\n") > + ans<-vertex(x,y) > + len<-length(ans$vx) > + #cat("vx=",ans$vx,"\n") > + #cat("vy=",ans$vy,"\n") > + } > + return(wx=x,wy=y)} >> ### >> Epower<-function(nc, d, ne, pc = (0:((1 - d) * nc))/nc, alpha = 0.05) > + { > + #------------------------------------- > + # nc???? sample size in historical control > + # d????? the increase of response rate between historical and experiment > + # ne???? sample size of corresonding rc of 0 to nc*(1-d) > + # pc???? the response rate of control group, where we compute the > + #??????? expected power > + # alpha? the size of test > + #------------------------------------- > + kk <- length(pc) > + rc <- 0:(nc * (1 - d)) > + pp <- rep(NA, kk) > + ppp <- rep(NA, kk) > + for(i in 1:(kk)) { > + pe <- pc[i] + d > + lhood <- dbinom(rc, nc, pc[i]) > + pp <- power1.f(rc, nc, ne, pe, alpha) > + ppp[i] <- sum(pp * lhood)/sum(lhood) > + } > + return(ppp) > + } >> >> # adapted from the old biss2 >> ss.rand<-function(rc,nc,d,alpha=.05,power=.8,tol=.01) > + { > + ne<-nc > + ne1<-nc+50 > + while(abs(ne-ne1)>tol & ne1<100000){ > + ne<-ne1 > + pe<-d+rc/nc > + ne1<-nef2(rc,nc,pe*ne,ne,alpha,power) > + > + ## if(is.na(ne1)) print(paste('rc=',rc,',nc=',nc,',pe=',pe,',ne=',ne)) > + } > + if (ne1>100000) return(NA) > + else return(ne1) > + } >> ### >> power1.f<-function(rc,nc,ne,pie,alpha=0.05){ > + #------------------------------------- > + # rcnumber of response in historical control > + # ncsample size in historical control > + # ne??? sample size in experitment group > + # pietrue response rate for experiment group > + # alphasize of the test > + #------------------------------------- > + > + za<-qnorm(1-alpha) > + re<-ne*pie > + xe<-asin(sqrt((re+0.375)/(ne+0.75))) > + xc<-asin(sqrt((rc+0.375)/(nc+0.75))) > + ans<-za*sqrt(1+(ne+0.5)/(nc+0.5))-(xe-xc)/sqrt(1/(4*(ne+0.5))) > + return(1-pnorm(ans)) > + } > > ? ? ? ?[[alternative HTML version deleted]] > > > ______________________________________________ > R-help at 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. > >-- Joshua Wiley Ph.D. Student, Health Psychology Programmer Analyst II, ATS Statistical Consulting Group University of California, Los Angeles https://joshuawiley.com/
I had thought that the problem might
be:
return(ne=ne1)
since R doesn't support that any more.
But when I tried it, I got results
(just without the name on the output).
Better would be to change that line to:
list(ne=ne1)
('return' is seldom necessary in either
R or S+.)
I'd suggest putting 'print' or 'cat' statements
in to try to figure out where things go
wrong.
You'll find other hints in 'S Poetry' and
'The R Inferno'. There might be positive
probability that at least one of those hints
will be useful. Both those are available
on www.burns-stat.com
On 05/10/2011 02:53, Scott Raynaud wrote:> I'm trying to convert an S-Plus program to R. Since I'm a SAS
programmer I'm not facile is either S-Plus or R, so I need some help. All I
did was convert the underscores in S-Plus to the assignment operator<-. Here
are the first few lines of the S-Plus file:
>
> sshc _ function(rc, nc, d, method, alpha=0.05, power=0.8,
> tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5)
> {
> ### for method 1
> if (method==1) {
> ne1 _ ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01)
> return(ne=ne1)
> }
>
>
> My translation looks like this:
>
> sshc<-function(rc, nc=500, d=.5, method=3, alpha=0.05, power=0.8,
> tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5)
> {
> ### for method 1
> if (method==1) {
> ne1<-ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01)
> return(ne=ne1)
> }
>
> The program runs without throwing errors, but I'm not getting any
ourput in the console. This is where it should be, right? I think I have this
set up correctly. I'm using method=3 which only requires nc and d to be
specified. Any ideas why I'm not seeing output?
>
> Here is the entire output:
>
>> ## sshc.ssc: sample size calculation for historical control studies
>> ## J. Jack Lee (jjlee at mdanderson.org) and Chi-hong Tseng
>> ## Department of Biostatistics, Univ. of Texas M.D. Anderson Cancer
Center
>> ##
>> ## 3/1/99
>> ## updated 6/7/00: add loess
>> ##------------------------------------------------------------------
>> ######## Required Input:
>> #
>> # rc number of response in historical control group
>> # nc sample size in historical control
>> # d target improvement = Pe - Pc
>> # method 1=method based on the randomized design
>> # 2=Makuch& Simon method (Makuch RW, Simon RM. Sample size
considerations
>> # for non-randomized comparative studies. J of Chron Dis 1980;
3:175-181.
>> # 3=uniform power method
>> ######## optional Input:
>> #
>> # alpha size of the test
>> # power desired power of the test
>> # tol convergence criterion for methods 1& 2 in terms of sample
size
>> # tol1 convergence criterion for method 3 at any given obs Rc in
terms of difference
>> # of expected power from target
>> # tol2 overall convergence criterion for method 3 as the max absolute
deviation
>> # of expected power from target for all Rc
>> # cc range of multiplicative constant applied to the initial values
ne
>> # l.span smoothing constant for loess
>> #
>> # Note: rc is required for methods 1 and 2 but not 3
>> # method 3 return the sample size need for rc=0 to (1-d)*nc
>> #
>> ######## Output
>> # for methdos 1& 2: return the sample size needed for the
experimental group (1 number)
>> # for given rc, nc, d, alpha, and power
>> # for method 3: return the profile of sample size needed for given
nc, d, alpha, and power
>> # vector $ne contains the sample size corresponding
to rc=0, 1, 2, ... nc*(1-d)
>> # vector $Ep contains the expected power
corresponding to
>> # the true pc = (0, 1, 2, ..., nc*(1-d)) / nc
>> #
>> #------------------------------------------------------------------
>> sshc<-function(rc, nc=500, d=.5, method=3, alpha=0.05, power=0.8,
> + tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5)
> + {
> + ### for method 1
> + if (method==1) {
> + ne1<-ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01)
> + return(ne=ne1)
> + }
> + ### for method 2
> + if (method==2) {
> + ne<-nc
> + ne1<-nc+50
> + while(abs(ne-ne1)>tol& ne1<100000){
> + ne<-ne1
> + pe<-d+rc/nc
> + ne1<-nef(rc,nc,pe*ne,ne,alpha,power)
> + ## if(is.na(ne1))
print(paste('rc=',rc,',nc=',nc,',pe=',pe,',ne=',ne))
> + }
> + if (ne1>100000) return(NA)
> + else return(ne=ne1)
> + }
> + ### for method 3
> + if (method==3) {
> + if (tol1> tol2/10) tol1<-tol2/10
> + ncstar<-(1-d)*nc
> + pc<-(0:ncstar)/nc
> + ne<-rep(NA,ncstar + 1)
> + for (i in (0:ncstar))
> + { ne[i+1]<-ss.rand(i,nc,d,alpha=.05,power=.8,tol=.01)
> + }
> + plot(pc,ne,type='l',ylim=c(0,max(ne)*1.5))
> + ans<-c.searchd(nc, d, ne, alpha, power, cc, tol1)
> + ### check overall absolute deviance
> + old.abs.dev<-sum(abs(ans$Ep-power))
> + ##bad<-0
> + print(round(ans$Ep,4))
> + print(round(ans$ne,2))
> + lines(pc,ans$ne,lty=1,col=8)
> + old.ne<-ans$ne
> + ##while(max(abs(ans$Ep-power))>tol2& bad==0){ #### unnecessary
##
> + while(max(abs(ans$Ep-power))>tol2){
> + ans<-c.searchd(nc, d, ans$ne, alpha, power, cc, tol1)
> + abs.dev<-sum(abs(ans$Ep-power))
> + print(paste(" old.abs.dev=",old.abs.dev))
> + print(paste(" abs.dev=",abs.dev))
> + ##if (abs.dev> old.abs.dev) { bad<-1}
> + old.abs.dev<-abs.dev
> + print(round(ans$Ep,4))
> + print(round(ans$ne,2))
> + lines(pc,old.ne,lty=1,col=1)
> + lines(pc,ans$ne,lty=1,col=8)
> + ### add convex
> + ans$ne<-convex(pc,ans$ne)$wy
> + ### add loess
> + ###old.ne<-ans$ne
> + loess.ne<-loess(ans$ne ~ pc, span=l.span)
> + lines(pc,loess.ne$fit,lty=1,col=4)
> + old.ne<-loess.ne$fit
> + ###readline()
> + }
> + return(ne=ans$ne, Ep=ans$Ep)
> + }
> + }
>>
>> ## needed for method 1
>> nef2<-function(rc,nc,re,ne,alpha,power){
> + za<-qnorm(1-alpha)
> + zb<-qnorm(power)
> + xe<-asin(sqrt((re+0.375)/(ne+0.75)))
> + xc<-asin(sqrt((rc+0.375)/(nc+0.75)))
> + ans<- 1/(4*(xc-xe)^2/(za+zb)^2-1/(nc+0.5)) - 0.5
> + return(ans)
> + }
>> ## needed for method 2
>> nef<-function(rc,nc,re,ne,alpha,power){
> + za<-qnorm(1-alpha)
> + zb<-qnorm(power)
> + xe<-asin(sqrt((re+0.375)/(ne+0.75)))
> + xc<-asin(sqrt((rc+0.375)/(nc+0.75)))
> + ans<-(za*sqrt(1+(ne+0.5)/(nc+0.5))+zb)^2/(2*(xe-xc))^2-0.5
> + return(ans)
> + }
>> ## needed for method 3
>> c.searchd<-function(nc, d, ne, alpha=0.05, power=0.8,
cc=c(0.1,2),tol1=0.0001){
> + #---------------------------
> + # nc sample size of control group
> + # d the differece to detect between control and experiment
> + # ne vector of starting sample size of experiment group
> + # corresonding to rc of 0 to nc*(1-d)
> + # alpha size of test
> + # power target power
> + # cc pre-screen vector of constant c, the range should cover the
> + # the value of cc that has expected power
> + # tol1 the allowance between the expceted power and target power
> + #---------------------------
> + pc<-(0:((1-d)*nc))/nc
> + ncl<-length(pc)
> + ne.old<-ne
> + ne.old1<-ne.old
> + ### sweeping forward
> + for(i in 1:ncl){
> + cmin<-cc[1]
> + cmax<-cc[2]
> + ### fixed cci<-cmax bug
> + cci<-1
> + lhood<-dbinom((i:ncl)-1,nc,pc[i])
> + ne[i:ncl]<-(1+(cci-1)*(lhood/lhood[1])) * ne.old1[i:ncl]
> + Ep0<-Epower(nc, d, ne, pc, alpha)
> + while(abs(Ep0[i]-power)>tol1){
> + if(Ep0[i]<power) cmin<-cci
> + else cmax<-cci
> + cci<-(cmax+cmin)/2
> + ne[i:ncl]<-(1+(cci-1)*(lhood/lhood[1])) * ne.old1[i:ncl]
> + Ep0<-Epower(nc, d, ne, pc, alpha)
> + }
> + ne.old1<-ne
> + }
> + ne1<-ne
> + ### sweeping backward -- ncl:i
> + ne.old2<-ne.old
> + ne<-ne.old
> + for(i in ncl:1){
> + cmin<-cc[1]
> + cmax<-cc[2]
> + ### fixed cci<-cmax bug
> + cci<-1
> + lhood<-dbinom((ncl:i)-1,nc,pc[i])
> + lenl<-length(lhood)
> + ne[ncl:i]<-(1+(cci-1)*(lhood/lhood[lenl]))*ne.old2[ncl:i]
> + Ep0<-Epower(nc, d, cci*ne, pc, alpha)
> + while(abs(Ep0[i]-power)>tol1){
> + if(Ep0[i]<power) cmin<-cci
> + else cmax<-cci
> + cci<-(cmax+cmin)/2
> + ne[ncl:i]<-(1+(cci-1)*(lhood/lhood[lenl]))*ne.old2[ncl:i]
> + Ep0<-Epower(nc, d, ne, pc, alpha)
> + }
> + ne.old2<-ne
> + }
> + ne2<-ne
> + ne<-(ne1+ne2)/2
> + #cat(ccc*ne)
> + Ep1<-Epower(nc, d, ne, pc, alpha)
> + return(ne=ne, Ep=Ep1)
> + }
>> ###
>> vertex<-function(x,y)
> + { n<-length(x)
> + vx<-x[1]
> + vy<-y[1]
> + vp<-1
> + up<-T
> + for (i in (2:n))
> + { if (up)
> + { if (y[i-1]> y[i])
> + {vx<-c(vx,x[i-1])
> + vy<-c(vy,y[i-1])
> + vp<-c(vp,i-1)
> + up<-F
> + }
> + }
> + else
> + { if (y[i-1]< y[i]) up<-T
> + }
> + }
> + vx<-c(vx,x[n])
> + vy<-c(vy,y[n])
> + vp<-c(vp,n)
> + return(vx=vx,vy=vy,vp=vp)
> + }
>> ###
>> convex<-function(x,y)
> + {
> + n<-length(x)
> + ans<-vertex(x,y)
> + len<-length(ans$vx)
> + while (len>3)
> + {
> + #cat("x=",x,"\n")
> + #cat("y=",y,"\n")
> + newx<-x[1:(ans$vp[2]-1)]
> + newy<-y[1:(ans$vp[2]-1)]
> + for (i in (2:(len-1)))
> + {
> + newx<-c(newx,x[ans$vp[i]])
> + newy<-c(newy,y[ans$vp[i]])
> + }
> + newx<-c(newx,x[(ans$vp[len-1]+1):n])
> + newy<-c(newy,y[(ans$vp[len-1]+1):n])
> + y<-approx(newx,newy,xout=x)$y
> + #cat("new y=",y,"\n")
> + ans<-vertex(x,y)
> + len<-length(ans$vx)
> + #cat("vx=",ans$vx,"\n")
> + #cat("vy=",ans$vy,"\n")
> + }
> + return(wx=x,wy=y)}
>> ###
>> Epower<-function(nc, d, ne, pc = (0:((1 - d) * nc))/nc, alpha =
0.05)
> + {
> + #-------------------------------------
> + # nc sample size in historical control
> + # d the increase of response rate between historical and experiment
> + # ne sample size of corresonding rc of 0 to nc*(1-d)
> + # pc the response rate of control group, where we compute the
> + # expected power
> + # alpha the size of test
> + #-------------------------------------
> + kk<- length(pc)
> + rc<- 0:(nc * (1 - d))
> + pp<- rep(NA, kk)
> + ppp<- rep(NA, kk)
> + for(i in 1:(kk)) {
> + pe<- pc[i] + d
> + lhood<- dbinom(rc, nc, pc[i])
> + pp<- power1.f(rc, nc, ne, pe, alpha)
> + ppp[i]<- sum(pp * lhood)/sum(lhood)
> + }
> + return(ppp)
> + }
>>
>> # adapted from the old biss2
>> ss.rand<-function(rc,nc,d,alpha=.05,power=.8,tol=.01)
> + {
> + ne<-nc
> + ne1<-nc+50
> + while(abs(ne-ne1)>tol& ne1<100000){
> + ne<-ne1
> + pe<-d+rc/nc
> + ne1<-nef2(rc,nc,pe*ne,ne,alpha,power)
> +
> + ## if(is.na(ne1))
print(paste('rc=',rc,',nc=',nc,',pe=',pe,',ne=',ne))
> + }
> + if (ne1>100000) return(NA)
> + else return(ne1)
> + }
>> ###
>> power1.f<-function(rc,nc,ne,pie,alpha=0.05){
> + #-------------------------------------
> + # rcnumber of response in historical control
> + # ncsample size in historical control
> + # ne sample size in experitment group
> + # pietrue response rate for experiment group
> + # alphasize of the test
> + #-------------------------------------
> +
> + za<-qnorm(1-alpha)
> + re<-ne*pie
> + xe<-asin(sqrt((re+0.375)/(ne+0.75)))
> + xc<-asin(sqrt((rc+0.375)/(nc+0.75)))
> + ans<-za*sqrt(1+(ne+0.5)/(nc+0.5))-(xe-xc)/sqrt(1/(4*(ne+0.5)))
> + return(1-pnorm(ans))
> + }
>
> [[alternative HTML version deleted]]
>
>
>
>
> ______________________________________________
> R-help at 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.
--
Patrick Burns
pburns at pburns.seanet.com
twitter: @portfolioprobe
http://www.portfolioprobe.com/blog
http://www.burns-stat.com
(home of 'Some hints for the R beginner'
and 'The R Inferno')
On Wed, Oct 5, 2011 at 2:53 AM, Scott Raynaud <scott.raynaud at yahoo.com> wrote:> I'm trying to convert an S-Plus program to R.? Since I'm a SAS programmer I'm not facile is either S-Plus or R, so I need some help.? All I did was convert the underscores in S-Plus to the assignment operator <-.? Here are the first few lines of the S-Plus file: > > sshc _ function(rc, nc, d, method, alpha=0.05, power=0.8, > ???????????? tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5) > { > ### for method 1 > if (method==1) { > ne1 _ ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01) > return(ne=ne1) > ?????????????? } > > > My?translation looks like this: > > sshc<-function(rc, nc=500, d=.5, method=3, alpha=0.05, power=0.8, > ????????????? tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5) > { > ### for method 1 > if (method==1) { > ?ne1<-ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01) > ?return(ne=ne1) > ?????????????? } > > The program runs without throwing errors, but I'm not getting any ourput in the console.? This is where it should be, right?? I think I have this set up correctly.? I'm using method=3 which only requires nc and d to be specified.? Any ideas why I'm not seeing output?Long shot: the code you posted looked like (and hard to tell without indentation) just a bunch of function definitions. R won't actually do anything unless you call those functions with some parameters. So, when you say you get no output when you 'run' the code, what exactly do you mean by 'run' the code? What I would do is: 1. Put the code in a file called 'whatever.R'. 2. Start R, and do source("whatever.R"). That defines the functions. do "ls()" and you should see them. 3. Call one of the functions: sshc(100,10) I'd call that, in R terms, "calling the sshc function" rather than running anything. Barry
It looks like this code was written for S+ 4.5 (aka '2000') or before, which was based on S version 3. Try changing return(name1=value1, name2=value2) to return(list(name1=value1, name2=value2)) In S+ from 5.0 onwards return(name=value) or return(name1=value1, name2=value2) throws away the names and in R return only takes a single object (and also ignores the name). The c.search function in your code ends with return(ne=ne, Ep=Ep1) and the code calling c.search() acts as though the writer expects that function to return list(ne=ne, Ep=Ep1) ans <- c.searchd(nc, d, ne, alpha, power, cc, tol1) ... old.ne <- ans$ne Bill Dunlap Spotfire, TIBCO Software wdunlap tibco.com> -----Original Message----- > From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Scott Raynaud > Sent: Tuesday, October 04, 2011 6:53 PM > To: r-help at r-project.org > Subject: [R] SPlus to R > > I'm trying to convert an S-Plus program to R.? Since I'm a SAS programmer I'm not facile is either S- > Plus or R, so I need some help.? All I did was convert the underscores in S-Plus to the assignment > operator <-.? Here are the first few lines of the S-Plus file: > > sshc _ function(rc, nc, d, method, alpha=0.05, power=0.8, > ???????????? tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5) > { > ### for method 1 > if (method==1) { > ne1 _ ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01) > return(ne=ne1) > ?????????????? } > > > My?translation looks like this: > > sshc<-function(rc, nc=500, d=.5, method=3, alpha=0.05, power=0.8, > ????????????? tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5) > { > ### for method 1 > if (method==1) { > ?ne1<-ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01) > ?return(ne=ne1) > ?????????????? } > > The program runs without throwing errors, but I'm not getting any ourput in the console.? This is > where it should be, right?? I think I have this set up correctly.? I'm using method=3 which only > requires nc and d to be specified.? Any ideas why I'm not seeing output? > > Here is the entire output: > > > ## sshc.ssc: sample size calculation for historical control studies > > ## J. Jack Lee (jjlee at mdanderson.org) and Chi-hong Tseng > > ## Department of Biostatistics, Univ. of Texas M.D. Anderson Cancer Center > > ## > > ## 3/1/99 > > ## updated 6/7/00: add loess > > ##------------------------------------------------------------------ > > ######## Required Input: > > # > > # rc???? number of response in historical control group > > # nc???? sample size in historical control > > # d????? target improvement = Pe - Pc > > # method 1=method based on the randomized design > > #??????? 2=Makuch & Simon method (Makuch RW, Simon RM. Sample size considerations > > #????????? for non-randomized comparative studies. J of Chron Dis 1980; 3:175-181. > > #??????? 3=uniform power method > > ######## optional Input: > > # > > # alpha? size of the test > > # power? desired power of the test > > # tol??? convergence criterion for methods 1 & 2 in terms of sample size > > # tol1?? convergence criterion for method 3 at any given obs Rc in terms of difference > > #????????? of expected power from target > > # tol2?? overall convergence criterion for method 3 as the max absolute deviation > > #????????? of expected power from target for all Rc > > # cc???? range of multiplicative constant applied to the initial values ne > > # l.span smoothing constant for loess > > # > > # Note:? rc is required for methods 1 and 2 but not 3 > > #??????? method 3 return the sample size need for rc=0 to (1-d)*nc > > # > > ######## Output > > # for methdos 1 & 2: return the sample size needed for the experimental group (1 number) > > #??????????????????? for given rc, nc, d, alpha, and power > > # for method 3:????? return the profile of sample size needed for given nc, d, alpha, and power > > #??????????????????? vector $ne contains the sample size corresponding to rc=0, 1, 2, ... nc*(1-d) > > #??????????????????? vector $Ep contains the expected power corresponding to > > #????????????????????? the true pc = (0, 1, 2, ..., nc*(1-d)) / nc > > # > > #------------------------------------------------------------------ > > sshc<-function(rc, nc=500, d=.5, method=3, alpha=0.05, power=0.8, > +????????????? tol=0.01, tol1=.0001, tol2=.005, cc=c(.1,2), l.span=.5) > + { > + ### for method 1 > + if (method==1) { > + ne1<-ss.rand(rc,nc,d,alpha=.05,power=.8,tol=.01) > + return(ne=ne1) > +??????????????? } > + ### for method 2 > + if (method==2) { > + ne<-nc > + ne1<-nc+50 > + while(abs(ne-ne1)>tol & ne1<100000){ > + ne<-ne1 > + pe<-d+rc/nc > + ne1<-nef(rc,nc,pe*ne,ne,alpha,power) > + ## if(is.na(ne1)) print(paste('rc=',rc,',nc=',nc,',pe=',pe,',ne=',ne)) > + } > + if (ne1>100000) return(NA) > + else return(ne=ne1) > + } > + ### for method 3 > + if (method==3) { > + if (tol1 > tol2/10) tol1<-tol2/10 > + ncstar<-(1-d)*nc > + pc<-(0:ncstar)/nc > + ne<-rep(NA,ncstar + 1) > + for (i in (0:ncstar)) > + { ne[i+1]<-ss.rand(i,nc,d,alpha=.05,power=.8,tol=.01) > + } > + plot(pc,ne,type='l',ylim=c(0,max(ne)*1.5)) > + ans<-c.searchd(nc, d, ne, alpha, power, cc, tol1) > + ### check overall absolute deviance > + old.abs.dev<-sum(abs(ans$Ep-power)) > + ##bad<-0 > + print(round(ans$Ep,4)) > + print(round(ans$ne,2)) > + lines(pc,ans$ne,lty=1,col=8) > + old.ne<-ans$ne > + ##while(max(abs(ans$Ep-power))>tol2 & bad==0){? #### unnecessary ## > + while(max(abs(ans$Ep-power))>tol2){ > + ans<-c.searchd(nc, d, ans$ne, alpha, power, cc, tol1) > + abs.dev<-sum(abs(ans$Ep-power)) > + print(paste(" old.abs.dev=",old.abs.dev)) > + print(paste("???? abs.dev=",abs.dev)) > + ##if (abs.dev > old.abs.dev) { bad<-1} > + old.abs.dev<-abs.dev > + print(round(ans$Ep,4)) > + print(round(ans$ne,2)) > + lines(pc,old.ne,lty=1,col=1) > + lines(pc,ans$ne,lty=1,col=8) > + ### add convex > + ans$ne<-convex(pc,ans$ne)$wy > + ### add loess > + ###old.ne<-ans$ne > + loess.ne<-loess(ans$ne ~ pc, span=l.span) > + lines(pc,loess.ne$fit,lty=1,col=4) > + old.ne<-loess.ne$fit > + ###readline() > + } > + return(ne=ans$ne, Ep=ans$Ep) > +??????????????? } > + } > > > > ## needed for method 1 > > nef2<-function(rc,nc,re,ne,alpha,power){ > + za<-qnorm(1-alpha) > + zb<-qnorm(power) > + xe<-asin(sqrt((re+0.375)/(ne+0.75))) > + xc<-asin(sqrt((rc+0.375)/(nc+0.75))) > + ans<- 1/(4*(xc-xe)^2/(za+zb)^2-1/(nc+0.5)) - 0.5 > + return(ans) > + } > > ## needed for method 2 > > nef<-function(rc,nc,re,ne,alpha,power){ > + za<-qnorm(1-alpha) > + zb<-qnorm(power) > + xe<-asin(sqrt((re+0.375)/(ne+0.75))) > + xc<-asin(sqrt((rc+0.375)/(nc+0.75))) > + ans<-(za*sqrt(1+(ne+0.5)/(nc+0.5))+zb)^2/(2*(xe-xc))^2-0.5 > + return(ans) > + } > > ## needed for method 3 > > c.searchd<-function(nc, d, ne, alpha=0.05, power=0.8, cc=c(0.1,2),tol1=0.0001){ > + #--------------------------- > + # nc???? sample size of control group > + # d????? the differece to detect between control and experiment > + # ne???? vector of starting sample size of experiment group > + #??? corresonding to rc of 0 to nc*(1-d) > + # alpha? size of test > + # power? target power > + # cc? pre-screen vector of constant c, the range should cover the > + #??? the value of cc that has expected power > + # tol1?? the allowance between the expceted power and target power > + #--------------------------- > + pc<-(0:((1-d)*nc))/nc > + ncl<-length(pc) > + ne.old<-ne > + ne.old1<-ne.old > + ### sweeping forward > + for(i in 1:ncl){ > + cmin<-cc[1] > + cmax<-cc[2] > + ### fixed cci<-cmax bug > + cci <-1 > + lhood<-dbinom((i:ncl)-1,nc,pc[i]) > + ne[i:ncl]<-(1+(cci-1)*(lhood/lhood[1])) * ne.old1[i:ncl] > + Ep0 <-Epower(nc, d, ne, pc, alpha) > + while(abs(Ep0[i]-power)>tol1){ > + if(Ep0[i]<power) cmin<-cci > + else cmax<-cci > + cci<-(cmax+cmin)/2 > + ne[i:ncl]<-(1+(cci-1)*(lhood/lhood[1])) * ne.old1[i:ncl] > + Ep0<-Epower(nc, d, ne, pc, alpha) > + } > +? ne.old1<-ne > + } > + ne1<-ne > + ### sweeping backward -- ncl:i > + ne.old2<-ne.old > + ne???? <-ne.old > + for(i in ncl:1){ > + cmin<-cc[1] > + cmax<-cc[2] > + ### fixed cci<-cmax bug > + cci <-1 > + lhood<-dbinom((ncl:i)-1,nc,pc[i]) > + lenl <-length(lhood) > + ne[ncl:i]<-(1+(cci-1)*(lhood/lhood[lenl]))*ne.old2[ncl:i] > + Ep0 <-Epower(nc, d, cci*ne, pc, alpha) > + while(abs(Ep0[i]-power)>tol1){ > + if(Ep0[i]<power) cmin<-cci > + else cmax<-cci > + cci<-(cmax+cmin)/2 > + ne[ncl:i]<-(1+(cci-1)*(lhood/lhood[lenl]))*ne.old2[ncl:i] > + Ep0<-Epower(nc, d, ne, pc, alpha) > + } > +? ne.old2<-ne > + } > + ne2<-ne > + ne<-(ne1+ne2)/2 > + #cat(ccc*ne) > + Ep1<-Epower(nc, d, ne, pc, alpha) > + return(ne=ne, Ep=Ep1) > + } > > ### > > vertex<-function(x,y) > + { n<-length(x) > + vx<-x[1] > + vy<-y[1] > + vp<-1 > + up<-T > + for (i in (2:n)) > + { if (up) > + { if (y[i-1] > y[i]) > + {vx<-c(vx,x[i-1]) > +? vy<-c(vy,y[i-1]) > +? vp<-c(vp,i-1) > +? up<-F > + } > + } > + else > + { if (y[i-1] < y[i]) up<-T > + } > + } > + vx<-c(vx,x[n]) > + vy<-c(vy,y[n]) > + vp<-c(vp,n) > + return(vx=vx,vy=vy,vp=vp) > + } > > ### > > convex<-function(x,y) > + { > + n<-length(x) > + ans<-vertex(x,y) > + len<-length(ans$vx) > + while (len>3) > + { > + #cat("x=",x,"\n") > + #cat("y=",y,"\n") > + newx<-x[1:(ans$vp[2]-1)] > + newy<-y[1:(ans$vp[2]-1)] > + for (i in (2:(len-1))) > + { > +? newx<-c(newx,x[ans$vp[i]]) > + newy<-c(newy,y[ans$vp[i]]) > + } > + newx<-c(newx,x[(ans$vp[len-1]+1):n]) > + newy<-c(newy,y[(ans$vp[len-1]+1):n]) > + y<-approx(newx,newy,xout=x)$y > + #cat("new y=",y,"\n") > + ans<-vertex(x,y) > + len<-length(ans$vx) > + #cat("vx=",ans$vx,"\n") > + #cat("vy=",ans$vy,"\n") > + } > + return(wx=x,wy=y)} > > ### > > Epower<-function(nc, d, ne, pc = (0:((1 - d) * nc))/nc, alpha = 0.05) > + { > + #------------------------------------- > + # nc???? sample size in historical control > + # d????? the increase of response rate between historical and experiment > + # ne???? sample size of corresonding rc of 0 to nc*(1-d) > + # pc???? the response rate of control group, where we compute the > + #??????? expected power > + # alpha? the size of test > + #------------------------------------- > + kk <- length(pc) > + rc <- 0:(nc * (1 - d)) > + pp <- rep(NA, kk) > + ppp <- rep(NA, kk) > + for(i in 1:(kk)) { > + pe <- pc[i] + d > + lhood <- dbinom(rc, nc, pc[i]) > + pp <- power1.f(rc, nc, ne, pe, alpha) > + ppp[i] <- sum(pp * lhood)/sum(lhood) > + } > + return(ppp) > + } > > > > # adapted from the old biss2 > > ss.rand<-function(rc,nc,d,alpha=.05,power=.8,tol=.01) > + { > + ne<-nc > + ne1<-nc+50 > + while(abs(ne-ne1)>tol & ne1<100000){ > + ne<-ne1 > + pe<-d+rc/nc > + ne1<-nef2(rc,nc,pe*ne,ne,alpha,power) > + > + ## if(is.na(ne1)) print(paste('rc=',rc,',nc=',nc,',pe=',pe,',ne=',ne)) > + } > + if (ne1>100000) return(NA) > + else return(ne1) > + } > > ### > > power1.f<-function(rc,nc,ne,pie,alpha=0.05){ > + #------------------------------------- > + # rcnumber of response in historical control > + # ncsample size in historical control > + # ne??? sample size in experitment group > + # pietrue response rate for experiment group > + # alphasize of the test > + #------------------------------------- > + > + za<-qnorm(1-alpha) > + re<-ne*pie > + xe<-asin(sqrt((re+0.375)/(ne+0.75))) > + xc<-asin(sqrt((rc+0.375)/(nc+0.75))) > + ans<-za*sqrt(1+(ne+0.5)/(nc+0.5))-(xe-xc)/sqrt(1/(4*(ne+0.5))) > + return(1-pnorm(ans)) > + } > > [[alternative HTML version deleted]]