-
-
Notifications
You must be signed in to change notification settings - Fork 22
Description
Hello,
First, thank you for what you're doing with this package. It's been very helpful and user-friendly!
I'm trying to estimate a random effects model for an ordinal outcome using the clmm2 function from the ordinal package. However, when I try to estimate marginal means, I'm getting either a warning or an error depending on the specific model that says "Could not calculate standard errors for contrasts. This can happen when random effects are involved. You may try following:".
I'm unsure whether the message is meant to say more about what I can try? More importantly, I'm unsure what the underlying issue is/why this happens, though I recognize the example I'm sharing below may not be enough to uncover the statistical challenge beneath the warning.
Example that produces the warning:
> library(tidyverse)
> library(ordinal)
> library(modelbased)
>
> mtcars$carb <- factor(mtcars$carb)
> mtcars$gear <- factor(mtcars$gear, levels=c(3,4,5))
>
> mtcars <- mtcars %>% mutate(interact = interaction(vs, am, sep=":"))
> contrasts(mtcars$interact) <- contr.sum(levels(mtcars$interact))
>
> mod = clmm2(gear ~ interact, data=mtcars, random=carb, link="logistic", Hess=T, method="nlminb", doFit=T, model=T)
> summary(mod)
Cumulative Link Mixed Model fitted with the Laplace approximation
Call:
clmm2(location = gear ~ interact, random = carb, data = mtcars,
Hess = T, model = T, link = "logistic", doFit = T, method = "nlminb")
Random effects:
Var Std.Dev
carb 6.030997 2.455809
Location coefficients:
Estimate Std. Error z value Pr(>|z|)
interact1 -21.3452 307.5453 -0.0694 0.94467
interact2 -3.9713 173.1364 -0.0229 0.98170
interact3 14.2173 173.3322 0.0820 0.93463
No scale coefficients
Threshold coefficients:
Estimate Std. Error z value
3|4 -5.8279 173.1493 -0.0337
4|5 13.4510 173.3312 0.0776
log-likelihood: -10.1403
AIC: 32.2806
Condition number of Hessian: 555652.30
>
> est = estimate_means(mod, by="auto", estimate="average")
We selected `by=c("interact")`.
Could not calculate standard errors for contrasts. This can happen when random effects are involved. You may try following:
> est
Average Predictions
interact | Probability
----------------------
0:0 | 1.00
1:0 | 0.55
0:1 | 0.56
1:1 | 0.80
Variable predicted: gear
Predictors modulated: interact
Predictions are on the response-scale.
Notably, my own model has one level of its interaction that is missing by design (e.g., 0:1, though I have an interaction between a 3-level variable and a 4-level variable). I've dropped that level of the interaction, which is why I created a variable via interaction() rather than using the var1*var2 syntax. I don't have the sense that this is the problem, since I've been able to model binary outcome variables via glmer() using the same overall syntax without getting this warning.
This seems like it might be related to disussion #548 , though the error/warning is slightly different so I started a new issue.