Testing a linear combination of coefficients in R

Question

Let's say I have the following regression:

mort_probit <- glm(
  mort ~ age + I(age*age) + hs_grad + 
         some_college + college + post_grad + 
         black + hisp + other + rich + middle_class, 
  family = binomial("probit"), 
  data=data
)

This is set up such that the base case is a poor white. I want to test if a wealthy black individual has a lower risk for mortality than a poor white.

I set up a linear combination of coefficients in the following way:

summary(glht(mort_probit, linfct = c("(black + rich) = 0 ")))

But I'm not sure if this answers my question. What exactly is this command doing? I'm guessing it sets "black" and "rich" to 1 and test if that minus the base case = 0.

What is the correct way of doing this?

Answer 1

But I'm not sure if this answers my question. What exactly is this command doing? I'm guessing it sets "black" and "rich" to 1 and test if that minus the base case = 0.

It does, it your model (model type, predictor terms) is correctly specified.

It depends on the coding of variables black and rich. If black = 1 for a black person and black = 0 for a white person and rich = 1 for a rich person and rich = 0 for a poor person, then the sum of coefficients of black and rich indeed measures the difference in probit (inverse probability function of a normal distribution, or normal quantile value, or linear predictor) of mort between a black rich person (black = rich = 1) and a white poor person (black = rich = 0). For example, if $\beta_\text{black} = 0.5$ and $\beta_\text{rich} = -0.7$, then $\beta_\text{black} + \beta_\text{rich} = -0.2$ which means a black rich person has a lower predicted risk of mortality than a white poor person, controlling for all other predictors in the model.

Nevertheless, risk is usually measured and best understood as a probability instead of quantile (i.e., probit or logit). Because normal and logistic cumulative probability curves are nonlinear, the effect of any predictor on probability is nonlinear and depends on all predictors. A -0.2 difference in probit can be between quantiles at 1 and 0.8, or -0.4 and -0.6, or something else, which correspond to probabilities at .8413 and .7881 (-.0532), or .3446 and .2743 (-.0703 difference). So we should translate the contrast in linear predictor to probability to facilitate communication. See tutorials at https://marginaleffects.com/chapters/categorical.html.

However, you have a few other highly correlated variables hisp, other, and middle_class which may not be treated constant when varying black and rich. Can I person still be in the middle class without being rich? It is useful to show a percentage breakout of different combinations of these categorical groups. Also, when categorical predictors are involved, it is very important to examine interaction effects: add terms like black:rich black:age college:rich. The combined effect of being a black AND rich person may be different than the sum of the effects of being black and being rich separately. See my answer to Frank Harrell's interpretation of interaction in regression results.

What is the correct way of doing this?

Mortality and disease risks are usually studied in survival analysis instead of binary regression. It seems your data report the age of a person at death. A relevant model will be

survival::coxph(Surv(age, mort) ~ black + rich + black:rich + ...)

This assumes that age is continuous (i.e. calculated as days instead of years). Note that for a coxph() model, the predicted quantities include type=c("lp", "risk", "expected", "terms", "survival"). The last is survival probability at a give time. None of them gives (median) survival time, which must be acquired from survfit(coxph(), newdata = , times =).

If the data only record integer age, discrete survival analysis will be more appropriate. See discussions at Discrete-Time Event History (Survival) Model in R and book Tutz, G., & Schmid, M. (2016). Modeling discrete time-to-event data. Springer Nature. https://doi.org/10.1007/978-3-319-28158-2. It essentially entails glm(mort ~ s(age) + black + rich + black:rich + ...), family = binomial("cloglog")) where s() is a smoothing curve and the data must be restructured so that each patient has several rows, each for a year until the death year where age increases by 1 for each year. For your data, probably age is the only predictor that varies by row within the same patient.

Note that both glht() and marginaleffects::comparisons() conduct Wald tests, an extension of the single-coefficient z test based on large-sample asymptotic. For generalized linear models, maximum-likelihood estimates are consistent but biased, so the test is meaningful only if the sample size is large enough. While the former function tests only coefficients, the latter can also test difference in predicted probabilities.

Wald tests require a well-behaving likelihood function (sufficiently approximating to a quadratic curve at its maximum). Basically it means the coefficient and its standard error should not be very large. For the Wald test caveats, see Likelihood ratio vs. score vs. Wald test: Different p values, which to use?. To conduct likelihood-ratio on linear combinations of coefficients, which many test routines do not provide directly, see Likelihood-ratio and score tests of a (non)linear combination of coefficients.