[R] SPlus to R
Joshua Wiley
jwiley.psych at gmail.com
Wed Oct 5 08:08:56 CEST 2011
Hi Scott,
I am not familiar with S-Plus (though many aspects are quite similar
to R). I will say that your function looks approximately correct. I
am not familiar with the ss.rand function. I searched, and found some
things that I suspect are similar in the packages MBESS, but without
knowing more about it from S-Plus, it is tough to make a testable
example.
Do you have access to S-Plus? Can you provide more information about
this function, what it does, what is like, etc.? There are some
active members of this list who are quite familiar with S-Plus so one
of them may be more insightful.
Cheers,
Josh
On Tue, Oct 4, 2011 at 6:53 PM, Scott Raynaud <scott.raynaud at yahoo.com> wrote:
> 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 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=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]]
>
>
> ______________________________________________
> 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.
>
>
--
Joshua Wiley
Ph.D. Student, Health Psychology
Programmer Analyst II, ATS Statistical Consulting Group
University of California, Los Angeles
https://joshuawiley.com/
More information about the R-help
mailing list