2026-01-21
Model Selection
Building a model
Selecting variables
Prediction vs interpretation
Comparing potential models
Model Fitting
Find best fit line
Using OLS in this class
Parameter estimation
Categorical covariates
Interactions
Model Evaluation
Model Use (Inference)
Bad option for visualization:
ggplot(gapm, aes(x = fct_reorder(freedom_status, fs_order), y = life_exp)) +
geom_point() +
labs(x = "Freedom status",
y = "Country life expectancy (years)",
title = "Life expectancy vs. Freedom status") +
theme(axis.title = element_text(size = 20),
axis.text = element_text(size = 20),
title = element_text(size = 20))
Good option for visualization:
ggplot(gapm, aes(x = fct_reorder(freedom_status, fs_order), y = life_exp)) +
geom_jitter(size = 1, alpha = .6, width = 0.2) +
stat_summary(fun = mean, geom = "point", size = 8, shape = 18) +
labs(x = "Freedom status",
y = "Country life expectancy (years)",
title = "Life expectancy vs. Freedom status",
caption = "Diamonds = FS averages") +
theme(axis.title = element_text(size = 20),
axis.text = element_text(size = 20),
title = element_text(size = 20))
geom_jitter()Good option for visualization:
ggplot(gapm, aes(x = fct_reorder(freedom_status, fs_order), y = life_exp)) +
geom_boxplot() +
labs(x = "Freedom status",
y = "Country life expectancy (years)",
title = "Life expectancy vs. Freedom status") +
theme(axis.title = element_text(size = 20),
axis.text = element_text(size = 20),
title = element_text(size = 20))
When using a categorical covariate/predictor,
Instead we model the means of the outcome
In 511/525, we used Kruskal-Wallis test and our ANOVA table to test if groups means were statistically different from one another
We can do this using linear models AND we can include other variables in the model
We want to create a function that can map each freedom status to life expectancy
If not free: \(\widehat{LE} = 68.99\) years
If partly free: \(\widehat{LE} = 70.39\) years
If free: \(\widehat{LE} = 74.13\) years
\[\begin{aligned} \widehat{\textrm{LE}} = & 68.99 \cdot I(\text{not free}) + 70.39 \cdot I(\text{partly free}) + \\ &74.13 \cdot I(\text{free}) \end{aligned}\]
However, a linear regression equation still requires an intercept!
\[\begin{aligned} \widehat{\textrm{LE}} = & 68.99 + 1.4 \cdot I(\text{partly free}) + \\ &5.14 \cdot I(\text{free}) \\ \widehat{\textrm{LE}} = & \widehat\beta_0 + \widehat\beta_1 \cdot I(\text{partly free}) + \\ & \widehat\beta_2 \cdot I(\text{free}) \end{aligned}\]
\[\widehat{\textrm{LE}} = \widehat\beta_0 + \widehat\beta_1 \cdot I(FS = \text{partly free}) + \widehat\beta_2 \cdot I(FS = \text{free})\]
| Freedom status | Regression equation for FS | Average Life Expectancy for FS |
|---|---|---|
| Not free | \(\begin{aligned} \widehat{\textrm{LE}} = &\widehat\beta_0 + \widehat\beta_1 \cdot 0 + \widehat\beta_2 \cdot 0 \end{aligned}\) | \(\widehat{\textrm{LE}} = \widehat\beta_0\) |
| Partly free | \(\begin{aligned} \widehat{\textrm{LE}} = &\widehat\beta_0 + \widehat\beta_1 \cdot 1+ \widehat\beta_2 \cdot 0\end{aligned}\) | \(\widehat{\textrm{LE}} = \widehat\beta_0 + \widehat\beta_1\) |
| Free | \(\begin{aligned} \widehat{\textrm{LE}} = &\widehat\beta_0 + \widehat\beta_1 \cdot 0 + \widehat\beta_2 \cdot 1\end{aligned}\) | \(\widehat{\textrm{LE}} = \widehat\beta_0 + \widehat\beta_2\) |
Similar to as before:
Observed values \(Y\) are the values in the dataset
Fitted values \(\widehat{Y}\) are the values that fall on the best-fit line for a specific value of x are the means of the outcome stratified by the categorical predictor’s levels
Residuals (\(\widehat{\epsilon} = Y - \widehat{Y}\)) are the differences between the two
A few changes needed:
Make the income levels readable
ggplot(gapm2, aes(x = income_level_4, y = life_exp)) +
geom_jitter(size = 1, alpha = .6, width = 0.2) +
stat_summary(fun = mean, geom = "point", size = 8, shape = 18) +
labs(x = "Income levels",
y = "Country life expectancy (years)",
title = "Life expectancy vs. income levels",
caption = "Diamonds = Income level averages") +
theme(axis.title = element_text(size = 20),
axis.text = element_text(size = 20),
title = element_text(size = 20),
axis.text.x=element_text(angle = 20, vjust = 1, hjust=1))
Using visual assessment, linearity holds for our income levels (more in next lessons)
We can use the ordinal encoding for income levels
