Training data

Contains measurements of all variables.

birthwt_train <- slice_sample(birthwt, prop = 0.9)

New data

Only contains measurements of covariates … but we’d like to know the outcome

birthwt_new <- birthwt %>%
  anti_join(birthwt_train) %>%
  select(-birth_weight)

Step 1: Fit/train/learn a model

Ingredients:

1. Type of model : How complicated can the relationships between covariates and outcome be?
2. Formula : Which outcome and covariates are we considering?
3. Training data : Which measurements are we using to inform our guesses?

fit <- lm(
  formula = birth_weight ~ m_age + m_weight + m_smoker + m_prem_labors,
  data    = birthwt_train
)

Step 2: Predict outcomes

new_obs <- birthwt_new[1:5, ]

What’s our best educated guess for the birth weight of these people’s children?

predict(fit, newdata = new_obs)

##        1        2        3        4        5 
## 2.985043 2.715692 3.037831 3.029939 2.666900

Significance: will dropping this variable change the predictions?

predict(fit, newdata = new_obs)

##        1        2        3        4        5 
## 2.985043 2.715692 3.037831 3.029939 2.666900

What happens if we remove the non-significant variable m_age?

fit_no_age <- lm(
  formula = birth_weight ~ m_weight + m_smoker + m_prem_labors,
  data    = birthwt_train
)

predict(fit_no_age, newdata = new_obs)

##        1        2        3        4        5 
## 2.997543 2.739253 3.038551 3.030349 2.628531

How about the significant m_smoker?

fit_no_smoker <- lm(
  formula = birth_weight ~ m_age + m_weight + m_prem_labors,
  data    = birthwt_train
)

predict(fit_no_smoker, newdata = new_obs)

##        1        2        3        4        5 
## 2.889462 2.880035 2.945557 2.937779 2.849588

Bonus: local linear regression

Bonus: gradient-boosted tree regression

What if the model is wrong?

fit <- lm(y ~ time, data = df)
summary(fit)

## 
## Call:
## lm(formula = y ~ time, data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.145 -2.803 -1.165  1.915 20.722 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -3.0997     0.7787  -3.981 0.000132 ***
## time          3.7385     0.2268  16.485  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.983 on 98 degrees of freedom
## Multiple R-squared:  0.735,  Adjusted R-squared:  0.7322 
## F-statistic: 271.7 on 1 and 98 DF,  p-value: < 2.2e-16

Let’s plot the data

df %>%
  ggplot(aes(x = time, y = y)) +
  geom_point() +
  stat_smooth(method = "lm", se = FALSE, color = g_accent) +
  theme_burnet

The residuals are the errors in our predictions.

residuals <- resid(fit)

df %>%
  ggplot(aes(x = time, y = residuals)) +
  geom_point() +
  theme_burnet

Unbiased residuals:

What if the noise is not Gaussian?

Try to estimate the number of centimetres in an inch by measuring a bunch of objects with a good centimetre ruler and a really terrible inch ruler.

n_obs <- 1000

df <- tibble(
  len_cm = rnorm(n_obs, mean = 10),
  len_in = 2.54 * len_cm + rcauchy(n_obs)
)

Fit the model and find the slope - we expect it to be around 2.54:

fit <- lm(len_in ~ len_cm, data = df)
summary(fit)

## 
## Call:
## lm(formula = len_in ~ len_cm, data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -92.48  -2.57  -1.26   0.00 371.75 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   10.920      6.256   1.745   0.0812 .
## len_cm         1.581      0.623   2.537   0.0113 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 20.38 on 998 degrees of freedom
## Multiple R-squared:  0.006409,   Adjusted R-squared:  0.005413 
## F-statistic: 6.437 on 1 and 998 DF,  p-value: 0.01133

Let’s have a look at the data

df %>%
  ggplot(aes(y = len_in, x = len_cm)) +
  geom_point() +
  theme_burnet

Check if the residuals are normally distributed

resid <- resid(fit)
qqnorm(resid)

Normally distributed residuals