R help - May 2005 - Confidence intervals for prediction based on the logistic equation

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]]