Submission: Knit this file to HTML, publish to RPubs with the title epi553_hw02_lastname_firstname, and paste the RPubs URL in the Brightspace submission comments. Also upload this .Rmd file to Brightspace.

AI Policy: AI tools are NOT permitted on this assignment. See the assignment description for full details.

Part 0: Data Preparation (10 points)

# Load required packages
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.4.3
## Warning: package 'ggplot2' was built under R version 4.4.3
## Warning: package 'tibble' was built under R version 4.4.3
## Warning: package 'tidyr' was built under R version 4.4.3
## Warning: package 'readr' was built under R version 4.4.3
## Warning: package 'purrr' was built under R version 4.4.3
## Warning: package 'dplyr' was built under R version 4.4.3
## Warning: package 'stringr' was built under R version 4.4.3
## Warning: package 'forcats' was built under R version 4.4.3
## Warning: package 'lubridate' was built under R version 4.4.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.6
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.1     ✔ tibble    3.3.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.2
## ✔ purrr     1.2.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.4.3
## 
## Attaching package: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(broom)
## Warning: package 'broom' was built under R version 4.4.3
# Import the dataset — update the path if needed
bmd <- read.csv("/Users/sarah/OneDrive/Documents/EPI 553/data/bmd.csv")

# Quick check
glimpse(bmd)
## Rows: 2,898
## Columns: 14
## $ SEQN     <int> 93705, 93708, 93709, 93711, 93713, 93714, 93715, 93716, 93721…
## $ RIAGENDR <int> 2, 2, 2, 1, 1, 2, 1, 1, 2, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2…
## $ RIDAGEYR <int> 66, 66, 75, 56, 67, 54, 71, 61, 60, 60, 64, 67, 70, 53, 57, 7…
## $ RIDRETH1 <int> 4, 5, 4, 5, 3, 4, 5, 5, 1, 3, 3, 1, 5, 4, 2, 3, 2, 4, 4, 3, 3…
## $ BMXBMI   <dbl> 31.7, 23.7, 38.9, 21.3, 23.5, 39.9, 22.5, 30.7, 35.9, 23.8, 2…
## $ smoker   <int> 2, 3, 1, 3, 1, 2, 1, 2, 3, 3, 3, 2, 2, 3, 3, 2, 3, 1, 2, 1, 1…
## $ totmet   <int> 240, 120, 720, 840, 360, NA, 6320, 2400, NA, NA, 1680, 240, 4…
## $ metcat   <int> 0, 0, 1, 1, 0, NA, 2, 2, NA, NA, 2, 0, 0, 0, 1, NA, 0, NA, 1,…
## $ DXXOFBMD <dbl> 1.058, 0.801, 0.880, 0.851, 0.778, 0.994, 0.952, 1.121, NA, 0…
## $ tbmdcat  <int> 0, 1, 0, 1, 1, 0, 0, 0, NA, 1, 0, 0, 1, 0, 0, 1, NA, NA, 0, N…
## $ calcium  <dbl> 503.5, 473.5, NA, 1248.5, 660.5, 776.0, 452.0, 853.5, 929.0, …
## $ vitd     <dbl> 1.85, 5.85, NA, 3.85, 2.35, 5.65, 3.75, 4.45, 6.05, 6.45, 3.3…
## $ DSQTVD   <dbl> 20.557, 25.000, NA, 25.000, NA, NA, NA, 100.000, 50.000, 46.6…
## $ DSQTCALC <dbl> 211.67, 820.00, NA, 35.00, 13.33, NA, 26.67, 1066.67, 35.00, …
# Recode RIDRETH1 as a labeled factor
bmd <- bmd %>%
  mutate(
    RIDRETH1 = factor(RIDRETH1,
                      levels = 1:5,
                      labels = c("Mexican American", "Other Hispanic",
                                 "Non-Hispanic White", "Non-Hispanic Black", "Other")),

    # Recode RIAGENDR as a labeled factor
    RIAGENDR = factor(RIAGENDR,
                      levels = c(1, 2),
                      labels = c("Male", "Female")),

    # Recode smoker as a labeled factor
    smoker = factor(smoker,
                    levels = c(1, 2, 3),
                    labels = c("Current", "Past", "Never"))
  )
# Report missing values for the key variables
cat("Total N:", nrow(bmd), "\n")
## Total N: 2898
cat("Missing DXXOFBMD:", sum(is.na(bmd$DXXOFBMD)), "\n")
## Missing DXXOFBMD: 612
cat("Missing calcium:", sum(is.na(bmd$calcium)), "\n")
## Missing calcium: 293
# Create the analytic dataset (exclude missing DXXOFBMD or calcium)
bmd_analytic <- bmd %>%
  filter(!is.na(DXXOFBMD), !is.na(calcium))

cat("Final analytic N:", nrow(bmd_analytic), "\n")
## Final analytic N: 2129

Part 1: Exploratory Visualization (15 points)

Research Question: Is there a linear association between dietary calcium intake (calcium, mg/day) and total femur bone mineral density (DXXOFBMD, g/cm²)?

# Create a scatterplot with a fitted regression line
ggplot(bmd_analytic, aes(x = calcium, y = DXXOFBMD)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", se = FALSE, color = "steelblue") +
  labs(
    title   = "Scatterplot of Total Calcium Intake and Bone Mineral Density",
    x       = "Calcium Intake (mg/day",
    y       = "Bone Mineral Density (g/cm²)"
  ) +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

Written interpretation (3–5 sentences):

The scatterplot shows a weak positive linear relationship. As calcium intake increases, the bone mineral density slightly increases. The plot has no notable outliers.

Part 2: Simple Linear Regression (40 points)

Step 1 — Fit the Model (5 points)

# Fit the simple linear regression model
model <- lm(DXXOFBMD ~ calcium, data = bmd_analytic)

# Display the full model summary
summary(model)
## 
## Call:
## lm(formula = DXXOFBMD ~ calcium, data = bmd_analytic)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.55653 -0.10570 -0.00561  0.10719  0.62624 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 8.992e-01  7.192e-03 125.037  < 2e-16 ***
## calcium     3.079e-05  7.453e-06   4.131 3.75e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1582 on 2127 degrees of freedom
## Multiple R-squared:  0.007959,   Adjusted R-squared:  0.007493 
## F-statistic: 17.07 on 1 and 2127 DF,  p-value: 3.751e-05

Step 2 — Interpret the Coefficients (10 points)

A. Intercept (β₀):

The intercept is 8.992e-01. This represents the bone mineral density when the calcium intake is zero. Numerically, it is the point that passes the y-axis. It is not a meaningful quantity, since no one consumes 0mg of calcium.

B. Slope (β₁):

The slope of 3.079e-05 indicates that for every 1-unit increase in calcium a day, bone mineral density increases by 0.0000308 g/cm². This is a very small positive effect.

Step 3 — Statistical Inference (15 points)

# 95% confidence interval for the slope
confint(model)
##                    2.5 %       97.5 %
## (Intercept) 8.851006e-01 9.133069e-01
## calcium     1.617334e-05 4.540649e-05

State your hypotheses:

  • H₀: B1 = 0 There is no linear association between dietary calcium intake and femur bone mineral density.

  • H₁: B1 != 0 There is a linear relationship between dietary calcium intake and femur bone mineral density. The changes in calcium intake are related to the changes in BMD.

Report the test results:

The t-statistic is 4.131, with 2127 degrees of freedom and a p-value of 3.75e-05. We reject the null hypothesis since the p-value is below .05. This indicates that there is statistically significant evidence of a linear association between calcium intake and bone mineral density.

Interpret the 95% confidence interval for β₁:

The 95% confidence interval is .0000162 to .0000454. This means each additional mg/day of calcium intake is associated with an increase in femur bone mineral density anywhere between .000016 and .000045 g/cm².

Step 4 — Model Fit: R² and Residual Standard Error (10 points)

R² (coefficient of determination):

The R-squared value is .007493. About .8% of variance in femur bone mineral density is explained by dietary calcium intake. This suggests that other predictors play a larger role in explaining the variation in femur bone mineral density than calcium intake alone. Residual Standard Error (RSE):

The residual standard error is 0.1582 g/cm². This represents the typical size of the predictor error the model makes when estimating femur bone mineral density from calcium intake alone. The RSE suggests that the model’s predictions are off my about 0.1582 g/cm².

Part 3: Prediction (20 points)

# Create a new data frame with the target predictor value
new_data <- data.frame(calcium = 1000)

# 95% confidence interval for the mean response at calcium = 1000
predict(model, newdata = new_data, interval = "confidence")
##         fit       lwr      upr
## 1 0.9299936 0.9229112 0.937076
# 95% prediction interval for a new individual at calcium = 1000
predict(model, newdata = new_data, interval = "prediction")
##         fit       lwr      upr
## 1 0.9299936 0.6195964 1.240391

Written interpretation (3–6 sentences):

[Answer all four questions from the assignment description: The predicted BMD for someone consuming 1,000 mg of calcium a day is 0.930g/cm². The 95% confidence interval of 0.9229 to 0.9371 g/cm² represents the range of values for the average bone mineral density among individuals consuming 1,000mg/day of calcium. The 95% prediction interval is 0.6196 to 1.2404 g/cm². This interval is wider since it includes the uncertainty in mean response estimate and the individual person variability. A calcium intake of 1,000mg/day is a meaningful value to predict at, as it is a realistic and common value in adults.

Part 4: Reflection (15 points)

Write 200–400 words in continuous prose (not bullet points) addressing all three areas below.

A. Statistical Insight (6 points)

[What does the regression model tell you about the calcium–BMD relationship? Were the results surprising? What are the key limitations of interpreting SLR from a cross-sectional survey as causal evidence? What confounders might explain the observed association?]

B. From ANOVA to Regression (5 points)

[Homework 1 used one-way ANOVA to compare mean BMD across ethnic groups. Now you have used SLR to model BMD as a function of a continuous predictor. Compare these two approaches: what kinds of questions does each method answer? What does regression give you that ANOVA does not? When would you prefer one over the other?]

C. R Programming Growth (4 points)

[What was the most challenging part of this assignment from a programming perspective? How did you work through it? What R skill do you feel more confident about after completing this homework?]

The simple linear regression model of dietary calcium intake on femur bone mineral density shows a statistically significant but small positive association. The slope indicates higher calcium intake is linked to a slightly higher bone mineral density. Calcium intake only accounts for less than 1% of variance in bone mineral density. This is not surprising considering the various factors that affect bone mineral density. Since the data comes from a cross-sectional survey, the model can’t support causal evidence. Calcium intake and femur bone mineral density were measured at the same time; therefore, there’s no way of knowing whether exposure or outcome occurred first. Factors such as age, sex, BMI, physical activity, vitamin D status, etc. could explain the observed association or mask stronger relationships. Calcium intake alone is not a major drive of bone mineral density. Comparing this simple linear regression model approach to the one-way ANOVA from Homework 1, highlights how the two methods answer different questions. The one-way ANOVA was testing whether groups differ. In Homework 1, we were comparing average bone mineral density across different ethnic groups. On the other hand, the regression model estimates the direction and magnitude of the association. This model allows prediction, confidence intervals, and the ability to incorporate multiple predictors. Regression is preferred when the predictor is continuous or when determining how the outcome changes per unit increase of the predictor. ANOVA is most appropriate when comparing groups defined by categories. The most challenging part of this assignment was understanding the prediction and confidence intervals differed slightly. Working with different examples and working through the predict function clarified distinctions between prediction of a new individual and the confidence interval.

End of Homework 2