regression_df =
read_csv("data/nhis_data01.csv") %>%
janitor::clean_names() %>%
filter(year == 2020) %>%
select(worrx, deprx, age, sex, poverty, marst, cvddiag) %>%
mutate(
sex = recode_factor(sex,
"1" = "_Male",
"2" = "_Female"),
marital_status = recode_factor(marst,
"10" = "_Married", "11" = "_Married","12" = "_Married","13" = "_Married",
"20" = "_Widowed","30" = "_Divorced","40" = "_Separated",
"50" = "_Never married"),
poverty = recode_factor(poverty,
"11" = "_Less than 1.0", "12" = "_Less than 1.0",
"13" = "_Less than 1.0", "14" = "_Less than 1.0",
"21" = "_1.0-2.0", "22" = "_1.0-2.0",
"23" = "_1.0-2.0", "24" = "_1.0-2.0",
"25" = "_1.0-2.0",
"31" = "_2.0 and above","32" = "_2.0 and above",
"33" = "_2.0 and above","34" = "_2.0 and above",
"35" = "_2.0 and above","36" = "_2.0 and above",
"37" = "_2.0 and above","38" = "_2.0 and above"),
age = ifelse(age>=85, NA, age),
worrx = ifelse(worrx>=3, NA, worrx),
deprx = ifelse(deprx>=3, NA, deprx),
worrx = recode_factor(worrx,
'1' = 0,
'2' = 1),
deprx = recode_factor(deprx, '1' = 0, '2' = 1),
cvddiag = recode_factor(cvddiag,
"1" = "_Never had COVID-19",
"2" = "_Had COVID-19")
) %>%
drop_na(age, worrx, deprx, sex, poverty, marital_status, cvddiag)mosi_regression_df =
read_csv("data/nhis_data01.csv") %>%
janitor::clean_names() %>%
filter(year == 2020) %>%
select(worrx, deprx, age, sex, poverty, marst, cvddiag) %>%
mutate(
sex = recode_factor(sex,
"1" = "Male",
"2" = "Female"),
marital_status = recode_factor(marst,
"10" = "Married", "11" = "Married","12" = "Married","13" = "Married",
"20" = "Widowed","30" = "Divorced","40" = "Separated",
"50" = "Never married"),
poverty = recode_factor(poverty,
"11" = "Less than 1.0", "12" = "Less than 1.0",
"13" = "Less than 1.0", "14" = "Less than 1.0",
"21" = "1.0-2.0", "22" = "1.0-2.0",
"23" = "1.0-2.0", "24" = "1.0-2.0",
"25" = "1.0-2.0",
"31" = "2.0 and above","32" = "2.0 and above",
"33" = "2.0 and above","34" = "2.0 and above",
"35" = "2.0 and above","36" = "2.0 and above",
"37" = "2.0 and above","38" = "2.0 and above"),
age = ifelse(age>=85, NA, age),
worrx = ifelse(worrx>=3, NA, worrx),
deprx = ifelse(deprx>=3, NA, deprx),
worrx = recode_factor(worrx,
'1' = 'No',
'2' = 'Yes'),
deprx = recode_factor(deprx, '1' = 'No', '2' = 'Yes'),
cvddiag = recode_factor(cvddiag,
"1" = "Never had COVID-19",
"2" = "Had COVID-19")
) %>%
drop_na(age, worrx, deprx, sex, poverty, marital_status, cvddiag)Whether taken medication for depression is associated with COVID-19 adjusting for age, sex, family income level and current marital status.
The Mosaic Plot included four categorical variables and it was used to visualize the proportional relationship between these variables and the outcome (whether or not taken medication for depression) in the population.
Based on the plot, we can observe that compared to the other three variables (sex, family income level, and current marital status), there was no obvious difference in the proportion of people who took medication for depression when comparing those who had COVID-19 with those who never had COVID-19. That is, having had COVID-19 or not had no significant effect on whether or not taking medication for depression in the population.
In truth, when comparing male with female, we can observe the largest difference in the proportion of people who took medication for depression. It means that among those four variables, the variable sex has the largest difference in the proportion of people who took medication for depression.
sex_deprx =
mosi_regression_df %>%
ggplot() +
geom_mosaic(
aes(
x = product(deprx, sex),
fill = deprx
),
offset = 0.01,
show.legend = FALSE
)+
labs(
x="",
y=""
)+
theme(
axis.text.y = element_blank(),
axis.ticks.y=element_blank(),
axis.text.x = element_text(angle = 90, hjust = 1)
)
poverty_deprx =
mosi_regression_df %>%
ggplot() +
geom_mosaic(
aes(
x = product(deprx, poverty),
fill = deprx
),
offset = 0.01
)+
labs(
x="",
y="",
fill = "Whether taken medicine for depression")+
theme(
axis.text.y = element_blank(),
axis.ticks.y=element_blank(),
axis.text.x = element_text(angle = 90, hjust = 1)
)
marital_status_deprx =
mosi_regression_df %>%
ggplot() +
geom_mosaic(
aes(
x = product(deprx, marital_status),
fill = deprx
),
offset = 0.01,
show.legend = FALSE
)+
labs(
x="",
y=""
)+
theme(
axis.text.y = element_blank(),
axis.ticks.y=element_blank(),
axis.text.x = element_text(angle = 90, hjust = 1)
)
cvddiag_deprx =
mosi_regression_df %>%
ggplot() +
geom_mosaic(
aes(
x = product(deprx, cvddiag),
fill = deprx
),
offset = 0.01,
show.legend = FALSE
)+
labs(
x="",
y=""
)+
theme(
axis.text.y = element_blank(),
axis.ticks.y=element_blank(),
axis.text.x = element_text(angle = 90, hjust = 1),
)
sex_deprx = ggplotly(sex_deprx)
poverty_deprx =
ggplotly(poverty_deprx) %>%
layout(legend = list(orientation = "h"))
marital_status_deprx = ggplotly(marital_status_deprx)
cvddiag_deprx = ggplotly(cvddiag_deprx)
subplot(
style(sex_deprx, showlegend = F),
style(cvddiag_deprx, showlegend = F),
style(marital_status_deprx, showlegend = F),
poverty_deprx, nrows=1) %>%
layout(legend = list(x = 0, y = -0.4))We seek to validate our models and assess its goodness of fit. In our motivating example, we have two models: (1) Crude model and (2) Adjusted model.
We can assess which of these two models fit the data better using the likelihood ratio test.
depre_crude_model =
glm(deprx ~ cvddiag,
family=binomial(link='logit'),
data = regression_df)depre_adjusted_model =
glm(deprx ~ sex + poverty + age + marital_status + cvddiag,
family=binomial(link='logit'),
data = regression_df)lrtest(depre_crude_model, depre_adjusted_model) %>%
kbl(caption = "Likelihood ratio test for crude model and adjusted model", col.names = c("Total df for each model", "LogLik", "difference in df", "Chisq-statistic", "p-value")) %>%
kable_paper("striped", full_width = F) %>%
column_spec(1, bold = T)| Total df for each model | LogLik | difference in df | Chisq-statistic | p-value |
|---|---|---|---|---|
| 2 | -6028.718 | NA | NA | NA |
| 10 | -5790.814 | 8 | 475.8081 | 0 |
Since four of these main effects were categorical variables, so we need to create dummy variables that indicate which levels of the predictors a given individual belonged. The outcome for this logistic regression is binary (Whether taking medication for depression = Yes/No).
Sex has one dummy variable and
Sex_Male is the reference group.poverty has two dummy variables and
poverty_less than 1.0 is the reference group.marital_status has four dummy variables
and marital_status_Married is the reference group.cvddiag(COVID-19 status)has one dummy
variable and cvddiag_Had COVID-19 is the reference
group.depre_adjusted_model %>%
broom::tidy() %>%
mutate(OR = exp(estimate),
Lower_CI = exp(estimate -1.96*std.error),
Upper_CI = exp(estimate +1.96*std.error)) %>%
select(term, OR, Lower_CI, Upper_CI, statistic, p.value) %>%
kbl(caption = "Effect of Selected Predictors on whether taking medication for depression"
, col.names = c("Predictors", "OR", "Lower bound of 95% CI","Upper bound of 95% CI", "t-statistic", "p-value"),
digits= 2) %>%
kable_paper("striped", full_width = F) %>%
column_spec(1, bold = T)| Predictors | OR | Lower bound of 95% CI | Upper bound of 95% CI | t-statistic | p-value |
|---|---|---|---|---|---|
| (Intercept) | 0.07 | 0.06 | 0.10 | -20.43 | 0.00 |
| sex_Female | 2.49 | 2.24 | 2.76 | 16.86 | 0.00 |
| poverty_1.0-2.0 | 0.76 | 0.64 | 0.90 | -3.13 | 0.00 |
| poverty_2.0 and above | 0.59 | 0.51 | 0.69 | -6.91 | 0.00 |
| age | 1.01 | 1.00 | 1.01 | 3.79 | 0.00 |
| marital_status_Widowed | 1.13 | 0.95 | 1.35 | 1.41 | 0.16 |
| marital_status_Divorced | 1.45 | 1.28 | 1.66 | 5.64 | 0.00 |
| marital_status_Separated | 1.31 | 0.92 | 1.87 | 1.47 | 0.14 |
| marital_status_Never married | 1.19 | 1.04 | 1.36 | 2.54 | 0.01 |
| cvddiag_Had COVID-19 | 1.13 | 0.89 | 1.45 | 1.02 | 0.31 |
After adjustment for sex, age, poverty(family income level), and
Current Marital Status, we obtained the p-value for the
variable cvddiage(COVID-19 status)is 0.31 > α=0.05.
Hence, we find no statistically significant association between COVID-19
status and taking medication for depression (aOR: 1.13, 95% CI: 0.89,
1.45). In addition, the findings in the statistical analysis matches
what we get in the depression trend section of our exploratory analysis.
Due to the Table 2 fallacy, we should avoid interpreting covariates
other than the exposure of interest, but they are included here for
completeness.