<br><font size=2 face="sans-serif">[my original message to s-news & r-help is attached ]</font>
<br>
<br><font size=2 face="sans-serif">No one possessed or knew of any S/R code for the sequential t-test. Also it doesn't appear in the SAS index.</font>
<br><font size=2 face="sans-serif">One or two suggested obtaining the S+ seqtrial software which may (or may not) cover this, but this seemed to be a bit of a "hammer to crack a nut". </font>
<br>
<br><font size=2 face="sans-serif">I have written a function based on the treatment in Wetheril and Glazebrook 'Sequential Methods in Statistics', which seems to reproduce the published tables OK (eg, table L in O.L.Davies Design and Analysis of Industrial Experiments).</font>
<br>
<br><font size=2 face="sans-serif">The function is neither elegant nor fast, but anyone is free to make use of it ... there being no guarantees or support with it, of course.</font>
<br>
<br><font size=2 face="sans-serif"> cheers Bob</font>
<br>
<br><font size=2 face="sans-serif">Code follows ...</font>
<br>
<br><font size=1 face="Courier New">function(delta = 1, alpha = 0.05, beta = 0.05, minn = 1, maxn = 20, bot = -10, top = 10)</font>
<br><font size=1 face="Courier New">{</font>
<br><font size=1 face="Courier New"> # Create tables and graphs for Barnard's SPRT ("sequential t-test") </font>
<br><font size=1 face="Courier New"> # (a) uses expressions given in Wetherill & Glazebrook 'Sequential Methods in Statistics', p 60</font>
<br><font size=1 face="Courier New"> # (b) succesfully reproduces table L in OL Davies Design & Analysis of Industrial Experiments'</font>
<br><font size=1 face="Courier New"> #</font>
<br><font size=1 face="Courier New"> # notation is a mixture of a and b.</font>
<br><font size=1 face="Courier New"> #</font>
<br><font size=1 face="Courier New"> # delta is the difference-to-detect ( in std devs )</font>
<br><font size=1 face="Courier New"> # alpha, beta are the reqd type1 & 2 errors</font>
<br><font size=1 face="Courier New"> # minn:maxn is the range of sample-numbers</font>
<br><font size=1 face="Courier New"> # bot:top is the range used in solving for critical values</font>
<br><font size=1 face="Courier New"> #</font>
<br><font size=1 face="Courier New"> #</font>
<br><font size=1 face="Courier New"> Hh <- function(N = 5, X = 0, lolim = 0, hilim = Inf)</font>
<br><font size=1 face="Courier New"> {</font>
<br><font size=1 face="Courier New"> # an integral required in the likelehood ratio</font>
<br><font size=1 face="Courier New"> #</font>
<br><font size=1 face="Courier New"> integrand <- function(y, n, x)</font>
<br><font size=1 face="Courier New"> {</font>
<br><font size=1 face="Courier New"> ans <- ((y^n) * exp(-0.5 * (x + y)^2))/gamma(n + 1)</font>
<br><font size=1 face="Courier New"> ans</font>
<br><font size=1 face="Courier New"> }</font>
<br><font size=1 face="Courier New"> a <- integrate(f = integrand, lower = lolim, upper = hilim, subdivisions = 500, n = N, x = X)</font>
<br><font size=1 face="Courier New"> a$integral</font>
<br><font size=1 face="Courier New"> }</font>
<br><font size=1 face="Courier New"> #</font>
<br><font size=1 face="Courier New"> #</font>
<br><font size=1 face="Courier New"> Lik <- function(U, N, DELTA, lhs)</font>
<br><font size=1 face="Courier New"> {</font>
<br><font size=1 face="Courier New"> # roots of this function are the critical values ( 'U0' and 'U1' in Davies )</font>
<br><font size=1 face="Courier New"> #</font>
<br><font size=1 face="Courier New"> lik <- lhs - (exp(-0.5 * (DELTA^2) * (N - U^2)) * Hh(N - 1, - DELTA * U))/Hh(N - 1)</font>
<br><font size=1 face="Courier New"> lik</font>
<br><font size=1 face="Courier New"> }</font>
<br><font size=1 face="Courier New"> #</font>
<br><font size=1 face="Courier New"> #</font>
<br><font size=1 face="Courier New"> # </font>
<br><font size=1 face="Courier New"> U0 <- U1 <- rep(NA, maxn)</font>
<br><font size=1 face="Courier New"> #</font>
<br><font size=1 face="Courier New"> SampNo <- 1:maxn</font>
<br><font size=1 face="Courier New"> #</font>
<br><font size=1 face="Courier New"> # for each n in the required sequence, get the critical values U0 and U1</font>
<br><font size=1 face="Courier New"> #</font>
<br><font size=1 face="Courier New"> for(n in minn:maxn) {</font>
<br><font size=1 face="Courier New"> U1[n] <- uniroot(Lik, lower = bot, upper = top, N = n, DELTA = delta, lhs = (1 - beta)/alpha)$root</font>
<br><font size=1 face="Courier New"> U0[n] <- uniroot(Lik, lower = bot, upper = top, N = n, DELTA = delta, lhs = beta/(1 - alpha))$root</font>
<br><font size=1 face="Courier New"> cat("\n", n, "\t", U0[n], "\t", U1[n])</font>
<br><font size=1 face="Courier New"> }</font>
<br><font size=1 face="Courier New"> #</font>
<br><font size=1 face="Courier New"> # make a plot to provide hard-copy for users</font>
<br><font size=1 face="Courier New"> #</font>
<br><font size=1 face="Courier New"> ttl <- paste("delta=", delta, "alpha=", alpha, "beta=", beta)</font>
<br><font size=1 face="Courier New"> yr <- range(na.omit(c(U0, U1)))</font>
<br><font size=1 face="Courier New"> # dyr <- max(yr) - min(yr)</font>
<br><font size=1 face="Courier New"> # yr <- c(yr[1] - 0.2 * dyr, yr[2] + 0.2 * dyr)</font>
<br><font size=1 face="Courier New"> pyr <- pretty(yr)</font>
<br><font size=1 face="Courier New"> plot(c(1, maxn), yr, type = "n", xlab = "Sample no.", ylab = "", main = ttl)</font>
<br><font size=1 face="Courier New"> lines(1:maxn, U0, lwd = 5, col = 13)</font>
<br><font size=1 face="Courier New"> lines(1:maxn, U1, lwd = 5, col = 13)</font>
<br><font size=1 face="Courier New"> text(1, U1[1], "Difference detected")</font>
<br><font size=1 face="Courier New"> text(1, U0[1], "Difference not detected")</font>
<br><font size=1 face="Courier New"> abline(v = seq(1:maxn), h = pyr)</font>
<br><font size=1 face="Courier New"> #</font>
<br><font size=1 face="Courier New"> # table as in Davies</font>
<br><font size=1 face="Courier New"> #</font>
<br><font size=1 face="Courier New"> ans <- data.frame(cbind(SampNo, U0, U1))</font>
<br><font size=1 face="Courier New"> ans</font>
<br><font size=1 face="Courier New">}</font>
<br>
<br><font size=2 face="sans-serif"> <br>
................................<br>
Robert D. Kinley<br>
Site Statistician<br>
Eli Lilly & co<br>
Speke, UK<br>
++151 448 6347<br>
</font>
<br>
<br>
<br>
<table width=100%>
<tr valign=top>
<td>
<td><font size=1 face="sans-serif"><b>KINLEY_ROBERT@lilly.com</b></font>
<br><font size=1 face="sans-serif">Sent by: s-news-owner@lists.biostat.wustl.edu</font>
<p><font size=1 face="sans-serif">08/03/2002 11:21</font>
<br>
<td><font size=1 face="Arial"> </font>
<br><font size=1 face="sans-serif"> To: "s-news (E-mail)" <s-news@lists.biostat.wustl.edu></font>
<br><font size=1 face="sans-serif"> cc: </font>
<br><font size=1 face="sans-serif"> Subject: [S] any S/R Code for Barnard's sequential t-test ?</font>
<br></table>
<br>
<br>
<br><font size=2 face="sans-serif"><br>
Does anyone have, or know of, some S or R code to implement <br>
the sequential t-test ?</font><font size=3> <br>
</font><font size=2 face="sans-serif"><br>
The test is outlined in Kendall & Stuart vol 2A, ch 24, and a 'how <br>
to do it' with a set of tables is given in 'Design & Analysis of <br>
Industrial Experiments' ( O.L.Davies).</font><font size=3> <br>
</font><font size=2 face="sans-serif"><br>
I'm interested in obtaining/creating an Splus or R application which</font><font size=3> </font><font size=2 face="sans-serif"><br>
calculates the upper and lower curves appropriate for arbitrary <br>
choices of alpha , Beta and mean-shift. </font><font size=3> <br>
</font><font size=2 face="sans-serif"><br>
Can anyone give me any pointers ?</font><font size=3> <br>
</font><font size=2 face="sans-serif"><br>
ta Bob</font><font size=3> </font><font size=2 face="sans-serif"><br>
<br>
...............................<br>
Robert D. Kinley<br>
Site Statistician<br>
Eli Lilly & co<br>
Speke, UK<br>
++151 448 6347</font>
<br>
<br>