1
$\begingroup$

Context

The purpose of this analysis is to examine sex differences in Bath Ankylosing Spondylitis Functional Index (BASFI) scores in Axial spondyloarthritis (axSpA) patients who are starting their first Tumor Necrosis Factor inhibitor (TNFi) treatment over a two-year follow-up period. The BASFI score consists of ten domains, so a secondary research question is to determine if there are sex-based differences in these domains over time.

Population: axSpA patients initiating first TNFi

Dependent variable: BASFI-score (scale 0-100, higher is worse)

Independent variables of interest: sex, sex:time, and sex:time:domains

Data-structure

The data is organized such that each patient's ID is recorded 40 times (10 times for each domain at four different time points). A 4-level structure analysis was conducted to account for correlations between components of the composite score, repeated measurements within subjects, and subjects within countries. Time is coded as a numeric variable with one unit representing one year.

Here is a peek inside the dataset to further clarify how the data is structured:

##    ID time domain gender basfi 
## 1   1  0.0      1   Male    NA     
## 2   1  0.0      2   Male    NA     
## 3   1  0.0      3   Male    NA     
## 4   1  0.0      4   Male    NA     
## 5   1  0.0      5   Male    NA     
## 6   1  0.0      6   Male    NA     
## 7   1  0.0      7   Male    NA     
## 8   1  0.0      8   Male    NA     
## 9   1  0.0      9   Male    NA     
## 10  1  0.0     10   Male    NA     
## 11  1  0.5      1   Male    NA     
## 12  1  0.5      2   Male    NA     
## 13  1  0.5      3   Male    NA     
## 14  1  0.5      4   Male    NA     
## 15  1  0.5      5   Male    NA     
## 16  1  0.5      6   Male    NA     
## 17  1  0.5      7   Male    NA     
## 18  1  0.5      8   Male    NA     
## 19  1  0.5      9   Male    NA     
## 20  1  0.5     10   Male    NA     
## 21  1  1.0      1   Male     0     
## 22  1  1.0      2   Male     0     
## 23  1  1.0      3   Male     0     
## 24  1  1.0      4   Male     0     
## 25  1  1.0      5   Male     0     
## 26  1  1.0      6   Male     0     
## 27  1  1.0      7   Male     0     
## 28  1  1.0      8   Male     0     
## 29  1  1.0      9   Male     0     
## 30  1  1.0     10   Male     0     

Structure of the covariates:

## 'data.frame':    391200 obs. of  5 variables:
##  $ country                   : Factor w/ 12 levels "CH","CZ","DK",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ gender                    : Factor w/ 2 levels "Male","Female": 1 1 1 1 1 1 1 1 1 1 ...
##  $ ID                        : Factor w/ 9780 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ basfi                     : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ time                      : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ domain                    : Factor w/ 10 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...

Analysis

This is the code used in R.

mixed_final_model = lmer(basfi ~ 1 + gender + time + domain + gender*time + gender*domain + time*domain + gender*time*domain + (1 | ID) + (1| time) + (1| country), data = dat, REML = F, control=lmerControl(optimizer="bobyqa"))
summary(mixed_final_model)

Here is the output:

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

     AIC      BIC   logLik deviance df.resid 
 1514557  1514999  -757235  1514469   168660 

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-4.861 -0.603 -0.086  0.507  4.604 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 428.1    20.69   
 country  (Intercept)  85.7     9.26   
 time     (Intercept)  36.8     6.06   
 Residual             412.1    20.30   
Number of obs: 168704, groups:  ID, 5848; country, 9; time, 4

Fixed effects:
                             Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)                    32.550      5.645      7.033    5.77  0.00068 ***
genderFemale                   -1.821      0.730  16155.170   -2.49  0.01270 *  
time                           -8.015      4.109      3.763   -1.95  0.12732    
domain2                         9.518      0.425 162836.923   22.41  < 2e-16 ***
domain3                        -3.180      0.425 162837.375   -7.49  7.2e-14 ***
domain4                        -0.186      0.425 162836.848   -0.44  0.66204    
domain5                         9.448      0.425 162837.060   22.24  < 2e-16 ***
domain6                         1.170      0.425 162837.119    2.75  0.00589 ** 
domain7                        -1.688      0.425 162837.135   -3.97  7.1e-05 ***
domain8                        10.697      0.425 162837.384   25.18  < 2e-16 ***
domain9                        14.395      0.425 162837.350   33.87  < 2e-16 ***
domain10                       11.434      0.425 162837.929   26.90  < 2e-16 ***
genderFemale:time               1.937      0.466 163245.044    4.16  3.2e-05 ***
genderFemale:domain2           -1.358      0.670 162837.475   -2.03  0.04264 *  
genderFemale:domain3            7.045      0.670 162837.750   10.52  < 2e-16 ***
genderFemale:domain4            5.786      0.670 162837.483    8.64  < 2e-16 ***
genderFemale:domain5            7.058      0.670 162837.703   10.54  < 2e-16 ***
genderFemale:domain6            8.085      0.670 162837.665   12.07  < 2e-16 ***
genderFemale:domain7           12.346      0.670 162837.622   18.43  < 2e-16 ***
genderFemale:domain8           -2.757      0.670 162837.749   -4.12  3.9e-05 ***
genderFemale:domain9            6.544      0.670 162837.883    9.77  < 2e-16 ***
genderFemale:domain10          10.741      0.670 162838.300   16.03  < 2e-16 ***
time:domain2                   -1.474      0.401 162837.098   -3.68  0.00024 ***
time:domain3                    1.856      0.401 162837.336    4.63  3.7e-06 ***
time:domain4                   -0.609      0.401 162837.018   -1.52  0.12862    
time:domain5                   -1.951      0.401 162837.307   -4.87  1.1e-06 ***
time:domain6                    1.694      0.401 162837.123    4.23  2.4e-05 ***
time:domain7                    0.723      0.401 162837.384    1.80  0.07142 .  
time:domain8                    1.186      0.401 162837.388    2.96  0.00309 ** 
time:domain9                   -2.589      0.401 162837.558   -6.46  1.1e-10 ***
time:domain10                  -1.231      0.401 162837.887   -3.07  0.00216 ** 
genderFemale:time:domain2      -0.371      0.651 162837.767   -0.57  0.56875    
genderFemale:time:domain3      -1.736      0.651 162837.755   -2.67  0.00764 ** 
genderFemale:time:domain4      -1.221      0.651 162837.600   -1.88  0.06068 .  
genderFemale:time:domain5      -1.016      0.651 162837.803   -1.56  0.11846    
genderFemale:time:domain6      -1.816      0.651 162837.583   -2.79  0.00524 ** 
genderFemale:time:domain7      -2.134      0.651 162837.761   -3.28  0.00104 ** 
genderFemale:time:domain8       0.242      0.651 162837.661    0.37  0.71010    
genderFemale:time:domain9      -0.574      0.651 162837.844   -0.88  0.37758    
genderFemale:time:domain10     -1.509      0.651 162838.456   -2.32  0.02053 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

My interpretation so far: according to the 4-level structure analysis, there is a significant but small main effect of sex on the BASFI scores (-1.8, p-value = 0.013). This indicates that females have lower overall BASFI scores than males. The effect of gender:time is also significant, strong, and positive (1.93, p-value < 0.001). This suggests that the impact of gender on the BASFI scores changes over time. Specifically, females tend to have lower mean BASFI scores than males at earlier time points, but the difference between males and females becomes smaller or even reverses at later time points (after 2+ years).

Question: The effect of gender:time:domain is significant for domains 3, 6, 7, and 10, with effect sizes ranging from -1.5 to -2.1. All of these effects have a negative direction, with domain 1 serving as the reference category. This result suggests that the magnitude of the sex differences in these domains changes over time. Based on this information, is it possible to determine if the magnitude of the sex differences in BASFI over time increases in domains 3,6,7, and 10 compared to domain 1? Or does it decrease? Is another interpretation more appropriate?

EDIT: The purpose of coding time as a numeric variable in our model is to examine the potential linear relationship between time and the outcome variable (BASFI). By considering time as a continuous variable, we aim to investigate if any sex differences in the BASFI exist over the entire 24-month follow-up period, rather than being limited to specific time points. Additionally, we included interaction terms between time and gender in the model to explore whether the effect of gender on the BASFI changes over time. To account for the correlation between the components/domains of the BASFI over time, time was also added with random intercept.

$\endgroup$
4
  • $\begingroup$ Please edit the question to say why you are modeling time both as a fixed effect (including some fixed-effect interaction terms) and with its own random intercept in (1|time), particularly when it seems only to have 4 distinct values. It also looks like the fixed effect is modeled as a linear numeric value (only 1 coefficient reported), while the random effect necessarily treats it as categorical. That will make it difficult to interpret any "differences over time." Please provide that information by editing the question, as comments are easy to overlook and can be deleted. $\endgroup$
    – EdM
    Jan 10, 2023 at 15:14
  • $\begingroup$ @EdM question has been edited as requested $\endgroup$
    – Pashtun
    Jan 10, 2023 at 19:10
  • $\begingroup$ I added what I think are the correct meanings of the abbreviations you used; please correct if I was wrong. Also, please explain the way that the basfi values and the domain values are coded. From the structure of the model it seems that each of the outcome basfi values is for just one of the 10 domain scores, with the domain specified as a categorical predictor in the data row. But the text of the question makes it seem that the outcome values might instead be the overall BASFI score. Please edit the question to clarify. $\endgroup$
    – EdM
    Jan 12, 2023 at 3:27
  • $\begingroup$ @EdM Thank you for your response and for editing the question for further clarification of the abbreviations. Each of the basfi outcomes is indeed just one of the 10 domain scores with the domain specified as a factor. I am not sure how the text of the question implies something else, my apologies for that. I will add the first 10 rows of the dataset and data structure for further clarification. $\endgroup$
    – Pashtun
    Jan 12, 2023 at 10:23

1 Answer 1

1
$\begingroup$

First, as the breakdown by domain is a secondary goal, I'm going to assume that you have a separate model for the total BASFI score as a function of time and gender. You should use that as a guide for how deeply to pursue the interactions you are evaluating in the model you present. With a large study you often will find "statistically significant" results that are of limited practical importance. Minor "significant" interactions need to be evaluated against how much influence they will have on the total BASFI score. The magnitudes of 3-way interactions in your model are typically less than 2 units, versus an Intercept of 32.55 units (the estimate for Male gender, domain1, time = 0).

Interactions

When predictors are involved in interaction terms, the lower-level coefficients for those predictors in the default R coding are for the situation where the interacting predictors are at a level of 0 (continuous interactors) or at the reference level (categorical interactors). That's a continuing source of confusion. Thus your interpretation of coefficients isn't correct.

For gender and time, both involved in interactions with domain and with each other, the individual coefficients only represent their associations with outcome for domain1, and when time = 0 (coded continuous) for the gender coefficient and gender is Male (reference category) for the time coefficient. The "significant but small main effect of sex (-1.8, p-value = 0.013)" is only for domain1. The apparent lack of a "significant main effect" of time is only for Male gender and domain1. The coefficients for the domain values are their differences from domain1 for Male gender at time = 0. The genderFemale:time interaction coefficient is only for domain1. The time:domain interaction coefficients are only for Male gender. The gender:domain interaction coefficients are only for time = 0.

To evaluate the associations of any predictor with outcome, when it's involved in an interaction, you need to do a combined analysis on all coefficients involving it. The Anova() (capital "A") function in the car package can do that via its default Type II tests.

With interactions, it's best to evaluate model estimates for particular combinations of predictor values rather than focusing on individual coefficients as you seem to be doing. That's the best way to answer the specific questions you have in mind. The emmeans package provides helpful tools for that, for predictions from fixed effects.

Modeling time

You might reconsider how you are handling time. In this case with a discrete-valued time variable (4 numeric levels), it's possible to include time both as a linear fixed effect and as a random effect, as you are doing. See this explanation of what such dual coding does. The fixed-effect component gives the linear association of the values with outcome, as you want, while the random intercept for time tries to account for differences from that linear trend. It's not clear how successful that will be with only 4 time values, and you would have to go an extra step to add in the correct random-intercept value for each time value to the fixed-effect predictions. It's also not clear (to me, at least) how successfully that random intercept for time will "account for the correlation between the components/domains of the BASFI over time."

Your current approach assumes linearity for all interactions of gender and domain with time. That can be a pretty strict assumption, and I wonder how well that will work here, where I suspect that there are rapid initial improvements that end up being maintained at late times. The random time intercept doesn't account for non-linearity of time in the interactions.

With only 4 time points and such a large study, you could model time as a 4-level factor including the interactions, and then evaluate the linearity in time after the initial modeling, based on the coefficient estimates. Alternatively, you could keep your linear time numeric coding and augment it with indicator variables for just 2 of your specific time values. Tests on the "significance" of the combined 2 indicator variables then is a test for non-linearity in time. Either method would allow separate evaluation of linearity in time on its own and in its interactions involving gender and domain.

Possibly useful for your study, you could recode time into 3 indicator variables to show whether time has exceeded each of time = 0, time = 0.5, and time = 1. Then each regression coefficient would represent the change from the prior time value. See Section 2.7.3 of Harrell's Regression Modeling Strategies about these different ways of handling ordered predictors with small numbers of values.

$\endgroup$
6
  • $\begingroup$ Thank you for your detailed response. I appreciate the insights provided. In regards to the evaluation of the overall" score, I utilized a three-level model in which there were four observations per patient - one for each time point. The observations were the overall scores, and time was not included as a random effect in this analysis. I only included the four-level structure in the analysis of the different domains. I understand the value of using the car package in analyzing interactions, and I can see its importance in addition to analyzing individual coefficients. $\endgroup$
    – Pashtun
    Jan 13, 2023 at 13:37
  • $\begingroup$ By using the emmeans package, I can estimate the marginal means for males and females over time for the different domains if I code time as a factor. I can then plot the results to see if there is a clear relationship. Is this the approach you were suggesting? I have done this before but I was wondering if this is the best way to address this issue. If you had something else in mind, please let me know. Additionally, I would like to thank you for further clarifying how to interpret the main effect of sex, which only concerned domain 1 in this case. This was a valuable insight for me. $\endgroup$
    – Pashtun
    Jan 13, 2023 at 13:47
  • $\begingroup$ "You are correct that the relationship between the outcome and time is not linear. Patients improve significantly in the first six months but progress slows down at month 12 and 24. I have also noticed significant heteroscedascity in the model, especially around the median. This could be due to the non-linear relationship of the BASFI with time. To address this, you suggest adding time as a factor. However, how can I then answer my research question if the difference is maintained over the entire period? Could you please elaborate further on this? $\endgroup$
    – Pashtun
    Jan 13, 2023 at 13:53
  • 1
    $\begingroup$ @Pashtun treating time as a factor and then plotting marginal means with standard errors against time is what I had in mind. That removes the troublesome linearity-in-time assumption. The consec comparisons in emmeans, for each time point against the prior one, might be of most interest: the first comparison shows initial response, the subsequent two how well that is maintained. Or try the trt.vs.ctrl comparisons in emmeans, with time=0 as the control, for response relative to initial value at each subsequent time. $\endgroup$
    – EdM
    Jan 13, 2023 at 14:54
  • 1
    $\begingroup$ @Pashtun with the size of your study I don't think you have to worry about loss of power or overfitting due to the increased number of coefficient estimates. That will make display of individual coefficient estimates even messier, but you can summarize with an Anova() table and plot predictions. If you treat time as a 4-level categorical predictor as I recommend, however, you can't also include it as a random effect. See this explanation of why that dual coding only works for a numeric treatment of a small number of time values. $\endgroup$
    – EdM
    Jan 14, 2023 at 20:47

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.