gamlss fails with random effects

[mattansb]

It seems like emmeans fails for models fit with random effects (via the re() function term):

# sim some data
n <- 10
k <- 2

mus <- rnorm(n, mean = 100, sd = 20)
mus_cond <- rnorm(k, mean = 100, sd = 20)
sigmas <- rexp(n, rate = 1/70)
taus <- rexp(n, rate = 1/50)

ddd <- expand.grid(Subj = factor(1:n),
                   cond = factor(1:k),
                   trial = 1:50)
ddd$mu <- mus[ddd$Subj] + mus_cond[ddd$cond]
ddd$sigma <- sigmas[ddd$Subj]
ddd$tau <- taus[ddd$Subj]
ddd$rt <- rexGAUS(nrow(ddd),
                  mu = ddd$mu,
                  sigma = ddd$sigma,
                  nu = ddd$tau)

m <- gamlss(rt ~ cond + re(random = ~1|Subj),
            ~ cond + re(random = ~1|Subj),
            ~ cond,
            family = exGAUS(),
            data = ddd)

library(emmeans)

emmeans(m, ~cond, what = "mu")
#> Error in model.frame.default(formula, data = data, ...) : 
#>   variable lengths differ (found for 're(random = ~1 | Subj)')

Same model without random effects:

m <- gamlss(rt ~ cond,
            ~ cond,
            ~ cond,
            family = exGAUS(),
            data = ddd)

emmeans(m, ~cond, what = "mu")
#>  cond emmean   SE  df asymp.LCL asymp.UCL
#>  1       140 6.53 Inf       127       153
#>  2       138 6.22 Inf       126       150
#> 
#> Confidence level used: 0.95

[rvlenth]

I'm not sure if this is fixable. Part of the issue is that in trying to recover the data, I still have the extra re() term in the formula. But even if I remove that, I still hit another issue in that apparently this invokes a smoother, and (as is stated in the "models" vignette), smoothers are as yet unsupported. Note the presence of mu.s and sigma.s among the slots; those are matrices for smoothers.

> names(mm)
 [1] "family"             "parameters"         "call"               "y"                 
 [5] "control"            "weights"            "G.deviance"         "N"                 
 [9] "rqres"              "iter"               "type"               "method"            
[13] "contrasts"          "converged"          "residuals"          "noObs"             
[17] "mu.fv"              "mu.lp"              "mu.wv"              "mu.wt"             
[21] "mu.link"            "mu.terms"           "mu.x"               "mu.qr"             
[25] "mu.coefficients"    "mu.offset"          "mu.xlevels"         "mu.formula"        
[29] "mu.df"              "mu.nl.df"           "mu.s"               "mu.var"            
[33] "mu.coefSmo"         "mu.lambda"          "mu.pen"             "df.fit"            
[37] "pen"                "df.residual"        "sigma.fv"           "sigma.lp"          
[41] "sigma.wv"           "sigma.wt"           "sigma.link"         "sigma.terms"       
[45] "sigma.x"            "sigma.qr"           "sigma.coefficients" "sigma.offset"      
[49] "sigma.xlevels"      "sigma.formula"      "sigma.df"           "sigma.nl.df"       
[53] "sigma.s"            "sigma.var"          "sigma.coefSmo"      "sigma.lambda"      
[57] "sigma.pen"          "nu.fv"              "nu.lp"              "nu.wv"             
[61] "nu.wt"              "nu.link"            "nu.terms"           "nu.x"              
[65] "nu.qr"              "nu.coefficients"    "nu.offset"          "nu.xlevels"        
[69] "nu.formula"         "nu.df"              "nu.nl.df"           "nu.pen"            
[73] "P.deviance"         "aic"                "sbc"

This is confirmed by looking at the documentation for random() and re(), that they return a function of class "smooth"; I currently just don't know how to handle this.

BTW, please don't reuse the name of an object in a bug report (in this case two different ms) -- that causes confusion. In this case, I named the first model mm so I can distinguish it from the second one.

[mattansb]

I see...

Does this also mean that I can't use qdrg() here? Because the SEs wont be correct?

qdrg(
  ~ cond,
  data = ddd,
  coef = coef(m, what = "mu")[1:2],
  vcov = vcov(m)[1:2, 1:2]
)
#> 'emmGrid' object with variables:
#>     cond = 1, 2

[rvlenth]

Yes, I kind of think that's the case. But one way to find out is to fit a simple mixed model (just location) using games and the same using glmer or something, and compare the results. You'd need to pick an example where there are non-negligible random effects.

[rvlenth]

Here's such a comparison, I think:

> m3 <- gamlss(rt ~ cond + re(random = ~1|Subj),
+              data = ddd)
GAMLSS-RS iteration 1: Global Deviance = 12402.62 
GAMLSS-RS iteration 2: Global Deviance = 12402.62 

> m3.lme <- nlme::lme(rt ~ cond, random = ~1|Subj,
+              data = ddd)

> m3.rg = qdrg(
+     ~ cond,
+     data = ddd,
+     coef = coef(m3, what = "mu")[1:2],
+     vcov = vcov(m3)[1:2, 1:2]
+ )
Warning message:
In vcov.gamlss(m3) : Additive terms exists in the  mu formula. 
  Standard errors for the linear terms maybe are not appropriate

> confint(m3.rg)
 cond prediction   SE  df asymp.LCL asymp.UCL
 1           248 5.34 Inf       237       258
 2           272 5.34 Inf       262       283

Confidence level used: 0.95 

> confint(ref_grid(m3.lme))
 cond prediction   SE df lower.CL upper.CL
 1           248 11.6  9      222      274
 2           272 11.6  9      246      298

Degrees-of-freedom method: containment 
Confidence level used: 0.95 

The SEs from the gamlss fit via qdrg() are quite a bit smaller. Are these comparable models? I think so, but may be wrong.

[mattansb]

That seems correct to me.

Well.... that sucks :/

Thanks Russell!

[rvlenth]

I agree, and I will try to understand this. It's hard to understand why accounting for a random effect can be implemented by a smoother. And maybe I can try again the get the vamps developer interested.

Russ

Sent from my iPad

On May 27, 2021, at 12:06 AM, Mattan S. Ben-Shachar @.***> wrote:



That seems correct to me.

Well.... that sucks :/

Thanks Russell!

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/rvlenth/emmeans/issues/277#issuecomment-849322625, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AGMJPL2ERAA3LFRGDP66ZR3TPXHLNANCNFSM45M7IBIA.