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.
