2
$\begingroup$

I am trying to fit a beta regression to my data using mixed models, as there are 4 repeated observations per subject.

Legend:

p = (time) point: t1...t4

ID = subject ID

When I try:

summary(glmmTMB(y ~ p +(1|ID), dat, family=beta_family(link = "logit")))

I get the following output

 Family: beta  ( logit )
Formula:          y ~ p + (1 | ID)
Data: dat

     AIC      BIC   logLik deviance df.resid 
  -664.5   -640.6    338.3   -676.5      394 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 ID     (Intercept) 0.01381  0.1175  
Number of obs: 400, groups:  ID, 100

Overdispersion parameter for beta family (): 15.4 

Conditional model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.62540    0.06498  -25.01   <2e-16 ***
pt2          1.42948    0.08081   17.69   <2e-16 ***
pt3          2.27841    0.08282   27.51   <2e-16 ***
pt4          3.04969    0.08894   34.29   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

But when I allow the effect of p to vary from subject to subject, adding p also as a random slope, I get different fixed effects. This surprises me.

summary(glmmTMB(y ~ p +(p|ID), dat, family=beta_family(link = "logit")))


Family: beta  ( logit )
Formula:          y ~ p + (p | ID)
Data: dat

     AIC      BIC   logLik deviance df.resid 
      NA       NA       NA       NA      385 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev. Corr              
 ID     (Intercept) 0.2296   0.4791                     
        pt2         0.2723   0.5218   -0.99             
        pt3         0.1975   0.4444   -0.96  0.95       
        pt4         0.4276   0.6539   -0.40  0.49  0.24 
Number of obs: 400, groups:  ID, 100

Overdispersion parameter for beta family (): 26.2 

Conditional model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.67415    0.07124  -23.50   <2e-16 ***
pt2          1.47345    0.08351   17.64   <2e-16 ***
pt3          2.34437    0.08028   29.20   <2e-16 ***
pt4          3.11739    0.09856   31.63   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The fixed estimates are very similar, yet different. I guess it should be exactly the same.

Also a warning is displayed:

1: In fitTMB(TMBStruc) : Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting') 2: In fitTMB(TMBStruc) : Model convergence problem; singular convergence (7). See vignette('troubleshooting')

Same happens, when I set REML=TRUE.

Can just adding random slopes alter the fixed part? Or is this just an optimization issue?

My data is quite large, so I didn't put it here.

EDIT: Is this possibly related to this issue? Why and how does the inclusion of random effects in mixed models influence the fixed-effect intercept term?

So, to answer your general question: if you have a nonlinear link, a random effect can change the intercept because the random term is applied on a nonlinear link function.

$\endgroup$
3
$\begingroup$

The time variable p is a factor. Hence, when you include it in the random effects you specify that you want a different random effect per time point. Also, the covariance matrix for these random effects is taken to be unstructured. This is a rather complex model that is asking too much from the data you have. This is probably why the second model you fitted has not actually converged.

You could try treating the time variable as numeric at least in the random-effects part of the model.

In addition, note that glmmTMB() fits the model using the Laplace approximation that is known to be less optimal than the adaptive Gaussian quadrature one. The latter is provided in the GLMMadaptive package. You can find an example of a Beta mixed model in this package here.

$\endgroup$
  • $\begingroup$ Thank you very much, Sir. I had no idea such a powerful tool exists. And you prepared a lot of examples. Very nice documentation. I am truly impressed. Is it possible, to specify the covariance structure for either the random effects or the residuals? Like AR1, for instance? I could not find any relevant option. $\endgroup$ – humbleasker Nov 7 '19 at 19:17
  • $\begingroup$ At the moment these options are not available in GLMMadaptive. But I do plan to include them. $\endgroup$ – Dimitris Rizopoulos Nov 8 '19 at 15:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.