Thank you very much for the code Bill! I only had to change very few things to
make it work for probit and lmer (instead of glmmPQL) and it works perfectly!
Here's my code (I had some trouble with the output style glm.dose, so I just
have it come out as an ugly list now, which isn't a problem since I only have
to put it into a table).
model4 <- lmer (y~time + (1|blc/instar),REML=FALSE,
family=binomial(link="probit"))
summary(model4)
dose.p.glmm <- function(model4, cf = 1:2, p = 0.5) {
eta <- probit(p)
b <- fixef(model4)[cf]
k <- (eta - b[1])/b[2]
names(k) <- paste("p = ", format(p), ":", sep = "")
pd <- -cbind(1, k)/b[2]
SE <- sqrt(((pd %*% vcov(model4)[cf,cf]) * pd) %*% c(1, 1))
list(k, SE)
}
dose.p.glmm (model4, cf=1:2, p=0.5)
With this, I don't even need the pnorm- probit transformation.
Thank you very much!
Linda
> Date: Mon, 11 Jan 2010 10:21:50 -0500
> Subject: Re: [R] LD50 and SE in GLMM (lmer)
> From: billpikou...@gmail.com
> To: patili_bue...@hotmail.com
> CC: r-help@r-project.org
>
> 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_bue...@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@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.
> >
_________________________________________________________________
Hotmail: Trusted email with Microsofts powerful SPAM protection.
[[alternative HTML version deleted]]
______________________________________________
R-help@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.
