Team 13 Lab Work
Qingyang Shi, Yiling Zhou, Jaroslav Krc, Yao Lin, and Nathan Ribeiro
2023-11-27
# load the following libraries:
# tidyverse, broom, table1, GGally, car, easystats, datawizard, and modelsummary
lapply(
list('tidyverse', 'broom', 'table1', 'GGally', 'car', 'easystats', 'datawizard', 'modelsummary'),
require,
character.only = TRUE
)## Loading required package: tidyverse## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.3 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.3 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── 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
## Loading required package: broom
##
## Loading required package: table1
##
##
## Attaching package: 'table1'
##
##
## The following objects are masked from 'package:base':
##
## units, units<-
##
##
## Loading required package: GGally
##
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
##
## Loading required package: car
##
## Loading required package: carData
##
##
## Attaching package: 'car'
##
##
## The following object is masked from 'package:dplyr':
##
## recode
##
##
## The following object is masked from 'package:purrr':
##
## some
##
##
## Loading required package: easystats
##
## # Attaching packages: easystats 0.6.0 (red = needs update)
## ✔ bayestestR 0.13.1 ✔ correlation 0.8.4
## ✔ datawizard 0.9.0 ✔ effectsize 0.8.6
## ✖ insight 0.19.6 ✔ modelbased 0.8.6
## ✖ performance 0.10.5 ✖ parameters 0.21.2
## ✔ report 0.5.7 ✖ see 0.8.0
##
## Restart the R-Session and update packages in red with `easystats::easystats_update()`.
##
##
## Loading required package: modelsummary
##
##
## Attaching package: 'modelsummary'
##
##
## The following object is masked from 'package:parameters':
##
## supported_models
##
##
## The following object is masked from 'package:insight':
##
## supported_models## [[1]]
## [1] TRUE
##
## [[2]]
## [1] TRUE
##
## [[3]]
## [1] TRUE
##
## [[4]]
## [1] TRUE
##
## [[5]]
## [1] TRUE
##
## [[6]]
## [1] TRUE
##
## [[7]]
## [1] TRUE
##
## [[8]]
## [1] TRUE# load your dataset and assign it to a variable using read_rds("data/natality.rds")
data <- read_rds("data/natality.rds")The lab specifications can be found [here]
Exercise 1
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.
# Write your code here
data %>%
mutate(CIG_REC = case_when(
CIG_REC == "No" ~ "Non smoker",
TRUE ~ "Smoker"
)) -> data
# Defining pretty names for the labels
label(data$DBWT) <- "Birth weigth (g)"
label(data$risks) <- "Risk factors reported"
label(data$MRACEHISP) <- "Mother race/Hispanic origin"
label(data$DMAR) <- "Marital status"
label(data$MEDUC) <- "Mother education"
label(data$PRIORLIVE) <- "Prior births (living)"
label(data$PRIORDEAD) <- "Prior births (dead)"
label(data$BMI) <- "Mother BMI"
caption <- "Maternal demographics stratified by smoking"
table1(~ DBWT + risks + MRACEHISP + DMAR + MEDUC + PRIORLIVE + PRIORDEAD + BMI | CIG_REC,
data = data)| Non smoker (N=1564) |
Smoker (N=124) |
Overall (N=1688) |
|
|---|---|---|---|
| Birth weigth (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 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] |
Exercise 2
Use the ggpairs function to explore the association
between the different variables. Try to replicate the figure below, or
create your own variation.
# Write your code here
p1 <- ggpairs(data,
columns = c('DBWT', 'risks', 'MAGER', 'BMI'),
aes(color = CIG_REC),
title = 'Natality in the United States (2018)',
columnLabels = c(label(data$DBWT), label(data$risks),
label(data$MAGER), label(data$BMI)),
legend = c(3,3),
lower = list(continuous = "density", combo = "box_no_facet")
) + scale_fill_discrete(name = "")
p1
Exercise 3
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 \(R^2\) and the AIC?
# Write your code here
fit_null <- lm(DBWT ~ 1, data)
fit_unadjusted <- lm(DBWT ~ CIG_REC, data)
fit_full <- lm(DBWT ~ ., data)
fit_optimal <- step(fit_null,
scope = list(upper = fit_full),
direction = "both",
data = data)## 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## Analysis of Variance Table
##
## Model 1: DBWT ~ 1
## Model 2: DBWT ~ CIG_REC
## Model 3: DBWT ~ DMAR + MRACEHISP + BMI + risks + PRIORLIVE + CIG_REC +
## PRIORDEAD
## Model 4: DBWT ~ CIG_REC + risks + MAGER + MRACEHISP + DMAR + MEDUC + PRIORLIVE +
## PRIORDEAD + BMI
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 1687 611839087
## 2 1686 608153901 1 3685187 11.2491 0.0008144 ***
## 3 1675 548005046 11 60148855 16.6914 < 2.2e-16 ***
## 4 1671 547416336 4 588710 0.4493 0.7730024
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Analysis of Variance Table
##
## Model 1: DBWT ~ CIG_REC
## Model 2: DBWT ~ CIG_REC + risks + MAGER + MRACEHISP + DMAR + MEDUC + PRIORLIVE +
## PRIORDEAD + BMI
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 1686 608153901
## 2 1671 547416336 15 60737564 12.36 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modelsummary(
list(
"null model" = fit_full,
"unadjusted model" = fit_unadjusted,
"optimal model" = fit_optimal,
"full model" = fit_full
),
stars = T
)| null model | unadjusted model | optimal model | full model | |
|---|---|---|---|---|
| (Intercept) | 3162.373*** | 3247.903*** | 3080.559*** | 3162.373*** |
| (107.021) | (15.187) | (64.625) | (107.021) | |
| CIG_RECSmoker | −156.593** | −179.096** | −160.083** | −156.593** |
| (56.325) | (56.032) | (55.632) | (56.325) | |
| risksYes | −152.902*** | −158.808*** | −152.902*** | |
| (32.055) | (31.655) | (32.055) | ||
| MAGER | −3.836 | −3.836 | ||
| (3.085) | (3.085) | |||
| MRACEHISPNH Black | −290.175*** | −292.140*** | −290.175*** | |
| (41.829) | (41.690) | (41.829) | ||
| MRACEHISPNH Other | −120.774* | −121.695* | −120.774* | |
| (52.960) | (52.796) | (52.960) | ||
| MRACEHISPHispanic | −72.464+ | −74.492+ | −72.464+ | |
| (39.120) | (38.075) | (39.120) | ||
| DMARUnmarried | −189.782*** | −184.275*** | −189.782*** | |
| (35.084) | (31.823) | (35.084) | ||
| MEDUCHS | 8.448 | 8.448 | ||
| (49.909) | (49.909) | |||
| MEDUCSome college | −0.240 | −0.240 | ||
| (50.838) | (50.838) | |||
| MEDUCUni | 31.610 | 31.610 | ||
| (56.902) | (56.902) | |||
| PRIORLIVE1 | 140.483*** | 130.686*** | 140.483*** | |
| (34.725) | (33.814) | (34.725) | ||
| PRIORLIVE2 | 123.385** | 107.170* | 123.385** | |
| (43.880) | (41.746) | (43.880) | ||
| PRIORLIVE3 | 161.298* | 138.776* | 161.298* | |
| (63.149) | (59.972) | (63.149) | ||
| PRIORLIVE4+ | 162.508* | 131.949* | 162.508* | |
| (67.562) | (62.333) | (67.562) | ||
| PRIORDEAD1+ | 265.467+ | 258.434+ | 265.467+ | |
| (154.184) | (153.905) | (154.184) | ||
| BMI | 10.493*** | 10.261*** | 10.493*** | |
| (2.277) | (2.259) | (2.277) | ||
| Num.Obs. | 1688 | 1688 | 1688 | 1688 |
| R2 | 0.105 | 0.006 | 0.104 | 0.105 |
| R2 Adj. | 0.097 | 0.005 | 0.098 | 0.097 |
| AIC | 26246.1 | 26393.7 | 26239.9 | 26246.1 |
| BIC | 26343.8 | 26410.0 | 26315.9 | 26343.8 |
| Log.Lik. | −13105.039 | −13193.843 | −13105.946 | −13105.039 |
| F | 12.291 | 10.217 | 16.259 | 12.291 |
| RMSE | 569.47 | 600.23 | 569.78 | 569.47 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
Answer: The conclusion according to ANOVA of model comparison is that unadjusted model is significantly better than the null model but worse then the optimal model, while there is no significant difference between optimal model and full model. If only relying on the AIC, then I can’t differentiate the random difference of AIC or true difference of AIC without a statistical test. And also sometime MSE can indicator the prediction performance but need splitting a test data set.
Exercise 4
Center and scale all the continuous variables, and create a forest plot for the optimal model.
# Write your code here
data %>%
standardize(two_sd = T) -> data.sta
data.sta %>%
lm(data = .,
formula = DBWT ~ DMAR + MRACEHISP + BMI + risks + PRIORLIVE + CIG_REC + PRIORDEAD) -> fit_optimal_sta
fit_optimal_sta %>%
model_parameters() %>%
plot()
Answer: The forest plot shows some variable are negatively associated with the birth weight and some do so positively. The race, smoking and married status seems influence the birth weight most and cause the low birth weight. But this is only a association without grounding by causality assumption. If we justice all other variables are the confounders of the causality between smoking and birth weight, then smoking would be the etiology of low birth weight.
Exercise 5
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.
# Write your code here
fit_unadjusted_sta <- lm(DBWT ~ CIG_REC, data = data.sta)
fit_full_sta <- lm(DBWT ~ ., data = data.sta)
fit_optimal_sta %>%
tidy(conf.int = T) %>%
mutate(
across(where(is.numeric), \(x) round(x, 2))
) %>%
filter(term == "CIG_RECSmoker") %>%
mutate(model = 'optimal') %>%
add_row(
{
fit_unadjusted_sta %>%
tidy(conf.int = T) %>%
mutate(
across(where(is.numeric), \(x) round(x, 2))
) %>%
filter(term == "CIG_RECSmoker") %>%
mutate(model = "unadjusted")
}
) %>%
add_row(
{
fit_full_sta %>%
tidy(conf.int = T) %>%
mutate(
across(where(is.numeric), \(x) round(x, 2))
) %>%
filter(term == "CIG_RECSmoker") %>%
mutate(model = "full")
}
) %>%
select(model, term, estimate, conf.low, conf.high, p.value)## # A tibble: 3 × 6
## model term estimate conf.low conf.high p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 optimal CIG_RECSmoker -0.13 -0.22 -0.04 0
## 2 unadjusted CIG_RECSmoker -0.15 -0.24 -0.06 0
## 3 full CIG_RECSmoker -0.13 -0.22 -0.04 0.01Answer: The results show no difference in terms of effect of smoking between optimal model and full model. But the p-value seems a little bit large in full model than that in optimal model. While unadjusted model shows difference compared with others. The reason seems there exist confounders which did not include in the model so than the effect of smoking is overestimated in this model by 0.02.
Exercise 6
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
a dumbell
plot 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: When creating the dataset for plotting,
include the value zero for each predictor not included in the
interaction (the reference level for all categorical variables). Picking
the value zero in a centered predictor is like choosing the mean of that
variable, which is a reasonable choice.
# Write your code here
variable_names <- names(data.sta)[-1]
formula <- paste0(variable_names, collapse = "+")
formula <- paste0("DBWT~", formula, "+CIG_REC:MEDUC:MRACEHISP")
fit_full_interaction <- lm(as.formula(formula), data = data)
anova(fit_full_interaction)## Analysis of Variance Table
##
## Response: DBWT
## Df Sum Sq Mean Sq F value Pr(>F)
## CIG_REC 1 3685187 3685187 11.2456 0.0008162 ***
## risks 1 3206912 3206912 9.7861 0.0017891 **
## MAGER 1 4173062 4173062 12.7344 0.0003693 ***
## MRACEHISP 3 27855950 9285317 28.3347 < 2.2e-16 ***
## DMAR 1 11099717 11099717 33.8715 7.060e-09 ***
## MEDUC 3 334058 111353 0.3398 0.7965591
## PRIORLIVE 4 6248946 1562236 4.7673 0.0007963 ***
## PRIORDEAD 1 860120 860120 2.6247 0.1054027
## BMI 1 6958800 6958800 21.2352 4.375e-06 ***
## CIG_REC:MRACEHISP:MEDUC 22 7037893 319904 0.9762 0.4920570
## Residuals 1649 540378443 327701
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(marginaleffects)
predictions(
fit_full_interaction,
newdata = datagrid(
MRACEHISP = c("NH White", "NH Black", "NH Other", "Hispanic"),
MEDUC = c("lower", "HS", "Some college", "Uni"),
CIG_REC = c("Non smoker", "Smoker")
)
) %>%
as_tibble() %>%
mutate(DBWT = estimate) %>%
select(DBWT, MRACEHISP, MEDUC, CIG_REC) %>%
ggplot() +
aes(x = CIG_REC, y = DBWT, color = MEDUC, group = MEDUC) +
geom_point() +
geom_line() +
scale_colour_viridis_d(name = "Education level") +
facet_wrap(~MRACEHISP) +
labs(x="Smoking status", y = "Birth weight (g)", title = "Effect of smoking on birth weight stratified by education and race") +
theme_bw()## Warning: Model matrix is rank deficient. Some variance-covariance parameters are
## missing.
Answer: The effect of smoking on birth weight not only modified by different education level of mother but also by different races. At high level of education (e.g., university), the large magnitude of positive effect was shown except for the white people which a reverse association was shown. In Hispanic people, the effects are determined almost by the different level of education. In black people, education by university seems to prevent the negative effect caused by smoking. Overall, there seems to exist heterogeneous effects of smoking on birth weight, however, according to ANOVA, the interaction is not significant, perhaps, due to lack of detection power.
Exercise 7
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.
# Write your code here
fit_optimal_interaction <- lm(
DBWT ~ DMAR + MRACEHISP + BMI + risks + PRIORLIVE + CIG_REC + PRIORDEAD + BMI:MRACEHISP,
data = data)
anova(fit_optimal_interaction)## Analysis of Variance Table
##
## Response: DBWT
## Df Sum Sq Mean Sq F value Pr(>F)
## DMAR 1 28672236 28672236 87.6460 < 2.2e-16 ***
## MRACEHISP 3 15359037 5119679 15.6500 4.857e-10 ***
## BMI 1 4532441 4532441 13.8549 0.000204 ***
## risks 1 6059709 6059709 18.5235 1.775e-05 ***
## PRIORLIVE 4 5580134 1395034 4.2644 0.001949 **
## CIG_REC 1 2707991 2707991 8.2779 0.004064 **
## PRIORDEAD 1 922493 922493 2.8199 0.093289 .
## MRACEHISP:BMI 3 1032376 344125 1.0519 0.368521
## Residuals 1672 546972670 327137
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1predictions(
fit_optimal_interaction,
newdata = datagrid(
MRACEHISP = c("NH White", "NH Black", "NH Other", "Hispanic"),
BMI = c(mean(data$BMI), mean(data$BMI)+sd(data$BMI))
)
) %>%
as_tibble() %>%
mutate(DBWT = estimate) %>%
mutate(
BMI =
case_when(
BMI == mean(data$BMI) ~ "mean BMI",
TRUE ~ "mean BMI + sd BMI"
)
) %>%
select(DBWT, MRACEHISP, BMI) %>%
ggplot() +
aes(x = as.factor(BMI), y = DBWT, color = MRACEHISP, group = MRACEHISP) +
geom_point() +
geom_line() +
scale_colour_viridis_d(name = "Race") +
labs(x="BMI", y = "Birth weight (g)", title = "Effect of BMI on birth weight stratified by race") +
theme_bw()
Answer: Although the figure seems to illustrate a interaction by race, the ANOVA showed no significant difference of the interaction term. So the interaction effect seems to be attributed to random error.
Exercise 8
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.
# Write your code here
variable_names <- names(data.sta)[-1]
formula <- paste0(variable_names, collapse = "+")
formula <- paste0("DBWT~", formula, "+BMI:MAGER")
fit_full_interaction_age <- lm(as.formula(formula), data = data)
anova(fit_full_interaction_age)## Analysis of Variance Table
##
## Response: DBWT
## Df Sum Sq Mean Sq F value Pr(>F)
## CIG_REC 1 3685187 3685187 11.2442 0.0008166 ***
## risks 1 3206912 3206912 9.7849 0.0017899 **
## MAGER 1 4173062 4173062 12.7328 0.0003695 ***
## MRACEHISP 3 27855950 9285317 28.3312 < 2.2e-16 ***
## DMAR 1 11099717 11099717 33.8673 7.059e-09 ***
## MEDUC 3 334058 111353 0.3398 0.7965900
## PRIORLIVE 4 6248946 1562236 4.7667 0.0007967 ***
## PRIORDEAD 1 860120 860120 2.6244 0.1054221
## BMI 1 6958800 6958800 21.2326 4.377e-06 ***
## MAGER:BMI 1 87600 87600 0.2673 0.6052290
## Residuals 1670 547328737 327742
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1predictions(
fit_full_interaction_age,
newdata = datagrid(
MAGER = c(20, 30, 40),
BMI = c(mean(data$BMI), mean(data$BMI)+sd(data$BMI))
)
) %>%
as_tibble() %>%
mutate(DBWT = estimate) %>%
mutate(
BMI =
case_when(
BMI == mean(data$BMI) ~ "mean BMI",
TRUE ~ "mean BMI + sd BMI"
)
) %>%
select(DBWT, MAGER, BMI) %>%
ggplot() +
aes(x = as.factor(BMI), y = DBWT, color = as.factor(MAGER), group = as.factor(MAGER)) +
geom_point() +
geom_line() +
scale_colour_viridis_d(name = "Mother's age") +
labs(x="BMI", y = "Birth weight (g)", title = "Effect of BMI on birth weight stratified by mother's age") +
theme_bw()
Answer: Well, both figure and ANOVA test showed the same thing that there is no interaction by mother’s age.