[my original message to s-news & r-help is attached ]
No one possessed or knew of any S/R code for the sequential t-test. Also
it doesn't appear in the SAS index.
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".
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).
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.
cheers Bob
Code follows ...
function(delta = 1, alpha = 0.05, beta = 0.05, minn = 1, maxn = 20, bot =
-10, top = 10)
{
# Create tables and graphs for Barnard's SPRT ("sequential
t-test")
# (a) uses expressions given in Wetherill & Glazebrook
'Sequential
Methods in Statistics', p 60
# (b) succesfully reproduces table L in OL Davies Design &
Analysis of Industrial Experiments'
#
# notation is a mixture of a and b.
#
# delta is the difference-to-detect ( in std devs )
# alpha, beta are the reqd type1 & 2 errors
# minn:maxn is the range of sample-numbers
# bot:top is the range used in solving for critical values
#
#
Hh <- function(N = 5, X = 0, lolim = 0, hilim = Inf)
{
# an integral required in the likelehood ratio
#
integrand <- function(y, n, x)
{
ans <- ((y^n) * exp(-0.5 * (x + y)^2))/gamma(n +
1)
ans
}
a <- integrate(f = integrand, lower = lolim, upper =
hilim, subdivisions = 500, n = N, x = X)
a$integral
}
#
#
Lik <- function(U, N, DELTA, lhs)
{
# roots of this function are the critical values ( 'U0'
and 'U1' in Davies )
#
lik <- lhs - (exp(-0.5 * (DELTA^2) * (N - U^2)) * Hh(N -
1, - DELTA * U))/Hh(N - 1)
lik
}
#
#
#
U0 <- U1 <- rep(NA, maxn)
#
SampNo <- 1:maxn
#
# for each n in the required sequence, get the critical values U0
and U1
#
for(n in minn:maxn) {
U1[n] <- uniroot(Lik, lower = bot, upper = top, N = n,
DELTA = delta, lhs = (1 - beta)/alpha)$root
U0[n] <- uniroot(Lik, lower = bot, upper = top, N = n,
DELTA = delta, lhs = beta/(1 - alpha))$root
cat("\n", n, "\t", U0[n], "\t",
U1[n])
}
#
# make a plot to provide hard-copy for users
#
ttl <- paste("delta=", delta, "alpha=", alpha,
"beta=", beta)
yr <- range(na.omit(c(U0, U1)))
# dyr <- max(yr) - min(yr)
# yr <- c(yr[1] - 0.2 * dyr, yr[2] + 0.2 * dyr)
pyr <- pretty(yr)
plot(c(1, maxn), yr, type = "n", xlab = "Sample
no.", ylab = "",
main = ttl)
lines(1:maxn, U0, lwd = 5, col = 13)
lines(1:maxn, U1, lwd = 5, col = 13)
text(1, U1[1], "Difference detected")
text(1, U0[1], "Difference not detected")
abline(v = seq(1:maxn), h = pyr)
#
# table as in Davies
#
ans <- data.frame(cbind(SampNo, U0, U1))
ans
}
................................
Robert D. Kinley
Site Statistician
Eli Lilly & co
Speke, UK
++151 448 6347
KINLEY_ROBERT at lilly.com
Sent by: s-news-owner at lists.biostat.wustl.edu
08/03/2002 11:21
To: "s-news (E-mail)" <s-news at
lists.biostat.wustl.edu>
cc:
Subject: [S] any S/R Code for Barnard's sequential t-test ?
Does anyone have, or know of, some S or R code to implement
the sequential t-test ?
The test is outlined in Kendall & Stuart vol 2A, ch 24, and a 'how
to do it' with a set of tables is given in 'Design & Analysis of
Industrial Experiments' ( O.L.Davies).
I'm interested in obtaining/creating an Splus or R application which
calculates the upper and lower curves appropriate for arbitrary
choices of alpha , Beta and mean-shift.
Can anyone give me any pointers ?
ta Bob
...............................
Robert D. Kinley
Site Statistician
Eli Lilly & co
Speke, UK
++151 448 6347
-------------- next part --------------
An HTML attachment was scrubbed...
URL:
https://stat.ethz.ch/pipermail/r-help/attachments/20020322/d8d24882/attachment.html
