5

I am fitting a linear mixed effect model in R (function lme), and I get a Var-Cov matrix with negative entries (Log-Cholesky). This does not allow me to compute confidence intervals on the standard deviations of the random effects. Shouldn't this matrix theoretically be postive-definite? Does this warn me of something specific?

The model has two fixed effects (factors) with interaction term, and nested random effects on intercept and slope.

[UPDATE]

I have uploaded the data at this link. The explanation of the experiment and the data-set are detailed in this question. To sanity-check what the problem might be, I first simplify the dataset by considering a single level of the factor spd_des:

> spdDes <- 's15'
> dat <- dat[dat$spd_des==spdDes,]

Nevertheless, when I try to fit a model with random slope on mPair, I get negative entries on the var-cov matrix. 

> lmer(cc_marg ~ mPair + (mPair|ratID), data = dat, REML=TRUE, na.action=na.omit)

This is true using lme, lmer, as well as stan_lmer. In particular, stan_lmer returns the warning There were 8 divergent transitions after warmup. When I plot the traces with stan_trace I do not see any convergence: all the traces look like random noise.

Robert Long
  • 53,316
  • 10
  • 84
  • 148
Cristiano
  • 507
  • 2
  • 11
  • What is the model formula ? – Robert Long May 15 '19 at 15:57
  • `lme(cc_marg ~ mPair*spd_des , random = ~mPair|ratID/spd_des/cycle, data=dat_trf, na.action=na.omit, method = "ML", control=lCtr )` – Cristiano May 15 '19 at 16:22
  • That's a lot of random effects. Do you really need random slopes over all 3 levels ?The data may not support such a structure. Is there aný special reason for using `lme` rather than, for example, `lmer` in the `lme4` package or `mixed_model` in `GLMMadaptive` ? – Robert Long May 15 '19 at 16:29
  • I was actually trying to figure out how to have random slope only on ratId/spd_des, and random intercept also on ratID/spd_des/cycle. How do I set this up? As for lme, I thought it was as good of a package as the other ones. Any specific reason to switch to lmer? – Cristiano May 15 '19 at 16:38
  • `lme4` is much more recent (same primary author) with better support to diagnose and identify problems like these and it allows much more flexibility with random effects provided that you dont want to model the residual covariance structure. – Robert Long May 15 '19 at 16:44
  • So, the reason my just be that the model is too complicated? You confirm that something is wrong if the var-cov matrix has negative entries? – Cristiano May 15 '19 at 17:42
  • Something is definitely wrong, and my guess is that the model is misspecified but there could be other reasons, it's very hard to say without much more detail about the data. Try lme4, there is much support to diagnose problems, or tell us more about the data. – Robert Long May 15 '19 at 17:54

1 Answers1

3

Very often this problem occurs in the optimization routine used by the software, and indicates that the the solution is on or very near the boundary of the parameter space, usually in mixed models, that is a variance component of zero (or a correlation between random effects of 1 or -1). It may be the case that the particular variance component actually is (or is very near) zero (in which case it can be omitted from the model), or that the variance component is meaningfully positive, but the optimization routine fails, often because the number of clusters/levels of the relevant factor is small. Sometimes this can be overcome by using a different optimiser. In R, it is often better, when fitting linear mixed models without a residual covariance structure to use one of the newer packages such as lme4, which offers much better support for diagnosing convergence problems. Another option is to use a Bayesian approach, such as the rstanarm package in R which again can be useful for diagnosing convergence issues. Another package that may be helpful is GLMMadaptive which implements Gaussian Adaptive Qundrature rather than the Laplace method.

Lastly, I would say that in this specific case, the random effects structure is quite complex as it is fitting nested random effects for 3 levels, and random slopes over all 3 levels, and their relevant correlations. This may be an over-specified model, particularly if any of the numbers of clusters are small. I would certainly start with a much simpler model in lme4 with just random intercepts, and add random slopes if warranted by theory to the relevant levels one at a time. Uncorrelated random effects can also be specified if necessary. Building a model up this way will often be a good way to identify and overcome these kinds of problem

Demidenko's book "Mixed Models. Theory and Applications in R" discusses these problems in depth in section 2.15.

The answers to this question and this question are also relevant.

Robert Long
  • 53,316
  • 10
  • 84
  • 148
  • I followed your suggestion of using `lme4` and simplifying the model. Unfortunately I have the same issue. First, I have reduced the dataset to a single level of the factor `spd_des`, so I could remove it from the model. Then, I fit a model with random intercept only without any issue. However, when I fit the following model `lmer(cc_marg ~ mPair + (mPair|ratID), data = dat_trf, REML=TRUE, na.action=na.omit)` I still get negative entry in the var-cov matrix. Knowing the underlying question, I do think that random effects can indeed be highly positively correlated. I am not sure how to proceed. – Cristiano May 24 '19 at 22:05
  • is `mpair` a factor ? How many levels does it have ? How many `ratID`s are there ? How many observations in total and are there any singleton clusters? – Robert Long May 25 '19 at 04:57
  • `mPair` is a factor with 6 levels. There are 10 levels of `ratID`. For each of those levels of `ratID` I have about 100 (different number across rats, with a minimum of 30) observations for each level of `mPair`. There are a few combinations of `mPair` and `ratID` with no data (not many), but there is no rats with no observations in any mPair. – Cristiano May 25 '19 at 15:16
  • Seems reasonable. Try plugging in `rstanarm` with identical syntax then look at trace plots to see if it is converging to reasonable values. – Robert Long May 25 '19 at 17:35
  • ok so...I fit the simplified model with `stan_lmer`. I still obtain negative entries in the var-cov matrix. These entries are all close to zero, and are similar to those obtained with `lmer`. `stan_lmer` returns the warning `There were 8 divergent transitions after warmup.` When I plot the traces with `stan_trace` I do not see any convergence: all the traces look like random noise. Any clue? – Cristiano May 29 '19 at 17:49
  • i realized that it is not true that all the negative entries are close to zero. There is also a -8 and a -9. – Cristiano May 29 '19 at 20:14
  • Can you post your data? – Robert Long May 30 '19 at 08:27
  • I posted the data as an update to the main question. Thanks a lot for your help! – Cristiano Jun 03 '19 at 16:48