Logistic regression: Interaction continuous/categorical variables and multiple comparisons

Question

I am fitting a logistic regression model with two independent variables, one continuous (length, here lun) and one categorical (Year = 2013, 2014, 2015). My dependent variable is the maturity stage of an individual (stage, 1= mature & 0= immature). The aim is to understand if the resulting logistic curves among the three years are statistically different.

Here the code:

confronto_mod <- glm(stage~lun*Year,data=confronto_macro,family=binomial) 
summary(confronto_mod)

    Call:
glm(formula = stage ~ lun + Year + lun:Year, family = binomial, 
    data = confronto_macro)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.8590  -0.3734  -0.0219   0.5665   2.3230  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -9.15307    2.13656  -4.284 1.84e-05 ***
lun            0.29703    0.06505   4.566 4.97e-06 ***
Year2014      -6.86345    2.95156  -2.325 0.020053 *  
Year2015     -11.58898    3.43635  -3.372 0.000745 ***
lun:Year2014   0.21789    0.09202   2.368 0.017887 *  
lun:Year2015   0.39219    0.10985   3.570 0.000357 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 957.81  on 695  degrees of freedom
Residual deviance: 459.78  on 690  degrees of freedom
AIC: 471.78

Number of Fisher Scoring iterations: 7

Now, as far I understand here we got the following.

• The increase of a unit in length in 2014 yields an increase in the odds ratio of a factor exp(0.21789)= 1.24345, with respect to the reference taken by R here of 2013; the same in 2015 with an increase in the odds ratio of a factor exp(0.39219)=1.480219. In 2013 a unit increase in length yields an increase of a factor exp(0.29703)= 1.345856
• So the three Years appear to yields a significantly different increase of the probability of being mature for every increase in the length unit, is this right?
• The test behind is the Wald test right?
• How can I compare the three Years to explore more these differences? Basically to understand if 2014 and 2015 are significantly different.

I tried to do the following but unsure if it is correct:

library('multcomp')    
glht(confronto_mod, linfct=c('lun:Year2014 = 0', 'lun:Year2015 = 0',
                                         'lun:Year2015 - lun:Year2014 = 0'))


    Simultaneous Tests for General Linear Hypotheses

    Fit: glm(formula = stage ~ lun * Year, family = binomial, data = confronto_macro)

    Linear Hypotheses:
                                     Estimate Std. Error z value Pr(>|z|)   
    lun:Year2014 == 0                 0.21789    0.09202   2.368  0.04674 * 
    lun:Year2015 == 0                 0.39219    0.10985   3.570  0.00103 **
    lun:Year2015 - lun:Year2014 == 0  0.17430    0.10987   1.586  0.25005   
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    (Adjusted p values reported -- single-step method)

Otherwise, it is sufficient to change for example the 2014 as a reference in the model and explore the resulting output?

Thank you so much for your feedbacks!

Answer 1

The interpretation you provided for your model is not correct.

Your fitted model is in effect a collection of 3 fitted sub-models. Each fitted sub-model relates the log odds of being mature rather than immature in a given year to length (lun).

Year 2013:

log odds of being mature = -9.15307 + 0.29703*lun

Year 2014:

log odds of being mature = (-9.15307 - 6.86345) + (0.29703 + 0.21789)*lun

Year 2015:

log odds of being mature = (-9.15307 - 11.58898) + (0.29703 + 0.39219)*lun

The first fitted sub-model indicates that, for the year 2013, a 1-unit increase in length is associated with an estimated 35% increase in the odds of being mature. Indeed, exp(0.29703) = 1.35 and (1.35-1)x100% = 35%. You can also state that, for the year 2013, the odds of being mature were estimated to increase by a multiplicative factor of 1.35 for each extra unit of length.

The second fitted sub-model indicates that, for the year 2014, a 1-unit increase in length is associated with a 67% increase in the odds of being mature. Indeed, exp(0.29703 + 0.21789) = 1.67 and (1.67-1)x100% = 67%. Another way of stating this is that, for the year 2014, the odds of being mature were estimated to increase by a multiplicative factor of 1.67 for each extra unit of length.

The third fitted sub-model indicates that, for the year 2015, a 1-unit increase in length is associated with a 99% increase in the odds of being mature. Indeed, exp(0.29703 + 0.39219) = 1.99 and (1.99 - 1)x100% = 99%. In other words, for the year 2015, the odds of being mature were estimated to increase by a multiplicative factor of 1.99 for each extra unit of length.

The multiplicative factors 1.35, 1.67 and 1.99 are referred to as odds ratios. They allow you to compare the odds of being mature in a given year between two lengths which differ by 1-unit.

Because these odds ratios are all positive, you can also conclude that length is positively associated with the probability of being mature on each of the three years (not just with the log odds of being mature). While length has a positive linear effect on the log odds of being mature for each year, it has a positive non-linear effect on the probability of being mature for that year.

To visualize the effect of length on the probability of being mature, you can use the effects package in R:

require(effects)
plot(allEffects(confronto_mod))

You are right that the p-values reported for each year's sub-model come from Wald tests.

The emmeans package will allow you to investigate your fitted model further. For example:

require(emmeans) 
emtrends(confronto_mod, ~ Year, var="lun")

should produce estimated values and confidence intervals for the effect of length (lun) on the log odds of being mature for each year. You can check that this is indeed the case by comparing the estimates provided by emtrends in the lun.trend column against the estimates 0.29703 for 2013, 0.29703 + 0.21789 for 2014 and 0.29703 + 0.39219 for 2015.

Furthermore, the command:

emtrends(confronto_mod, pairwise ~ Year, var="lun")

should allow you to compare the effect of length (lun) on the log odds of being mature between pairs of year. Check whether the values reported in the Estimate column for pairs of years are obtained as differences between the estimates 0.29703 for 2013, 0.29703 + 0.21789 for 2014 and 0.29703 + 0.39219 for 2015.

The reason these checks are important is because it is possible that the emtrends might perform a back-transformation of results behind the scenes (e.g., a transformation from the log odds scale to the probability scale). I haven't used this function much so please do your due dilligence to understand how it works. You can see here for a related example of emtrends but in a linear regression context: https://stats.idre.ucla.edu/r/seminars/interactions-r/#s4a.