Incorrect estimates, SE, and CI with nlme::lme() #295
Labels
Comments
|
For Problem 1, it happens whenever weights = ~ involves multiple variables, because weighting_var is expected to have length 1. See a simple case below. Without mean structure predictor data("Orthodont", package = "nlme")
summary(Model <- lme( # a model of variance only
distance ~ 1, data = Orthodont, # grand mean
weights = varConstPower(form = ~ age | Sex))) # var = b0 + |age|^b1
" AIC BIC logLik
514.5723 533.2821 -250.2861
Random effects:
Formula: ~1 | Subject
(Intercept) Residual
StdDev: 1.84895 0.109692
Variance function:
Structure: Constant plus power of variance covariate, different strata
Formula: ~age | Sex
Parameter estimates:
Male Female
const 0.0001144912 13.0138472
power 1.3237960288 -0.4748515
Fixed effects: distance ~ 1
Value Std.Error DF t-value p-value
(Intercept) 23.42394 0.4046777 81 57.88296 0"
find_predictors(Model) # NULL
find_predictors(Model, effects = "all", component = "all") # NULL
find_variables(Model)
'$response
[1] "distance"'
find_weights(Model) # correct
'"age" "Sex"'
formula(Model$modelStruct$varStruct) # ~age | Sex
names(get_data(Model)) # okay
ggpredict(Model) # empty list!!! Should give a grand mean
'list()
attr(,"class")
[1] "ggalleffects" "list"
attr(,"model.name")
[1] "Model"'
ggemmeans(Model) # error, but might be as expected
'"Error: `terms` needs to be a character vector with at least one
predictor name: one term used for the x-axis, more optional
terms as grouping factors.
In addition: Warning message:
In !is.null(weighting_var) && !weighting_var %in% colnames(mf) :
length(x) = 2 > 1 in coercion to logical(1)'
ggeffect(Model) # NULL, but might be as expectedWith mean structure predictor summary(Model <- lme( # a model of variance of multiple terms
distance ~ age, data = Orthodont, # mean change with age
weights = varConstPower(form = ~ age | Sex))) # var = b0 + |age|^b1
" AIC BIC logLik
439.3828 466.0172 -209.6914
Random effects:
Formula: ~age | Subject
Structure: General positive-definite
StdDev Corr
(Intercept) 1.6329350 (Intr)
age 0.1974838 -0.387
Residual 10.8517576
Variance function:
Structure: Constant plus power of variance covariate, different strata
Formula: ~age | Sex
Parameter estimates:
Male Female
const 0.1168654 9.637857e-05
power -1.4312273 -1.194815e+00
Fixed effects: distance ~ age
Value Std.Error DF t-value p-value
(Intercept) 17.268816 0.6113480 80 28.24711 0
age 0.610175 0.0589014 80 10.35926 0"
find_predictors(Model)
'$conditional
[1] "age"'
find_predictors(Model, effects = "all", component = "all") # same
find_variables(Model)
'$response
[1] "distance"
$conditional
[1] "age"'
find_weights(Model) # correct
'"age" "Sex"'
names(get_data(Model)) # okay
ggpredict(Model) # no CI!!!
'Error: Confidence intervals could not be computed.
$age
# Predicted values of distance
age | Predicted
---------------
8 | 22.15
10 | 23.37
12 | 24.59
14 | 25.81
attr(,"class")
[1] "ggalleffects" "list"
attr(,"model.name")
[1] "Model"
Warning messages:
1: In !is.null(weighting_var) && !weighting_var %in% colnames(mf) :
length(x) = 2 > 1 in coercion to logical(1)
2: In colnames(datlist)[ncol(datlist)] <- w :
number of items to replace is not a multiple of replacement length'
ggemmeans(Model) # same error as previous model, but might be as expected
"Error: `terms` needs to be a character vector with at least one
predictor name: one term used for the x-axis, more optional
terms as grouping factors.
In addition: Warning message:
In !is.null(weighting_var) && !weighting_var %in% colnames(mf) :
'length(x) = 2 > 1' in coercion to 'logical(1)'"
ggeffect(Model) # has CI but warning
'age | Predicted | 95% CI
--------------------------------
8 | 22.15 | [21.37, 22.93]
10 | 23.37 | [22.55, 24.19]
12 | 24.59 | [23.68, 25.50]
14 | 25.81 | [24.76, 26.87]
attr(,"class")
[1] "ggalleffects" "list"
attr(,"model.name")
[1] "Model"
Warning message:
In !is.null(weighting_var) && !weighting_var %in% colnames(mf) :
length(x) = 2 > 1 in coercion to logical(1)'SE and CI in the last |
|
Further note: in a different example, ggemmeans(lme(), type = "random") calculated SE much much larger than ggpredict() had (before recent updates, not ggpredict does not give any SE). |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Hi @strengejacke
There are a few problems with ggeffects() of a lme() model
Some similar issues were reported in {marginaleffects} at vincentarelbundock/marginaleffects#99 (comment)
Model set up. {insight} find_ and get_ seem to work well.
predict.lme() as a reference. It has a level = c() argument to include or not random effects by clustering level.
ggeffects results. See comments inside
The text was updated successfully, but these errors were encountered: