[R-sig-ME] lmer nonconvergent: care to run and explain?
Paul Johnson
pauljohn32 at gmail.com
Wed Oct 14 19:02:13 CEST 2015
Every now and then, I come back to my old project that simulates
multi-level data. I've asked about this from time to time here.
Today I noticed some unfamiliar convergence problems in lme4_1.1-9. I
updated to lme4-1.1-10 and got new and different problems. Rescaling
variables did solve convergence issues in penultimate version, but not
any more.
Here's my data simulation code. Appears to me this should work. I'm
giving lots of groups, plenty of observations in each one. Maybe my
random effect variance is too small, but it seems to me that should
not drive convergence problems. I've got Stata output for comparison
below.
Am I setting conditions that make this illconditioned?
library(MASS) ## for rmvnorm
set.seed(1234)
## Clusters
M <- 100
## Minimum per cluster
Mmin <- 134
## Variance in number of members per cluster
Mvar <- 20
N <- Mmin + rpois(M, lambda = Mvar)
Mcolors <- rainbow(M)
## Create Mind and Iind indexes for book keeping.
Mind <- unlist(lapply(1:M, function(x) rep(x, each = N[x])))
Iind <- unlist(lapply(1:M, function(x) seq(1, N[x], by = 1)))
##Create 3 variables from a multivariate Normal with these
## characteristics
Xmeans <- c(100, 200, 150)
Xsds <- c(20, 30, 40)
Xrho12 <- 0.4
Xrho13 <- 0.2
Xrho23 <- 0.0
Xcorr.mat <- diag(1, 3)
Xcorr.mat[1,2] <- Xcorr.mat[2,1] <- Xrho12
Xcorr.mat[1,3] <- Xcorr.mat[3,1] <- Xrho13
Xcorr.mat[2,3] <- Xcorr.mat[3,2] <- Xrho23
if(!isSymmetric(Xcorr.mat))stop("Xcorr.mat is not symmetric")
Xsigma <- diag(Xsds) %*% Xcorr.mat %*% diag(Xsds)
X.mat <- mvrnorm(n = sum(N), mu = Xmeans, Sigma = Xsigma)
dimnames(X.mat)[[2]] <- c("x1", "x2", "x3")
## The true fixed effects of beta0, beta1, beta2, beta3 in
## y = beta0 + beta1*x1 + beta2*x2 + beta3*x3
beta <- c(1.5, 0.25, -0.2, 0.05)
## Standard deviation of error term at individual level
STDE <- 25
## Create a dependent variable that has no clustering effects.
## FIXED Part + RANDOM PART
y1 <- cbind(1, X.mat) %*% beta + rnorm(sum(N), m = 0, s = STDE)
dat <- cbind(data.frame(Mind, Iind, y1), as.data.frame(X.mat))
rownames(dat) <- paste(dat$Mind, dat$Iind, sep=".") ##may help bookkeeping later
rm(Mind, Iind, y1, X.mat) ## cleanup workspace
summary(lm(y1 ~ x1 + x2 + x3, data=dat))
## Layer on additive group level error in y2 and y3
## Parameters for random effects:
## STDEb0: standard deviation of clustered intercepts.
STDEb0 <- 0.5
## STDEb1: standard deviation of slopes across cluster units
STDEb1 <- 0.25
Mrho <- 0.2
## I'm tempted to get fancy here with a matrix for STDEb0, STDEb1, and
## Mrho. It would pay off, I expect. But be harder to teach.
Msigma <- diag(c(STDEb0, STDEb1)) %*% matrix(c(1, Mrho, Mrho, 1), 2,
2) %*% diag(c(STDEb0, STDEb1))
Mreffects <- mvrnorm(M, mu = c(0, 0), Sigma = Msigma)
colnames(Mreffects) <- c("b0", "b1")
## In y2, include random intercept
dat$y2 <- dat$y1 + Mreffects[ ,"b0"][dat$Mind]
## In y3, add in random slope effect
dat$y3 <- dat$y2 + dat$x1 * Mreffects[ ,"b1"][dat$Mind]
## Do some diagnostics on Mreffects
plot(Mreffects[, 1], Mreffects[,2], xlab = "Intercept Random Effect",
ylab = "Slope Random Effect", main = "True Random Effects")
## library(rockchalk)
## summarize(Mreffects)
summary(Mreffects)
apply(Mreffects, 2, sd)
cor(Mreffects)
plot(Mreffects[ ,"b0"][dat$Mind], dat$x1 * Mreffects[
,"b1"][dat$Mind], col=dat$Mind,
main = "Random Slopes and Intercepts",
xlab = "Intercept Random Effect (b0)",
ylab = "Slope Random Effect (b1*x)")
## The "grand" regression?
m1 <- lm(y3 ~ x1 + x2 + x3, data=dat)
summary(m1)
## The "dummy variable" regression?
m2 <- lm(y3 ~ x1 + x2 + x3 + as.factor(Mind), data=dat)
summary(m2)
library(lme4)
## M separate regressions, with 3 predictors
m3list <- lmList(y3 ~ x1 + x2 + x3 | Mind, data = dat, pool = FALSE)
## Predicted values set x2 and x3 at their cluster-specific means.
plot(y3 ~ x1, data = dat, col = Mcolors[Mind], lwd = 0.6, main = "lm
on Separate Clusters")
for( i in seq_along(m3list)){
m3mf <- model.frame(m3list[[i]]) #data set for group i
x1range <- range(m3mf$x1) ## use group-specific ranges this time
pgroup <- predict( m3list[[i]], newdata = data.frame(x1=x1range,
x2=mean(m3mf$x2), x3=mean(m3mf$x3)))
lines(x1range, pgroup, col=Mcolors[i])
}
## If a cluster has 3 or fewer cases, this warning will appear
## Warning message:
## In predict.lm(m3list[[i]], newdata = data.frame(x1 = x1range, x2 =
mean(m3mf$x2), :
## prediction from a rank-deficient fit may be misleading
## Record keeping
m3newdat <- lapply(m3list, function(x) {
m3mf <- model.frame(x)
ndf = data.frame(x1 = range(m3mf$x1), x2 = mean(m3mf$x2), x3 =
mean(m3mf$x3))
ndf$m3pred <- predict(x, newdata = ndf)
ndf} )
m3newdat <- do.call("rbind", m3newdat)
## Recover the Mind value by butchering the row names
m3newdat$Mind <- as.integer(gsub("\\.[0-9].*", "", row.names(m3newdat)))
## Draw new graphs on a new device, so we can compare
dev.new()
##
## Estimate this as a mixed-effects model with lme4
## mm2: Just the random intercept, no random slope
mm2 <- lmer(y2 ~ x1 + x2 + x3 + (1 | Mind), data=dat)
summary(mm2)
mm2VarCorr <- VarCorr(mm2)
mm2VarCorr
cor(fitted(mm2), dat$y2)
## Both random intercept and random slope, not correlated with each other
mm3 <- lmer( y3 ~ x1 + x2 + x3 + (1|Mind) + (0 + x1 | Mind), data=dat,
verbose=3)
## That produces a non-convergent model.
## At return
## eval: 479 fn: 153252.65 par: 0.0237800 0.00122944
## Warning messages:
## 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 1.54298 (tol = 0.002,
component 1)
## 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Ok, lets scale them
dat[ , c("x1s", "x2s", "x3s")] <- lapply(dat[, c("x1", "x2", "x3")], scale)
dat$y3s <- scale(dat$y3)
mm3s <- lmer( y3s ~ x1s + x2s + x3s + (1 | Mind) + (0 + x1 | Mind),
data=dat, verbose = 3)
mm4 <- lmer( y3 ~ x1 + x2 + x3 + (x1 | Mind), data = dat, verbose = 2)
summary(mm4)
mm4s <- lmer( y3s ~ x1s + x2s + x3s + (x1s | Mind), data = dat, verbose = 2)
summary(mm4)
Here's what I see at the end.
At return
eval: 453 fn: 34253.131 par: 1.45320e-05 0.00891596
Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.00204447 (tol = 0.002,
component 1)
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?
> sessionInfo()
R version 3.2.2 (2015-08-14)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 15.04
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lme4_1.1-10 Matrix_1.2-2 MASS_7.3-44
loaded via a namespace (and not attached):
[1] compiler_3.2.2 minqa_1.2.4 tools_3.2.2 Rcpp_0.12.1
[5] splines_3.2.2 nlme_3.1-122 grid_3.2.2 nloptr_1.0.4
[9] lattice_0.20-33
Earlier today, with the previous lme4, the scaling of the variables
produced a convergent model.
Upon exporting the data to Stata 14, we observe:
. mixed y1 x1 x2 x3 || Mind: x1, cov(unstruct)
Performing EM optimization:
Performing gradient-based optimization:
Iteration 0: log likelihood = -71367.808
Iteration 1: log likelihood = -71366.333
Iteration 2: log likelihood = -71366.311
Iteration 3: log likelihood = -71366.308
Iteration 4: log likelihood = -71366.304
Iteration 5: log likelihood = -71366.304
Computing standard errors:
Mixed-effects ML regression Number of obs = 15419
Group variable: Mind Number of groups = 100
Obs per group: min = 146
avg = 154.2
max = 169
Wald chi2(3) = 1095.11
Log likelihood = -71366.304 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
y1 | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
x1 | .2534893 .0110537 22.93 0.000 .2318245 .2751541
x2 | -.1914348 .0072503 -26.40 0.000 -.2056452 -.1772245
x3 | .050516 .0051398 9.83 0.000 .0404421 .0605899
_cons | -.8820319 1.575815 -0.56 0.576 -3.970573 2.206509
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval]
-----------------------------+------------------------------------------------
Mind: Unstructured |
var(x1) | 7.47e-07 . . .
var(_cons) | .2079643 . . .
cov(x1,_cons) | -.0003941 . . .
-----------------------------+------------------------------------------------
var(Residual) | 613.2876 . . .
------------------------------------------------------------------------------
LR test vs. linear regression: chi2(3) = 0.06 Prob > chi2 = 0.9964
Note: LR test is conservative and provided only for reference.
Fixed parameter estimates are good, and the residual variance is not
far wrong. The estimated variance of the slope random effect is
miniscule. But there's no convergence error.
Maybe I'm not generating the random effect for the regression on x1
correctly. In the simulated data,
> rockchalk::summarize(Mreffects)
$numerics
Mreffects_1 Mreffects_2
0% -1.4122 -0.8380
25% -0.3799 -0.1887
50% -0.0438 -0.0244
75% 0.2899 0.0976
100% 0.9951 0.5505
mean -0.0615 -0.0495
sd 0.4650 0.2455
var 0.2162 0.0602
NA's 0.0000 0.0000
N 100.0000 100.0000
$factors
NULL
> cor(Mreffects)
b0 b1
b0 1.0000000 0.1445065
b1 0.1445065 1.0000000
--
Paul E. Johnson
Professor, Political Science Director
1541 Lilac Lane, Room 504 Center for Research Methods
University of Kansas University of Kansas
http://pj.freefaculty.org http://crmda.ku.edu
More information about the R-sig-mixed-models
mailing list