[R] Confidence intervals for prediction based on the logistic equation

From: Bock, Michael <MBock_at_arcadis-us.com>
Date: Fri 27 May 2005 - 02:01:17 EST


Greetings,

We are performing a meta-analysis of mink pup survival data versus chemical concentration. We have modeled percent survival successfully using nls as shown below and the plot. What we need to do is construct a confidence interval on the concentration at which we get 50% survival (aka the EC50, although we may want other percent survivals in the future). My first question is, what seems like the best method for determining such a confidence interval? Is the model we have chosen being estimated in the best way? We could use the number of survival and the sample size and use the binomial glm model but we only have access to summaries, as such any estimate would require the use of an estimate of sample size which seems inappropriate to me.

In the literature a boot strap approach has been suggested. I have tried this approach but I get model convergence errors. Examples of the errors can be found in the comments of the code below.

Also, I have tried the predict and profile procedure to help determine the EC50 but have been unsuccessful. What I am hoping is either someone can suggest a better method to get what we want or some way to get the bootstraping to work that is statistically valid (it has been suggested that we just ignore bootstrap runs in which convergence is not attained and I have an issue with that idea).

Thanks in advance,
Mike
FYI R 2.1.0 on Win XP SP2

library(boot)
library(MASS)
#Dataset

BC <-NULL
Mink <- structure(list(LogWBCTEQa = c(-0.6477, -1.0592, -0.7454, -0.629,

-3.8486, -2.4847, -1.5025, -0.0256, -1.9292, -1.6282, -1.3272,
-1.0264, -2.1327, 0.165, -0.2807, -0.8197, -1.2969, -1.4345,
-1.664, -0.6525, -0.8067, -1.961, -1.4804, -1.4764, -1.5813,
-1.453, -1.2584, -2.102, -1.8009, -1.4999, -2.1008, -1.7998,
-1.4988, -2.1013, -1.791, -1.4681, -3.3718, -2.8219, -2.6779,
-2.5255, -2.3224, -1.6907), PCTSurv = c(0, 0, 0, 0, 100, 100,
0, 0, 0, 0, 0, 0, 99, 0, 0, 0, 0, 0, 100, 0, 0, 56.78, 0, 73, 25, 12, 0, 100, 99, 37, 92, 4, 7, 71, 4, 0, 100, 92, 70, 100,
100, 52), PCTFail = c(100, 100, 100, 100, 0, 0, 100, 100, 100,
100, 100, 100, 1, 100, 100, 100, 100, 100, 0, 100, 100, 43.22,
100, 27, 75, 88, 100, 0, 1, 63, 8, 96, 93, 29, 96, 100, 0, 8,
30, 0, 0, 48), Psurv = c(0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0,
0.99, 0, 0, 0, 0, 0, 1, 0, 0, 0.567764706, 0, 0.73, 0.25, 0.12,
0, 1, 0.99, 0.37, 0.92, 0.04, 0.07, 0.71, 0.04, 0, 1, 0.92, 0.7,
1, 1, 0.52)), .Names = c("LogWBCTEQa", "PCTSurv", "PCTFail",
"Psurv"), class = "data.frame", row.names = c("1", "2", "3",
"4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15",
"16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26",
"27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "40", "41", "42"))
#function to bootstrap the model

EC50C <- function (xdata,i) {
nlsdata <- xdata[i,]
Sur.nls <- nls( PCTSurv ~ 100/(100 + exp(( xmid - (LogWBCTEQa) )/scal ) ),
                       data = nlsdata,
                       start = list( xmid = 0, scal = 1 ),
                       alg = "plinear", trace = TRUE )

#output xmid, modify to calc EC50 (if PCTSurv = 50 what is LogWBCTEQa)
C <- coef(Sur.nls)[1]
EC50C <- C
}
#bootstrap

BC<- boot(Mink, EC50C, R=10, stype = "i") boot.ci(boot.out = BC, conf = c(0.10, 0.5 , 0.9, 0.95), type = c("norm"))
#Try again without using bootstrap for testing
#sometimes BC fails due to:
#Error in nls(PCTSurv ~ 100/(100 + exp((xmid - (LogWBCTEQa))/scal)),
data = nlsdata, :
# singular gradient
#or
#Error in qr.solve(QR.B, cc) : singular matrix 'a' in solve
#or
#Error in nls(PCTSurv ~ 100/(100 + exp((xmid - (LogWBCTEQa))/scal)),
data = nlsdata, :
# number of iterations exceeded maximum of 50
#or
#Error in nls(PCTSurv ~ 100/(100 + exp((xmid - (LogWBCTEQa))/scal)),
data = nlsdata, :
# step factor 0.000488281 reduced below 'minFactor' of
0.000976563
TS <- nls( PCTSurv ~ 100/(100 + exp(( xmid - (LogWBCTEQa) )/scal ) ),
                       data = Mink,
                       start = list( xmid = 0, scal = 1 ),
                       alg = "plinear", trace = TRUE )
plot(PCTSurv ~ LogWBCTEQa, Mink)
tt <- seq(-4,0,length = 100)
lines(tt,predict(TS,list (LogWBCTEQa = tt)))
profile(TS)
confint(TS,"xmid")
predict(TS,


Michael J. Bock, PhD.
ARCADIS
24 Preble St. Suite 100
Portland, ME 04101
207.828.0046
fax 207.828.0062

        [[alternative HTML version deleted]]



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Fri May 27 02:19:55 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:32:07 EST