Added interaction terms to logistic regression - how to interpret from here?

Question

I'm using logistic regression in r to examine correlates of experiencing headache (yes(1)/no (0)) when taking a medication. My correlates (all categorical) include:

• number of tablets of the medication (tabs)
• year
• gender
• age
• health district (HD)
• location (home, work, outside)

I have two interaction terms that seem to be significant - Health district (HD) x year, and HD x # of tablets

I ran the model with the interaction terms (see below), and I'm just not sure where to go from here. I know that for the variables involved in the interaction, I can't interpret the ORs as I normally would. Instead, I plotted the tabs by health district and year.. is that enough? When showing my regression table, should I show all of the ORs, but only interpret the ones not involved in interaction?

Thanks in advance for any advice.

Here are the regression results with coefficients exponentiated:

                               exp(Est.) 2.5%   97.5%    P
(Intercept)                         0.20 0.01    4.17 0.30
tabs2                               2.15 0.36   12.82 0.40
tabs3                               0.70 0.05    9.37 0.79
tabs>3                              6.66 0.70   63.57 0.10
year2016                            0.17 0.01    4.57 0.29
year2017                            0.08 0.00    1.35 0.08
year2018                            0.06 0.00    1.18 0.06
genderFemale                        0.80 0.57    1.11 0.19
genderOther                         0.00 0.00     Inf 0.99
genderUnknown                       1.16 0.58    2.32 0.67
age19 - 30                          2.88 0.65   12.73 0.16
age31 - 60                          3.22 0.73   14.10 0.12
ageOver 60                          1.57 0.23   10.71 0.65
ageUnknown                          2.48 0.53   11.67 0.25
HDVIHA                              0.11 0.00    3.11 0.19
HDFHA                               0.42 0.03    6.81 0.54
HDIHA                               0.49 0.03    8.02 0.61
HDVCH                               0.11 0.01    2.12 0.14
HDUnknown                           2.36 0.30   18.52 0.41
locationhome                        0.63 0.45    0.90 0.01
locationwork                        0.67 0.36    1.26 0.21
locationOther                       0.57 0.39    0.83 0.00
tabs2:HDVIHA                        0.52 0.05    5.19 0.58
tabs3:HDVIHA                        0.60 0.02   19.77 0.77
tabs>3:HDVIHA                       0.00 0.00     Inf 0.98
tabs2:HDFHA                         0.78 0.12    5.20 0.80
tabs3:HDFHA                         1.68 0.11   24.87 0.71
tabs>3:HDFHA                        0.27 0.02    3.21 0.30
tabs2:HDIHA                         0.44 0.06    3.40 0.43
tabs3:HDIHA                         6.12 0.38   97.95 0.20
tabs>3:HDIHA                        0.69 0.06    8.68 0.78
tabs2:HDVCH                         0.26 0.03    2.21 0.21
tabs3:HDVCH                         0.64 0.03   13.12 0.77
tabs>3:HDVCH                        0.45 0.03    6.41 0.55
tabs2:HDUnknown                     0.61 0.09    4.19 0.62
tabs3:HDUnknown                     3.36 0.23   49.51 0.38
tabs>3:HDUnknown                    0.40 0.04    4.52 0.46
year2016:HAVIHA                    51.05 0.88 2966.51 0.06
year2017:HDVIHA                    48.41 1.23 1913.35 0.04
year2018:HDVIHA                     0.00 0.00     Inf 0.99
year2016:HDFHA                      3.26 0.10  104.55 0.50
year2017:HDFHA                      5.80 0.30  111.77 0.24
year2018:HDFHA                      7.48 0.32  174.24 0.21
year2016:HDIHA                      4.34 0.13  148.79 0.42
year2017:HDIHA                      8.45 0.40  177.91 0.17
year2018:HDIHA                      6.23 0.25  153.68 0.26
year2016:HDVCH                     36.08 0.97 1336.63 0.05
year2017:HDVCH                     81.30 3.04 2173.58 0.01
year2018:HDVCH                     32.33 0.93 1123.52 0.05
year2016:HDUnknown                  0.00 0.00     Inf 0.98
year2017:HDUnknown                  2.35 0.48   11.49 0.29
year2018:HDUnknown                  1.00 1.00    1.00  NaN                                          

Answer 1

The usual rule of thumb to avoid overfitting in logistic regression is to have about 15 cases in the minority class per predictor you are evaluating, unless you are using some type of penalization like with ridge regression. That's about 17 total predictors including interactions for your data set with 255 headache cases. Your model, with 3 times that many predictors, is probably very badly overfit and the p-values are unreliable.

To cut down on the number of predictors, try to use variables like age and year as continuous rather than categorical predictors. Consider modeling the health districts HD as random effects in a mixed model instead of fixed effects.

Instead of setting aside separate "unknown" categories, use multiple imputation to get multiple complete data sets with the missing data estimated in a principled way. Besides cutting down several categories, that avoids bias resulting from the "missingness" pattern and lets you get estimates that take the variability both in the model and in the data imputation into account.

The situation you face is very common in biomedical studies. See Frank Harrell's course notes and book on Regression Modeling Strategies. Besides what Harrell has to say about data imputation, Stef van Buuren has a helpful and freely available book specifically on that topic.

It's worth your while to analyze your data well. Take advantage of this as an opportunity to learn more about this type of modeling, if possible collaborating with an experienced statistician.