[R] LD50 and SE in GLMM (lmer)

[Bill Pikounis]

Sorry for the delay in response. I had a somewhat similar need
recently with the difference that I used a logit link for a bioassay.
The design had different dose-response "replicates" that I modeled as
"blocks". It looks like you are concentrating on estimation of fixed
effects and thus the population / marginal LD50 estimate. If so, then
there is a function called dose.p in the MASS package, courtesy of
Venables and Ripley, which is used in the context of an example on 190
- 194 of the 4th edition of their book (2002), 4th ediiion, that I
think would be very helpful to study. The example code can also be
found in the ch07.R file in the scripts sub-directory/folder of the
MASS package directory/folder. The example illustrates the use of GLM
with a logit link. To adapt it for use with a GLMM, I came up with the
following, which is nearly identical to how dose.p is defined in R
2.10.0

dose.p.glmm <-  function(obj, cf = 1:2, p = 0.5) {
  eta <- obj$family$linkfun(p)
  b <- fixef(obj)[cf]
  x.p <- (eta - b[1L])/b[2L]
  names(x.p) <- paste("p = ", format(p), ":", sep = "")
  pd <- -cbind(1, x.p)/b[2L]
    SE <- sqrt(((pd %*% vcov(obj)[cf, cf]) * pd) %*% c(1, 1))
  res <- structure(x.p, SE = SE, p = p)
  class(res) <- "glm.dose"
  res
}

Essentially only the fixef() call in the 2nd line of the body was
needed to replace the coef() call. Please also note that I used this
for a glmmPQL() call from the MASS package, not lmer().

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

Unfortunately I don't know offhand, and do not have a reference handy
to check to be sure, so perhaps you can find a local statistician to
help?  I myself  always have a preference to use the logit / logistic
over probit, as they are both symmetric around 0.5 and are often
reported to provide similar results.

Hope that helps,
Bill

###############

Bill Pikounis
Statistician



2010/1/7 Linda Bürgi <patili_buergi at hotmail.com>:
>
> 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
>
>
>
> _________________________________________________________________
> Hotmail: Trusted email with powerful SPAM protection.
>
>        [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>