Bock, Michael
2005-May-26 16:01 UTC
[R] 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]]