How do I work out the significance of main effects in negative binomial models with more than two factors?

Question

I am trying to look at how the number of events X is affected by the three factors A (4 levels), B (2 levels) and C (2 levels) using a negative binomial model as follows:

model.count.nb<-glm.nb(X.total ~ A*B*C, link = "log", data = count.X)
summary(model.count.nb)
Call:
glm.nb(formula = X.total ~ A * B * C, 
    data = count.X, link = "log", init.theta = 4.864199324)
Deviance Residuals: 
    Min       1Q   Median       3Q      Max

-3.7430  -0.5601   0.0705   0.4939   2.0461
Coefficients:
                                      Estimate Std. Error z value Pr(>|z|)

(Intercept)                           5.529870   0.152589  36.240   <2e-16 ***
A2                                   -0.011079   0.193019  -0.057    0.954

A3                                    0.060288   0.215733   0.279    0.780

A4                                   -0.270010   0.206022  -1.311    0.190

B2                                    0.098550   0.215697   0.457    0.648

C2                                   -0.249439   0.205999  -1.211    0.226

A2:B2                                 0.095448   0.300565   0.318    0.751

A3:B2                                -0.001116   0.281570  -0.004    0.997

A4:B2                                 0.330682   0.302998   1.091    0.275

A2:C2                                 0.244077   0.282310   0.865    0.387

A3:C2                                -0.215565   0.294696  -0.731    0.464

A4:C2                                 0.172165   0.276917   0.622    0.534

B2:C2                                -0.137896   0.303376  -0.455    0.649

A2:B2:C2                             -0.025005   0.440519  -0.057    0.955

A3:B2:C2                              0.401948   0.396986   1.012    0.311

A4:B2:C2                              0.027745   0.403595   0.069    0.945

Signif. codes:  0 ‘*’ 0.001 ‘’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for Negative Binomial(4.8642) family taken to be 1)
Null deviance: 221.54  on 185  degrees of freedom

Residual deviance: 193.99  on 170  degrees of freedom
  (1 observation deleted due to missingness)
AIC: 2280.9
Number of Fisher Scoring iterations: 1
          Theta:  4.864 
      Std. Err.:  0.505 


2 x log-likelihood:  -2246.898 

In order to determine significance of each term in the model I remove each term from the model and then perform a likelihood ratio test to obtain a p-value for each term. For example:

#To determine significance of the 3-way interaction
model.count.nb.no3way<-glm.nb(X.total ~A+B+C+A:B+A:C+B:C, link = "log", data = count.X)
lrtest(model.count.nb, model.count.nb.no3way)#Interaction insignificant
Likelihood ratio test
Model 1: X.total ~ A * B * C
Model 2: X.total ~ A + B + C + A:B + A:C + B:C
  #Df  LogLik Df  Chisq Pr(>Chisq)
1  17 -1123.5

2  14 -1124.3 -3 1.6887     0.6394
#To determine significance of interaction between A and B
model.count.nb.noAB<-glm.nb(X.total ~ A+B+C+A:C+B:C, link = "log", data = count.X)
summary(model.count.nb.noAB)
lrtest(model.count.nb.no3way, model.count.nb.noAB)#Interaction insignificant
Likelihood ratio test
Model 1: X.total ~ A + B + C + A:B + A:C + B:C
Model 2: X.total ~ A + B + C + A:C + B:C
  #Df  LogLik Df  Chisq Pr(>Chisq)
1  14 -1124.3

2  11 -1125.8 -3 3.0111     0.3899

When determining the significance of either the 2-way or 3-way interactions only one term needs to be dropped from the model in order to use a likelihood ratio test however I am not sure which terms to drop if looking for the significance of the main effects and which model to compare this to. My first thought was to compare the model with no-3way interaction with the model without any terms involving a particular main effect i.e.:

model.count.nb.noC<-glm.nb(X.total ~ A+B+A:B, link = "log", data = count.X)
summary(model.count.nb.noC)
lrtest(model.count.nb.no3way, model.count.nb.noC)

The likelihood ratio test does work however, would this not represent the significance of C+A:C+B:C rather than just the main effect C.

Any help anyone can provide with my query would be greatly appreciated.

Answer 1

First, your stepwise approach, starting with the full model and then removing "insignificant" interaction terms, will make subsequent inference incorrect. Once you start using the results of a model to alter the model, the assumptions used for calculating p-values no longer hold. See this answer and this answer. I'll continue based on the full model with all pre-specified interactions.

Second, there's a problem defining a "main effect" for a predictor that's involved in an interaction. That's the case for any regression model. There's nothing specific here to a negative binomial model, except for the interpretation of coefficient estimates and some details of which versions of tests to perform. Some possibilities follow.

If you want to assess the overall importance of C in your full model, then what you want to evaluate is C + A:C + B:C + A:B:C. You could do that via a likelihood-ratio test. It can be more convenient to do a Wald "chunk test" on all coefficients involving C.

At the other extreme, evaluating the "main effect" in terms of the individual coefficient for C in your model is a fool's errand. Changing the reference level (the level internally coded as 0) of either A or B would change the individual coefficient estimate for C, as they both interact with C. Changing the reference (0) level is the same as centering a continuous predictor involved in an interaction. The fundamental model and its predictions would be the same, but the value reported for the individual coefficient of C would be different and thus so would the apparent "significance" of its difference from 0.

You can perform some type of analysis of variance. If you have a perfectly balanced model, you could call anova() directly on your model.count.nb object to get sequential likelihood ratio tests on all predictors and their interactions, via the anova.negbin() function that the MASS package invokes for such models. If the model isn't fully balanced, however, I think that will give you problems.

The Anova() function (note the capital "A") in the car package can avoid the problems of imbalance. Applied to your model.count.nb object, its default "Type-II" test would evaluate, for example, the significance of C after taking into account all terms (including interactions) not involving C. From the help page:

Type-II tests are calculated according to the principle of marginality, testing each term after all others, except ignoring the term's higher-order relatives.

It's also possible to perform "Type-III" tests, "testing each term in the model after all of the others." (Emphasis added.) This page outlines the different type of ANOVA calculations. To do that with Anova() in the car package you would have to re-code your predictors to use contr.sum() or another set of orthogonal contrasts (not the default "treatment" contrasts), as explained in the help page linked above. Alternatively, the emmeans package can apparently perform Type-III tests without predictor recoding, according to the package's author in an answer illustrating functions in its predecessor lsmeans package.

This answer explains a way to get an effect measure for each of A, B and C in your model, evaluated by "averaging [predictor values] over the N cases that were actually observed." If you choose that approach, you can get standard errors and confidence intervals by combining the formula for the variance of a weighted sum of variables with the coefficient variance-covariance matrix.

Answer 2

Assuming your primary goal is inferential, you should not base your model building strategy solely around statistical significance testing but also on what is theoretically or 'clinically' meaningful.

Is B an important potential effect moderator for the relationship of A and X.total? What about C? Then they should be included whether or not statistically significant.

There are also other measures you can use to help with model selection such as AIC/BIC and pseudo R-squared, but these come with their own limitations as well.