LD50 and SE in GLMM (lmer)

And one more question: is it correct to use pnorm (where John Maindonald used exp(hat)/(1+exp(hat)))?

Hi All!



I am desperately needing some help figuring out how to calculate LD50 with a GLMM (probit link) or, more importantly, the standard error of the LD50.



I conducted a cold temperature experiment and am trying to assess after how long 50% of the insects had died (I had 3 different instars (non significant fixed effect) and several different blocks (I did 4 replicates at a time)= random effect).



Since there is no predict function for lmer, I used the following to get predicted values (thanks to a post by John Maindonald (I'll attach his post below)):




model4 <- lmer (y~time + (1|blc/instar),family=binomial(link="probit"))
summary(model4)



?b <- fixef(model4)
?X <- (model.matrix(terms(model4),zerotest))
?hat <- X%*%b
?pxal <- pnorm(hat) ? ?# probit link, for logit it would be: pval <- exp(hat)/(1+exp(hat))
?pval


Once I get the pval, I see where the 0.5 predicted value lies and I adjust the x's in zerotest to be more detailed in that range, eg. x: 1-420hours, I see that 0.5 is in the 320hours area, so I adjust x to be 320.1, 320.2, 320.3, etc. to get the precise 0.500. Very clumsy but I guess it's correct?



Now my biggest problem: how do I get the SE?



John Maindonald goes on to do this:

U <- chol(as.matrix(summary(model4)@vcov))

se <- sqrt(apply(X%*%t(U), 1, function(x)sum(x^2)))

list(hat=hat, se=se, x=X[,xcol])



Unfortunately, I could not figure out what the chol(as.matrix...) part is about (chol does what?) and therefore I have no idea, how to use this code to get my LD50 SE (I would need the SE to be expressed in terms of x).

Could anybody help me with this?



And one more question: is it correct to use pnorm (where John Maindonald used exp(hat)/(1+exp(hat)))?



Thanks so much in advance!



Linda





Previous post by John Maindonald:




ciplot <- ?function(obj=model4, data=zerotest,xcol=2,nam="hours"){
? ?cilim <- function(obj, xcol){
? ? ?b <- fixef(obj)
? ? ?vcov <- summary(obj)@vcov
? ? ?X <- unique(model.matrix(obj))
? ? ?hat <- X%*%b
?pval <- exp(hat)/(1+exp(hat)) ? # NB, designed for logit link
? ? U <- chol(as.matrix(summary(obj)@vcov))
? ? ?se <- sqrt(apply(X%*%t(U), 1, function(x)sum(x^2)))
? ? ?list(hat=pval, se=se, x=X[,xcol])
? ? ? ? ? }
? ?limfo <- cilim(obj, xcol)
? ?hat <- limfo$hat
? ?se <- limfo$se
? ?x <- limfo$x
? ?upper <- hat+se
? ?lower <- hat-se
? ?ord <- order(x)
? ?plot(x, hat, yaxt="n", type="l", xlab=nam, ylab="")
? ?rug(x)
? ?lines(x[ord], lower[ord])
? ?lines(x[ord], upper[ord])
? ?ploc <- c(0.01, 0.05, 0.1, 0.2, 0.5, 0.8, 0.9)
? ?axis(2, at=log(ploc/(1-ploc)), labels=paste(ploc), las=2)
? ?}

## Usage
model4 <- lmer (y~time + (1|blc/instar),family=binomial)

ciplot(obj=model4)





Linda Buergi
Environmental Science, Policy and Management
UC Berkeley, Berkeley, California