Estimates and Integration with lme4 package
Closed this issue · 4 comments
I have been using the mitml package to calculate variance explained for a set of multilevel models, but noticed a few issues:
- When specifying the model using lme4, the estimates are different from when I specify the model using nlme. Specifically, the RB2 calculation changes depending on which package I use.
- When I run the function "multilevelR2" for my model that was specified using lme4 I get the following error: Error in multilevelR2(h1.1):Calculation of R-squared statistics not supported for models of class.
Thanks for reporting this!
Regarding issue 1
I can't reproduce this behavior. In the example below, the results are exactly the same. Generally, smaller differences can occur between nlme
and lme4
, but these are usually very small. Can you send me a specific example that I can use to reproduce this (e.g., via email)?
library(mitml)
library(lme4)
library(nlme)
data(studentratings)
# impute
imp <- panImpute(formula = ReadDis + ReadAchiev ~ 1 + (1|ID), data = studentratings, seed = 1234)
implist <- mitmlComplete(imp)
# fit models
fit1 <- with(implist, lmer(ReadDis ~ ReadAchiev + (1|ID)))
fit2 <- with(implist, lme(fixed = ReadDis ~ ReadAchiev, random = ~ 1 | ID,
data = data.frame(ReadDis, ReadAchiev, ID)))
# results as in output
testEstimates(fit1)
# Estimate Std.Error t.value df P(>|t|) RIV FMI
# (Intercept) 3.549 0.144 24.563 626.039 0.000 0.136 0.123
# ReadAchiev -0.002 0.000 -7.105 998.386 0.000 0.105 0.097
testEstimates(fit2)
# Estimate Std.Error t.value df P(>|t|) RIV FMI
# (Intercept) 3.549 0.144 24.563 626.039 0.000 0.136 0.123
# ReadAchiev -0.002 0.000 -7.105 998.385 0.000 0.105 0.097
# results with higher precision
testEstimates(fit1)$estimates
# Estimate Std.Error t.value df P(>|t|) RIV FMI
# (Intercept) 3.548747357 0.1444768667 24.562738 626.0392 0.000000e+00 0.1362350 0.12269860
# ReadAchiev -0.001957248 0.0002754753 -7.104983 998.3855 2.286393e-12 0.1049052 0.09675261
testEstimates(fit2)$estimates
# Estimate Std.Error t.value df P(>|t|) RIV FMI
# (Intercept) 3.548747357 0.1444768668 24.562738 626.0392 0.000000e+00 0.1362350 0.12269860
# ReadAchiev -0.001957248 0.0002754753 -7.104983 998.3855 2.286393e-12 0.1049052 0.09675261
# R²
multilevelR2(fit1)
# RB1 RB2 SB MVP
# 0.06090104 0.19417320 0.08411631 0.06834806
multilevelR2(fit2)
# RB1 RB2 SB MVP
# 0.06090105 0.19417318 0.08411631 0.06834806
Regarding issue 2
The error message indicates that the class attribute of the fitted models is changed from nlme
to something else. This can happen, for example, when using the lmerTest
package, which overrides many functions in lme4
. Do you use any additional packages that may cause this?
At the present time, lmerTest
is not supported by mitml
. Therefore, the only workaround is to fit the models without loading lmerTest
.
Thanks for the additional information. Unfortunately, I still can't reproduce this behavior. For example, with the studentratings
data set provided with mitml
(see below), the results are all fine (tested on Linux and Windows with the most recent versions of R and all relevant packages).
I would like to investigate this further, but I need a reproducible example that shows this behavior. Could you provide me with (1) a reproducible example with both data and code, and (2) the output of your sessionInfo()
after running that example? You can send them to my email address listed here.
library(mitml)
library(lme4)
library(nlme)
data(studentratings)
studentratings <- na.omit(studentratings[,c("ID","ReadDis","ReadAchiev")])
# fit models
fit1 <- lmer(ReadDis ~ ReadAchiev + (1|ID), data = studentratings)
fit1.1 <- lmer(ReadDis ~ ReadAchiev + (1|studentratings$ID), data = studentratings)
fit2 <- lme(fixed = ReadDis ~ ReadAchiev, random = ~ 1|ID, data = studentratings)
# R²
multilevelR2(fit1)
# RB1 RB2 SB MVP
# 0.06224511 0.21019639 0.08648722 0.07231612
multilevelR2(fit1.1)
# RB1 RB2 SB MVP
# 0.06224511 0.21019639 0.08648722 0.07231612
mitml:::.getRsquared(fit2, print=c("RB1", "RB2", "SB", "MVP"), method="nlme")
# RB1 RB2 SB MVP
# 0.06224532 0.21019331 0.08648691 0.07231607
Closing this for now, because the problem is still not reproducible with no more response.