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.