[R] SPlus script

Scott Raynaud scott.raynaud at yahoo.com
Wed Jun 5 15:20:35 CEST 2013


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 (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=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 (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, 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))
}



More information about the R-help mailing list