# Interpreting the Explanatory Effect of Dropping Non-Significant Variables

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)

dat_loc <- "https://github.com/b-shelton/stack_questions/blob/main/explanatory_effect_change_example.csv?raw=true"

# 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")