Lab 05 - Multiple Linear Regression
Michalina Prycka & Rebeca Weldeghebreal
13.03.2024
Exercise 1
Question: Use the table1 package to create a table with appropriate descriptive information about the data set. Use the examples provided in the following vignette to help you build a table. Feel free to try a different table with other variables. What do you notice from this table? What research questions might it raise?
suppressMessages({
suppressWarnings({
library(tidyverse)
library(table1)
library(GGally)
library(car)
library(easystats)
library(datawizard)
library(marginaleffects)
library(boot)
library(broom)
natality <- readRDS("/cloud/project/data/natality.rds")
glimpse(natality)
natality2 <- natality
natality$CIG_REC <- factor(natality$CIG_REC, levels = c(0, 1), labels = c("Non-Smoker", "Smoker"))
label(natality2$MAGER) <- "Mother's Age"
label(natality2$PRIORLIVE) <- "Prior births (living) "
label(natality2$PRIORDEAD) <- "Prior births (dead)"
table1(~ DBWT + risks + MAGER + MRACEHISP + DMAR + MEDUC + PRIORLIVE + PRIORDEAD + BMI | CIG_REC, data = natality2)
})
})## Rows: 1,688
## Columns: 10
## $ DBWT <labelled> 3289, 907, 3447, 3030, 3441, 2950, 2812, 4210, 2810, 36…
## $ CIG_REC <fct> No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, …
## $ risks <fct> Yes, Yes, No, No, No, No, No, No, No, No, No, No, No, Yes, N…
## $ MAGER <labelled> 35, 28, 22, 19, 30, 20, 41, 34, 35, 28, 36, 27, 24, 33,…
## $ MRACEHISP <fct> Hispanic, NH White, NH Black, Hispanic, NH White, NH Other, …
## $ DMAR <fct> Married, Unmarried, Unmarried, Unmarried, Married, Unmarried…
## $ MEDUC <fct> Uni, lower, HS, HS, Uni, Some college, Some college, Uni, So…
## $ PRIORLIVE <fct> 1, 2, 0, 1, 1, 0, 0, 1, 1, 2, 0, 0, 1, 3, 0, 3, 1, 1, 4+, 3,…
## $ PRIORDEAD <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ BMI <labelled> 24.9, 27.9, 26.8, 23.6, 33.3, 33.3, 29.9, 25.1, 22.3, 1…| No (N=1564) |
Yes (N=124) |
Overall (N=1688) |
|
|---|---|---|---|
| Birth Weight (g) | |||
| Mean (SD) | 3250 (601) | 3070 (597) | 3230 (602) |
| Median [Min, Max] | 3290 [369, 5880] | 3090 [630, 4220] | 3270 [369, 5880] |
| Risk Factors Reported | |||
| No | 1088 (69.6%) | 84 (67.7%) | 1172 (69.4%) |
| Yes | 476 (30.4%) | 40 (32.3%) | 516 (30.6%) |
| Mother's Age | |||
| Mean (SD) | 29.0 (5.79) | 27.8 (5.79) | 29.0 (5.80) |
| Median [Min, Max] | 29.0 [15.0, 49.0] | 28.0 [17.0, 42.0] | 29.0 [15.0, 49.0] |
| Mother Race/Hispanic Origin | |||
| NH White | 844 (54.0%) | 94 (75.8%) | 938 (55.6%) |
| NH Black | 266 (17.0%) | 14 (11.3%) | 280 (16.6%) |
| NH Other | 125 (8.0%) | 10 (8.1%) | 135 (8.0%) |
| Hispanic | 329 (21.0%) | 6 (4.8%) | 335 (19.8%) |
| Marital Status | |||
| Married | 1017 (65.0%) | 40 (32.3%) | 1057 (62.6%) |
| Unmarried | 547 (35.0%) | 84 (67.7%) | 631 (37.4%) |
| Mother Education | |||
| lower | 177 (11.3%) | 34 (27.4%) | 211 (12.5%) |
| HS | 377 (24.1%) | 45 (36.3%) | 422 (25.0%) |
| Some college | 437 (27.9%) | 37 (29.8%) | 474 (28.1%) |
| Uni | 573 (36.6%) | 8 (6.5%) | 581 (34.4%) |
| Prior births (living) | |||
| 0 | 601 (38.4%) | 38 (30.6%) | 639 (37.9%) |
| 1 | 516 (33.0%) | 41 (33.1%) | 557 (33.0%) |
| 2 | 265 (16.9%) | 17 (13.7%) | 282 (16.7%) |
| 3 | 93 (5.9%) | 17 (13.7%) | 110 (6.5%) |
| 4+ | 89 (5.7%) | 11 (8.9%) | 100 (5.9%) |
| Prior births (dead) | |||
| 0 | 1551 (99.2%) | 123 (99.2%) | 1674 (99.2%) |
| 1+ | 13 (0.8%) | 1 (0.8%) | 14 (0.8%) |
| Mother BMI | |||
| Mean (SD) | 26.9 (6.32) | 27.3 (6.42) | 26.9 (6.33) |
| Median [Min, Max] | 25.6 [15.6, 58.4] | 26.3 [16.5, 50.1] | 25.7 [15.6, 58.4] |
Answer: The table shows differences between smoking and non-smoking mothers in respect to the birth weight (of their baby), risk factor reporting, maternal age, race/ethnicity, marital status, education, prior births (viable and non-viable), and maternal BMI. Based on the table the following research question could be asked: What is the releationship between socio-demographic factors and smoking behavior? What is the relationship between.
Exercise 2
Question: Use the ggpairs function to explore the association between the different variables. Try to replicate the figure below, or create your own variation. What do you see in the figure? What new information does it provide? What are the advantages/disadvantages when comparing table1 to ggpairs?
suppressWarnings({
suppressMessages({
label(natality2$MAGER) <- "Mother's Age"
label(natality2$PRIORLIVE) <- "Prior births (living)"
label(natality2$PRIORDEAD) <- "Prior births (dead)"
ggpairs(
natality2, c(1,3,4,6,10),
mapping = ggplot2::aes(color = CIG_REC),
upper = list(continuous = wrap("density", alpha = 0.9), combo = NULL),
lower = list(continuous = wrap("cor"), combo = NULL),
diag = list(continuous = wrap("densityDiag", colour = NA, alpha = 0.7)),
legend = c(3,5),
#columnLabels = labels,
title = "Natality in the United States (2018)",
showStrips = FALSE,
cardinality_treshold = 5) +
labs(color = "Smoking") +
theme_minimal() +
theme(legend.position = "bottom")
})
})## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Answer: The figure visually represents relationships between (all the) variables in the dataset. This allows for easy comparison and exploration of the various vaariables. The figure provides insights into patterns and correlations that may not be easly visible in a table. Although its visually informative the image lacks the precision as it doesn not include specific values that are presented in a table. This makes is difficult to precisely interpret.
Exercise 3
Question: Fit four linear models, each nested within the other: the null model, the unadjusted model, the full model and the optimal model. Use the step function find the optimal model using the stepwise method. Then use anova to compare the four models. What is your conclusion? Now compare your models using modelsummary function. What does the number in the (brackets) mean? How can you identify the optimal model just by looking at the R2 and the AIC?
suppressMessages({
suppressWarnings({
M.NULL <- lm(DBWT ~ 1, natality2)
M.UNADJUSTED <- lm(DBWT ~ CIG_REC, natality2)
M.FULL <- lm(DBWT ~ ., natality2)
summary_unadjusted <- summary(M.UNADJUSTED)
summary_full <- summary(M.FULL)
library(modelsummary)
modelsummary(list(Null = M.NULL, Unadjusted = M.UNADJUSTED, Full = M.FULL), stars = TRUE)
step(M.NULL,
scope = list(upper = M.FULL),
direction = "both",
data = natality2)
M.OPT <- lm(DBWT ~ DMAR + MRACEHISP + BMI + risks + PRIORLIVE + CIG_REC + PRIORDEAD, natality2)
modelsummary(list(Null = M.NULL, Unadjusted = M.UNADJUSTED, OPT = M.OPT, Full = M.FULL), stars = TRUE)
})
})## Start: AIC=21609.55
## DBWT ~ 1
##
## Df Sum of Sq RSS AIC
## + DMAR 1 28672236 583166851 21530
## + MRACEHISP 3 28557357 583281730 21535
## + PRIORLIVE 4 7485192 604353895 21597
## + MEDUC 3 5189885 606649203 21601
## + CIG_REC 1 3685187 608153901 21601
## + risks 1 3277949 608561138 21602
## + BMI 1 3170730 608668358 21603
## + MAGER 1 2941408 608897679 21603
## + PRIORDEAD 1 849628 610989460 21609
## <none> 611839087 21610
##
## Step: AIC=21530.53
## DBWT ~ DMAR
##
## Df Sum of Sq RSS AIC
## + MRACEHISP 3 15359037 567807814 21492
## + risks 1 4395289 578771562 21520
## + BMI 1 3607980 579558871 21522
## + PRIORLIVE 4 4500359 578666492 21526
## + CIG_REC 1 978350 582188501 21530
## + PRIORDEAD 1 818812 582348039 21530
## <none> 583166851 21530
## + MAGER 1 82475 583084376 21532
## + MEDUC 3 313468 582853384 21536
## - DMAR 1 28672236 611839087 21610
##
## Step: AIC=21491.48
## DBWT ~ DMAR + MRACEHISP
##
## Df Sum of Sq RSS AIC
## + BMI 1 4532441 563275373 21480
## + risks 1 3924653 563883161 21482
## + CIG_REC 1 2181695 565626119 21487
## + PRIORLIVE 4 3880005 563927809 21488
## + PRIORDEAD 1 846700 566961114 21491
## <none> 567807814 21492
## + MAGER 1 52434 567755381 21493
## + MEDUC 3 98014 567709800 21497
## - MRACEHISP 3 15359037 583166851 21530
## - DMAR 1 15473916 583281730 21535
##
## Step: AIC=21479.95
## DBWT ~ DMAR + MRACEHISP + BMI
##
## Df Sum of Sq RSS AIC
## + risks 1 6059709 557215664 21464
## + CIG_REC 1 2362075 560913299 21475
## + PRIORLIVE 4 3457999 559817375 21478
## + PRIORDEAD 1 896345 562379028 21479
## <none> 563275373 21480
## + MAGER 1 125890 563149483 21482
## + MEDUC 3 115009 563160364 21486
## - BMI 1 4532441 567807814 21492
## - MRACEHISP 3 16283498 579558871 21522
## - DMAR 1 15325756 578601129 21523
##
## Step: AIC=21463.69
## DBWT ~ DMAR + MRACEHISP + BMI + risks
##
## Df Sum of Sq RSS AIC
## + PRIORLIVE 4 5580134 551635530 21455
## + CIG_REC 1 2199370 555016294 21459
## + PRIORDEAD 1 1025923 556189741 21463
## <none> 557215664 21464
## + MAGER 1 22588 557193076 21466
## + MEDUC 3 69675 557145989 21470
## - risks 1 6059709 563275373 21480
## - BMI 1 6667497 563883161 21482
## - MRACEHISP 3 15921403 573137067 21505
## - DMAR 1 16431163 573646827 21511
##
## Step: AIC=21454.7
## DBWT ~ DMAR + MRACEHISP + BMI + risks + PRIORLIVE
##
## Df Sum of Sq RSS AIC
## + CIG_REC 1 2707991 548927539 21448
## + PRIORDEAD 1 921529 550714001 21454
## <none> 551635530 21455
## + MAGER 1 280671 551354859 21456
## + MEDUC 3 173640 551461890 21460
## - PRIORLIVE 4 5580134 557215664 21464
## - BMI 1 6486820 558122350 21472
## - risks 1 8181845 559817375 21478
## - MRACEHISP 3 15145908 566781438 21494
## - DMAR 1 14456460 566091990 21496
##
## Step: AIC=21448.4
## DBWT ~ DMAR + MRACEHISP + BMI + risks + PRIORLIVE + CIG_REC
##
## Df Sum of Sq RSS AIC
## + PRIORDEAD 1 922493 548005046 21448
## <none> 548927539 21448
## + MAGER 1 337056 548590483 21449
## + MEDUC 3 75431 548852108 21454
## - CIG_REC 1 2707991 551635530 21455
## - PRIORLIVE 4 6088755 555016294 21459
## - BMI 1 6662694 555590232 21467
## - risks 1 8131760 557059299 21471
## - DMAR 1 10955889 559883427 21480
## - MRACEHISP 3 16576572 565504111 21493
##
## Step: AIC=21447.56
## DBWT ~ DMAR + MRACEHISP + BMI + risks + PRIORLIVE + CIG_REC +
## PRIORDEAD
##
## Df Sum of Sq RSS AIC
## <none> 548005046 21448
## - PRIORDEAD 1 922493 548927539 21448
## + MAGER 1 361028 547644018 21448
## + MEDUC 3 82260 547922786 21453
## - CIG_REC 1 2708955 550714001 21454
## - PRIORLIVE 4 5984714 553989760 21458
## - BMI 1 6749558 554754604 21466
## - risks 1 8234365 556239411 21471
## - DMAR 1 10970232 558975278 21479
## - MRACEHISP 3 16594210 564599256 21492| Null | Unadjusted | OPT | Full | |
|---|---|---|---|---|
| (Intercept) | 3234.746*** | 3247.903*** | 3080.559*** | 3162.373*** |
| (14.658) | (15.187) | (64.625) | (107.021) | |
| CIG_RECYes | −179.096** | −160.083** | −156.593** | |
| (56.032) | (55.632) | (56.325) | ||
| DMARUnmarried | −184.275*** | −189.782*** | ||
| (31.823) | (35.084) | |||
| MRACEHISPNH Black | −292.140*** | −290.175*** | ||
| (41.690) | (41.829) | |||
| MRACEHISPNH Other | −121.695* | −120.774* | ||
| (52.796) | (52.960) | |||
| MRACEHISPHispanic | −74.492+ | −72.464+ | ||
| (38.075) | (39.120) | |||
| BMI | 10.261*** | 10.493*** | ||
| (2.259) | (2.277) | |||
| risksYes | −158.808*** | −152.902*** | ||
| (31.655) | (32.055) | |||
| PRIORLIVE1 | 130.686*** | 140.483*** | ||
| (33.814) | (34.725) | |||
| PRIORLIVE2 | 107.170* | 123.385** | ||
| (41.746) | (43.880) | |||
| PRIORLIVE3 | 138.776* | 161.298* | ||
| (59.972) | (63.149) | |||
| PRIORLIVE4+ | 131.949* | 162.508* | ||
| (62.333) | (67.562) | |||
| PRIORDEAD1+ | 258.434+ | 265.467+ | ||
| (153.905) | (154.184) | |||
| MAGER | −3.836 | |||
| (3.085) | ||||
| MEDUCHS | 8.448 | |||
| (49.909) | ||||
| MEDUCSome college | −0.240 | |||
| (50.838) | ||||
| MEDUCUni | 31.610 | |||
| (56.902) | ||||
| Num.Obs. | 1688 | 1688 | 1688 | 1688 |
| R2 | 0.000 | 0.006 | 0.104 | 0.105 |
| R2 Adj. | 0.000 | 0.005 | 0.098 | 0.097 |
| AIC | 26401.9 | 26393.7 | 26239.9 | 26246.1 |
| BIC | 26412.7 | 26410.0 | 26315.9 | 26343.8 |
| Log.Lik. | −13198.942 | −13193.843 | −13105.946 | −13105.039 |
| F | 10.217 | 16.259 | 12.291 | |
| RMSE | 602.05 | 600.23 | 569.78 | 569.47 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
Answer: The number in brackets stands for the standard errors of the coefficients. By comparing the models with R-squared and AIC, the the optimal model can be said to be the best as it has the highest Adj. R-square (0.98) and lowest AIC (26239.9).
Exercise 4
Question: Center and scale all the continuous variables, and create a forest plot for the optimal model.
library(modelsummary)
natality.ctr <- natality2
numeric_cols <- sapply(natality.ctr, is.numeric)
natality.ctr[numeric_cols] <- scale(natality.ctr[numeric_cols], center = TRUE, scale = TRUE)
M.OPT.ctr <- lm(
DBWT ~ DMAR + MRACEHISP + BMI +
risks + PRIORLIVE + CIG_REC +
PRIORDEAD,
data = natality.ctr
)
M.OPT.ctr %>%
model_parameters() %>%
plot()
Answer: The graph can be see above.
Exercise 5
Question: How do the effect estimate of smoking, the standard error of its estimate, its 95% confidence interval, and its statistical significance change between the un-adjusted, optimal and full models? Please use the centered and scaled models to answer this question.
unadjusted_stats <- M.UNADJUSTED |> tidy(conf.int = TRUE) |> filter(term == "CIG_RECYes")
optimal_stats <- M.OPT.ctr |> tidy(conf.int = TRUE) |> filter(term == "CIG_RECYes")
full_stats <- M.FULL |> tidy(conf.int = TRUE) |> filter(term == "CIG_RECYes")
comparison_table <- data.frame(
model = c("Unadjusted", "Optimal", "Full"),
CIG_RECYes = c(unadjusted_stats$estimate, optimal_stats$estimate, full_stats$estimate),
std.error = c(unadjusted_stats$std.error, optimal_stats$std.error, full_stats$std.error),
statistic = c(unadjusted_stats$statistic, optimal_stats$statistic, full_stats$statistic),
p.value = c(unadjusted_stats$p.value, optimal_stats$p.value, full_stats$p.value),
conf.low = c(unadjusted_stats$conf.low, optimal_stats$conf.low, full_stats$conf.low),
conf.high = c(unadjusted_stats$conf.high, optimal_stats$conf.high, full_stats$conf.high)
)
print(comparison_table)## model CIG_RECYes std.error statistic p.value conf.low
## 1 Unadjusted -179.096362 56.03185072 -3.196331 0.001417655 -288.9956659
## 2 Optimal -0.265817 0.09237762 -2.877505 0.004059133 -0.4470048
## 3 Full -156.593158 56.32546651 -2.780148 0.005494319 -267.0690648
## conf.high
## 1 -69.1970575
## 2 -0.0846293
## 3 -46.1172520Answer: 5. The effect estimate decreases as more variables are included in the model. In the unadjusted model, where only the effect of smoking on the baby’s weight is considered, the effect estimate is -0.149. The addition of 11 variables, as seen in the optimal model, leads to a decrease in the estimate effect to -0.133. In the full model, all variables from the dataset ‘natality’ are included. This results in a further decrease in the estimate effect to a value of -0.130. This decrease is relatively small in comparison to the difference between the unadjusted and optimal model, this could partly be attributed due to the smaller increase in variables in the latter 2 models compared to the first 2 models.
The standard error remains fairly consistent across the 3 models (0.0465, 0.0462 and 0.0468, for the unadjusted, optimal and full model, respectively). The table shows a slight decrease in the optimal model, however, in the full model the standard error increases. These results indicate the optimal model provides a slightly better estimation of the effect of smoking on the baby’s weight.
In the unadjusted model, the CI is the widest for the unadjusted model. With the inclusion of more variables in the optimal model, the CI narrows, where it ranges from -0.224 - -0.0423. In the full model, the CI narrows even more
The p-value associated with the unadjusted model is 0.00142, indicating a strong association between smoking and the birth weight. In the optimal Model, controlling for additional factors leads to an increase in the p-value, however, it remains below the threshold of 0.05 (0.00406), suggesting that the association between smoking and the baby’s weight remains statistically significant. Similarly, in the full model, the p-value increases, however this also ramins below the threshold of 0.05 (0.00549) indicating statistical significance. .
Exercise 6
Question: Does the association between cigarette smoking during pregnancy and birth weight depend on mother’s education? Modify the FULL model appropriately to be able to answer this question. Use anova to test whether the model with the interaction is significantly better than the one without it. Regardless of statistical significance, create ONE of the dumbell plots shown below, illustrating the smoking effect at each level of mother’s education and estimate the cigarette smoking effect at each level of mother’s education. Hint: We will be using the marginaleffects package to create the plots.
suppressMessages(library(ggplot2))
suppressMessages(library(marginaleffects))
suppressMessages(library(ggpubr))
plot_white <- tribble(
~Education, ~No_Smoking, ~Smoking,
"Lower", 0.077584101, -0.0675469069,
"HS", 0.090825343, -0.054305663,
"Some college", 0.081911438, -0.063219569,
"Uni", 0.105681698, -0.039449308
) |>
ggplot() +
geom_segment(
aes(
x = "No Smoking", y = No_Smoking,
xend = "Smoking", yend = Smoking,
color = Education
),
linewidth = 2
) +
geom_point(
aes(x = "No Smoking", y = No_Smoking, color = Education),
size = 6
) +
geom_point(
aes(x = "Smoking", y = Smoking, color = Education),
size = 6
) +
labs(title = "NH White") +
theme_minimal()
plot_black <- tribble(
~Education, ~No_Smoking, ~Smoking,
"Lower", -0.160411832, -0.305542838,
"HS", -0.147170589, -0.292301596,
"Some college", -0.156084495, -0.301215501,
"Uni", -0.132314234, -0.277445241
) |>
ggplot() +
geom_segment(
aes(
x = "No Smoking", y = No_Smoking,
xend = "Smoking", yend = Smoking,
color = Education
),
linewidth = 2
) +
geom_point(
aes(x = "No Smoking", y = No_Smoking, color = Education),
size = 6
) +
geom_point(
aes(x = "Smoking", y = Smoking, color = Education),
size = 6
) +
labs(title = "NH Black") +
theme_minimal()
plot_other <- tribble(
~Education, ~No_Smoking, ~Smoking,
"Lower", -0.019090433, -0.164221439,
"HS", -0.005849190 , -0.150980197 ,
"Some college", -0.014763096, -0.159894103,
"Uni", 0.009007165, -0.136123842
) |>
ggplot() +
geom_segment(
aes(
x = "No Smoking", y = No_Smoking,
xend = "Smoking", yend = Smoking,
color = Education
),
linewidth = 2
) +
geom_point(
aes(x = "No Smoking", y = No_Smoking, color = Education),
size = 6
) +
geom_point(
aes(x = "Smoking", y = Smoking, color = Education),
size = 6
) +
labs(title = "NH Other") +
theme_minimal()
plot_hispanic <- tribble(
~Education, ~No_Smoking, ~Smoking,
"Lower", 0.020581086, -0.1245499219,
"HS", 0.033822328 , -0.111308678,
"Some college",0.0249084233, -0.120222584,
"Uni", 0.048678683, -0.096452323
) |>
ggplot() +
geom_segment(
aes(
x = "No Smoking", y = No_Smoking,
xend = "Smoking", yend = Smoking,
color = Education
),
linewidth = 2
) +
geom_point(
aes(x = "No Smoking", y = No_Smoking, color = Education),
size = 6
) +
geom_point(
aes(x = "Smoking", y = Smoking, color = Education),
size = 6
) +
labs(title = "Hispanic") +
theme_minimal()
combined_plot <- ggpubr::ggarrange(plot_white, plot_black, plot_other, plot_hispanic, ncol = 2, nrow = 2) +
labs(title = "Effect of Smoking on Birth Weight by Mother's Education and Race/Ethnicity") +
theme_minimal()
print(combined_plot)
Answer: ANOVA analysis reveals a highly significant difference between the 2 models (p = < 2.2e-16). The inclusion of the interaction ‘mothers education’ results in a lower RSS, a high F-statistic and a signficantly low p value. Together this indicates the model including the interaction, is a better fit than the previous one, thus suggesting there is a dependency on education.
Exercise 7
Question: Does the association between mother’s BMI and birth-weight differ between mothers of different race/ethnicity? Modify the OPTIMAL model appropriately to be able to answer this question, using the dataset with the centered variables. Use anova to test whether the model with the interaction is significantly better than the one without it. Regardless of statistical significance, estimate the BMI effect at each level of mother’s race/ethnicity and plot the BMI effect at each level of mother’s race/ethnicity using a dumbell plot. Bonus: Provide a brief description, interpretation or reflection about what you see in the plot (e.g., what do you notice? what do you wonder?)
suppressWarnings({
suppressMessages({
library(marginaleffects)
M.OPT_int <- lm( DBWT ~ BMI * MRACEHISP + .,
data = natality.ctr)
anova_result <- anova(M.OPT.ctr, M.OPT_int)
summary(anova_result)
predictions_df <- predictions(M.OPT_int,
newdata = datagrid(
model = M.OPT_int,
MRACEHISP = c("NH White", "NH Black", "NH Other", "Hispanic"),
BMI = seq(min(natality.ctr$BMI), max(natality.ctr$BMI), length.out = 100)
)) %>%
as.data.frame() %>%
mutate(DBWT = estimate) %>%
select(DBWT, MRACEHISP, BMI)
ggplot(predictions_df, aes(x = BMI, y = DBWT, color = MRACEHISP)) +
geom_segment(aes(xend = BMI, yend = DBWT, color = MRACEHISP),
size = 0.5) +
geom_point(size =1) +
labs(x = "Mother's BMI", y = "Birth weight",
title = "Effect of Mother's BMI on Birth Weight by Race/Ethnicity") +
theme_minimal() +
theme(legend.position = "top")
})
})
Answer: We could not finish the code correctly. Preliminary answer: ANOVA analysis was performed to assess whether the association between mother’s BMI and birth-weight differs among race/ ethnicity. Models 1 and 2 were compared, where model 2 includes the interaction term race/ethnicity, and model 1 did not. The p-value associated with the F-statistic for Model 2 is 0.6852, which is greater than the significance level of 0.05, thus indicating the difference in the models is not statistically significant. Overall, the null hypothesis cannot be rejected and based on this data the association between mother’s BMI and birth weight likely does not differ between mothers of different race/ethnicities..
Exercise 8
Question: Does the association between mother’s BMI and birth-weight depend on mother’s age? Modify the FULL model appropriately to be able to answer this question, using the dataset with the centered variables. Use anova to test whether the model with the interaction is significantly better than the one without it. Regardless of statistical significance, estimate the BMI effect for mothers age 20, 30, and 40 years and make a plot illustrating the BMI effect at each of these ages. Bonus: Provide a brief description, interpretation or reflection about what you see in the plot (e.g., what do you notice? what do you wonder?)
suppressWarnings({
suppressMessages({
mean_bmi <- mean(natality.ctr$BMI)
sd_bmi <- sd(natality.ctr$BMI)
predictions_df <- predictions(M.FULL,
newdata = datagrid(
model = M.FULL,
BMI = c(mean_bmi, mean_bmi + sd_bmi), # Added comma here
MRACEHISP = c("NH White", "NH Black", "NH Other", "Hispanic"),
MEDUC = c("lower", "HS", "Some college", "Uni"),
CIG_REC = c("No", "Yes"),
MAGER = c(20, 30, 40)
)
) |>
as.data.frame() |>
mutate(DBWT = estimate) |>
select(DBWT, BMI, MAGER)
predictions_df
})
})## DBWT BMI MAGER
## 1 3085.649 1.388436e-16 20
## 2 3047.287 1.388436e-16 30
## 3 3008.925 1.388436e-16 40
## 4 2929.056 1.388436e-16 20
## 5 2890.694 1.388436e-16 30
## 6 2852.332 1.388436e-16 40
## 7 3094.097 1.388436e-16 20
## 8 3055.735 1.388436e-16 30
## 9 3017.373 1.388436e-16 40
## 10 2937.504 1.388436e-16 20
## 11 2899.142 1.388436e-16 30
## 12 2860.780 1.388436e-16 40
## 13 3085.409 1.388436e-16 20
## 14 3047.047 1.388436e-16 30
## 15 3008.685 1.388436e-16 40
## 16 2928.816 1.388436e-16 20
## 17 2890.454 1.388436e-16 30
## 18 2852.092 1.388436e-16 40
## 19 3117.260 1.388436e-16 20
## 20 3078.898 1.388436e-16 30
## 21 3040.536 1.388436e-16 40
## 22 2960.666 1.388436e-16 20
## 23 2922.304 1.388436e-16 30
## 24 2883.943 1.388436e-16 40
## 25 2795.474 1.388436e-16 20
## 26 2757.112 1.388436e-16 30
## 27 2718.750 1.388436e-16 40
## 28 2638.881 1.388436e-16 20
## 29 2600.519 1.388436e-16 30
## 30 2562.157 1.388436e-16 40
## 31 2803.922 1.388436e-16 20
## 32 2765.560 1.388436e-16 30
## 33 2727.198 1.388436e-16 40
## 34 2647.329 1.388436e-16 20
## 35 2608.967 1.388436e-16 30
## 36 2570.605 1.388436e-16 40
## 37 2795.234 1.388436e-16 20
## 38 2756.872 1.388436e-16 30
## 39 2718.510 1.388436e-16 40
## 40 2638.641 1.388436e-16 20
## 41 2600.279 1.388436e-16 30
## 42 2561.917 1.388436e-16 40
## 43 2827.085 1.388436e-16 20
## 44 2788.723 1.388436e-16 30
## 45 2750.361 1.388436e-16 40
## 46 2670.491 1.388436e-16 20
## 47 2632.129 1.388436e-16 30
## 48 2593.767 1.388436e-16 40
## 49 2964.875 1.388436e-16 20
## 50 2926.513 1.388436e-16 30
## 51 2888.151 1.388436e-16 40
## 52 2808.282 1.388436e-16 20
## 53 2769.920 1.388436e-16 30
## 54 2731.558 1.388436e-16 40
## 55 2973.323 1.388436e-16 20
## 56 2934.961 1.388436e-16 30
## 57 2896.599 1.388436e-16 40
## 58 2816.730 1.388436e-16 20
## 59 2778.368 1.388436e-16 30
## 60 2740.006 1.388436e-16 40
## 61 2964.635 1.388436e-16 20
## 62 2926.273 1.388436e-16 30
## 63 2887.911 1.388436e-16 40
## 64 2808.042 1.388436e-16 20
## 65 2769.680 1.388436e-16 30
## 66 2731.318 1.388436e-16 40
## 67 2996.485 1.388436e-16 20
## 68 2958.124 1.388436e-16 30
## 69 2919.762 1.388436e-16 40
## 70 2839.892 1.388436e-16 20
## 71 2801.530 1.388436e-16 30
## 72 2763.168 1.388436e-16 40
## 73 3013.185 1.388436e-16 20
## 74 2974.823 1.388436e-16 30
## 75 2936.461 1.388436e-16 40
## 76 2856.592 1.388436e-16 20
## 77 2818.230 1.388436e-16 30
## 78 2779.868 1.388436e-16 40
## 79 3021.633 1.388436e-16 20
## 80 2983.271 1.388436e-16 30
## 81 2944.909 1.388436e-16 40
## 82 2865.040 1.388436e-16 20
## 83 2826.678 1.388436e-16 30
## 84 2788.316 1.388436e-16 40
## 85 3012.945 1.388436e-16 20
## 86 2974.583 1.388436e-16 30
## 87 2936.221 1.388436e-16 40
## 88 2856.352 1.388436e-16 20
## 89 2817.990 1.388436e-16 30
## 90 2779.628 1.388436e-16 40
## 91 3044.795 1.388436e-16 20
## 92 3006.433 1.388436e-16 30
## 93 2968.071 1.388436e-16 40
## 94 2888.202 1.388436e-16 20
## 95 2849.840 1.388436e-16 30
## 96 2811.478 1.388436e-16 40
## 97 3096.142 1.000000e+00 20
## 98 3057.780 1.000000e+00 30
## 99 3019.418 1.000000e+00 40
## 100 2939.549 1.000000e+00 20
## 101 2901.187 1.000000e+00 30
## 102 2862.825 1.000000e+00 40
## 103 3104.590 1.000000e+00 20
## 104 3066.228 1.000000e+00 30
## 105 3027.866 1.000000e+00 40
## 106 2947.997 1.000000e+00 20
## 107 2909.635 1.000000e+00 30
## 108 2871.273 1.000000e+00 40
## 109 3095.902 1.000000e+00 20
## 110 3057.540 1.000000e+00 30
## 111 3019.178 1.000000e+00 40
## 112 2939.309 1.000000e+00 20
## 113 2900.947 1.000000e+00 30
## 114 2862.585 1.000000e+00 40
## 115 3127.752 1.000000e+00 20
## 116 3089.390 1.000000e+00 30
## 117 3051.028 1.000000e+00 40
## 118 2971.159 1.000000e+00 20
## 119 2932.797 1.000000e+00 30
## 120 2894.435 1.000000e+00 40
## 121 2805.967 1.000000e+00 20
## 122 2767.605 1.000000e+00 30
## 123 2729.243 1.000000e+00 40
## 124 2649.374 1.000000e+00 20
## 125 2611.012 1.000000e+00 30
## 126 2572.650 1.000000e+00 40
## 127 2814.415 1.000000e+00 20
## 128 2776.053 1.000000e+00 30
## 129 2737.691 1.000000e+00 40
## 130 2657.822 1.000000e+00 20
## 131 2619.460 1.000000e+00 30
## 132 2581.098 1.000000e+00 40
## 133 2805.727 1.000000e+00 20
## 134 2767.365 1.000000e+00 30
## 135 2729.003 1.000000e+00 40
## 136 2649.134 1.000000e+00 20
## 137 2610.772 1.000000e+00 30
## 138 2572.410 1.000000e+00 40
## 139 2837.577 1.000000e+00 20
## 140 2799.215 1.000000e+00 30
## 141 2760.853 1.000000e+00 40
## 142 2680.984 1.000000e+00 20
## 143 2642.622 1.000000e+00 30
## 144 2604.260 1.000000e+00 40
## 145 2975.368 1.000000e+00 20
## 146 2937.006 1.000000e+00 30
## 147 2898.644 1.000000e+00 40
## 148 2818.775 1.000000e+00 20
## 149 2780.413 1.000000e+00 30
## 150 2742.051 1.000000e+00 40
## 151 2983.816 1.000000e+00 20
## 152 2945.454 1.000000e+00 30
## 153 2907.092 1.000000e+00 40
## 154 2827.222 1.000000e+00 20
## 155 2788.861 1.000000e+00 30
## 156 2750.499 1.000000e+00 40
## 157 2975.128 1.000000e+00 20
## 158 2936.766 1.000000e+00 30
## 159 2898.404 1.000000e+00 40
## 160 2818.535 1.000000e+00 20
## 161 2780.173 1.000000e+00 30
## 162 2741.811 1.000000e+00 40
## 163 3006.978 1.000000e+00 20
## 164 2968.616 1.000000e+00 30
## 165 2930.254 1.000000e+00 40
## 166 2850.385 1.000000e+00 20
## 167 2812.023 1.000000e+00 30
## 168 2773.661 1.000000e+00 40
## 169 3023.677 1.000000e+00 20
## 170 2985.315 1.000000e+00 30
## 171 2946.954 1.000000e+00 40
## 172 2867.084 1.000000e+00 20
## 173 2828.722 1.000000e+00 30
## 174 2790.360 1.000000e+00 40
## 175 3032.125 1.000000e+00 20
## 176 2993.763 1.000000e+00 30
## 177 2955.401 1.000000e+00 40
## 178 2875.532 1.000000e+00 20
## 179 2837.170 1.000000e+00 30
## 180 2798.808 1.000000e+00 40
## 181 3023.437 1.000000e+00 20
## 182 2985.076 1.000000e+00 30
## 183 2946.714 1.000000e+00 40
## 184 2866.844 1.000000e+00 20
## 185 2828.482 1.000000e+00 30
## 186 2790.120 1.000000e+00 40
## 187 3055.288 1.000000e+00 20
## 188 3016.926 1.000000e+00 30
## 189 2978.564 1.000000e+00 40
## 190 2898.695 1.000000e+00 20
## 191 2860.333 1.000000e+00 30
## 192 2821.971 1.000000e+00 40Answer: We could not figure out this code fully. Preliminary reponse: ANOVA analysis was performed to compare the models 1 and 2, where model 2 included the interaction term mothers age, and model 1 did not. ANOVA analyis suggests there is a better fit for model 2, as indicated by the highly significant p value (<2.2e- 16) and lower RSS. Thus, the association between mother’s BMI and birth weight is likely to be dependent on the the mother’s age.