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.
