[R] SPlus script

Ista Zahn istazahn at gmail.com
Thu Jun 6 16:15:45 CEST 2013


Presumably something like

r <- sshc(50)
print(r)

But if you were getting output before than you already have a script
that does something like this. It would be better to find it...

Best,
Ista

On Thu, Jun 6, 2013 at 9:02 AM, Scott Raynaud <scott.raynaud at yahoo.com> wrote:
> Ok.  Now I see that the sshc function is not being called.  Thanks for pointing that out.
> I'm not certain about the solution, however.  I tried putting call("sshc") at the end of the
> program, but nothing happened.  My memory about all of this is fuzzy.  Suggestions
> on how to call the function appreciated.
>
> ----- Original Message -----
> From: William Dunlap <wdunlap at tibco.com>
> To: Scott Raynaud <scott.raynaud at yahoo.com>; "r-help at r-project.org" <r-help at r-project.org>
> Cc:
> Sent: Wednesday, June 5, 2013 2:17 PM
> Subject: RE: [R] SPlus script
>
> 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.
>
> ______________________________________________
> 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.



More information about the R-help mailing list