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))
}
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.raynaud@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 (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=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@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]]
See my comments below. ________________________________ From: Pascal Oettli <kridox at ymail.com> To: Scott Raynaud <scott.raynaud at yahoo.com> Sent: Wednesday, June 5, 2013 10:02 AM Subject: Re: [R] SPlus script Hello, 1) It is always nice to say something as "Hello", Tried, but could make it come out in anything other than four letters. 2) What do you want us to do with that script, without the required "commented, minimal, self-contained, reproducible code"? I want it to run as it did before.? I cannot explain why the exact same code used to run and now doesn't.? I've tried everything I can think of, so maybe brighter minds than mine can shed some light. 3) The lastest version of R is 3.0.1. Yes, I know.? I was trying to keep things as constant as possible.? Besides, do you really think the version wouold make a difference?? The code runs without error.? It just doesn't produce output.? Regards, Pascal 2013/6/5 Scott Raynaud <scott.raynaud at 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 (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)) >} > >______________________________________________ >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. >
Both the R and S+ versions (which seem to differ only in the use of _ for assignment in the S+ version) do nothing but define some functions. You would not expect any printed output unless you used those functions on some data. Is there another script that does that? 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: Wednesday, June 05, 2013 6:21 AM > To: r-help at r-project.org > Subject: [R] SPlus script > > 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)) > } > > ______________________________________________ > 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.