Error with nested random effect

[brian-lau]

Hi,

I am having trouble fitting a model with a nested random effect. The formula looks something like this:

y ~ fixedEffects + (CONDITION | SUBJECT/SIDE)

This fits fine using lmer, but I get the following error with rlmer

mr = rlmer(f,df,init=lmerNoFit(f,df))

Error in validObject(.Object) : 
  invalid class “rlmerMod” object: parameters are not allowed to be in multiple blocks simultaneously
In addition: Warning message:
In unlist(lapply(object@blocks, function(x) unique(as.vector(x)))) !=  :
  longer object length is not a multiple of shorter object length

Removing the nesting by removing /SIDE produces a warning:

Warning message:
In rlmer.fit.DAS.nondiag(lobj, verbose, max.iter, rel.tol) :
  algorithm converged, but estimating equations are not satisfied.

Any ideas?

Thanks in advance.

[BERENZ]

As @kollerma writes in the description "Both support only diagonal and unstructured random effect covariance structures." which means only random slopes are possible, not random coefficients.

[brian-lau]

Thanks for the pointer. However, the preceding sentence says "The two functions have the same limitations in the support of different random effect and residual error covariance structures." where the reference is to lmer and rlmer, suggesting that if the model can be fit using lmer it can also be fit using rlmer. And the random effects covariance is unstructured here.

[BERENZ]

Your right. I was wrong it is possible to fit a random coefficient models (they are even in examples).

[kollerma]

The first error message is from a test that checks that a random effects coefficient only appears in a single block type and not in multiple ones. It's hard to say whether the test is buggy or this is indeed the case here. Probably the former - these types of models should work.

But in any case, such a complicated model will require quite a big dataset or good starting values. This is what the second error message for the simplified model hints at. Try to start with a model as simple as possible and then use the parameter estimates of the simple model as initial parameters of the more complicated one (use the list form for the init argument). Iterate the process.

Another approach to get better initial coefficients is this. You should also be able to see from the robustness weights of a simpler model what observations are problematic in order to get an lmer fit on the rest of the data only. You could then use those parameter estimates as initial parameters for the full data fit with rlmer.

[kollerma]

In order to see whether this is a bug or not, it would be good if you could produce a simple example. Or maybe send me the data in private.

[brian-lau]

Thanks @kollerma for the suggestions. I started with simpler models using the same dataset with initial parameters from lmer. The robust models converge with a random intercept and even with nesting ((1|ID/S). This works even with the most complex fixed effects specification. However, adding a random intercept fails with the check mentioned initially, even with the simplest fixed effects specification:

f = as.formula("y ~ C + (1 + C|ID/S)")
mr = rlmer(f,data,verbose=1,init=lmerNoFit(f,data))

Error in validObject(.Object) : 
  invalid class “rlmerMod” object: parameters are not allowed to be in multiple blocks simultaneously
In addition: Warning message:
In unlist(lapply(object@blocks, function(x) unique(as.vector(x)))) !=  :
  longer object length is not a multiple of shorter object length

sessionInfo()
R version 3.3.3 (2017-03-06)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X El Capitan 10.11.6

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

other attached packages:
 [1] bindrcpp_0.2    ggeffects_0.2.1 lmerTest_2.0-33 robustlmm_2.1-3 perturb_2.05    sjPlot_2.3.3    MuMIn_1.15.6    r2glmm_0.1.2    effects_4.0-0  
[10] carData_3.0-0   car_2.1-5       ggplot2_2.2.1   lme4_1.1-13     Matrix_1.2-11  

I can send you data necessary to replicate this if you send me your email (brian.lau@upmc.fr).

[kollerma]

Hi Brian,

Sorry - this fell under the table. There is another branch that disables that check. I haven't had time to evaluate precisely whether the results will make sense, but plan to do so in the near future.

You can install the version of the branch using the following commands:

install.packages("robustbase", repos="http://R-Forge.R-project.org")
library(devtools)
install_github("kollerma/fastGHQuad")
install_github("kollerma/robustlmm", ref = "multi-grouping")

Best, Manuel

[kollerma]

This was fixed with the release of version 2.2-1.

I've disabled the test. Too complicated models with result in variance component estimates equal to zero. Note that estimating classic linear mixed effects models is challenging already and does not always work (hence the option of lme4 to use a different optimizer). Doing everything robustly does not make this easier.