Interpreting the Explanatory Effect of Dropping Non-Significant Variables

Question

I have a data set comprised of an overdispersed Poisson response variable, two standardized continuous features, and a categorical feature. I am using deviation coding on my categorical feature and I am manually applying the deviation coding, because I want to visualize how the removal of non-significant categories from effect the coefficients of all the other model inputs.

As you can see in the reproducible example below, removing the non-significant categories from the model makes the 1 WHITE category go from significant to not significant, using p-value < .05 as the standard for variable significance.

Should I consider the 1 WHITE category as significant or not, since it is significant when including non-significant variables, but becomes non-significant when the other non-significant variables are removed.

My goal at this point is explanation, not prediction.

library(dplyr)
library(MASS)
library(ggplot2)

# Download data
dat_loc <- "https://github.com/b-shelton/stack_questions/blob/main/explanatory_effect_change_example.csv?raw=true"
dat <- read.csv(dat_loc)


# Manually create deviation coding features
races <- array(sort(unique(dat$race)))
contrasts(dat$race) <- contr.sum(length(races))
race_contrasts <- contrasts(dat$race)

for (i in c(1:(length(races)-1))) {
  focal_race <- races[i]
  contrast_df <- data.frame("race" = row.names(race_contrasts))
  contrast_df[, focal_race] <- race_contrasts[, i]
  dat <- dplyr::inner_join(dat, contrast_df, by = "race")
  print(paste0("Created contrast-encoded variable: ", focal_race))
}



# Negative Binomial model using all variables
mod1 <- glm.nb(response ~ conditions
                 + doses_per_month
                 +`1 WHITE`
                 +`2 HISPANIC`
                 +`3 BLACK`
                 +`4 ASIAN`
                 +`5 NATIVE HAWAIIAN`
                 +`6 AMERICAN INDIAN`
                 +`7 OTHER`
                 +`8 DECLINED` 
               , data = dat)

summary(mod1)

#Call:
#glm.nb(formula = response ~ conditions + doses_per_month + `1 WHITE` + 
#    `2 HISPANIC` + `3 BLACK` + `4 ASIAN` + `5 NATIVE HAWAIIAN` + 
#    `6 AMERICAN INDIAN` + `7 OTHER` + `8 DECLINED`, data = dat, 
#    init.theta = 0.3100583261, link = log)
#
#Deviance Residuals: 
#    Min       1Q   Median       3Q      Max  
#-1.9391  -0.6595  -0.5773  -0.4190   5.4887  
#
#Coefficients:
#                    Estimate Std. Error z value Pr(>|z|)    
#(Intercept)         -1.50614    0.16753  -8.990  < 2e-16 ***
#conditions           0.29339    0.02585  11.351  < 2e-16 ***
#doses_per_month      0.30618    0.02554  11.987  < 2e-16 ***
#`1 WHITE`            0.34367    0.17429   1.972 0.048627 *  
#`2 HISPANIC`         0.31243    0.17062   1.831 0.067076 .  
#`3 BLACK`            0.70055    0.17746   3.948 7.89e-05 ***
#`4 ASIAN`           -0.74274    0.19497  -3.810 0.000139 ***
#`5 NATIVE HAWAIIAN`  0.24217    0.63144   0.384 0.701334    
#`6 AMERICAN INDIAN` -1.38039    0.99884  -1.382 0.166976    
#`7 OTHER`            0.75020    0.18799   3.991 6.59e-05 ***
#`8 DECLINED`        -0.34630    0.63741  -0.543 0.586926    
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#(Dispersion parameter for Negative Binomial(0.3101) family taken to be 1)
#
#    Null deviance: 6117.7  on 9999  degrees of freedom
#Residual deviance: 5322.6  on 9989  degrees of freedom
#AIC: 14154
#
#Number of Fisher Scoring iterations: 1
#
#
#              Theta:  0.3101 
#          Std. Err.:  0.0138 
#
# 2 x log-likelihood:  -14130.0440 

Using the summary above, I create a second model that only includes variables that were significant to the first model, and then I compare the Incidence Rate Ratio 95% ranges.

# Creating a second model, 
# removing variables proving non-significant from the first model
mod2 <- glm.nb(response ~ conditions
                 + doses_per_month
                 +`1 WHITE`
                 #+`2 HISPANIC`
                 +`3 BLACK`
                 +`4 ASIAN`
                 #+`5 NATIVE HAWAIIAN`
                 #+`6 AMERICAN INDIAN`
                 +`7 OTHER`
                 #+`8 DECLINED`
               , data = dat)


# Create a table that compares the Incidence Rate Ratios between models 1 and 2
conf_mod1 <- data.frame(exp(confint(mod1)))
conf_mod1$coefficient <- exp(coefficients(mod1))
conf_mod1$feature <- row.names(exp(confint(mod1)))
conf_mod1$lower2.5 <- conf_mod1[,1]-1
conf_mod1$upper2.5 <- conf_mod1[,2]-1
conf_mod1$coeff_impact <- conf_mod1$coefficient-1
conf_mod1 <- conf_mod1[c("feature", "lower2.5", "upper2.5", "coeff_impact")]
conf_mod1$model <- "mod1"

conf_mod2 <- data.frame(exp(confint(mod2)))
conf_mod2$coefficient <- exp(coefficients(mod2))
conf_mod2$feature <- row.names(exp(confint(mod2)))
conf_mod2$lower2.5 <- conf_mod2[,1]-1
conf_mod2$upper2.5 <- conf_mod2[,2]-1
conf_mod2$coeff_impact <- conf_mod2$coefficient-1
conf_mod2 <- conf_mod2[c("feature", "lower2.5", "upper2.5", "coeff_impact")]
conf_mod2$model <- "mod2"

# Combine IRRs for models 1 and 2
conf_mods <- rbind(conf_mod1, conf_mod2)


# Order variables for plot
feature_order <- conf_mod1 %>%
  mutate(sig = ifelse((lower2.5 < 0 & upper2.5 < 0) | (lower2.5 > 0 & upper2.5 > 0),
                        1, 0)) %>%
  arrange(desc(sig), desc(upper2.5))

conf_mods$feature <- factor(conf_mods$feature, levels = feature_order$feature)
conf_mod1$feature <- factor(conf_mod1$feature, levels = feature_order$feature)


# Plot comparisons/outcomes
ggplot(conf_mods) + 
  geom_linerange(aes(x=feature, 
                     ymin=lower2.5, 
                     ymax=upper2.5, 
                     colour=model, 
                     group=model), size=6, position = position_dodge(.8)) +
  #geom_point(aes(x=feature, y=coeff_impact)) +
  geom_hline(yintercept=0, color="red", linetype="dashed") +
  theme_bw() +
  theme(axis.text.x=element_text(angle=45, hjust=1)) +
  labs(title = "Negative Binomial Model Incident Rate Ratio Range (95%) Comparisons\nComplete Model (mod1) vs Reduced Model (mod2)",
       x="", y="Incidence Rate Ratio Range")