Hello, 1) It is always nice to say something as "Hello", 2) What do you want us to do with that script, without the required "commented, minimal, self-contained, reproducible code"? 3) The lastest version of R is 3.0.1.
Regards, Pascal 2013/6/5 Scott Raynaud <scott.rayn...@yahoo.com> > This originally was an SPlus script that I modifeid about a > year-and-a-half ago. It worked perfectly then. Now I can't get any output > despite not receiving an error message. I'm providing the SPLUS script as > a reference. I'm running R15.2.2. Any help appreciated. > > ************************************MY > MODIFICATION********************************************************************* > ## sshc.ssc: sample size calculation for historical control studies > ## J. Jack Lee (jj...@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=1092, d=.085779816, 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(list(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(list(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(list(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(list(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){ > #------------------------------------- > # rc number of response in historical control > # nc sample size in historical control > # ne sample size in experitment group > # pie true response rate for experiment group > # alpha size 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)) > } > > > > *************************************ORIGINAL SPLUS > SCRIPT************************************************************ > ## sshc.ssc: sample size calculation for historical control studies > ## J. Jack Lee (jj...@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, 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) > } > ### 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){ > #------------------------------------- > # rc number of response in historical control > # nc sample size in historical control > # ne sample size in experitment group > # pie true response rate for experiment group > # alpha size 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)) > } > > ______________________________________________ > R-help@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. > [[alternative HTML version deleted]] ______________________________________________ R-help@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.