[R] sequential t-test - replies

[KINLEY_ROBERT@lilly.com]

[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