[R-sig-ME] Why lme4 fails on schizo dataset, but glmmTMB succeeds?
Sun, John
j@un20 @end|ng |rom @|b@ny@edu
Tue Oct 11 22:17:30 CEST 2022
Dear All,
I am writing to ask why lme4 fails to work with random-slopes+intercepts on schizo dataset, but glmmTMB works.
>## data from http://www.uic.edu/~hedeker/SCHIZX1.DAT.txt
## inspiration from http://www.uic.edu/~hedeker/long.html
#may have to manually download from website.
fn <- "schizx1.dat"
if (!file.exists(fn)) {
download.file("http://www.uic.edu/~hedeker/SCHIZX1.DAT.txt",dest=fn)
}
## almost works: get extra ^Z, need to strip it ...
dd0 <- read.table("schizx1.dat")
names(dd0) <- c("id","imps79","imps79b","imps79o","int","tx",
"week","sweek","txswk")
apply(dd0,2,function(x) sum(x==-9,na.rm=TRUE))
## dangerous to use na.strings==c("-9","NA") in general?
## what if there are legitimate -9 values in *some* columns?
na.vals <- function(x,crit) {
x[crit] <- NA
return(x)
}
library(plyr) ## for mutate()
library(MASS) ## for glmmPQL
library(lme4)
dd <- mutate(dd0,
imps79=na.vals(imps79,imps79<0),
imps79b=na.vals(imps79b,imps79b<0),
imps79b=factor(imps79b,labels=c("le mild","ge moderate")),
tx=factor(tx,labels=c("placebo","drug")))
with(dd,table(tx,sweek))
with(na.omit(dd),table(tx,sweek))
## reduce this to sequences: most data are observed at only 4 time periods
dd$rweek <- round(dd$sweek^2)
dd2 <- na.omit(subset(dd,rweek %in% c(0,1,3,6)))
nrow(na.omit(dd))-nrow(dd2) ## lose only 34 cases (out of 1600)
dd3 <- ddply(dd2,"id",summarise,
tx=tx[1],
trans=paste(imps79b,collapse=""))
ttab <- with(dd3,table(trans,tx))
ttab[rowSums(ttab)>10,]
m0 <- glm(imps79b~tx*sweek,dd,family=binomial)
m3 <- glmer(imps79b~tx*sweek+(1|id),dd,family=binomial)
## extremely slow: lots of 'step-halving' action
m4 <- glmer(imps79b~tx*sweek+(1+sweek|id),dd,family=binomial,
verbose=100)
## eventually fails to converge with maxgrad=180.4 (!!)
library(glmmTMB)
m5<-glmmTMB(imps79b~tx*sweek+us(1+sweek|id),dd,family=binomial(link="logit"))
summary(m5)
Family: binomial ( logit )
Formula:
imps79b ~ tx * sweek + us(1 + sweek | id)
Data: dd
AIC BIC logLik deviance df.resid
1173.7 1211.4 -579.9 1159.7 1596
Random effects:
Conditional model:
Groups Name Variance Std.Dev. Corr
id (Intercept) 2184.9 46.74
sweek 365.9 19.13 -0.81
Number of obs: 1603, groups: id, 437
Conditional model:
Estimate Std. Error z value
(Intercept) 45.2732 7.6836 5.892
txdrug 0.8314 6.9096 0.120
sweek -14.6736 3.1210 -4.702
txdrug:sweek -7.3960 2.9027 -2.548
Pr(>|z|)
(Intercept) 3.81e-09 ***
txdrug 0.9042
sweek 2.58e-06 ***
txdrug:sweek 0.0108 *
---
Signif. codes:
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
https://github.com/lme4/lme4/blob/master/misc/issues/schizo.R
Best regards,
John
More information about the R-sig-mixed-models
mailing list