3
$\begingroup$

Background

I investigated if sex differences over time (gender:time) in treatment response (PPA) were dependent on smoking status (gender:time:smoking_status), accounting for correlations between repeated measurements within subjects (1|ID) and between subjects within countries (1|country) using linear mixed model analyses (lme4/lmerTest). PPA runs from 0 to 100. Gender is "male" (reference) or "female". Smoking status is "never" (reference), "current", or "past". Time is categorical (0, 0.5, 1, and 2 years) with baseline as reference. The data is organized such that each patient's ID is recorded 4 times for each time point.

Model

library(lme4)     
library(lmerTest)  
mixed_smoking_interaction = lmer(pga ~ 1 + gender + time + smoking_status +
                                   gender*time + time*smoking_status + gender*smoking_status + 
                                   gender*time*smoking_status + (1|ID) + (1|country), data = dat, REML = F, control=lmerControl(optimizer="bobyqa"))

Output

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: ppa ~ 1 + gender + time + smoking_status + gender * time + time *  
    smoking_status + gender * smoking_status + gender * time *      smoking_status + (1 | ID) + (1 | country)
   Data: dat
Control: lmerControl(optimizer = "bobyqa")

     AIC      BIC   logLik deviance df.resid 
  316997   317226  -158472   316943    35064 

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-4.260 -0.617 -0.117  0.582  3.605 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 214.6    14.65   
 country  (Intercept)  14.8     3.84   
 Residual             343.3    18.53   
Number of obs: 35091, groups:  ID, 12424; country, 13

Fixed effects:
                                            Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)                                   57.016      1.198    14.736   47.61  < 2e-16 ***
genderFemale                                   3.498      0.633 29240.607    5.53  3.3e-08 ***
time0.5                                      -32.421      0.469 23415.545  -69.17  < 2e-16 ***
time1                                        -34.523      0.513 24947.778  -67.25  < 2e-16 ***
time2                                        -35.219      0.600 25694.302  -58.65  < 2e-16 ***
smoking_statuscurrent                          4.845      0.697 28780.332    6.95  3.6e-12 ***
smoking_statuspast                             1.039      0.720 28992.122    1.44  0.14890    
genderFemale:time0.5                           5.565      0.714 23402.390    7.79  6.8e-15 ***
genderFemale:time1                             5.002      0.803 25356.535    6.23  4.8e-10 ***
genderFemale:time2                             4.145      0.965 26460.845    4.29  1.8e-05 ***
time0.5:smoking_statuscurrent                  1.592      0.780 23361.254    2.04  0.04113 *  
time1:smoking_statuscurrent                    0.200      0.848 24803.146    0.24  0.81331    
time2:smoking_statuscurrent                    0.645      0.954 25627.537    0.68  0.49940    
time0.5:smoking_statuspast                     3.438      0.808 23371.499    4.25  2.1e-05 ***
time1:smoking_statuspast                       4.101      0.888 25033.836    4.62  3.9e-06 ***
time2:smoking_statuspast                       4.159      1.080 26146.881    3.85  0.00012 ***
genderFemale:smoking_statuscurrent             0.572      1.157 28977.151    0.49  0.62131    
genderFemale:smoking_statuspast                0.922      1.144 29036.548    0.81  0.42005    
genderFemale:time0.5:smoking_statuscurrent    -1.032      1.314 23286.470   -0.79  0.43209    
genderFemale:time1:smoking_statuscurrent       2.543      1.471 25213.547    1.73  0.08383 .  
genderFemale:time2:smoking_statuscurrent       2.245      1.737 26206.302    1.29  0.19616    
genderFemale:time0.5:smoking_statuspast       -2.575      1.294 23362.436   -1.99  0.04666 *  
genderFemale:time1:smoking_statuspast          0.309      1.477 25528.520    0.21  0.83441    
genderFemale:time2:smoking_statuspast         -0.536      1.897 26747.673   -0.28  0.77757

The output of the model states that genderFemale:time0.5:smoking_statuspast is significantly different (p=0.047). I always assumed that this meant that the sex difference at t=0.5 (Male - Female) in past smokers was significantly different compared never smokers at t=0.5 and that the -2.575 indicated the mean difference in units. However, I suspect that my interpretation is incorrect because when I tested the marginal means, I discovered that it was not smoking status past, but smoking status current, that had a significantly different mean sex difference (M - F) at 1 year, and I am struggling to understand the difference in these results.

Marginal means

gender time smoking_status emmean   SE   df lower.CL upper.CL
 Male   0    never            57.0 1.20 14.7     54.5     59.6
 Female 0    never            60.5 1.22 15.9     57.9     63.1
 Male   0.5  never            24.6 1.19 14.5     22.0     27.1
 Female 0.5  never            33.7 1.22 15.8     31.1     36.2
 Male   1    never            22.5 1.21 15.3     19.9     25.1
 Female 1    never            31.0 1.25 17.6     28.4     33.6
 Male   2    never            21.8 1.24 17.2     19.2     24.4
 Female 2    never            29.4 1.32 22.0     26.7     32.2
 Male   0    current          61.9 1.25 17.7     59.2     64.5
 Female 0    current          65.9 1.38 26.1     63.1     68.8
 Male   0.5  current          31.0 1.25 17.7     28.4     33.7
 Female 0.5  current          39.6 1.38 26.3     36.8     42.5
 Male   1    current          27.5 1.28 19.1     24.9     30.2
 Female 1    current          39.1 1.46 32.5     36.2     42.1
 Male   2    current          27.3 1.31 21.3     24.6     30.0
 Female 2    current          37.8 1.61 48.4     34.5     41.0
 Male   0    past             58.0 1.27 18.8     55.4     60.7
 Female 0    past             62.5 1.36 24.7     59.7     65.3
 Male   0.5  past             29.1 1.27 18.6     26.4     31.7
 Female 0.5  past             36.5 1.36 24.6     33.7     39.3
 Male   1    past             27.6 1.30 20.5     24.9     30.3
 Female 1    past             37.4 1.45 31.9     34.4     40.3
 Male   2    past             27.0 1.40 27.7     24.1     29.9
 Female 2    past             35.0 1.72 62.0     31.6     38.5

The plot of the three-way interaction enter image description here

Code used to test if the marginal means (Male - Female) for current and past smokers differ from never smokers for each time point.

library(emmeans)
emm_model1  <- emmeans(mixed_smoking_interaction, ~gender*time*smoking_status)

#Marginal estimated means
emm_model1

#Create a matrix to be used for a custom contrast later

#Sex differences (M - F) in never smokers
Diff_S0_t0 <- rep(0, 24)
Diff_S0_t0[1] <- 1
Diff_S0_t0[2] <- -1

Diff_S0_t0.5 <- rep(0, 24)
Diff_S0_t0.5[3] <- 1
Diff_S0_t0.5[4] <- -1

Diff_S0_t1 <- rep(0, 24)
Diff_S0_t1[5] <- 1
Diff_S0_t1[6] <- -1

Diff_S0_t2 <- rep(0, 24)
Diff_S0_t2[7] <- 1
Diff_S0_t2[8] <- -1

#Sex differences (M - F) in current smokers
Diff_S1_t0 <- rep(0, 24)
Diff_S1_t0[9] <- 1
Diff_S1_t0[10] <- -1

Diff_S1_t0.5 <- rep(0, 24)
Diff_S1_t0.5[11] <- 1
Diff_S1_t0.5[12] <- -1

Diff_S1_t1 <- rep(0, 24)
Diff_S1_t1[13] <- 1
Diff_S1_t1[14] <- -1

Diff_S1_t2 <- rep(0, 24)
Diff_S1_t2[15] <- 1
Diff_S1_t2[16] <- -1

#Sex differences (M - F) in past smokers
Diff_S2_t0 <- rep(0, 24)
Diff_S2_t0[17] <- 1
Diff_S2_t0[18] <- -1

Diff_S2_t0.5 <- rep(0, 24)
Diff_S2_t0.5[19] <- 1
Diff_S2_t0.5[20] <- -1

Diff_S2_t1 <- rep(0, 24)
Diff_S2_t1[21] <- 1
Diff_S2_t1[22] <- -1

Diff_S2_t2 <- rep(0, 24)
Diff_S2_t2[23] <- 1
Diff_S2_t2[24] <- -1

#Do the sex differences in current smokers and past smokers differ from patients with never smokers, separated for every time point?
contrast(emm_model1, method = list("T0_Diff_S0 - Diff_S1" = Diff_S0_t0 - Diff_S1_t0,
                                   "T0.5_Diff_S0 - Diff_S1" = Diff_S0_t0.5 - Diff_S1_t0.5,
                                   "T1_Diff_S0 - Diff_S1" = Diff_S0_t1 - Diff_S1_t1,
                                   "T2_Diff_S0 - Diff_S1" = Diff_S0_t2 - Diff_S1_t2,
                                   "T0_Diff_S0 - Diff_S2" = Diff_S0_t0 - Diff_S2_t0,
                                   "T0.5_Diff_S0 - Diff_S2" = Diff_S0_t0.5 - Diff_S2_t0.5,
                                   "T1_Diff_S0 - Diff_S2" = Diff_S0_t1 - Diff_S2_t1,
                                   "T2_Diff_S0 - Diff_S2" = Diff_S0_t2 - Diff_S2_t2))

Final Output

 contrast               estimate   SE    df t.ratio p.value
 T0_Diff_S0 - Diff_S1      0.572 1.16 28977   0.494  0.6213
 T0.5_Diff_S0 - Diff_S1   -0.461 1.15 28847  -0.399  0.6897
 T1_Diff_S0 - Diff_S1      3.115 1.32 32821   2.366  0.0180
 T2_Diff_S0 - Diff_S1      2.816 1.61 35073   1.754  0.0794
 T0_Diff_S0 - Diff_S2      0.922 1.14 29037   0.806  0.4200
 T0.5_Diff_S0 - Diff_S2   -1.653 1.13 28616  -1.460  0.1444
 T1_Diff_S0 - Diff_S2      1.231 1.32 33185   0.930  0.3524
 T2_Diff_S0 - Diff_S2      0.387 1.78 34849   0.217  0.8279

In this output, the capital T indicates the time point, Diff indicates the sex difference (male - female), and S indicates the smoking status (S0 = never smokers, S1 = current smokers, and S2 = past smokers). From this output, it can be deduced that the mean sex difference at one year between current smokers and never smokers is 3.1 units, statistically different from 0 (p=0.018).

Question: I do not understand how the output of the marginal means tells us a different story (current smoker t=1 significant) from the output of the linear mixed model (past smoker t=0.5 significant) and would appreciate help unravelling this.

EDIT: Adding three plots, status x sex by time, time x status by gender, and time x gender by status.

Status x sex by time

Time x Status by gender [Time x Gender by status 4

Improve this question
$\endgroup$

2 Answers 2

Reset to default
2
$\begingroup$

You have not actually tested the 3-way interaction as you claim. Rather, you looked at the 8 separate coefficients, chose the one that was most significant, and compared it to the 0.05 cutoff which ignores the issue of multiplicity.

To formally test the presence of a 3-way interaction, I would suggest a likelihood ratio test from library(lmtest).

mod1 = lmer(pga ~ 1 + gender + time + smoking_status +
                                   gender*time + time*smoking_status + gender*smoking_status + 
                                   gender*time*smoking_status + (1|ID) + (1|country), data = dat, REML = F, control=lmerControl(optimizer="bobyqa")) 
mod2 = lmer(pga ~ 1 + gender + time + smoking_status +
                                   gender*time + time*smoking_status + gender*smoking_status + 
                                   + (1|ID) + (1|country), data = dat, REML = F, control=lmerControl(optimizer="bobyqa")) 

lrtest(mod1, mod2)
Improve this answer
$\endgroup$
6
  • $\begingroup$ Thanks for the response. I know there are different ways to test if a three-way interaction is significant, such as the Anova type 3 test and the one you described. We aimed to further investigate any three-way interaction with a significant regression coefficient (e.g., is it clinically relevant?). Could you please instead answer the question? $\endgroup$
    – Pashtun
    Commented Mar 7, 2023 at 18:16
  • $\begingroup$ @Pashtun since we don't have a reproducible example, guesswork is completely out of the question. HOWEVER it should be noted that inference from a fixed effects model applies to the conditional mean, whereas emmeans models the marginal means. They are not going to give the same result. Try replacing the mixed effects model with an OLS. $\endgroup$
    – AdamO
    Commented Mar 8, 2023 at 17:28
  • $\begingroup$ Thanks for your response. Cant share the data unfortunately because of patiënt confidentiality. Perhaps you could write out steps for me to undertake to investigate the issue? Or ask what other information you would need? $\endgroup$
    – Pashtun
    Commented Mar 8, 2023 at 17:53
  • $\begingroup$ Did you mean the function lrtest? $\endgroup$
    – COOLSerdash
    Commented Mar 9, 2023 at 16:32
  • 1
    $\begingroup$ Thanks for confirming. I made the edit. Feel free to change if necessary. $\endgroup$
    – COOLSerdash
    Commented Mar 9, 2023 at 17:11
1
$\begingroup$

A 3-way interaction indicates that any of the three composite 2-way interactions differ. As you said, Sex x Status could differ by Time. Or, Sex x Time could differ by Status, or Time x Status could differ by Sex. Note that the effect you found with emmeans is nominally significant in your original regression. More generally, do you have a significant interaction? A single p=0.046 out of 6 tests is not very convincing... have you tested whether adding 3-way interactions significantly improves model fit over a model without them?

Improve this answer
$\endgroup$
5
  • $\begingroup$ Two way interactions show that smoking:time and gender:time are significant but not gender:smoking. How to interpret? What do you mean with the effect you found in emmeans is nominally significant in the original regression? Yes, you are right a single significant p-value is not convincing but we decided to report any model with any significant p-value as opposed to doing an overall Anova type 3 test for example. I did not check the models fit as this is an association model and not prediction model. The plan was made a priori in study protocol to analyze this model regardless of fit. $\endgroup$
    – Pashtun
    Commented Mar 1, 2023 at 18:21
  • 2
    $\begingroup$ @Pashtun - 2-way interactions do not need to be significant for a 3-way interaction to be significant. I listed 3 possible effects driving what you're seeing, so far you've only explored the first. Regarding my comment about the nominal effect, in your results table it reads: genderFemale:time1:smoking_statuscurrent 2.543 1.471 25213.547 1.73 0.08383 . Is this not the interaction driving what you're seeing with emmeans - a nominal effect? $\endgroup$
    – David B
    Commented Mar 1, 2023 at 18:31
  • $\begingroup$ Thanks for the explanation. Yes now that you mention this, it does make sense that what I found could be the nominal effect in the results of the lmer(). I will test the other two scenarios and report the results back tomorrow. $\endgroup$
    – Pashtun
    Commented Mar 1, 2023 at 19:18
  • $\begingroup$ I added the three graphs that you requested. Please explain how this can help with the interpretation. $\endgroup$
    – Pashtun
    Commented Mar 2, 2023 at 20:44
  • $\begingroup$ Not saying that you necessarily need to provide the plots, but that you could repeat your emmeans procedure the other two ways that I described. $\endgroup$
    – David B
    Commented Mar 2, 2023 at 20:48

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.