As you can see, my original modification was to replace all _ with
<-. It worked after I
did that. This is a simulation that generates its own data based on
the given parameters.
I should obtain a plot of the number of experimental subjects as a
function of the repsonse
rate in the historical group. Maybe I'm forgetting something
here... It's been a while since I
ran this. Thanks for your help.
----- Original Message -----
From: William Dunlap <wdun...@tibco.com>
To: Scott Raynaud <scott.rayn...@yahoo.com>; "r-help@r-project.org"
<r-help@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-boun...@r-project.org
[mailto:r-help-boun...@r-project.org] On Behalf
> Of Scott Raynaud
> Sent: Wednesday, June 05, 2013 6:21 AM
> To: r-help@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 (jj...@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 (jj...@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.
______________________________________________
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.