[R-sig-ME] Struggle with glmer().

[Rolf Turner]

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)
}