[R-sig-ME] Struggle with glmer().
Rolf Turner
r.turner at auckland.ac.nz
Sun Jul 16 11:57:07 CEST 2017
To all and sundry, near and far,
Ben Bolker in particular :-) .... [1]
For a considerable time I have been struggling with a problem that would
appear to require application of the glmer() function from lme4. The
main thrust of the problem involves estimation of lethal dose levels
(e.g. LD99) for a certain pest elimination treatment. It seems to be
impossible to get glmer() to return a fit of the relevant model without
producing a number of dire and ominous warnings.
The principal investigator on the problem has managed to get a
linearisation of the model to produce an almost-satisfactory fit, but
that's an "almost", but not quite. I won't go into details about this
now, because it's mostly off the point.
The data involved are confidential, so I can't show you those; instead I
have written a function to generate artificial data (which are
superficially similar to the real data but are much tidier and without
the quirks and peculiarities of the real data). The artificial data are
nicely balanced (unlike the real data).
The simulated data set has 1080 records [2]. I would not have thought
this to be an overwhelmingly large data set.
I have attached the code of my data generating function "artSim()".
The coefficients of the fixed effects ("beta0" and "beta1") and the
variance components ("Sigma") are designed to be roughly similar to
estimates obtained from the "linearisation" mentioned above.
Having written "artSim()" I tried:
set.seed(42)
X <- artSim()
library(lme4)
fit <- glmer(cbind(Dead,Alive) ~ (Trt + 0)/x + (x | Rep),
family=binomial(link="cloglog"),data=X)
After a fairly lengthy wait, I was not too surprised to get the same
ominous warnings that I got with the real data:
> Warning messages:
> 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
> unable to evaluate scaled gradient
> 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
> Model failed to converge: degenerate Hessian with 1 negative eigenvalues
Could I prevail upon you wise and knowledgeable people to have a look at
my simulation function and the fitting procedure and let me know if I am
doing something silly?
One last comment --- with the real data, I have tried using estimates
obtained from the "linearised" model as starting values for glmer().
This simply made matters worse! In addition to the foregoing warnings I
got a warning about "failure to converge in 10000 evaluations".
I haven't bothered to try this device with the simulated data ....
Thanks for any avuncular (or materteral!) advice that anyone can provide.
cheers,
Rolf Turner
--
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276
[1] Cf. "King John's Christmas" by A. A. Milne
[2] The pest elimination treatment under study is *not* 1080!!! The
presence of that number is purely coincidental!!! :-)
-------------- next part --------------
artSim <- function(){
#
# Function to simulate "artificial" data which is at least superficially
# similar to some real data.
#
require(MASS)
link <- "cloglog"
B <- binomial(link=link)
linkfun <- B$linkfun
linkinv <- B$linkinv
# Construct (artificial) treatment factor, covariate, and
# (random) replicate factor.
x <- seq(0,28,by=2)
Trt <- LETTERS[1:24]
Rep <- 1:3 # Three reps per treatment.
Xdat <- expand.grid(x=x,Trt=Trt,Rep=Rep)
uRep <- with(Xdat,factor(paste0(Rep,Trt)))
Xdat$Rep <- with(Xdat,factor(as.numeric(uRep)))
beta0 <- seq(-3,0.45,by=0.15)
beta1 <- rep(seq(0.05,0.3,by=0.05),4)
names(beta0) <- Trt
names(beta1) <- Trt
Sigma <- matrix(c(0.06,-0.001,-0.001,0.0001),nrow=2)
lb0 <- beta0[match(Xdat$Trt,names(beta0))]
lb1 <- beta1[match(Xdat$Trt,names(beta1))]
nrep <- 72
imat <- match(Xdat$Rep,1:nrep)
Z <- mvrnorm(nrep,c(0,0),Sigma)[imat,]
linpr <- lb0 + Z[,1] + (lb1 + Z[,2])*Xdat$x
p <- linkinv(linpr)
nsize <- 25
Dead <- rbinom(nrow(Xdat),nsize,p)
Alive <- nsize - Dead
x0 <- (linkfun(0.99) - beta0)/beta1
Xdat$Dead <- Dead
Xdat$Alive <- Alive
attr(Xdat,"trueLD99") <- x0
return(Xdat)
}
More information about the R-sig-mixed-models
mailing list