HW #6: Linear Regression
Dr. Teresa Vasquez, EdD, MBA, MSDS (in progress)
October 10, 2026

Linear regression is an important statistical model for data analysis in various research including biomedical research. You will recall that the linear model relies on a number of assumptions:

  1. Independent observations
  2. Identically distributed observations
  3. Homoscedasticity
  4. \(\epsilon \sim N(0, \sigma^2)\)

Using the data defined as d1n (NHANES, first 300 complete rows). Note that SCr is the outcome of interest.

# Setup: run this first
suppressPackageStartupMessages(library(tidyverse))
dark <- function() par(bg = "black", col.axis = "white", col.lab = "white", col.main = "white", fg = "white")

# The homework's exact code (downloads NHANES, needs internet)
library(Hmisc, warn.conflicts = FALSE, quietly = TRUE)
got <- tryCatch({ Hmisc::getHdata(nhgh); TRUE }, error = function(e) FALSE)
if (got) {
  d1 <- nhgh %>%
    select(wt, bmi, SCr, age, sex, re, income) %>%
    filter(complete.cases(wt, bmi, SCr, age, sex, re, income))

  # use only n = 300
  d1n = d1[1:300,]
} else {
  d1n <- readRDS("d1n.rds")   # same 300 rows, saved from the code above
}
dim(d1n)
## [1] 300   7

0. Look at the data

d1n has 7 variables: weight wt (kg), body mass index bmi (kg/m²), serum creatinine SCr (mg/dL), age (years), sex, race/ethnicity re (5 groups), and family income (14 brackets). SCr is a blood test for kidney function. Higher SCr usually means the kidneys are filtering less.

0a. What does SCr look like?

How many people are in d1n? Find the mean, median, smallest, and largest SCr, and draw a histogram. Is anything unusual?

Rows: 300 people and 7 variables.

Mean, median, min, max of SCr: mean 0.9061, median 0.83, min 0.44, max 9.13 mg/dL.

Anything unusual? Two people sit far above everyone else: SCr = 9.13 and 7.15 mg/dL (both men, both Non-Hispanic Black, ages 51 and 55). Everyone else is under 3. A typical adult SCr is about 0.6 to 1.2, so these two values are consistent with kidney disease. They pull the mean (0.9061) above the median (0.83) and give SCr a long right tail.

dim(d1n)
## [1] 300   7
summary(d1n$SCr)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.4400  0.7100  0.8300  0.9061  0.9625  9.1300
d1n %>% filter(SCr > 3) %>% select(SCr, age, sex, re, bmi)
##    SCr      age  sex                 re   bmi
## 1 9.13 51.33333 male Non-Hispanic Black 26.52
## 2 7.15 55.33333 male Non-Hispanic Black 36.74
dark()
hist(d1n$SCr, breaks = 60, main = "SCr (mg/dL)", xlab = "SCr", col = "gray30", border = "gray60")

1. SCr and bmi

Question 1. Please assess linear association between serum creatinine concentration (SCr) and body mass index (bmi).

Model 1: \(Y = \beta_0 + \beta_1 x_1 + \epsilon_1\), with \(Y\) = SCr and \(x_1\) = bmi.

1a. Fit model 1 and test the slope

Fit SCr ~ bmi. Write the fitted line. Test \(H_0: \beta_1 = 0\) vs \(H_1: \beta_1 \ne 0\) at \(\alpha = 0.05\), give the 95% CI for \(\beta_1\), and report \(R^2\). Is there a linear association?

Fitted line: \(\widehat{SCr} = 0.736282 + 0.005932\,\text{bmi}\)

Test statistic and decision: \(\dfrac{\hat\beta_1}{SE(\hat\beta_1)} = \dfrac{0.005932}{0.005795} = 1.024\). Since \(|1.024| < T_{1-(0.05/2),\,300-2} = 1.968\), do not reject \(H_0: \beta_1 = 0\) (\(p = 0.307\)).

95% CI: \(0.005932 \pm 1.968 \times 0.005795 = (-0.005471,\ 0.017336)\). It contains 0.

\(R^2\): 0.003505. bmi explains about 0.35% of the variance in SCr.

Conclusion: there is no significant linear association between SCr and bmi. Each extra kg/m² of bmi goes with 0.0059 mg/dL more SCr, but a slope of 0 fits the data just as well.

m1 <- lm(SCr ~ bmi, data = d1n)
summary(m1)
## 
## Call:
## lm(formula = SCr ~ bmi, data = d1n)
## 
## Residuals:
## Creatinine [mg/dLs] 
##     Min      1Q  Median      3Q     Max 
## -0.4699 -0.2018 -0.0705  0.0579  8.2364 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.736282   0.170001   4.331 2.03e-05 ***
## bmi         0.005932   0.005795   1.024    0.307    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6443 on 298 degrees of freedom
## Multiple R-squared:  0.003505,   Adjusted R-squared:  0.0001606 
## F-statistic: 1.048 on 1 and 298 DF,  p-value: 0.3068
confint(m1)
##                    2.5 %     97.5 %
## (Intercept)  0.401726518 1.07083713
## bmi         -0.005471399 0.01733569
qt(0.975, 300 - 2)
## [1] 1.967957

1b. Check the four assumptions

For model 1, plot residuals vs fitted values and a normal Q-Q plot of the residuals. Go through the four assumptions one at a time. Which ones hold?

1. Independent: reasonable. Each row is a different person in the survey, so one person’s SCr doesn’t affect another’s.

2. Identically distributed: reasonable, apart from the two kidney-failure cases. They look like they come from a different population (people with kidney disease).

3. Homoscedasticity: mostly holds. The residuals form an even band across the fitted values, except for the two points far above it.

4. \(\epsilon \sim N(0, \sigma^2)\): violated. The Q-Q plot bends sharply upward at the right end, and the Shapiro-Wilk test rejects normality (\(W = 0.299\), \(p < 2.2 \times 10^{-16}\)). The two extreme SCr values drive this. With \(n = 300\) the slope’s sampling distribution is still close to normal (CLT), but those two residuals inflate \(SE(\hat\beta_1)\) (see 1c).

dark(); par(mfrow = c(1, 2))
plot(fitted(m1), resid(m1), pch = 19, col = "deepskyblue",
     xlab = "Fitted SCr", ylab = "Residual", main = "Residuals vs fitted")
abline(h = 0, col = "orange", lwd = 2)
qqnorm(resid(m1), pch = 19, col = "deepskyblue", main = "Normal Q-Q of residuals")
qqline(resid(m1), col = "orange", lwd = 2)

par(mfrow = c(1, 1))
shapiro.test(resid(m1))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m1)
## W = 0.29938, p-value < 2.2e-16

1c. (extra check) Refit without the two outliers

Refit model 1 without the two people with SCr above 3 (\(n = 298\)). Does the conclusion change? What happens to \(SE(\hat\beta_1)\)?

Slope, SE, p-value: without the two people with SCr above 3 (\(n = 298\)), \(\dfrac{\hat\beta_1}{SE(\hat\beta_1)} = \dfrac{0.003231}{0.002143} = 1.508 < T_{0.975,\,296} = 1.968\), \(p = 0.133\).

Does the conclusion change? No. bmi is still not significantly associated with SCr. What changes is precision: \(SE(\hat\beta_1)\) drops from 0.005795 to 0.002143, because those two points carried most of the residual variance. I report the full-data model as the answer and keep this as a sensitivity check, not a reason to drop data.

keep <- subset(d1n, SCr < 3)
m1_keep <- lm(SCr ~ bmi, data = keep)
c(n = nrow(keep))
##   n 
## 298
coef(summary(m1_keep))
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept) 0.765115357 0.062826385 12.178249 6.237235e-28
## bmi         0.003231245 0.002142866  1.507908 1.326446e-01
qt(0.975, 298 - 2)
## [1] 1.968011

2. SCr and bmi, adjusted for confounders

Question 2. Please assess linear association between SCr and bmi while adjusting for potential confounders in biomedical research.

Model 2: \(Y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \epsilon\), with bmi plus the confounders.

2a. Pick the confounders

Which variables in d1n would you adjust for, and why? Which would you leave out, and why?

Adjust for: - age and sex, the two most common confounders in biomedical research. Both are tied to bmi and to SCr: SCr rises with age and is higher in men (more muscle mass). - race/ethnicity (re): SCr differs by race, and bmi differs by group.

Leave out: - wt: bmi = weight / height², and \(\text{cor}(\text{wt}, \text{bmi}) = 0.867\). Both in one model would measure nearly the same thing (collinearity), and the bmi slope would turn into “bmi at fixed weight,” a different question. - income: Question 4 tests it on its own.

A confounder is linked to both \(x\) (bmi) and \(Y\) (SCr). Leaving one out lets its effect leak into the bmi slope.

cor(d1n$wt, d1n$bmi)
## [1] 0.8670951

2b. Fit model 2 and test bmi

Fit SCr ~ bmi + age + sex + re. Test the bmi slope and give its 95% CI. How does it compare with model 1? Which covariates matter?

bmi slope, test, CI: model 2 is SCr ~ bmi + age + sex + re, with \(p = 7\) predictors (bmi, age, the female dummy, 4 race dummies), so \(df = 300 - 7 - 1 = 292\). \(\dfrac{\hat\beta_1}{SE(\hat\beta_1)} = \dfrac{0.0008701}{0.0056147} = 0.155 < T_{0.975,\,292} = 1.968\), \(p = 0.877\). 95% CI: \((-0.010180,\ 0.011920)\), which contains 0.

Compared with model 1: the bmi slope drops from 0.005932 to 0.000870, about 85% smaller. Part of the crude slope came from age, sex, and race. After adjusting, there is still no significant linear association between SCr and bmi.

Covariates that matter: age (+0.005076 mg/dL per year, \(p = 0.006\)), female (−0.3037 mg/dL vs male, \(p < 0.001\)), and Non-Hispanic Black (+0.4984 mg/dL vs Mexican American, \(p < 0.001\)). Both extreme SCr values (9.13 and 7.15) are in the Non-Hispanic Black group; without them (SCr < 3, n = 298), the Non-Hispanic Black coefficient drops to +0.2355 mg/dL (\(p = 1.5 \times 10^{-9}\)), so it stays significant.

m2 <- lm(SCr ~ bmi + age + sex + re, data = d1n)
summary(m2)
## 
## Call:
## lm(formula = SCr ~ bmi + age + sex + re, data = d1n)
## 
## Residuals:
## Creatinine [mg/dLs] 
##     Min      1Q  Median      3Q     Max 
## -0.7760 -0.1588 -0.0395  0.0854  7.6908 
## 
## Coefficients:
##                                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                          0.6571081  0.1878159   3.499  0.00054 ***
## bmi                                  0.0008701  0.0056147   0.155  0.87695    
## age                                  0.0050759  0.0018353   2.766  0.00604 ** 
## sexfemale                           -0.3036787  0.0705014  -4.307 2.26e-05 ***
## reOther Hispanic                     0.0746771  0.1351053   0.553  0.58087    
## reNon-Hispanic White                 0.1171904  0.0958266   1.223  0.22234    
## reNon-Hispanic Black                 0.4984289  0.1142977   4.361 1.80e-05 ***
## reOther Race Including Multi-Racial  0.0081367  0.1482445   0.055  0.95627    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6059 on 292 degrees of freedom
## Multiple R-squared:  0.1367, Adjusted R-squared:  0.116 
## F-statistic: 6.607 on 7 and 292 DF,  p-value: 2.924e-07
confint(m2)["bmi", ]
##       2.5 %      97.5 % 
## -0.01018020  0.01192049
qt(0.975, 300 - 7 - 1)
## [1] 1.968121
m2_keep <- lm(SCr ~ bmi + age + sex + re, data = subset(d1n, SCr < 3))
coef(summary(m2_keep))["reNon-Hispanic Black", ]
##     Estimate   Std. Error      t value     Pr(>|t|) 
## 2.354890e-01 3.768325e-02 6.249170e+00 1.469205e-09

3. Compare model 1 and model 2

Question 3. Please compare the first model with the second one. Which model fits the data better?

Model 1 is the reduced model (bmi only). Model 2 is the full model (bmi + age + sex + re). Model 1 is model 2 with the 6 confounder \(\beta\)s set to 0, so the F-test applies.

3a. F-test, adjusted \(R^2\), and AIC

Run the F-test of model 1 vs model 2 by hand and check it with anova(). Compare adjusted \(R^2\) and AIC. Which model fits better?

F by hand: model 1 (bmi only) is model 2 with the 6 confounder \(\beta\)s set to 0. \(df_1 = 6\), \(df_2 = 300 - 7 - 1 = 292\).

\[F\text{-}val = \frac{(SSR_{full} - SSR_{reduced})/df_1}{SSE_{full}/df_2} = \frac{(SSE_{reduced} - SSE_{full})/df_1}{SSE_{full}/df_2} = \frac{(123.7252 - 107.1832)/6}{107.1832/292} = 7.511\]

After looking at this again, I realized the numbers I plugged in are the error sums of squares (SSE), not the regression sums of squares (SSR), so I’m giving both forms. They are equal, because SST is the same for both models: \(SSR_{full} - SSR_{reduced} = SSE_{reduced} - SSE_{full}\).

Decision: \(F(6, 292, 0.05) = 2.130\). Since \(7.511 > 2.130\), reject \(H_0\) (\(p = 1.7 \times 10^{-7}\)). anova(m1, m2) gives the same F.

Adjusted \(R^2\) and AIC: adjusted \(R^2\) is 0.0002 for model 1 and 0.1160 for model 2. AIC is 591.65 for model 1 and 560.59 for model 2 (lower is better).

Better model: model 2 fits the data better. All three measures agree. The improvement comes from age, sex, and race, not from bmi.

anova(m1, m2)
## Analysis of Variance Table
## 
## Model 1: SCr ~ bmi
## Model 2: SCr ~ bmi + age + sex + re
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    298 123.72                                  
## 2    292 107.18  6    16.542 7.5109 1.685e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sse_red <- sum(resid(m1)^2); sse_full <- sum(resid(m2)^2)
F_val <- ((sse_red - sse_full) / 6) / (sse_full / 292)
c(F_hand = F_val, F_crit = qf(0.95, 6, 292),
  adjR2_m1 = summary(m1)$adj.r.squared, adjR2_m2 = summary(m2)$adj.r.squared)
##       F_hand       F_crit     adjR2_m1     adjR2_m2 
## 7.5109052812 2.1296863842 0.0001606191 0.1160406631
AIC(m1, m2)
##    df      AIC
## m1  3 591.6473
## m2  9 560.5903

4. SCr and income

Question 4. Please assess linear association between SCr and income. Is this association statistically significant?

income is categorical with 14 brackets, so R makes 13 dummy variables. “Is income associated with SCr?” means “are all 13 income \(\beta\)s equal to 0?” That is one F-test, not 13 t-tests.

Professor’s note (Oct 7): add two potential confounders, age and sex, to the full model. So the full model is SCr ~ income + age + sex and the reduced model is SCr ~ age + sex.

4a. Fit SCr ~ income + age + sex and test income

Fit SCr ~ income + age + sex. State \(H_0\) and \(H_1\) for income, compute the F statistic against the reduced model SCr ~ age + sex, and decide at \(\alpha = 0.05\). One income bracket has \(p < 0.05\) in the unadjusted model. Does that mean income matters?

\(H_0\) / \(H_1\): income has 14 brackets, so the model has 13 income dummies. \(H_0\): all 13 income \(\beta\)s are 0 (after adjusting for age and sex). \(H_1\): at least one is not 0.

F and decision: full = SCr ~ income + age + sex (\(p = 15\): 13 income dummies, age, female), reduced = SCr ~ age + sex. \(df_1 = 13\), \(df_2 = 300 - 15 - 1 = 284\).

\[F\text{-}val = \frac{(SSR_{full} - SSR_{reduced})/df_1}{SSE_{full}/df_2} = \frac{(SSE_{reduced} - SSE_{full})/df_1}{SSE_{full}/df_2} = \frac{(116.0020 - 111.4649)/13}{111.4649/284} = 0.889\]

After looking at this again, I realized the numbers I plugged in are the error sums of squares (SSE), not the regression sums of squares (SSR), so I’m giving both forms. They are equal, because SST is the same for both models: \(SSR_{full} - SSR_{reduced} = SSE_{reduced} - SSE_{full}\).

\(F(13, 284, 0.05) = 1.755\). Since \(0.889 < 1.755\), do not reject \(H_0\) (\(p = 0.565\)). After adjusting for age and sex, the association between SCr and income is not statistically significant.

About the one significant bracket: without the confounders, \([65000, 75000)\) had \(p = 0.024\) on its own. After adjusting for age and sex it is \(p = 0.059\). Either way, one bracket only compares that group with the baseline \([0, 5000)\), and with 13 dummies tested at 0.05, I’d expect \(13 \times 0.05 = 0.65\) false positives by chance. The question is about income as a whole, so the joint F-test is the right test. In the adjusted model, age (\(p = 0.014\)) and sex (\(p < 0.001\)) are the terms that matter.

m4_red <- lm(SCr ~ age + sex, data = d1n)            # reduced: confounders only
m4 <- lm(SCr ~ income + age + sex, data = d1n)       # full: income + confounders
summary(m4)
## 
## Call:
## lm(formula = SCr ~ income + age + sex, data = d1n)
## 
## Residuals:
## Creatinine [mg/dLs] 
##     Min      1Q  Median      3Q     Max 
## -0.6649 -0.1711 -0.0438  0.0815  7.8946 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           0.749011   0.190575   3.930 0.000107 ***
## income[5000,10000)    0.015965   0.234570   0.068 0.945785    
## income[10000,15000)   0.078250   0.231116   0.339 0.735180    
## income[15000,20000)   0.005745   0.221753   0.026 0.979350    
## income[20000,25000)   0.074879   0.214443   0.349 0.727214    
## income[25000,35000)   0.049977   0.196476   0.254 0.799396    
## income[35000,45000)   0.246516   0.203728   1.210 0.227275    
## income[45000,55000)   0.052690   0.231345   0.228 0.820001    
## income[55000,65000)  -0.044229   0.257116  -0.172 0.863546    
## income[65000,75000)   0.422064   0.222257   1.899 0.058579 .  
## income> 20000        -0.083637   0.251646  -0.332 0.739862    
## income< 20000         0.002183   0.654441   0.003 0.997340    
## income[75000,100000) -0.028058   0.210672  -0.133 0.894142    
## income>= 100000       0.057957   0.203022   0.285 0.775489    
## age                   0.004673   0.001893   2.469 0.014128 *  
## sexfemale            -0.259862   0.074283  -3.498 0.000543 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6265 on 284 degrees of freedom
## Multiple R-squared:  0.1023, Adjusted R-squared:  0.05483 
## F-statistic: 2.156 on 15 and 284 DF,  p-value: 0.007946
anova(m4_red, m4)                                    # F-test: are all 13 income betas 0?
## Analysis of Variance Table
## 
## Model 1: SCr ~ age + sex
## Model 2: SCr ~ income + age + sex
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1    297 116.00                           
## 2    284 111.47 13    4.5371 0.8892 0.5649
c(SSE_red = sum(resid(m4_red)^2), SSE_full = sum(resid(m4)^2),
  F_crit = qf(0.95, 13, 300 - 15 - 1))
##    SSE_red   SSE_full     F_crit 
## 116.002025 111.464882   1.754733

5. SCr and age, with a curve

Question 5. Often, age is nonlinearly associated with SCr. Please test if age and \(age^2\) are associated with SCr.

Professor’s note (Oct 7): add two potential confounders, age and sex, to the full model. Age is already in the model, so sex is added.

Model 5: \(Y = \beta_0 + \beta_1 x_1 + \beta_2 x_1^2 + \beta_3 x_2 + \epsilon\), with \(x_1\) = age and \(x_2\) = 1 if female, 0 if male.

5a. Fit age + age² + sex and test the age terms

Fit SCr ~ age + I(age^2) + sex. (1) Test whether age and \(age^2\) together are associated with SCr, adjusting for sex (reduced model SCr ~ sex). (2) Test whether the \(age^2\) term is needed beyond a straight line (reduced model SCr ~ age + sex). Plot the fitted curves over the data.

Fitted model: \(\widehat{SCr} = 0.5864 + 0.01864\,\text{age} - 0.0001531\,\text{age}^2 - 0.2927\,\text{female}\)

(1) Joint test, \(H_0: \beta_1 = \beta_2 = 0\): reduced = SCr ~ sex. \(df_1 = 2\), \(df_2 = 300 - 3 - 1 = 296\).

\[F\text{-}val = \frac{(118.5730 - 115.0030)/2}{115.0030/296} = 4.594 > F(2, 296, 0.05) = 3.026\]

\(p = 0.011\), so reject \(H_0\). After adjusting for sex, age and \(age^2\) together are associated with SCr.

(2) Is \(age^2\) needed, \(H_0: \beta_2 = 0\): \(\dfrac{\hat\beta_2}{SE(\hat\beta_2)} = \dfrac{-0.0001531}{0.00009548} = -1.604\), and \(|-1.604| < T_{0.975,\,296} = 1.968\), so do not reject (\(p = 0.110\)). The same test as an F (age + sex vs age + age² + sex) gives \(F = 2.571 = (-1.604)^2\), \(p = 0.110\).

Conclusion: adjusting for sex, age is significantly associated with SCr (\(p = 0.011\)), but there’s no evidence the relationship is curved (\(p = 0.11\)). In these 300 people, a straight line in age is enough. Women’s SCr is 0.29 mg/dL lower than men’s at the same age (\(p < 0.001\)).

m5 <- lm(SCr ~ age + I(age^2) + sex, data = d1n)
summary(m5)
## 
## Call:
## lm(formula = SCr ~ age + I(age^2) + sex, data = d1n)
## 
## Residuals:
## Creatinine [mg/dLs] 
##     Min      1Q  Median      3Q     Max 
## -0.4937 -0.1823 -0.0687  0.0521  7.9900 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.864e-01  1.833e-01   3.199  0.00153 ** 
## age          1.864e-02  8.871e-03   2.102  0.03642 *  
## I(age^2)    -1.531e-04  9.548e-05  -1.604  0.10988    
## sexfemale   -2.927e-01  7.228e-02  -4.049 6.57e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6233 on 296 degrees of freedom
## Multiple R-squared:  0.07375,    Adjusted R-squared:  0.06437 
## F-statistic: 7.857 on 3 and 296 DF,  p-value: 4.627e-05
anova(lm(SCr ~ sex, data = d1n), m5)                 # (1) joint test of age and age^2
## Analysis of Variance Table
## 
## Model 1: SCr ~ sex
## Model 2: SCr ~ age + I(age^2) + sex
##   Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
## 1    298 118.57                              
## 2    296 115.00  2    3.5699 4.5942 0.01084 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(SCr ~ age + sex, data = d1n), m5)           # (2) is age^2 needed?
## Analysis of Variance Table
## 
## Model 1: SCr ~ age + sex
## Model 2: SCr ~ age + I(age^2) + sex
##   Res.Df RSS Df Sum of Sq      F Pr(>F)
## 1    297 116                           
## 2    296 115  1   0.99901 2.5713 0.1099
c(SSE_red = sum(resid(lm(SCr ~ sex, data = d1n))^2), SSE_full = sum(resid(m5)^2),
  F_crit2 = qf(0.95, 2, 296), T_crit = qt(0.975, 296))
##    SSE_red   SSE_full    F_crit2     T_crit 
## 118.572958 115.003012   3.026257   1.968011
ages <- seq(12, 80, by = 0.5)
dark()
plot(d1n$age, d1n$SCr, pch = ifelse(d1n$sex == "male", 19, 1), col = "deepskyblue",
     xlab = "Age (years)", ylab = "SCr (mg/dL)", ylim = c(0.3, 2.5),
     main = "SCr ~ age + age^2 + sex (2 points above 2.5 not shown)")
lines(ages, predict(m5, data.frame(age = ages, sex = "male")), col = "orange", lwd = 3)
lines(ages, predict(m5, data.frame(age = ages, sex = "female")), col = "white", lwd = 3, lty = 2)
legend("topleft", c("male (filled dots, solid line)", "female (open dots, dashed line)"),
       col = c("orange", "white"), lwd = 3, lty = c(1, 2), bty = "n", text.col = "white")

5b. (extra check) Age without the two outliers

Repeat 5a (with sex in the model) without the two people with SCr above 3. What changes?

Joint test: without the two people with SCr above 3, \(F = 15.43\) on 2 and 294 df, \(p = 4.3 \times 10^{-7}\). Much stronger than \(F = 4.594\) (\(p = 0.011\)) with all 300.

\(age^2\) term: \(F = 1.289\), \(p = 0.257\). Still no curve.

What changed: the two kidney-failure cases are middle-aged men, so they add a lot of noise without following the age trend. Removing them shrinks \(SSE_{full}\) from 115.0030 to 13.0809, and the age signal stands out. The conclusion (age matters, no curve) is the same either way. This is a sensitivity check; my answer to Question 5 uses all 300 rows.

keep <- subset(d1n, SCr < 3)
m5_keep <- lm(SCr ~ age + I(age^2) + sex, data = keep)
anova(lm(SCr ~ sex, data = keep), m5_keep)           # joint test of age and age^2
## Analysis of Variance Table
## 
## Model 1: SCr ~ sex
## Model 2: SCr ~ age + I(age^2) + sex
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    296 14.454                                  
## 2    294 13.081  2    1.3727 15.426 4.257e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(SCr ~ age + sex, data = keep), m5_keep)     # is age^2 needed?
## Analysis of Variance Table
## 
## Model 1: SCr ~ age + sex
## Model 2: SCr ~ age + I(age^2) + sex
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1    295 13.138                           
## 2    294 13.081  1  0.057335 1.2886 0.2572
c(SSE_full = sum(resid(m5_keep)^2))
## SSE_full 
## 13.08091

↑ Back to top