# Focal Predictions from a linear model: How to test for difference between factor levels (pairwise) instead of comparing errorbars?

Sorry for the cluelessness, I know this is a topic that arises often in different variations. Still, I couldn't find an answer for my situation.

In my work, I sample people and try to generalize the results to say something meaningful about the population from which the sample was taken. To do so, I typically need to account for potential biases in my sample that are due to unrepresentativeness. For example, if I want to study driving habits and ask people how often they speed over the limit, I need to verify that variables such as age or gender -- in the sample I work on -- are distributed in a similar way to the entire population I wish to say something about. Otherwise, if I remain oblivious to influences of confounding variables (e.g., risky driving habits are more likely in young male adults), and take my sample as-is without accounting for those variables, it would damage the validity of any conclusions I would make regarding the population.

My solution for addressing the influence of confounding variables is to run a multiple regression. The predicted variable is my variable of interest, whereas the predictor variables are potential variables that (I have a reason to believe) could confound the results. After fitting the model, I predict. For the predictors, I specify values that are of the "average person in the population". That is, for age I enter the average as it is in the population. For gender, I do something else. If my sample data has gender coded as a factor with 0/1 for male/female, then gender for the average person in the population would be 0.5. The prediction I get reflects the average response in the population, while controlling for confounding variables.

If I want to test whether any of those predictors/confounding variables has a significant influence on the predicted variable, I check the model summary and see which predictors are significant (by their p-value).

Often, my predictors are factors with multiple levels. For example, when testing driving habits, variable such as education might confound the results. Typically, education is a factor variable with the following levels 1 = no high school diploma / 2 = high school diploma / 3 = 2-yrs at College / 4 = 4-yrs at College or Bachelor degree / 5 = Master's degree / 6 = PhD.

When I embark on my investigation and go to sample people, I collect a lot of demographic information. My primary goal when analyzing people's responses is to control for confounding variables while I'm on my mission to say something general about the population. So I throw into the regression any demographic variables that might confound the results. I then predict the "average response in the population" using the method I described above.

However, sometimes, after I get my result, I'm curious to see whether a certain predictor happened to be significant in the model. Furthermore, in the case of a factor variable, I'm curious whether different levels of that predictor were different from each other. My go-to solution is to compare errorbars. Do they overlap or don't? But I know this method is inaccurate.

So my question is, given the procedure detailed above, how can I compare different pairs of factor levels and say whether the difference I see between those levels is indeed significant, rather than just random error?

# Demonstrating with an example

## Multiple linear regression - Measuring Driving habits

### Data

• "How often do you speed over the limit?"
• "How often do you text while driving?"

Choosing one value on a scale of 0 ("Never") to 5 ("Always"). In addition, I measure each person's age, gender, and level of education.

   person_id q_speeding q_texting   age is_female education
<int>      <int>     <int> <int>     <dbl>     <int>
1         1          3         4    31         1         3
2         2          2         1    41         1         6
3         3          5         4    25         1         2
4         4          2         1    25         1         1
5         5          3         0    78         1         3
6         6          5         4    64         1         1
7         7          5         5    46         1         1
8         8          3         3    21         1         5
9         9          2         2    56         1         6
10        10          0         1    61         0         1
# ... with 990 more rows


### Fitting a model

lm(person_response ~ age + is_female + education)


### Predicting the average response in the population (i.e., dealing with sample unrepresentativeness)

prediction_data <- data.frame(age = 45,        ## average age in the population is 45

is_female = 0.5, ## across gender, hence 0.5

education = 2)   ## the average level of education in the population
## is w/ high school diploma, hence 2


### Now I wonder... are any predictors significant?

• Speeding question
Call:
lm(formula = subject_response ~ age + is_female + education,
data = .x)

Residuals:
Min      1Q  Median      3Q     Max
-2.7880 -1.4289 -0.1604  1.4433  3.4144

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  2.622815   0.196877  13.322  < 2e-16 ***
age          0.001013   0.002979   0.340   0.7339
is_female    0.262707   0.116062   2.264   0.0238 *
education   -0.177603   0.040060  -4.433 1.03e-05 ***        ## p-value indicates that gender and
## education are significant!
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.712 on 996 degrees of freedom
Multiple R-squared:  0.02353,   Adjusted R-squared:  0.02059
F-statistic: 8.001 on 3 and 996 DF,  p-value: 2.844e-05

• texting question
Call:
lm(formula = subject_response ~ age + is_female + education,
data = .x)

Residuals:
Min      1Q  Median      3Q     Max
-2.9521 -1.5222  0.2574  1.4463  2.6586

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  2.169711   0.193563  11.209   <2e-16 ***
age          0.004719   0.002929   1.611   0.1074
is_female   -0.019391   0.114108  -0.170   0.8651
education    0.082511   0.039385   2.095   0.0364 *     ## p-value < 0.05 suggests that education is
## significant in this question too.
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.683 on 996 degrees of freedom
Multiple R-squared:  0.006677,  Adjusted R-squared:  0.003685
F-statistic: 2.232 on 3 and 996 DF,  p-value: 0.08298


### Given that education seems to be a significant predictor, I want to split up the results by levels of education

Diving into the levels of education within each question, I'd like to be able to say which levels differ indeed ("significantly") from each other, and which differences might simply reflect "noise" in the data. So far, I've used the malpractice of comparing errorbars. I know it's not the right way, and I want to do better.

## In summary

1. I will be thankful for advice on how to achieve the testing for differences between factor levels given the focal predictions I got.
2. It will be helpful if someone could frame my problem in context. I'm familiar with relevant terminology: "setting contrasts", "ANOVA", "pairwise comparisons", "post-hoc comparisons", etc. But I cannot organize all of this in relation to my problem: dealing with predicted data.
3. Any hint on utilizing the relevant code or functions in R would be very helpful too. I'm posting this question here and not on Stack Overflow because I don't know how to begin asking there. This is first and foremost a question about statistics and only afterwards it's a question about implementing it with programming.

### Appendix -- Code for generating the data, models, and plots

## Data ##

library(tidyverse)
library(magrittr)

set.seed (2022)

n <- 1000

df <-
data.frame(person_id = 1:n,
q_speeding = sample(c(0:5), size = n, replace = TRUE),
q_texting = sample(c(0:5), size = n, replace = TRUE),
age = sample(18:80, size = n, replace = TRUE),
is_female = sample(c(0, 1), prob = c(0.3, 0.7), size = n, replace = TRUE)
)

df %<>%
mutate(education = case_when(q_speeding > 3 ~ sample(c(1:6), prob = c(0.5, 0.3, 0.05, 0.05, 0.05, 0.05), size = n, replace = TRUE),
q_speeding <= 3 ~ sample(c(1:6), prob = c(0.5, 0.2, 0.6, 0.05, 0.05, 0.05), size = n, replace = TRUE)))

as_tibble(df)

## First Plot (Only Means) and Model Underlying It ##

## prediction data ##
prediction_data <- data.frame(age = 45,
is_female = 0.5,
education = 2)

## model ##
model_fits_and_ci <-
df %>%
pivot_longer(starts_with("q_"), values_to = "subject_response") %>%
group_by(name) %>%
tidyr::nest() %>%
mutate(model_fit = map(data, ~ lm(data = .x, subject_response ~ age + is_female + education
)),
predicted_values = map(model_fit, ~ bind_cols(prediction_data,
as.data.frame(predict(newdata = prediction_data, .x,
type = "response", se.fit = T))) %>%
rowwise() %>%
mutate(estimate =  fit,
conf.low =  fit - qt(.975, df) * se.fit,
conf.high = fit + qt(.975, df) * se.fit)))

## plot ##
model_fits_and_ci %>%
unnest(predicted_values) %>%
ggplot(aes(x = name, y = estimate)) +
geom_bar(stat = "identity", position = position_dodge(width = .95), fill = "royalblue", width = 0.4) +
geom_errorbar(stat = "identity", aes(ymin = conf.low, ymax = conf.high),
colour = "black", width = .15, position = position_dodge(width = .95)) +
geom_text(aes(label = round(estimate, 2)), vjust=1.6, color="white", size=4.5, position = position_dodge(width = 1)) +
ylab("average response") +
xlab("question")


## Second Plot (Split by education levels) ##

## prediction data ##
prediction_data <- expand_grid(age = 45,
is_female = 0.5,
education = 1:6)

## model ##
model_fits_and_ci_by_education <-
df %>%
pivot_longer(starts_with("q_"), values_to = "subject_response") %>%
group_by(name) %>%
tidyr::nest() %>%
mutate(model_fit = map(data, ~ lm(data = .x, subject_response ~ age + is_female + as.factor(education)
)),
predicted_values = map(model_fit, ~ bind_cols(prediction_data,
as.data.frame(predict(newdata = prediction_data, .x,
type = "response", se.fit = T))) %>%
rowwise() %>%
mutate(estimate =  fit,
conf.low =  fit - qt(.975, df) * se.fit,
conf.high = fit + qt(.975, df) * se.fit)))

## plot ##
model_fits_and_ci_by_education %>%
unnest(predicted_values) %>%
ggplot(aes(x = name, y = estimate, group = education, fill = as.factor(education))) +
geom_bar(stat = "identity", position = position_dodge(width = .95), width = 0.8) +
geom_errorbar(stat = "identity", aes(ymin = conf.low, ymax = conf.high),
colour = "black", width = .15, position = position_dodge(width = .95)) +
geom_text(aes(label = round(estimate, 1)), vjust=1.6, color="white", size=3.5, position = position_dodge(width = 1)) +
scale_fill_hue(labels = c("no high school diploma",
"high school diploma",
"2-yrs at college",
"4-yrs at college or bachelor degree",
"master's degree",
"phd")) +
labs(fill = "education level") +
ylab("average response") +
xlab("question") +
facet_wrap(~ name, scales = "free")

• I think you have used education as a single predictor in your model. That is why you see only a single coefficient for this predictor. However, if you use as.factor(education), you will have all education levels as separate groups.
– Gijs
Commented Oct 12, 2020 at 16:04
• That's true, but the question about comparing different levels of education to each other remains... Commented Oct 12, 2020 at 16:12
• cran.r-project.org/web/packages/emmeans/vignettes/… Commented Oct 12, 2020 at 16:59
• @BenBolker, I've seen this before. However, it's unclear to me how I can utilize these functions based on the predicted values. As far as I understand, this is good for testing pairs in the model itself, regardless of specific predictions (as in my question). Commented Oct 12, 2020 at 18:26

Very interesting question! I'll provide a couple of hints to hopefully get you started.

First, please refer to the 2019 article A General Framework for Comparing Predictions and Marginal Effects across Models by Mize et al.: https://journals.sagepub.com/doi/full/10.1177/0081175019852763. While this article covers a more general setting than yours (i.e., comparing predictions across multiple models), it does include a nice discussion on comparing predictions within the same model in its section 4.1. Testing the Equality of Predictions and Marginal Effects.

The discussion is simplified, in the sense that it considers a simpler setting where one would be interesting in comparing two predictions. This would be akin to you having only two education levels - no highschool versus highschool or beyond - and wanting to compare the true value of your response variable between the average subject having no highschool education and the average subject having highschool education or beyond. You could accomplish this comparison using a Wald test whose test statistic is computed as $$z = num/denum$$ where:

$$num = pred1 - pred2$$

$$denum = \sqrt{SE(pred1)^2 + SE(pred2)^2 -2Cov(pred1,pred2)}$$

The null hypothesis of equality in the true mean values of the response variable for the two average individuals is rejected if 𝑧$$z$$ exceeds the critical value. (I presume we are dealing with the critical value of the standard normal distrubution, but am not sure this is the case and the article is vague about it.)

Here, pred1 is the predicted value of the response variable corresponding to the average subject with no high scool education and pred2 is the predicted value of the response variable corresponding to the average subject with high scool education or beyond. SE denotes the standard error of the prediction and Cov denotes the covariance between the two predictions.

This link explains in detail how you can get the variance-covariance matrix for your predictions:

https://stackoverflow.com/questions/39337862/linear-model-with-lm-how-to-get-prediction-variance-of-sum-of-predicted-value

In essence, what you would need from this link is something like this:

## compute predictions from lm model
oo <- lm_predict(model, newdat, FALSE)
oo

# extract variance-covariance matrix for
# estimated mean response values at specified
# predictor values
var.fit <- oo$var.fit # extract resudual variance residual.var <- oo$residual.var

## compute variance-covariance matrix for
## predicted individual response values
## at specified predictor values
VCOV_adj <- with(oo, var.fit + diag(residual.var, nrow(var.fit)))


You will need to define your newdata so that it includes all the predictor variables in your model, set to the values you need. An easy way to do that is to use the expand.grid() function. You will also need to copy the definition of the lm_predict() function from the link provided and paste it in your R working space.

The article I mentioned above states that it is possible to expand this Wald testing approach to situations where you would need to compare more than 2 predicted values. The article suggests that "this test can be generalized to test more complex hypotheses, such as the equality of more than two effects (Cameron and Trivedi 2005:135–39)".

I have not yet checked the Cameron and Trivedi reference, but here it is:

Cameron, A. Colin, Trivedi, Pravin K. 2005. Microeconometrics: Methods and Applications. New York: Cambridge University Press.

• @BenBolker: I provided some hints for approaching this question, but I must confess that the Cameron and Trivedi reference is too technical for me. The reference covers Wald Hypothesis Tests of Linear Restrictions in its section 5.5.1 and I suspect that is what is needed (?). Apparently, there is also something like post-prediction inference, though that doesn’t seem to be what is going on here? Commented Oct 13, 2020 at 3:39
• Thank you Isabella. Your answer advances my understanding greatly. But I want to verify that I understood correctly -- the method you proposed cannot allow me to pick 2 specific predictions out of the 6 education levels (e.g., 2-yr college VS. 4-yr college). Rather, I must recode my data so that the factor variable (education in this example) has 2 levels only, and then fit the model + predictions, then conduct the Wald test. Did I understand correctly? Commented Oct 13, 2020 at 7:45
• The method I described in my answer is a simplified case of the method you need. In principle, you don’t need to artificially consider just 2 categories of education - that simplification was done for pedagogical reasons and also because that was the case covered explicitly in the paper I cited. I hope someone who is a theoretical statistician from this forum can chime in and help explain how we should interpret that Cameroni and Trivedi reference in the context of your problem. Commented Oct 13, 2020 at 16:08
• My superficial understanding of the Cameron and Trivedi reference is that you would need to conduct some multivariate form of the Wald test because we are testing multiple hypotheses, but I can’t quite tell from the reference how to implement that. So we are stuck for the moment unless someone else here can help us get unstuck. Commented Oct 13, 2020 at 16:10