[R] SPlus to R
Spencer Graves
spencer.graves at structuremonitoring.com
Wed Oct 5 09:00:06 CEST 2011
When R was invented, nearly all of the core R functions were
written to produce exactly the same answers as those returned by S-Plus.
Some very minor exceptions were made for time series functions, for
example, where better algorithms in R produced slightly better fits.
There may be other, less benign examples, but none come to mind at the
moment. Some of the core R functions have arguments that were not
available in S-Plus at the time they were written. For example,
probability functions included log and lower.tail arguments that, if my
memory is correct, were not available in S-Plus at that time. More
recent versions of S-Plus can reportedly run R scripts and packages
without change.
One set of examples comes to mind where you can find subtle
differences between S-Plus and R: If you get Pinheiro and Bates (2000)
Mixed-Effects Models in S and S-Plus (Springer) and try to run the code
exactly as written in the book, you may get the answers in the book
under S-Plus but not in a few cases in R. However, the nlme package
companion to the book includes a "scripts" directory containing files
required to obtain all the results in that book. In most cases, the
functions and syntax, etc., are exactly the same in S-Plus and R.
However, there are a very few cases where the answers are NOT the same,
because some defaults have changed. I don't remember the details now.
If you have an S-Plus script in which you replaced "_" by "<-",
that should take care of most if not all of the required conversion.
However, you need to check. If it breaks when you try to run it, find
where it breaks, read the documentation and try to understand why it
breaks and how to fix it.
The "dubug" function in R can help immensely with this: If you
say debug(fn), then each time fn is called, it will stop and allow you
to walk through fn line by line evaluating what it does, even changing
the code and values of variables on the fly if you like. This can help
you fix anything that breaks AND understand whether you are still
getting the correct answer.
Hope this helps.
Spencer
On 10/4/2011 11:08 PM, Joshua Wiley wrote:
> 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.
>>
--
Spencer Graves, PE, PhD
President and Chief Technology Officer
Structure Inspection and Monitoring, Inc.
751 Emerson Ct.
San José, CA 95126
ph: 408-655-4567
web: www.structuremonitoring.com
More information about the R-help
mailing list