0

I am trying to fit a logistic random intercept model using glmer function from package lme4. Unfortunately I am getting the following warning messages and clearly wrong results (for the coefficients).

Warning messages:
1: In vcov.merMod(object, use.hessian = use.hessian) :
  variance-covariance matrix computed from finite-difference Hessian is
not positive definite: falling back to var-cov estimated from RX
2: In vcov.merMod(object, correlation = correlation, sigm = sig) :
  variance-covariance matrix computed from finite-difference Hessian is
not positive definite: falling back to var-cov estimated from RX

Doing some research, I found that glmer produces merMod objects. Relying on this : http://jaredknowles.com/journal/2014/5/17/mixed-effects-tutorial-2-fun-with-mermod-objects and reproducing the "Exploring the Internals of a merMod Object" section I had the following results as far as my model is concerned:

### [1] "standardGeneric"
###attr(,"package")
###[1] "methods"

which is clearly different from

### [1] "lmerMod"
### attr(,"package")
### [1] "lme4"

as indicated on the tutorial.

My question why my object is not of merMod class? Are the above warnings related to this and how could I fix them?

This is the code I used to create the data frame and run the model

diagn00<- rep(0,240)

drug00<- rep(0,240)

time00<- c(rep(0,80),rep(1,80),rep(2,80))

response00<- c(rep(0,39),rep(1,41),rep(0,33),rep(1,47),rep(0,26),rep(1,54))

patients00<- rep(1:80,3)

test<- data.frame(patients00,diagn00,drug00,time00,response00)

diagn01<- rep(0,210)

drug01<- rep(1,210)

time01<- c(rep(0,70),rep(1,70),rep(2,70))

response01<- c(rep(0,33),rep(1,37),rep(0,15),rep(1,55),rep(0,2),rep(1,68))

patients01<- rep(81:150,3)

diagn10<- rep(1,300)

drug10<- rep(0,300)

time10<- c(rep(0,100),rep(1,100),rep(2,100))

response10<- c(rep(0,79),rep(1,21),rep(0,72),rep(1,28),rep(0,54),rep(1,46))

patients10<- rep(151:250,3)

diagn11<- rep(1,270)

drug11<- rep(1,270)

time11<- c(rep(0,90),rep(1,90),rep(2,90))

response11<-  c(rep(0,74),rep(1,16),rep(0,45),rep(1,45),rep(0,15),rep(1,75))

patients11<- rep(251:340,3)

diagnosis<- c(diagn00,diagn01,diagn10,diagn11)
diagnosis<- as.factor(diagnosis)

id<- c(patients00,patients01,patients10,patients11)
id<- as.factor(id)

drug<- c(drug00,drug01,drug10,drug11)
drug<- as.factor(drug)
time<- c(time00,time01,time10,time11)

response<- c(response00,response01,response10,response11)

id<- c(patients00,patients01,patients10,patients11)

data<- data.frame(id, response, diagnosis, drug, time)
e<- order(data$id)
d<- data[e,]

library(lme4)

d2<- data.frame(d)
d2$response<- as.factor(d2$response)
d2$time<- as.factor(d2$time)
d2$id<- as.factor(d2$id)

t<- glmer(response ~ diagnosis + drug + time + time:drug +
(1 | id), family=binomial, data=d2)

session info

e: [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages: [1] stats graphics grDevices utils datasets methods base

other attached packages: [1] glmmML_1.0 nlme_3.1-111 MASS_7.3-29 lme4_1.1-6 Rcpp_0.11.1 Matrix_1.1-3 MuMIn_1.10.0 [8] gee_4.13-18 geepack_1.1-6

loaded via a namespace (and not attached): [1] grid_3.0.2 lattice_0.20-29 minqa_1.2.3 RcppEigen_0.3.2.1.1 splines_3.0.2
[6] tools_3.0.2

package info

packageVersion("lme4")
[1] ‘1.1.6’
swihart
  • 2,648
  • 2
  • 18
  • 42
Manos
  • 103
  • 1
  • 1
  • 6
  • I can't reproduce this. If I put in your `control` argument I get an error ("fails with unused argument (maxIter=500)"). If I comment it out the model runs fine and produces apparently sensible results. Results of `sessionInfo()`/`packageVersion("lme4")` ... ? – Ben Bolker Jul 10 '14 at 10:01
  • You are right, I edited my code and now it produces the following warnings 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 266.714 (tol = 0.001) 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge: degenerate Hessian with 1 negative eigenvalues Is there anything I can do? – Manos Jul 10 '14 at 10:21
  • I would like to add that calling summary(t) produces the first two warnings I mentioned in the beginning. – Manos Jul 10 '14 at 10:25
  • I have looked at this problem some more, and it's *very* unstable -- very large parameter estimates suggest complete separation. I will keep digging, but this may not be easy to solve. – Ben Bolker Jul 11 '14 at 06:15
  • @Ben Bolker. I solved the warning issues by increasing the number of quadrature points, but, I can't figure out whether the random intercept is significant. By complete separation you mean that the random intercept is insignificant? – Manos Jul 11 '14 at 06:24

0 Answers0