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