I'm re-posting this since it did not appear at the end of the thread. Sorry
for the inconvenience. Not sure why it's giving the message: An embedded and
charset-unspecified text was scrubbed... As far as I know my replies are set
up as text and not html.
I did some testing and I believe the program is operating properly. It takes
some time to finish, especially as sample sizes get larger, but I seem to be
able to reproduce the results from the original paper. Right now I'm most
interested in method 3. I set nc=40 and d=.2 as in the paper. The results are
below. The program iterates to find a solution and the more subjects (nc), the
more iterations it needs. The initial pass uses estimates from the randomized
method as starting values and this is plotted in black on a graph. The next
pass uses the algorithm to generate sample size estimates which are plotted in
gray on the first pass. Each new estimate is plotted in dark blue as the
previous one turns black. This reults in several curves being plotted but one
can get an idea of the relevant one by looking for the blue plot and examining
the $ne vector below.
Reading the graph is more troublesome as nc increases. It takes many more
iterations to complete and the graph gets quite crowded. I'm not sure what the
best solution is. Which would be easiest?
1) Plot ne against pc for the final iteration. The problem here is pc is not
saved in a vector.
2) Cut and paste $ne into Excel and re-create the corresponding pc's.
It's not clear to me how pc is incremented. Seems like it's in the code block
below:
### 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)
> sshc()
[1] 0.8000 0.7889 0.7825 0.7786 0.7763 0.7751 0.7747 0.7750 0.7759 0.7772
[11] 0.7790 0.7815 0.7852 0.7899 0.7943 0.7968 0.7979 0.7989 0.8002 0.8021
[21] 0.8045 0.8072 0.8100 0.8123 0.8137 0.8137 0.8124 0.8105 0.8095 0.8103
[31] 0.8134 0.8180 0.8244
[1] 11.64 19.62 27.50 36.44 46.54 58.19 71.81 87.93 107.30 130.38
[11] 158.47 192.34 235.47 291.93 370.48 471.26 569.17 626.36 626.18 588.94
[21] 529.27 459.67 389.30 323.70 263.50 210.58 164.26 124.25 90.49 63.23
[31] 42.13 26.23 8.13
[1] " old.abs.dev= 0.434571909427754"
[1] " abs.dev= 0.259398042835904"
[1] 0.8000 0.8031 0.8063 0.8094 0.8120 0.8140 0.8152 0.8156 0.8152 0.8141
[11] 0.8122 0.8097 0.8073 0.8055 0.8045 0.8037 0.8025 0.8015 0.8006 0.7995
[21] 0.7982 0.7969 0.7956 0.7946 0.7939 0.7930 0.7916 0.7897 0.7879 0.7872
[31] 0.7876 0.7886 0.7887
[1] 11.64 20.91 30.51 42.13 56.19 73.71 95.59 123.29 159.20 204.83
[11] 263.66 337.90 424.41 520.76 623.55 727.31 804.98 810.60 733.27 621.58
[21] 504.83 397.60 310.92 242.78 190.93 151.72 120.81 95.24 72.67 52.84
[31] 36.01 22.42 6.32
[1] " old.abs.dev= 0.259398042835904"
[1] " abs.dev= 0.160699574601494"
[1] 0.8000 0.7985 0.7972 0.7958 0.7944 0.7930 0.7917 0.7907 0.7903 0.7907
[11] 0.7917 0.7929 0.7941 0.7951 0.7964 0.7977 0.7982 0.7984 0.7989 0.7995
[21] 0.8000 0.8006 0.8013 0.8023 0.8036 0.8048 0.8054 0.8055 0.8055 0.8063
[31] 0.8081 0.8103 0.8125
[1] 11.64 20.56 29.46 39.64 51.27 64.86 80.67 99.46 122.30 149.90
[11] 184.62 230.35 286.04 352.93 429.63 508.43 588.14 635.70 616.86 561.50
[21] 489.16 408.17 333.63 269.21 215.86 173.66 139.24 110.43 84.50 61.28
[31] 41.38 25.49 7.42
[1] " old.abs.dev= 0.160699574601494"
[1] " abs.dev= 0.0893104146067729"
[1] 0.8000 0.8006 0.8013 0.8021 0.8029 0.8036 0.8043 0.8048 0.8050 0.8051
[11] 0.8049 0.8045 0.8037 0.8024 0.8019 0.8019 0.8015 0.8012 0.8010 0.8009
[21] 0.8007 0.8002 0.7995 0.7990 0.7987 0.7985 0.7980 0.7970 0.7958 0.7950
[31] 0.7950 0.7950 0.7939
[1] 11.64 20.72 29.92 40.75 53.42 68.99 87.79 111.53 141.92 179.12
[11] 226.54 290.16 365.88 456.60 556.96 640.04 722.83 757.33 711.48 624.30
[21] 522.26 423.34 335.91 264.10 206.62 162.86 128.85 101.60 77.86 56.72
[31] 38.30 23.41 6.49
[1] " old.abs.dev= 0.0893104146067729"
[1] " abs.dev= 0.0506515395975703"
[1] 0.8000 0.7996 0.7994 0.7990 0.7986 0.7982 0.7978 0.7974 0.7970 0.7969
[11] 0.7970 0.7975 0.7981 0.7979 0.7982 0.7990 0.7991 0.7990 0.7990 0.7992
[21] 0.7994 0.7996 0.7996 0.7998 0.8002 0.8009 0.8015 0.8015 0.8014 0.8016
[31] 0.8025 0.8035 0.8041
[1] 11.64 20.65 29.71 40.22 52.31 66.93 84.02 104.94 131.25 161.98
[11] 200.06 251.89 310.71 386.59 476.45 545.33 622.06 664.45 636.92 570.09
[21] 486.12 402.82 326.45 262.37 208.03 165.78 131.88 104.87 81.05 59.38
[31] 40.18 24.54 6.99
$ne
[1] 11.637555 20.651336 29.710258 40.223468 52.314568 66.928887
[7] 84.022673 104.936774 131.254491 161.981416 200.060785 251.890694
[13] 310.710449 386.591379 476.449601 545.327456 622.056417 664.447078
[19] 636.923065 570.091413 486.117837 402.817857 326.451672 262.370265
[25] 208.029592 165.777369 131.875077 104.873469 81.051410 59.379349
[31] 40.182140 24.544626 6.990464
$Ep
[1] 0.7999941 0.7996387 0.7993605 0.7990384 0.7986411 0.7982104 0.7978227
[8] 0.7974117 0.7970255 0.7969174 0.7970263 0.7974569 0.7980592 0.7979356
[15] 0.7982330 0.7989812 0.7991050 0.7990024 0.7990236 0.7992147 0.7994437
[22] 0.7995969 0.7996266 0.7997703 0.8002355 0.8009229 0.8014502 0.8015201
[29] 0.8013688 0.8016032 0.8025078 0.8035148 0.8040642
________________________________
From: Bert Gunter <gunter.ber...@gene.com>
Cc: r-help@r-project.org
Sent: Wednesday, October 5, 2011 2:25 PM
Subject: Re: [R] SPlus to R
Consider:
> f <- function(x){ x<- 10;x^2}
> f()
[1] 100
If the argument is not needed, there is no error in omitting it.
R uses "lazy evaluation" -- arguments are not evaluated until needed.
-- Bert
te:
It seems I have things set up correctly. I suspect that the arguments
sshc(100,10) are the isuue. It seems that the 100,10 is not necessary since
the code itself specifies the arguments. It runs and produces a power curve if
I simply type sshc() but it also seems to try to keep running somethng as I
have to click stop to get back to a prompt in the console.
>
>Why specify 100,10? There are 9 arguments, 3 which are required and the rest
>optional. Shouldn't I have to specify the 3 required arguments, nc, d and
>method at a minimum? It would look like sshc(nc=500, d=.5, method=3), right?
>I;m still not sure, however, why that would be necessary since it's hard coded.
>
>
>________________________________
>From: Barry Rowlingson <b.rowling...@lancaster.ac.uk>
>
>Cc: "r-help@r-project.org" <r-help@r-project.org>
>Sent: Wednesday, October 5, 2011 9:27 AM
>Subject: Re: [R] SPlus to R
>
>
>ote:
>> Hope I did this right. I repeated what I'd done before:
>>
>> 1) Opened script
>> 2) Selected run all (this produced my inital post
>>
>> Then as suggested I:
>>
>> 3) Typed ls()
>> 4) Saw that the function was present and issued sshc(100,10)
>>
>> Here's what I got:
>>
>>> ls()
>> [1] "c.searchd" "convex" "Epower" "nef" "nef2" "power1.f"
>> [7] "ss.rand" "sshc" "vertex"
>>> sshc(100,10)
>> Error in return(ne = ne, Ep = Ep1) :
>> multi-argument returns are not permitted
>> So it looks like I need to change the return(ne = ne, Ep = Ep1) to two
>> separate lines, correct?
>>
>> On a brighter note, I did get a power curve as expected. One thing I don't
>> understand is the meaning of the arguments in sshc(100,10).
>
>There are some comments in the function code that tell you:
>
># 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:
>
>- and so on.
>
>Beyond that, I'll have to defer to people who know what this is
>actually trying to compute...
>
>Also, its highly possible that this code has already been ported to R
>- lots of things have. If you know what its meant to compute then a
>quick search might get you running quicker.
>
>Barry
> [[alternative HTML version deleted]]
>
>
>______________________________________________
>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.
>
>
--
"Men by nature long to get on to the ultimate truths, and will often be
impatient with elementary studies or fight shy of them. If it were possible to
reach the ultimate truths without the elementary studies usually prefixed to
them, these would not be preparatory studies but superfluous diversions."
-- Maimonides (1135-1204)
Bert Gunter
Genentech Nonclinical Biostatistics
________________________________
From: William Dunlap <wdun...@tibco.com>
@r-project.org>
Sent: Wednesday, October 5, 2011 11:25 AM
Subject: RE: [R] SPlus to R
I think you only have to change the multi-argument
returns to call list. You can remove the name from
the single argument return, as it will be ignore
return(name=value) -> return(value)
return(n1=v1, n2=v2) -> return(list(n1=v1, n2=v2))
(I say "I think" because I don't have easy access to
S+ 4.5, from 1999 or so, to check this out.)
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 William Dunlap
> Sent: Wednesday, October 05, 2011 9:09 AM
> To: Scott Raynaud; r-help@r-project.org
> Subject: Re: [R] SPlus to R
>
> It looks like this code was written for S+ 4.5 (aka '2000')
> or before, which was based on S version 3. Try changing
> return(name1=value1, name2=value2)
> to
> return(list(name1=value1, name2=value2))
> In S+ from 5.0 onwards return(name=value) or return(name1=value1,
> name2=value2) throws away the names and in R return only takes
> a single object (and also ignores the name).
>
> The c.search function in your code ends with
> return(ne=ne, Ep=Ep1)
> and the code calling c.search() acts as though the writer
> expects that function to return list(ne=ne, Ep=Ep1)
> ans <- c.searchd(nc, d, ne, alpha, power, cc, tol1)
> ...
> old.ne <- ans$ne
>
> 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: Tuesday, October 04, 2011 6:53 PM
> > To: r-help@r-project.org
> > Subject: [R] SPlus to R
> >
> > 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 (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=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@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.
[[alternative HTML version deleted]]
______________________________________________
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.
