I'm trying to convert an S-Plus program to R.  Since I'm a SAS programmer I'm not facile is either S-Plus or R, so I need some help.  All I did was convert the underscores in S-Plus to the assignment operator <-.  Here are the first few lines of the S-Plus file:
 
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) 
               }

 
My translation looks like this:
 
sshc<-function(rc, nc=500, d=.5, 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) 
               }

The program runs without throwing errors, but I'm not getting any ourput in the console.  This is where it should be, right?  I think I have this set up correctly.  I'm using method=3 which only requires nc and d to be specified.  Any ideas why I'm not seeing output?
 
Here is the entire output:
 
> ## 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=500, d=.5, 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(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){
+ #-------------------------------------
+ # rcnumber of response in historical control
+ # ncsample size in historical control
+ # ne    sample size in experitment group
+ # pietrue response rate for experiment group
+ # alphasize 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))
+ }

	[[alternative HTML version deleted]]

