0
$\begingroup$

I selected the best fitted gamlss model for my data :

gamlss(formula = Ratio ~ PROTECTION * JOUR, nu.formula = ~JOUR, family = BEINF, data = D_E1, trace = FALSE)

I want to interpret the output of summary() to conclude which level of PROTECTION, JOUR and their INTERACTION are significantly different from each other.

Summary for mu gives :

Mu link function:  logit
Mu Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)         1.4235     0.7227   1.970 0.054443 .  
PROTECTIONG        -1.9529     0.9998  -1.953 0.056405 .  
PROTECTIONR        -3.2066     0.8776  -3.654 0.000619 ***
PROTECTIONT        -2.9102     0.8705  -3.343 0.001575 ** 
JOUR2              -1.8951     0.8462  -2.239 0.029605 *  
JOUR3              -1.4190     0.9103  -1.559 0.125322    
PROTECTIONG:JOUR2   2.6029     1.2884   2.020 0.048726 *  
PROTECTIONR:JOUR2   1.7200     1.2066   1.426 0.160223    
PROTECTIONT:JOUR2   2.4595     1.1220   2.192 0.033061 *  
PROTECTIONG:JOUR3   0.6028     1.2225   0.493 0.624094    
PROTECTIONR:JOUR3   1.7602     1.1407   1.543 0.129117    
PROTECTIONT:JOUR3   3.5884     1.1042   3.250 0.002068 **  

I understand that because there are 2 categorial variables, we can't interpretate the overall estimate for PROTECTION_G without considering the other variable. The log odds for PROTECTION_G and JOUR_1 (base) would be (intercept + estimate PROTECTION_G). The log odds for PROTECTION_G and JOUR_2 would be (intercept + estimate PROTECTION_G+ estimate JOUR_2)

Thus, I think I should rather calculate and compare the estimated means of the interaction, but it seems that the formula EXP(intercept + estimate)/(1+EXP(intercept+estimate) does not applicate to the interaction levels as it gives me :

                           Estimated
                             mean
PROTECTIONR:JOUR2   1.7200 -> 0.96   1.2066   1.426 0.160223    
PROTECTIONT:JOUR2   2.4595 -> 0.98   1.1220   2.192 0.033061 *  
PROTECTIONG:JOUR3   0.6028 -> 0.88   1.2225   0.493 0.624094    
PROTECTIONR:JOUR3   1.7602 -> 0.96   1.1407   1.543 0.129117    
PROTECTIONT:JOUR3   3.5884 -> 0.99   1.1042   3.250 0.002068 ** 

Indeed, the estimated mean are too high, especially looking at my actual data.

However, If I calculate the estimated means as EXP(estimate)/(1+EXP(estimate)), without the intercept, I get :


                           Estimated
                             mean
PROTECTIONR:JOUR2   1.7200 -> 0.84   1.2066   1.426 0.160223    
PROTECTIONT:JOUR2   2.4595 -> 0.92   1.1220   2.192 0.033061 *  
PROTECTIONG:JOUR3   0.6028 -> 0.65   1.2225   0.493 0.624094    
PROTECTIONR:JOUR3   1.7602 -> 0.85   1.1407   1.543 0.129117    
PROTECTIONT:JOUR3   3.5884 -> 0.97   1.1042   3.250 0.002068 ** 

Which is more probable but still not very consistent with my actual data :

Boxplot of the response variable, wrapped by PROTECTION , x = JOUR]1

My questions are :

-How to calculated the estimated mean of the interactions from the log odds ?

-Why do all the interaction including PROTECTIONF and JOUR1 do not appear in the summary ?

-Is there a way to perform pairwise comparison on a gamlss output ? Like emmeans and mulcomp ?

Thank you in advance Best regards

$\endgroup$

1 Answer 1

1
$\begingroup$

Let's tackle your second/title question first: the main and interaction effects for PROTECTIONF and JOUR1 are not shown because they are the reference levels, that is, the overall intercept (for PROTECTIONF at JOUR1) and JOUR2-3 (for PROTECTIONF) predict the response at that level. Doing otherwise would overparametrize the model, meaning it would no longer have a unique solution.

You were on the right track for recovering the predicted probability from the parameters through the expit $\pi=1/(1+e^{-\mu})$, but it is important to recall that your interaction parameter expresses only the additive difference in log-odds from the reference level(s). If you want to recover the prediction for PROTECTIONR at JOUR2 you need to include more than just that parameter, namely the overall intercept (prediction for PROTECTIONF at JOUR1), the PROTECTIONR parameter (difference between PROTECTIONF and PROTECTIONR at JOUR1), the JOUR2 effect (difference in PROTECTIONF between JOUR1 and JOUR2), and finally the interaction (difference in PROTECTIONR between JOUR1 and JOUR2).

From the above it's clear that a main effect certainly still has an interpretation in presence of an interaction: it's simply the prediction at the reference level(s) of that interaction. For example, PROTECTIONG is the (difference vs. PROTECTIONF in) log-odds at reference level JOUR1.

Plugging in the numbers you'll get some pretty funky predictions: the predicted response for PROTECTIONF at JOUR1 is $expit(1.4235)\approx 0.806$ whereas at JOUR2 it is $expit(1.4235-1.8951)\approx 0.384$, which clearly does not match the descriptive plot. The reason for that can be found in your choice of distribution though: you are using a beta inflated distribution, separately modelling the probability of response = 0 using JOUR. So, the beta model only sees responses that are >0, i.e. (I'm guessing) only the top two points for PROTECTIONF at JOUR1. This is much more consistent with what the model is predicting.

Finally, I've never done so myself but it should be possible to do (multiplicity-corrected) pairwise comparisons via ghlt() in multcomp, see also here.

$\endgroup$
1
  • $\begingroup$ Thank,you so much ! It is very helpful. I get an error message "dimensions of coefficients and covariance matrix don't match" when I run glht(Mymodel, linfct = c("PROTECTION= 0","JOUR=0", "PROTECTION:JOUR =0")). I do not understand this error message, I think I lack understanding of the matrix construction. But it would be enough for me to be able to change the intercept. However, the code '''data$PROTECTION <- relevel(data$PROTECTION, ref= "F") ''' does not change anything, do you know how to fix it ? $\endgroup$ Commented Oct 30, 2023 at 11:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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