GAD and PHQ9

Load Packages

library(tidyverse)
library(codebookr)
library(summarytools)
library(broom)
library(performance)
library(gt)
library(gtsummary)
library(janitor)
library(forcats)
library(here)
library(margins)

Import Data

load(file = "~/Desktop/R-Code/SDOH_ALL/mh_vax_co_sub_mm.rda")
View(mh_vax_co_sub_mm)

GAD data exploration

GAD7 dichotomous

mh_vax_co_sub_mm$gad_2 <- cut(mh_vax_co_sub_mm$gad_score,
                              breaks=c(0, 9, 21),
                              labels=c('None_mild', 'mod_severe'))

Baseline Characteristics by GAD7 score

baseline %>% tbl_summary(by = gad_fct,
         label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "English Speaking", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", depression_2 ~ "Depression", anxiety_2 ~ "Anxiety", any_psych_dx_2 ~ "Any Psychiatric Diagnosis", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartile", act_tob ~ "Active Tobacco Use"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p()

14415 observations missing `gad_fct` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `gad_fct` column before passing to `tbl_summary()`.
There was an error in 'add_p()/add_difference()' for variable 'race_5', p-value omitted:
Error in stats::fisher.test(structure(c(2L, 1L, 3L, 1L, 1L, 1L, 1L, 1L, : FEXACT error 501.
The hash table key cannot be computed because the largest key
is larger than the largest representable int.
The algorithm cannot proceed.
Reduce the workspace, consider using 'simulate.p.value=TRUE' or another algorithm.
There was an error in 'add_p()/add_difference()' for variable 'mstat_5', p-value omitted:
Error in stats::fisher.test(structure(c(1L, 1L, 2L, 3L, 2L, 3L, 4L, 3L, : FEXACT error 501.
The hash table key cannot be computed because the largest key
is larger than the largest representable int.
The algorithm cannot proceed.
Reduce the workspace, consider using 'simulate.p.value=TRUE' or another algorithm.

Characteristic	None, N = 2561	Mild, N = 2701	Moderate, N = 2881	Severe, N = 3421	p-value2
Age	50 (19)	48 (17)	43 (17)	42 (15)	<0.001
Gender					0.4
Male	94 (37%)	85 (31%)	88 (31%)	115 (34%)
Female	162 (63%)	185 (69%)	200 (69%)	227 (66%)
Race
White	223 (87%)	230 (85%)	245 (85%)	293 (86%)
Black	16 (6.2%)	24 (8.9%)	19 (6.6%)	27 (7.9%)
Other	6 (2.3%)	7 (2.6%)	15 (5.2%)	12 (3.5%)
Asian	11 (4.3%)	8 (3.0%)	7 (2.4%)	9 (2.6%)
Native	0 (0%)	1 (0.4%)	2 (0.7%)	1 (0.3%)
Ethnicity					0.3
NonHispanic	244 (99%)	258 (97%)	277 (99%)	332 (98%)
UNKNOWN	0 (0%)	0 (0%)	0 (0%)	0 (0%)
CHOOSE NOT TO DISCLOSE	0 (0%)	0 (0%)	0 (0%)	0 (0%)
Hispanic	3 (1.2%)	7 (2.6%)	2 (0.7%)	7 (2.1%)
(Missing)	9	5	9	3
English Speaking					>0.9
English	254 (99%)	269 (100%)	286 (99%)	340 (99%)
Other	2 (0.8%)	1 (0.4%)	2 (0.7%)	2 (0.6%)
Any Religious Affiliation					0.5
Yes	150 (60%)	150 (57%)	163 (58%)	180 (54%)
No	102 (40%)	113 (43%)	119 (42%)	156 (46%)
PATIENT REFUSED	0 (0%)	0 (0%)	0 (0%)	0 (0%)
UNKNOWN	0 (0%)	0 (0%)	0 (0%)	0 (0%)
(Missing)	4	7	6	6
Marital Status
Married	125 (49%)	121 (45%)	98 (34%)	97 (28%)
Unknown	30 (12%)	45 (17%)	34 (12%)	58 (17%)
Unmarried	85 (33%)	83 (31%)	128 (44%)	148 (43%)
DivorcedSeparated	11 (4.3%)	15 (5.6%)	25 (8.7%)	33 (9.6%)
Widow	5 (2.0%)	6 (2.2%)	3 (1.0%)	6 (1.8%)
Active Tobacco Use					<0.001
No	242 (95%)	241 (89%)	252 (88%)	248 (73%)
Yes	13 (5.1%)	29 (11%)	35 (12%)	91 (27%)
NOT ASKED	0 (0%)	0 (0%)	0 (0%)	0 (0%)
(Missing)	1	0	1	3
Depression	97 (38%)	161 (60%)	196 (68%)	247 (72%)	<0.001
Anxiety	113 (44%)	168 (62%)	236 (82%)	284 (83%)	<0.001
ptsd_2	1 (0.4%)	9 (3.3%)	24 (8.3%)	18 (5.3%)	<0.001
Any Psychiatric Diagnosis	160 (62%)	228 (84%)	265 (92%)	323 (94%)	<0.001
gad_score	2 (1)	7 (1)	12 (1)	18 (2)	<0.001
gad_2					<0.001
None_mild	256 (100%)	270 (100%)	0 (0%)	0 (0%)
mod_severe	0 (0%)	0 (0%)	288 (100%)	342 (100%)
phq_fct					<0.001
None	108 (52%)	44 (18%)	20 (7.7%)	8 (2.5%)
Mild	61 (29%)	90 (38%)	50 (19%)	24 (7.5%)
Moderate	24 (12%)	68 (28%)	79 (30%)	63 (20%)
Moderately Severe	11 (5.3%)	21 (8.8%)	77 (30%)	112 (35%)
Severe	3 (1.4%)	16 (6.7%)	35 (13%)	111 (35%)
(Missing)	49	31	27	24
Total SVI	0.28 (0.22)	0.33 (0.26)	0.35 (0.25)	0.34 (0.24)	0.002
(Missing)	3	1	2	5
Soceioeconomic Status	0.24 (0.22)	0.31 (0.26)	0.33 (0.24)	0.32 (0.24)	<0.001
(Missing)	5	3	3	10
Household Composition	0.27 (0.22)	0.31 (0.25)	0.35 (0.27)	0.32 (0.27)	0.018
(Missing)	3	1	2	5
Minority Status and Language	0.56 (0.29)	0.53 (0.29)	0.53 (0.30)	0.56 (0.28)	0.4
(Missing)	3	1	2	5
Housing and Transportation	0.37 (0.28)	0.42 (0.29)	0.42 (0.28)	0.42 (0.28)	0.067
(Missing)	3	2	2	5
SVI Quartile					0.032
First	139 (55%)	116 (43%)	115 (40%)	140 (42%)
Second	71 (28%)	83 (31%)	99 (35%)	118 (35%)
Third	34 (13%)	47 (17%)	49 (17%)	55 (16%)
Fourth	9 (3.6%)	23 (8.6%)	23 (8.0%)	24 (7.1%)
(Missing)	3	1	2	5
1 Mean (SD); n (%)
2 Kruskal-Wallis rank sum test; Pearson's Chi-squared test; Fisher's exact test

Percent anxiety diagnosis by GAD7 score

mh_vax_co_sub_mm %>% 
  dplyr::select(gad_fct, anxiety_2) -> anx_gad
anx_gad %>% 
  tbl_summary(by = gad_fct,
         label = list(anxiety_2 ~ "Anxiety Diagnosis"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p

14415 observations missing `gad_fct` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `gad_fct` column before passing to `tbl_summary()`.

Characteristic	None, N = 2561	Mild, N = 2701	Moderate, N = 2881	Severe, N = 3421	p-value2
Anxiety Diagnosis	113 (44%)	168 (62%)	236 (82%)	284 (83%)	<0.001
1 n (%)
2 Pearson's Chi-squared test

Gad7 by race

mh_vax_co_sub_mm %>% 
  dplyr::select(race_5, gad_fct, anxiety_2) -> gad_race
gad_race %>% 
  tbl_summary(by = race_5,
         label = list(gad_fct ~ "GAD7 Score", anxiety_2 ~ "Anxiety Diagnosis"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p()

There was an error in 'add_p()/add_difference()' for variable 'gad_fct', p-value omitted:
Error in stats::fisher.test(structure(c(NA, NA, NA, NA, NA, NA, NA, NA, : FEXACT error 501.
The hash table key cannot be computed because the largest key
is larger than the largest representable int.
The algorithm cannot proceed.
Reduce the workspace, consider using 'simulate.p.value=TRUE' or another algorithm.

Characteristic	White, N = 13,5021	Black, N = 9871	Other, N = 6511	Asian, N = 3751	Native, N = 561	p-value2
GAD7 Score
None	223 (23%)	16 (19%)	6 (15%)	11 (31%)	0 (0%)
Mild	230 (23%)	24 (28%)	7 (18%)	8 (23%)	1 (25%)
Moderate	245 (25%)	19 (22%)	15 (38%)	7 (20%)	2 (50%)
Severe	293 (30%)	27 (31%)	12 (30%)	9 (26%)	1 (25%)
(Missing)	12,511	901	611	340	52
Anxiety Diagnosis	2,665 (20%)	172 (17%)	111 (17%)	40 (11%)	14 (25%)	<0.001
1 n (%)
2 Pearson's Chi-squared test

Percent with anxiety diagnosis based on GAD7, stratified by race

# White 
white.anx <- subset(mh_vax_co_sub_mm, race_5=="White", select= c(gad_fct, anxiety_2))
  white.anx %>% 
    tbl_summary(by = gad_fct,
         label = list(anxiety_2 ~ "Anxiety Diagnosis"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p

12511 observations missing `gad_fct` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `gad_fct` column before passing to `tbl_summary()`.

Characteristic	None, N = 2231	Mild, N = 2301	Moderate, N = 2451	Severe, N = 2931	p-value2
Anxiety Diagnosis	106 (48%)	151 (66%)	197 (80%)	243 (83%)	<0.001
1 n (%)
2 Pearson's Chi-squared test

  
black.anx <- subset(mh_vax_co_sub_mm, race_5=="Black", select= c(gad_fct, anxiety_2))
  black.anx %>% 
    tbl_summary(by = gad_fct,
         label = list(anxiety_2 ~ "Anxiety Diagnosis"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p

901 observations missing `gad_fct` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `gad_fct` column before passing to `tbl_summary()`.

Characteristic	None, N = 161	Mild, N = 241	Moderate, N = 191	Severe, N = 271	p-value2
Anxiety Diagnosis	3 (19%)	12 (50%)	19 (100%)	23 (85%)	<0.001
1 n (%)
2 Pearson's Chi-squared test

asian.anx <- subset(mh_vax_co_sub_mm, race_5=="Asian", select= c(gad_fct, anxiety_2))
  asian.anx %>% 
    tbl_summary(by = gad_fct,
         label = list(anxiety_2 ~ "Anxiety Diagnosis"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p

340 observations missing `gad_fct` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `gad_fct` column before passing to `tbl_summary()`.

Characteristic	None, N = 111	Mild, N = 81	Moderate, N = 71	Severe, N = 91	p-value2
Anxiety Diagnosis	1 (9.1%)	1 (12%)	5 (71%)	6 (67%)	0.004
1 n (%)
2 Fisher's exact test

Percent with anxiety diagnosis based on GAD7, stratified by SVI Quartile

# First
First.anx <- subset(mh_vax_co_sub_mm, RPL_4=="First", select= c(gad_fct, anxiety_2))
  First.anx %>% 
    tbl_summary(by = gad_fct,
         label = list(anxiety_2 ~ "Anxiety Diagnosis"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p

5392 observations missing `gad_fct` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `gad_fct` column before passing to `tbl_summary()`.

Characteristic	None, N = 1391	Mild, N = 1161	Moderate, N = 1151	Severe, N = 1401	p-value2
Anxiety Diagnosis	62 (45%)	81 (70%)	90 (78%)	117 (84%)	<0.001
1 n (%)
2 Pearson's Chi-squared test

# Second  
second.anx <- subset(mh_vax_co_sub_mm, RPL_4=="Second", select= c(gad_fct, anxiety_2))
  second.anx %>% 
    tbl_summary(by = gad_fct,
         label = list(anxiety_2 ~ "Anxiety Diagnosis"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p

4249 observations missing `gad_fct` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `gad_fct` column before passing to `tbl_summary()`.

Characteristic	None, N = 711	Mild, N = 831	Moderate, N = 991	Severe, N = 1181	p-value2
Anxiety Diagnosis	27 (38%)	53 (64%)	88 (89%)	93 (79%)	<0.001
1 n (%)
2 Pearson's Chi-squared test

# Third 
third.anx <- subset(mh_vax_co_sub_mm, RPL_4=="Third", select= c(gad_fct, anxiety_2))
  third.anx %>% 
    tbl_summary(by = gad_fct,
         label = list(anxiety_2 ~ "Anxiety Diagnosis"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p

2953 observations missing `gad_fct` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `gad_fct` column before passing to `tbl_summary()`.

Characteristic	None, N = 341	Mild, N = 471	Moderate, N = 491	Severe, N = 551	p-value2
Anxiety Diagnosis	21 (62%)	25 (53%)	38 (78%)	49 (89%)	<0.001
1 n (%)
2 Pearson's Chi-squared test

  
# Fourth 
fourth.anx <- subset(mh_vax_co_sub_mm, RPL_4=="Fourth", select= c(gad_fct, anxiety_2))
  fourth.anx %>% 
    tbl_summary(by = gad_fct,
         label = list(anxiety_2 ~ "Anxiety Diagnosis"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p

1544 observations missing `gad_fct` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `gad_fct` column before passing to `tbl_summary()`.

Characteristic	None, N = 91	Mild, N = 231	Moderate, N = 231	Severe, N = 241	p-value2
Anxiety Diagnosis	2 (22%)	8 (35%)	19 (83%)	22 (92%)	<0.001
1 n (%)
2 Fisher's exact test

GAD bivariate analysis

tbl_uv_gad7 <-
  tbl_uvregression(
    mh_vax_co_sub_mm[c("gad_score", "age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "relig_affil", "mstat_5", "RPL_THEMES",  "RPL_THEME1", "RPL_THEME2",  "RPL_THEME3", "RPL_THEME4", "act_tob")],
    method = lm,
    y = gad_score,
    label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "English Speaking", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", act_tob ~ "Active Tobacco Use"))
    
print(tbl_uv_gad7, method = render)

`...` must be empty.
✖ Problematic argument:
• method = render

Characteristic	N	Beta	95% CI1	p-value
Age	1,326	-0.08	-0.10, -0.06	<0.001
Gender	1,326
Male		—	—
Female		1.2	0.49, 2.0	0.001
Race	1,326
White		—	—
Black		0.24	-1.1, 1.6	0.7
Other		-0.19	-2.1, 1.7	0.8
Asian		-1.5	-3.5, 0.55	0.2
Native		-0.96	-6.3, 4.3	0.7
Ethnicity	1,296
NonHispanic		—	—
Hispanic		-1.3	-3.8, 1.2	0.3
English Speaking	1,326
English		—	—
Other		1.9	-3.0, 6.8	0.4
Any Religious Affiliation	1,299
Yes		—	—
No		0.69	-0.04, 1.4	0.064
Marital Status	1,326
Married		—	—
Unknown		2.0	0.89, 3.1	<0.001
Unmarried		2.5	1.8, 3.3	<0.001
DivorcedSeparated		3.8	2.4, 5.3	<0.001
Widow		0.72	-1.9, 3.3	0.6
Total SVI	1,311	2.5	1.0, 3.9	<0.001
Soceioeconomic Status	1,300	3.3	1.8, 4.7	<0.001
Household Composition	1,311	1.8	0.38, 3.2	0.013
Minority Status and Language	1,311	0.15	-1.1, 1.4	0.8
Housing and Transportation	1,310	1.0	-0.23, 2.3	0.11
Active Tobacco Use	1,321
No		—	—
Yes		3.7	2.7, 4.7	<0.001
1 CI = Confidence Interval

NULL

GAD Multivariable models

GAD + RPL_THEMES

gad1 <- lm(gad_score ~ mstat_5 + lang_3 + age_yrs + race_5 + relig_affil
               + gender + ethnic_3 + RPL_THEMES,
              data = mh_vax_co_sub_mm)
summary(gad1)

Call:
lm(formula = gad_score ~ mstat_5 + lang_3 + age_yrs + race_5 + 
    relig_affil + gender + ethnic_3 + RPL_THEMES, data = mh_vax_co_sub_mm)

Residuals:
    Min      1Q  Median      3Q     Max 
-12.643  -5.401  -0.421   4.898  16.189 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)              10.43408    0.82274  12.682  < 2e-16 ***
mstat_5Unknown            1.17863    0.57797   2.039 0.041636 *  
mstat_5Unmarried          1.16402    0.46422   2.507 0.012287 *  
mstat_5DivorcedSeparated  4.03849    0.74988   5.386 8.64e-08 ***
mstat_5Widow              2.32634    1.35514   1.717 0.086287 .  
lang_3Other               2.92267    2.68158   1.090 0.275965    
age_yrs                  -0.07442    0.01208  -6.159 9.85e-10 ***
race_5Black              -0.86606    0.71209  -1.216 0.224131    
race_5Other               0.31029    1.11518   0.278 0.780872    
race_5Asian              -1.32379    1.05072  -1.260 0.207949    
race_5Native             -0.72793    2.64791  -0.275 0.783433    
relig_affilNo             0.25606    0.38133   0.671 0.502038    
genderFemale              1.28597    0.38349   3.353 0.000822 ***
ethnic_3Hispanic         -1.82718    1.35512  -1.348 0.177792    
RPL_THEMES                1.37734    0.77764   1.771 0.076779 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.364 on 1239 degrees of freedom
  (14317 observations deleted due to missingness)
Multiple R-squared:  0.08462,   Adjusted R-squared:  0.07428 
F-statistic: 8.181 on 14 and 1239 DF,  p-value: < 2.2e-16

broom::glance(gad1)

broom::tidy(gad1)

tbl_regression(gad1, label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "English Speaking", relig_affil ~ "Any Religious Affiliation", RPL_THEMES ~ "Total SVI", mstat_5 ~ "Marital Status"))

Characteristic	Beta	95% CI1	p-value
Marital Status
Married	—	—
Unknown	1.2	0.04, 2.3	0.042
Unmarried	1.2	0.25, 2.1	0.012
DivorcedSeparated	4.0	2.6, 5.5	<0.001
Widow	2.3	-0.33, 5.0	0.086
English Speaking
English	—	—
Other	2.9	-2.3, 8.2	0.3
Age	-0.07	-0.10, -0.05	<0.001
Race
White	—	—
Black	-0.87	-2.3, 0.53	0.2
Other	0.31	-1.9, 2.5	0.8
Asian	-1.3	-3.4, 0.74	0.2
Native	-0.73	-5.9, 4.5	0.8
Any Religious Affiliation
Yes	—	—
No	0.26	-0.49, 1.0	0.5
Gender
Male	—	—
Female	1.3	0.53, 2.0	<0.001
Ethnicity
NonHispanic	—	—
Hispanic	-1.8	-4.5, 0.83	0.2
Total SVI	1.4	-0.15, 2.9	0.077
1 CI = Confidence Interval

## Model performance
model_performance(gad1)

# Indices of model performance

AIC      |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
-------------------------------------------------------
8217.049 | 8299.195 | 0.085 |     0.074 | 6.326 | 6.364

performance::check_model(gad1, panel = TRUE)

## Margins
cplot(gad1, "RPL_THEMES", what = "prediction", main = "Predicted GAD7 Given SVI")

GAD + All Themes

gad2 <- lm(gad_score ~ mstat_5 + lang_3 + age_yrs + race_5 + relig_affil
               + gender + ethnic_3 + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 +
              RPL_THEME4,
              data = mh_vax_co_sub_mm)
summary(gad2)

Call:
lm(formula = gad_score ~ mstat_5 + lang_3 + age_yrs + race_5 + 
    relig_affil + gender + ethnic_3 + RPL_THEME1 + RPL_THEME2 + 
    RPL_THEME3 + RPL_THEME4, data = mh_vax_co_sub_mm)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.8995  -5.3670  -0.3562   4.8486  16.1453 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)              10.35355    0.90660  11.420  < 2e-16 ***
mstat_5Unknown            1.08157    0.58294   1.855 0.063783 .  
mstat_5Unmarried          1.11700    0.46629   2.395 0.016748 *  
mstat_5DivorcedSeparated  3.95798    0.75262   5.259 1.71e-07 ***
mstat_5Widow              2.20476    1.35753   1.624 0.104612    
lang_3Other               2.87694    2.68291   1.072 0.283789    
age_yrs                  -0.07304    0.01218  -5.999 2.60e-09 ***
race_5Black              -0.88481    0.71699  -1.234 0.217412    
race_5Other               0.39751    1.11616   0.356 0.721796    
race_5Asian              -1.18046    1.07010  -1.103 0.270185    
race_5Native             -0.59596    2.65727  -0.224 0.822580    
relig_affilNo             0.24074    0.38394   0.627 0.530762    
genderFemale              1.30129    0.38578   3.373 0.000766 ***
ethnic_3Hispanic         -1.84544    1.35714  -1.360 0.174144    
RPL_THEME1                1.93869    1.03646   1.870 0.061655 .  
RPL_THEME2                0.28442    0.87940   0.323 0.746426    
RPL_THEME3               -0.05825    0.67516  -0.086 0.931259    
RPL_THEME4               -0.38301    0.75179  -0.509 0.610516    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.361 on 1224 degrees of freedom
  (14329 observations deleted due to missingness)
Multiple R-squared:  0.08638,   Adjusted R-squared:  0.0737 
F-statistic: 6.808 on 17 and 1224 DF,  p-value: 7.007e-16

broom::glance(gad2)

broom::tidy(gad2)

tbl_regression(gad2, label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "English Speaking", relig_affil ~ "Any Religious Affiliation", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", mstat_5 ~ "Marital Status"))

Characteristic	Beta	95% CI1	p-value
Marital Status
Married	—	—
Unknown	1.1	-0.06, 2.2	0.064
Unmarried	1.1	0.20, 2.0	0.017
DivorcedSeparated	4.0	2.5, 5.4	<0.001
Widow	2.2	-0.46, 4.9	0.10
English Speaking
English	—	—
Other	2.9	-2.4, 8.1	0.3
Age	-0.07	-0.10, -0.05	<0.001
Race
White	—	—
Black	-0.88	-2.3, 0.52	0.2
Other	0.40	-1.8, 2.6	0.7
Asian	-1.2	-3.3, 0.92	0.3
Native	-0.60	-5.8, 4.6	0.8
Any Religious Affiliation
Yes	—	—
No	0.24	-0.51, 0.99	0.5
Gender
Male	—	—
Female	1.3	0.54, 2.1	<0.001
Ethnicity
NonHispanic	—	—
Hispanic	-1.8	-4.5, 0.82	0.2
Soceioeconomic Status	1.9	-0.09, 4.0	0.062
Household Composition	0.28	-1.4, 2.0	0.7
Minority Status and Language	-0.06	-1.4, 1.3	>0.9
Housing and Transportation	-0.38	-1.9, 1.1	0.6
1 CI = Confidence Interval

## Model performance
model_performance(gad2)

# Indices of model performance

AIC      |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
-------------------------------------------------------
8140.326 | 8237.691 | 0.086 |     0.074 | 6.315 | 6.361

performance::check_model(gad2, panel = TRUE)

## Margins
cplot(gad2, "RPL_THEME1", what = "prediction", main = "Predicted GAD7 Given THEME1")

cplot(gad2, "RPL_THEME2", what = "prediction", main = "Predicted GAD7 Given THEME2")

cplot(gad2, "RPL_THEME3", what = "prediction", main = "Predicted GAD7 Given THEME3")

cplot(gad2, "RPL_THEME4", what = "prediction", main = "Predicted GAD7 Given THEME4")

GAD7 by SVI Quartile

gad3 <- lm(gad_score ~ mstat_5 + lang_3 + age_yrs + race_5 + relig_affil
               + gender + ethnic_3 + RPL_4,
              data = mh_vax_co_sub_mm)
summary(gad3)

Call:
lm(formula = gad_score ~ mstat_5 + lang_3 + age_yrs + race_5 + 
    relig_affil + gender + ethnic_3 + RPL_4, data = mh_vax_co_sub_mm)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.8125  -5.3914  -0.4082   4.8865  17.0775 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)              10.48752    0.81305  12.899  < 2e-16 ***
mstat_5Unknown            1.20583    0.57855   2.084 0.037343 *  
mstat_5Unmarried          1.21350    0.46229   2.625 0.008772 ** 
mstat_5DivorcedSeparated  4.00318    0.74993   5.338 1.12e-07 ***
mstat_5Widow              2.29659    1.35354   1.697 0.090000 .  
lang_3Other               2.91446    2.68198   1.087 0.277388    
age_yrs                  -0.07483    0.01208  -6.196 7.85e-10 ***
race_5Black              -0.57632    0.71133  -0.810 0.417977    
race_5Other               0.28849    1.11422   0.259 0.795742    
race_5Asian              -1.34840    1.05001  -1.284 0.199320    
race_5Native             -1.09918    2.65302  -0.414 0.678715    
relig_affilNo             0.30741    0.38124   0.806 0.420201    
genderFemale              1.28494    0.38316   3.354 0.000822 ***
ethnic_3Hispanic         -1.84638    1.35637  -1.361 0.173678    
RPL_4Second               0.92308    0.42090   2.193 0.028484 *  
RPL_4Third               -0.16925    0.53745  -0.315 0.752883    
RPL_4Fourth               1.32837    0.76580   1.735 0.083058 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.356 on 1237 degrees of freedom
  (14317 observations deleted due to missingness)
Multiple R-squared:  0.08832,   Adjusted R-squared:  0.07653 
F-statistic:  7.49 on 16 and 1237 DF,  p-value: < 2.2e-16

broom::glance(gad3)

broom::tidy(gad3)

tbl_regression(gad3, label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "English Speaking", relig_affil ~ "Any Religious Affiliation", RPL_4 ~ "SVI Quartile", mstat_5 ~ "Marital Status"))

Characteristic	Beta	95% CI1	p-value
Marital Status
Married	—	—
Unknown	1.2	0.07, 2.3	0.037
Unmarried	1.2	0.31, 2.1	0.009
DivorcedSeparated	4.0	2.5, 5.5	<0.001
Widow	2.3	-0.36, 5.0	0.090
English Speaking
English	—	—
Other	2.9	-2.3, 8.2	0.3
Age	-0.07	-0.10, -0.05	<0.001
Race
White	—	—
Black	-0.58	-2.0, 0.82	0.4
Other	0.29	-1.9, 2.5	0.8
Asian	-1.3	-3.4, 0.71	0.2
Native	-1.1	-6.3, 4.1	0.7
Any Religious Affiliation
Yes	—	—
No	0.31	-0.44, 1.1	0.4
Gender
Male	—	—
Female	1.3	0.53, 2.0	<0.001
Ethnicity
NonHispanic	—	—
Hispanic	-1.8	-4.5, 0.81	0.2
SVI Quartile
First	—	—
Second	0.92	0.10, 1.7	0.028
Third	-0.17	-1.2, 0.89	0.8
Fourth	1.3	-0.17, 2.8	0.083
1 CI = Confidence Interval

## Model performance
model_performance(gad3)

# Indices of model performance

AIC      |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
-------------------------------------------------------
8215.975 | 8308.388 | 0.088 |     0.077 | 6.313 | 6.356

performance::check_model(gad3, panel = TRUE)

## Margins
cplot(gad3, "RPL_4", what = "prediction", main = "Predicted GAD7 Given SVI Quartile")

Anxiety model including gad

gad4 <- glm(anxiety_2 ~ mstat_5 + lang_3 + age_yrs + race_5 + relig_affil
               + gender + ethnic_3 + gad_score + RPL_THEMES,
              family = "binomial",
              data = mh_vax_co_sub_mm)
summary(gad4)

Call:
glm(formula = anxiety_2 ~ mstat_5 + lang_3 + age_yrs + race_5 + 
    relig_affil + gender + ethnic_3 + gad_score + RPL_THEMES, 
    family = "binomial", data = mh_vax_co_sub_mm)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.5144  -0.9171   0.4860   0.8266   2.1438  

Coefficients:
                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)              -0.439972   0.313089  -1.405 0.159943    
mstat_5Unknown           -0.176554   0.210421  -0.839 0.401441    
mstat_5Unmarried          0.066967   0.170822   0.392 0.695037    
mstat_5DivorcedSeparated  0.411338   0.301518   1.364 0.172498    
mstat_5Widow             -0.062526   0.488567  -0.128 0.898166    
lang_3Other              -0.339495   0.965871  -0.351 0.725220    
age_yrs                  -0.009383   0.004405  -2.130 0.033165 *  
race_5Black              -0.418127   0.261339  -1.600 0.109611    
race_5Other               0.619118   0.448672   1.380 0.167620    
race_5Asian              -1.532449   0.410458  -3.734 0.000189 ***
race_5Native              1.771725   1.258916   1.407 0.159326    
relig_affilNo            -0.136040   0.141475  -0.962 0.336260    
genderFemale              0.381743   0.141470   2.698 0.006967 ** 
ethnic_3Hispanic         -0.974035   0.507322  -1.920 0.054864 .  
gad_score                 0.161382   0.011898  13.564  < 2e-16 ***
RPL_THEMES               -0.123549   0.290671  -0.425 0.670802    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1652.7  on 1253  degrees of freedom
Residual deviance: 1339.8  on 1238  degrees of freedom
  (14317 observations deleted due to missingness)
AIC: 1371.8

Number of Fisher Scoring iterations: 4

broom::glance(gad4)

broom::tidy(gad4, exponentiate = TRUE)

model_performance(gad4)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
1371.823 | 1453.969 |     0.235 | 0.421 | 1.040 |    0.534 |      -Inf |       8.126e-04 | 0.643

tbl_regression(gad4, label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "English Speaking", relig_affil ~ "Any Religious Affiliation", RPL_THEMES ~ "Total SVI", mstat_5 ~ "Marital Status", gad_score ~ "GAD7 Score"), exponentiate = TRUE)  

Characteristic	OR1	95% CI1	p-value
Marital Status
Married	—	—
Unknown	0.84	0.55, 1.27	0.4
Unmarried	1.07	0.76, 1.49	0.7
DivorcedSeparated	1.51	0.85, 2.77	0.2
Widow	0.94	0.36, 2.49	0.9
English Speaking
English	—	—
Other	0.71	0.10, 4.93	0.7
Age	0.99	0.98, 1.00	0.033
Race
White	—	—
Black	0.66	0.39, 1.10	0.11
Other	1.86	0.79, 4.64	0.2
Asian	0.22	0.09, 0.47	<0.001
Native	5.88	0.67, 137	0.2
Any Religious Affiliation
Yes	—	—
No	0.87	0.66, 1.15	0.3
Gender
Male	—	—
Female	1.46	1.11, 1.93	0.007
Ethnicity
NonHispanic	—	—
Hispanic	0.38	0.14, 1.00	0.055
GAD7 Score	1.18	1.15, 1.20	<0.001
Total SVI	0.88	0.50, 1.56	0.7
1 OR = Odds Ratio, CI = Confidence Interval

## Performance 
model_performance(gad4)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
1371.823 | 1453.969 |     0.235 | 0.421 | 1.040 |    0.534 |      -Inf |       8.126e-04 | 0.643

performance::check_model(gad4, panel = TRUE)

## Margins 
cplot(gad4, "gad_score", what = "prediction", main = "Percent Likelihood of Anxiety Given Gad")

GAD7 logistic (non-mild, mod-severe)

gad5 <- glm(gad_2 ~age_yrs + gender + race_5 + ethnic_3 + lang_3 + mstat_5 + relig_affil + RPL_THEMES,
              family = "binomial",
              data = mh_vax_co_sub_mm)
summary(gad5)
broom::glance(gad5)
broom::tidy(gad5, exponentiate = TRUE)
model_performance(gad5)
tbl_regression(gad5, label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", RPL_THEMES ~ "Total SVI", mstat_5 ~ "Marital Status"), exponentiate = TRUE)  

## Performance 
model_performance(gad5)
performance::check_model(gad5, panel = TRUE)

PHQ9 Data Exploration

PHQ9 dichotomous

mh_vax_co_sub_mm$phq_2 <- cut(mh_vax_co_sub_mm$phq_score,
                              breaks=c(0, 9, 27),
                              labels=c('None_mild', 'mod_severe'))

Baseline Characteristics by PHQ9 score

baseline %>% tbl_summary(by = phq_fct,
         label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "English Speaking", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", depression_2 ~ "Depression", anxiety_2 ~ "Anxiety", any_psych_dx_2 ~ "Any Psychiatric Diagnosis", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", RPL_4 ~ "SVI Quartile", act_tob ~ "Active Tobacco Use"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p()

13568 observations missing `phq_fct` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `phq_fct` column before passing to `tbl_summary()`.
There was an error in 'add_p()/add_difference()' for variable 'race_5', p-value omitted:
Error in stats::fisher.test(structure(c(2L, 1L, 3L, 1L, 1L, 1L, 3L, 1L, : FEXACT error 501.
The hash table key cannot be computed because the largest key
is larger than the largest representable int.
The algorithm cannot proceed.
Reduce the workspace, consider using 'simulate.p.value=TRUE' or another algorithm.
There was an error in 'add_p()/add_difference()' for variable 'act_tob', p-value omitted:
Error in stats::fisher.test(structure(c(1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, : FEXACT error 7(location). LDSTP=18300 is too small for this problem,
  (pastp=280.516, ipn_0:=ipoin[itp=436]=4486, stp[ipn_0]=253.981).
Increase workspace or consider using 'simulate.p.value=TRUE'

Characteristic	None, N = 6701	Mild, N = 4941	Moderate, N = 3591	Moderately Severe, N = 2731	Severe, N = 2071	p-value2
Age	56 (22)	55 (22)	50 (21)	48 (18)	46 (18)	<0.001
Gender						0.011
Male	283 (42%)	187 (38%)	123 (34%)	92 (34%)	65 (31%)
Female	387 (58%)	307 (62%)	236 (66%)	181 (66%)	142 (69%)
Race
White	570 (85%)	414 (84%)	307 (86%)	238 (87%)	173 (84%)
Black	38 (5.7%)	38 (7.7%)	31 (8.6%)	18 (6.6%)	15 (7.2%)
Other	29 (4.3%)	19 (3.8%)	10 (2.8%)	8 (2.9%)	11 (5.3%)
Asian	31 (4.6%)	18 (3.6%)	10 (2.8%)	8 (2.9%)	6 (2.9%)
Native	2 (0.3%)	5 (1.0%)	1 (0.3%)	1 (0.4%)	2 (1.0%)
Ethnicity						0.5
NonHispanic	633 (98%)	474 (99%)	346 (98%)	265 (99%)	195 (97%)
UNKNOWN	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
CHOOSE NOT TO DISCLOSE	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
Hispanic	15 (2.3%)	7 (1.5%)	7 (2.0%)	4 (1.5%)	7 (3.5%)
(Missing)	22	13	6	4	5
English Speaking						0.6
English	660 (99%)	490 (99%)	355 (99%)	272 (100%)	204 (99%)
Other	10 (1.5%)	4 (0.8%)	4 (1.1%)	1 (0.4%)	3 (1.4%)
Any Religious Affiliation						0.4
Yes	409 (62%)	292 (60%)	202 (58%)	152 (57%)	116 (57%)
No	248 (38%)	194 (40%)	149 (42%)	116 (43%)	89 (43%)
PATIENT REFUSED	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
UNKNOWN	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
(Missing)	13	8	8	5	2
Marital Status						<0.001
Married	334 (50%)	201 (41%)	128 (36%)	87 (32%)	60 (29%)
Unknown	105 (16%)	70 (14%)	54 (15%)	48 (18%)	30 (14%)
Unmarried	184 (27%)	173 (35%)	144 (40%)	102 (37%)	86 (42%)
DivorcedSeparated	19 (2.8%)	29 (5.9%)	16 (4.5%)	29 (11%)	22 (11%)
Widow	28 (4.2%)	21 (4.3%)	17 (4.7%)	7 (2.6%)	9 (4.3%)
Active Tobacco Use
No	623 (93%)	448 (91%)	307 (86%)	203 (75%)	162 (79%)
Yes	45 (6.7%)	46 (9.3%)	51 (14%)	68 (25%)	44 (21%)
NOT ASKED	0 (0%)	0 (0%)	0 (0%)	0 (0%)	0 (0%)
(Missing)	2	0	1	2	1
Depression	141 (21%)	221 (45%)	236 (66%)	213 (78%)	187 (90%)	<0.001
Anxiety	187 (28%)	232 (47%)	207 (58%)	200 (73%)	165 (80%)	<0.001
ptsd_2	6 (0.9%)	12 (2.4%)	10 (2.8%)	22 (8.1%)	16 (7.7%)	<0.001
Any Psychiatric Diagnosis	268 (40%)	336 (68%)	298 (83%)	254 (93%)	201 (97%)	<0.001
gad_score	4 (4)	7 (5)	11 (5)	14 (5)	16 (5)	<0.001
(Missing)	427	242	122	51	42
gad_2						<0.001
None_mild	152 (84%)	151 (67%)	92 (39%)	32 (14%)	19 (12%)
mod_severe	28 (16%)	74 (33%)	142 (61%)	189 (86%)	146 (88%)
(Missing)	490	269	125	52	42
gad_fct						<0.001
None	108 (60%)	61 (27%)	24 (10%)	11 (5.0%)	3 (1.8%)
Mild	44 (24%)	90 (40%)	68 (29%)	21 (9.5%)	16 (9.7%)
Moderate	20 (11%)	50 (22%)	79 (34%)	77 (35%)	35 (21%)
Severe	8 (4.4%)	24 (11%)	63 (27%)	112 (51%)	111 (67%)
(Missing)	490	269	125	52	42
Total SVI	0.29 (0.23)	0.32 (0.25)	0.35 (0.24)	0.33 (0.24)	0.42 (0.26)	<0.001
(Missing)	16	5	2	5	1
Soceioeconomic Status	0.25 (0.22)	0.29 (0.25)	0.33 (0.25)	0.31 (0.24)	0.38 (0.26)	<0.001
(Missing)	23	9	6	5	2
Household Composition	0.29 (0.23)	0.30 (0.25)	0.33 (0.25)	0.32 (0.26)	0.41 (0.30)	<0.001
(Missing)	16	5	2	5	1
Minority Status and Language	0.54 (0.30)	0.54 (0.30)	0.53 (0.29)	0.56 (0.28)	0.56 (0.29)	0.5
(Missing)	16	5	2	5	1
Housing and Transportation	0.38 (0.29)	0.40 (0.29)	0.42 (0.28)	0.39 (0.27)	0.49 (0.27)	<0.001
(Missing)	20	5	2	5	1
SVI Quartile						<0.001
First	339 (52%)	234 (48%)	146 (41%)	112 (42%)	66 (32%)
Second	185 (28%)	140 (29%)	120 (34%)	96 (36%)	72 (35%)
Third	105 (16%)	82 (17%)	65 (18%)	42 (16%)	42 (20%)
Fourth	25 (3.8%)	33 (6.7%)	26 (7.3%)	18 (6.7%)	26 (13%)
(Missing)	16	5	2	5	1
1 Mean (SD); n (%)
2 Kruskal-Wallis rank sum test; Pearson's Chi-squared test; Fisher's exact test

Percent depression diagnosis by phq9 score

mh_vax_co_sub_mm %>% 
  dplyr::select(phq_fct, depression_2) -> dep_gad
dep_gad %>% 
  tbl_summary(by = phq_fct,
         label = list(depression_2 ~ "Depression Diagnosis"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p

13568 observations missing `phq_fct` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `phq_fct` column before passing to `tbl_summary()`.

Characteristic	None, N = 6701	Mild, N = 4941	Moderate, N = 3591	Moderately Severe, N = 2731	Severe, N = 2071	p-value2
Depression Diagnosis	141 (21%)	221 (45%)	236 (66%)	213 (78%)	187 (90%)	<0.001
1 n (%)
2 Pearson's Chi-squared test

PHQ9 and depression diagnosis by race

mh_vax_co_sub_mm %>% 
  dplyr::select(race_5, phq_fct, depression_2) -> phq_race
phq_race %>% 
  tbl_summary(by = race_5,
         label = list(phq_fct ~ "PHQ Score Score", depression_2 ~ "Depression Diagnosis"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p()

There was an error in 'add_p()/add_difference()' for variable 'phq_fct', p-value omitted:
Error in stats::fisher.test(structure(c(NA, NA, NA, NA, NA, NA, NA, NA, : FEXACT error 5.
The hash table key cannot be computed because the largest key
is larger than the largest representable int.
The algorithm cannot proceed.
Reduce the workspace, consider using 'simulate.p.value=TRUE' or another algorithm.

Characteristic	White, N = 13,5021	Black, N = 9871	Other, N = 6511	Asian, N = 3751	Native, N = 561	p-value2
PHQ Score Score
None	570 (33%)	38 (27%)	29 (38%)	31 (42%)	2 (18%)
Mild	414 (24%)	38 (27%)	19 (25%)	18 (25%)	5 (45%)
Moderate	307 (18%)	31 (22%)	10 (13%)	10 (14%)	1 (9.1%)
Moderately Severe	238 (14%)	18 (13%)	8 (10%)	8 (11%)	1 (9.1%)
Severe	173 (10%)	15 (11%)	11 (14%)	6 (8.2%)	2 (18%)
(Missing)	11,800	847	574	302	45
Depression Diagnosis	2,463 (18%)	184 (19%)	93 (14%)	36 (9.6%)	13 (23%)	<0.001
1 n (%)
2 Pearson's Chi-squared test

Percent with depression diagnosis based on PHQ9, stratified by race

# White 
white.dep <- subset(mh_vax_co_sub_mm, race_5=="White", select= c(phq_fct, depression_2))
  white.dep %>% 
    tbl_summary(by = phq_fct,
         label = list(depression_2 ~ "Depression Diagnosis"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p

11800 observations missing `phq_fct` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `phq_fct` column before passing to `tbl_summary()`.

Characteristic	None, N = 5701	Mild, N = 4141	Moderate, N = 3071	Moderately Severe, N = 2381	Severe, N = 1731	p-value2
Depression Diagnosis	126 (22%)	198 (48%)	205 (67%)	188 (79%)	158 (91%)	<0.001
1 n (%)
2 Pearson's Chi-squared test

  
# Black
black.dep <- subset(mh_vax_co_sub_mm, race_5=="Black", select= c(phq_fct, depression_2))
  black.dep %>% 
    tbl_summary(by = phq_fct,
         label = list(depression_2 ~ "Depression Diagnosis"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p

847 observations missing `phq_fct` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `phq_fct` column before passing to `tbl_summary()`.

Characteristic	None, N = 381	Mild, N = 381	Moderate, N = 311	Moderately Severe, N = 181	Severe, N = 151	p-value2
Depression Diagnosis	6 (16%)	15 (39%)	19 (61%)	16 (89%)	13 (87%)	<0.001
1 n (%)
2 Pearson's Chi-squared test

# Asian
asian.dep <- subset(mh_vax_co_sub_mm, race_5=="Asian", select= c(phq_fct, depression_2))
  asian.dep %>% 
    tbl_summary(by = phq_fct,
         label = list(depression_2 ~ "Depression Diagnosis"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p

302 observations missing `phq_fct` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `phq_fct` column before passing to `tbl_summary()`.

Characteristic	None, N = 311	Mild, N = 181	Moderate, N = 101	Moderately Severe, N = 81	Severe, N = 61	p-value2
Depression Diagnosis	1 (3.2%)	4 (22%)	6 (60%)	2 (25%)	4 (67%)	<0.001
1 n (%)
2 Fisher's exact test

factor(mh_vax_co_sub_mm$RPL_4)

   [1] Second Second Second Third  Third  First  Second First  Third  Fourth Third  Third  First 
  [14] Second First  Third  Second Fourth Second First  Second First  Second Second Fourth Third 
  [27] Second First  Third  Fourth Second First  First  First  First  Second Third  First  Third 
  [40] Second Second Fourth Second Second Third  Second Third  Fourth Third  Second Third  Second
  [53] Second Second <NA>   Third  Third  First  Second Second Second Second Fourth Second First 
  [66] Third  First  Fourth Second Third  Second Third  First  <NA>   Third  Third  Fourth Second
  [79] First  Second First  First  First  Second Third  First  First  First  First  Second Third 
  [92] Third  <NA>   First  First  Fourth First  First  First  First  First  First  First  Second
 [105] Second Second Second Second Third  Second Second Second First  Third  Third  First  First 
 [118] Second Third  Second Second Third  Second First  Fourth First  Third  <NA>   First  Second
 [131] First  First  Third  Third  Fourth First  First  Fourth First  Second Second Fourth Fourth
 [144] First  Second Third  First  Second <NA>   First  Third  Second First  Third  First  <NA>  
 [157] Third  First  First  First  Second First  Third  First  Second First  First  Third  Third 
 [170] <NA>   Second First  Fourth First  Third  Second Fourth Third  First  Third  Third  Second
 [183] <NA>   Third  Second Third  Fourth Third  <NA>   Third  Second First  Third  Second Second
 [196] Third  Second Fourth First  First  Third  Second Second Second First  Third  Third  First 
 [209] Fourth First  First  Third  First  First  Second First  Third  Third  Fourth First  First 
 [222] Second Second Fourth First  Fourth Third  Second First  Third  Third  First  Third  First 
 [235] First  Second <NA>   First  First  Second Third  First  First  First  First  Third  First 
 [248] First  Second Second First  Third  Third  Second Third  Second First  Second Second First 
 [261] First  Third  Second First  Fourth Second Fourth First  Fourth First  First  Second First 
 [274] First  First  <NA>   Second Second First  First  Second <NA>   First  First  Second Second
 [287] First  Fourth <NA>   Second First  First  Second Second Fourth Third  Third  First  Second
 [300] Fourth Fourth Fourth Second Third  Second First  First  Second Third  First  Third  Second
 [313] Third  Second Third  Third  First  Second Second Second Second Second First  Third  First 
 [326] First  Second First  First  Second First  Third  Third  First  First  First  Second Fourth
 [339] Second First  Fourth Third  Third  Second Fourth First  Second Second First  Third  First 
 [352] Fourth Third  First  <NA>   Second Second Third  First  Second Second Third  First  First 
 [365] Third  Second First  Fourth First  Second Second Second Fourth First  First  Fourth First 
 [378] Fourth Second Fourth Second First  Third  Third  Third  Third  Second First  <NA>   Second
 [391] First  Third  First  Second Fourth Second Third  First  Third  Third  <NA>   First  First 
 [404] First  Fourth Third  Second First  First  Third  First  Second <NA>   First  Third  Second
 [417] Third  Second <NA>   Fourth Fourth Fourth Fourth Second First  First  Third  First  First 
 [430] Second First  Fourth Second Third  First  Second Third  Second Third  Third  Third  First 
 [443] Third  First  Fourth Fourth Second Fourth Second First  Second First  Third  Fourth First 
 [456] Third  Fourth Fourth First  Third  Second Third  First  Third  Third  Third  Third  Second
 [469] Second Fourth Fourth First  <NA>   First  Fourth First  Third  First  Second Fourth Third 
 [482] First  First  First  Second First  Second Fourth First  Second First  Second First  Third 
 [495] First  First  Fourth Second Fourth Second First  First  Fourth Second First  Second Third 
 [508] Fourth Fourth First  Second Third  Fourth <NA>   Second First  First  First  First  Second
 [521] Second Second Second Third  First  Fourth First  Third  Second Third  Third  First  Second
 [534] Second First  Third  Second Fourth Third  Second First  First  First  Second First  First 
 [547] First  First  First  Second First  First  First  First  Second Second First  First  Second
 [560] Third  Fourth First  Second Second Second First  Second Second Second Fourth Fourth First 
 [573] First  <NA>   First  Second First  First  First  Third  Second Fourth First  First  First 
 [586] Second Fourth First  Fourth Second Second Second Fourth Fourth First  Second First  First 
 [599] First  Third  First  Third  Second Third  First  First  Third  First  First  First  Fourth
 [612] Second First  Fourth Fourth Fourth Fourth Fourth Fourth Fourth Second Second First  Second
 [625] Third  First  First  First  First  First  First  Fourth Second Second Third  <NA>   Fourth
 [638] Fourth Fourth Fourth Fourth Fourth Fourth Second Third  Third  First  Second First  Third 
 [651] First  First  Fourth Second First  First  Fourth First  First  First  Fourth Fourth First 
 [664] Second Second First  Third  Third  Second First  First  First  First  Fourth First  First 
 [677] First  Third  Second Third  Second First  Fourth Second First  Third  Second First  Third 
 [690] Third  First  Third  Third  Second Second Second Second Third  Second First  First  <NA>  
 [703] Second Second First  Fourth Third  First  Third  First  Second Fourth Second First  Second
 [716] First  Third  Second Fourth First  First  Fourth Second First  Second First  First  Third 
 [729] Third  First  Second First  Second First  Third  First  Fourth Second First  Fourth Second
 [742] First  Third  Third  Second Second Second Second First  First  First  Second Third  Fourth
 [755] First  Second Third  First  First  Fourth First  Third  First  Second First  First  <NA>  
 [768] Second Third  Second Second First  Fourth First  Second First  First  Third  First  Second
 [781] Third  Third  Third  First  Second First  Second Third  First  Third  Second Second Fourth
 [794] Third  First  First  Second First  Third  First  Third  Second First  Third  Second First 
 [807] Second First  Third  First  First  First  First  First  Second First  First  Third  Third 
 [820] Second Third  Third  Second Second First  Fourth First  First  First  Second First  <NA>  
 [833] Third  First  Second Second First  Second Second Third  Third  Second Second Second Second
 [846] Third  First  Second Second First  Second Third  Fourth First  First  Third  First  First 
 [859] Third  Second First  First  Third  Fourth Second First  Second Second First  First  Second
 [872] Third  Fourth First  Second First  First  Fourth First  Fourth Fourth Third  First  First 
 [885] Fourth First  First  Third  Second First  Third  Fourth First  Third  Third  Second First 
 [898] First  First  <NA>   Third  First  Third  First  Third  First  Second Third  Fourth Second
 [911] Second Third  Second Third  First  Third  Second Second First  Second Fourth <NA>   Fourth
 [924] Fourth Second Second First  Second Second Fourth First  First  Second Third  First  Second
 [937] First  First  First  First  First  First  Third  Second Second Third  Third  Third  Third 
 [950] First  First  Fourth First  First  Third  Third  First  Third  Second Third  Fourth Third 
 [963] First  Third  Fourth Fourth First  Fourth First  Third  First  Third  Third  First  Fourth
 [976] Second First  First  Third  Second Second First  Third  Second Third  Second Fourth Second
 [989] Third  Third  Second Third  Fourth First  Fourth First  First  Second First  Third 
 [ reached getOption("max.print") -- omitted 14571 entries ]
Levels: First Second Third Fourth

Percent with depression diagnosis based on PHQ9, stratified by SVI Quartile

# First
First.dep <- subset(mh_vax_co_sub_mm, RPL_4=="First", select= c(phq_fct, depression_2))
  First.dep %>% 
    tbl_summary(by = phq_fct,
         label = list(depression_2 ~ "Depression Diagnosis"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p

5005 observations missing `phq_fct` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `phq_fct` column before passing to `tbl_summary()`.

Characteristic	None, N = 3391	Mild, N = 2341	Moderate, N = 1461	Moderately Severe, N = 1121	Severe, N = 661	p-value2
Depression Diagnosis	70 (21%)	92 (39%)	91 (62%)	89 (79%)	58 (88%)	<0.001
1 n (%)
2 Pearson's Chi-squared test

# Second  
second.dep <- subset(mh_vax_co_sub_mm, RPL_4=="Second", select= c(phq_fct, depression_2))
  second.dep %>% 
    tbl_summary(by = phq_fct,
         label = list(depression_2 ~ "Depression Diagnosis"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p

4007 observations missing `phq_fct` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `phq_fct` column before passing to `tbl_summary()`.

Characteristic	None, N = 1851	Mild, N = 1401	Moderate, N = 1201	Moderately Severe, N = 961	Severe, N = 721	p-value2
Depression Diagnosis	42 (23%)	70 (50%)	84 (70%)	76 (79%)	66 (92%)	<0.001
1 n (%)
2 Pearson's Chi-squared test

# Third 
third.dep <- subset(mh_vax_co_sub_mm, RPL_4=="Third", select= c(phq_fct, depression_2))
  third.dep %>% 
    tbl_summary(by = phq_fct,
         label = list(depression_2 ~ "Depression Diagnosis"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p

2802 observations missing `phq_fct` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `phq_fct` column before passing to `tbl_summary()`.

Characteristic	None, N = 1051	Mild, N = 821	Moderate, N = 651	Moderately Severe, N = 421	Severe, N = 421	p-value2
Depression Diagnosis	23 (22%)	39 (48%)	44 (68%)	29 (69%)	39 (93%)	<0.001
1 n (%)
2 Pearson's Chi-squared test

  
# Fourth 
fourth.dep <- subset(mh_vax_co_sub_mm, RPL_4=="Fourth", select= c(phq_fct, depression_2))
  fourth.dep %>% 
    tbl_summary(by = phq_fct,
         label = list(depression_2 ~ "Depression Diagnosis"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p

1495 observations missing `phq_fct` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `phq_fct` column before passing to `tbl_summary()`.

Characteristic	None, N = 251	Mild, N = 331	Moderate, N = 261	Moderately Severe, N = 181	Severe, N = 261	p-value2
Depression Diagnosis	3 (12%)	17 (52%)	17 (65%)	16 (89%)	24 (92%)	<0.001
1 n (%)
2 Pearson's Chi-squared test

PHQ9 bivariate analysis

tbl_uv_phq9 <-
  tbl_uvregression(
    mh_vax_co_sub_mm[c("phq_score", "age_yrs", "gender", "race_5", "ethnic_3", "lang_3", "relig_affil", "mstat_5", "RPL_THEMES",  "RPL_THEME1", "RPL_THEME2",  "RPL_THEME3", "RPL_THEME4", "act_tob")],
    method = lm,
    y = phq_score,
    label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "English Speaking", relig_affil ~ "Any Religious Affiliation", mstat_5 ~ "Marital Status", RPL_THEMES ~ "Total SVI", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", act_tob ~ "Active Tobacco Use"))
    
print(tbl_uv_phq9, method = render)

`...` must be empty.
✖ Problematic argument:
• method = render

Characteristic	N	Beta	95% CI1	p-value
Age	2,340	-0.05	-0.06, -0.04	<0.001
Gender	2,340
Male		—	—
Female		1.8	1.2, 2.4	<0.001
Race	2,340
White		—	—
Black		1.2	0.00, 2.3	0.049
Other		0.19	-1.3, 1.7	0.8
Asian		-1.7	-3.2, -0.23	0.023
Native		1.6	-2.4, 5.7	0.4
Ethnicity	2,283
NonHispanic		—	—
Hispanic		-0.07	-2.1, 1.9	>0.9
English Speaking	2,340
English		—	—
Other		-0.77	-3.6, 2.0	0.6
Any Religious Affiliation	2,286
Yes		—	—
No		0.48	-0.11, 1.1	0.11
Marital Status	2,340
Married		—	—
Unknown		1.4	0.52, 2.2	0.001
Unmarried		2.5	1.8, 3.1	<0.001
DivorcedSeparated		5.1	3.8, 6.4	<0.001
Widow		1.5	0.06, 3.0	0.041
Total SVI	2,298	5.1	3.9, 6.3	<0.001
Soceioeconomic Status	2,278	5.3	4.1, 6.5	<0.001
Household Composition	2,298	3.5	2.3, 4.6	<0.001
Minority Status and Language	2,298	0.90	-0.10, 1.9	0.077
Housing and Transportation	2,293	2.6	1.6, 3.6	<0.001
Active Tobacco Use	2,333
No		—	—
Yes		4.3	3.4, 5.2	<0.001
1 CI = Confidence Interval

NULL

PHQ9 Multivariable models

PHQ9 + RPL_THEMES

phq1 <- lm(phq_score ~ mstat_5 + lang_3 + age_yrs + race_5 + relig_affil
               + gender + ethnic_3 + RPL_THEMES,
              data = mh_vax_co_sub_mm)
summary(phq1)

Call:
lm(formula = phq_score ~ mstat_5 + lang_3 + age_yrs + race_5 + 
    relig_affil + gender + ethnic_3 + RPL_THEMES, data = mh_vax_co_sub_mm)

Residuals:
    Min      1Q  Median      3Q     Max 
-13.163  -5.228  -1.716   4.547  20.668 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)               7.346296   0.653750  11.237  < 2e-16 ***
mstat_5Unknown            0.729094   0.451408   1.615  0.10642    
mstat_5Unmarried          1.081466   0.387857   2.788  0.00534 ** 
mstat_5DivorcedSeparated  4.461473   0.666677   6.692 2.79e-11 ***
mstat_5Widow              1.701267   0.765843   2.221  0.02642 *  
lang_3Other               0.184141   1.495478   0.123  0.90201    
age_yrs                  -0.044568   0.008403  -5.304 1.25e-07 ***
race_5Black              -0.755796   0.609172  -1.241  0.21485    
race_5Other              -0.168926   0.880318  -0.192  0.84784    
race_5Asian              -1.866076   0.787756  -2.369  0.01793 *  
race_5Native              0.982270   1.995110   0.492  0.62253    
relig_affilNo             0.302583   0.309042   0.979  0.32764    
genderFemale              1.543198   0.299749   5.148 2.86e-07 ***
ethnic_3Hispanic         -0.169237   1.074401  -0.158  0.87485    
RPL_THEMES                4.217378   0.632962   6.663 3.39e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.767 on 2174 degrees of freedom
  (13382 observations deleted due to missingness)
Multiple R-squared:  0.08711,   Adjusted R-squared:  0.08123 
F-statistic: 14.82 on 14 and 2174 DF,  p-value: < 2.2e-16

broom::glance(phq1)

broom::tidy(phq1)

tbl_regression(phq1, label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "English Speaking", relig_affil ~ "Any Religious Affiliation", RPL_THEMES ~ "Total SVI", mstat_5 ~ "Marital Status"))

Characteristic	Beta	95% CI1	p-value
Marital Status
Married	—	—
Unknown	0.73	-0.16, 1.6	0.11
Unmarried	1.1	0.32, 1.8	0.005
DivorcedSeparated	4.5	3.2, 5.8	<0.001
Widow	1.7	0.20, 3.2	0.026
English Speaking
English	—	—
Other	0.18	-2.7, 3.1	>0.9
Age	-0.04	-0.06, -0.03	<0.001
Race
White	—	—
Black	-0.76	-2.0, 0.44	0.2
Other	-0.17	-1.9, 1.6	0.8
Asian	-1.9	-3.4, -0.32	0.018
Native	0.98	-2.9, 4.9	0.6
Any Religious Affiliation
Yes	—	—
No	0.30	-0.30, 0.91	0.3
Gender
Male	—	—
Female	1.5	0.96, 2.1	<0.001
Ethnicity
NonHispanic	—	—
Hispanic	-0.17	-2.3, 1.9	0.9
Total SVI	4.2	3.0, 5.5	<0.001
1 CI = Confidence Interval

## Model performance
model_performance(phq1)

# Indices of model performance

AIC       |       BIC |    R2 | R2 (adj.) |  RMSE | Sigma
---------------------------------------------------------
14600.180 | 14691.239 | 0.087 |     0.081 | 6.744 | 6.767

performance::check_model(phq1, panel = TRUE)

## Margins
cplot(phq1, "RPL_THEMES", what = "prediction", main = "Predicted PHQ9 Given SVI")

PHQ + All themes

phq2 <- lm(phq_score ~ mstat_5 + lang_3 + age_yrs + race_5 + relig_affil
               + gender + ethnic_3 + RPL_THEME1 + RPL_THEME2 + RPL_THEME3 +
              RPL_THEME4,
              data = mh_vax_co_sub_mm)
summary(phq2)

Call:
lm(formula = phq_score ~ mstat_5 + lang_3 + age_yrs + race_5 + 
    relig_affil + gender + ethnic_3 + RPL_THEME1 + RPL_THEME2 + 
    RPL_THEME3 + RPL_THEME4, data = mh_vax_co_sub_mm)

Residuals:
    Min      1Q  Median      3Q     Max 
-13.158  -5.214  -1.693   4.570  20.442 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)               6.95264    0.72137   9.638  < 2e-16 ***
mstat_5Unknown            0.67366    0.45482   1.481  0.13872    
mstat_5Unmarried          1.08987    0.39090   2.788  0.00535 ** 
mstat_5DivorcedSeparated  4.48719    0.67153   6.682 2.99e-11 ***
mstat_5Widow              1.61148    0.76760   2.099  0.03590 *  
lang_3Other               0.17111    1.53329   0.112  0.91115    
age_yrs                  -0.04322    0.00847  -5.103 3.63e-07 ***
race_5Black              -0.83281    0.61564  -1.353  0.17627    
race_5Other              -0.01742    0.88498  -0.020  0.98430    
race_5Asian              -1.67504    0.80658  -2.077  0.03795 *  
race_5Native              0.89377    2.00000   0.447  0.65500    
relig_affilNo             0.33111    0.31214   1.061  0.28892    
genderFemale              1.52787    0.30192   5.060 4.54e-07 ***
ethnic_3Hispanic         -0.11019    1.07664  -0.102  0.91850    
RPL_THEME1                3.40973    0.84813   4.020 6.01e-05 ***
RPL_THEME2                1.39362    0.71640   1.945  0.05187 .  
RPL_THEME3                0.40477    0.54076   0.749  0.45422    
RPL_THEME4                0.13711    0.60698   0.226  0.82131    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.774 on 2146 degrees of freedom
  (13407 observations deleted due to missingness)
Multiple R-squared:  0.08861,   Adjusted R-squared:  0.08139 
F-statistic: 12.27 on 17 and 2146 DF,  p-value: < 2.2e-16

broom::glance(phq2)

broom::tidy(phq2)

tbl_regression(phq2, label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "English Speaking", relig_affil ~ "Any Religious Affiliation", RPL_THEME1 ~ "Soceioeconomic Status", RPL_THEME2 ~ "Household Composition", RPL_THEME3 ~ "Minority Status and Language", RPL_THEME4 ~ "Housing and Transportation", mstat_5 ~ "Marital Status"))

Characteristic	Beta	95% CI1	p-value
Marital Status
Married	—	—
Unknown	0.67	-0.22, 1.6	0.14
Unmarried	1.1	0.32, 1.9	0.005
DivorcedSeparated	4.5	3.2, 5.8	<0.001
Widow	1.6	0.11, 3.1	0.036
English Speaking
English	—	—
Other	0.17	-2.8, 3.2	>0.9
Age	-0.04	-0.06, -0.03	<0.001
Race
White	—	—
Black	-0.83	-2.0, 0.37	0.2
Other	-0.02	-1.8, 1.7	>0.9
Asian	-1.7	-3.3, -0.09	0.038
Native	0.89	-3.0, 4.8	0.7
Any Religious Affiliation
Yes	—	—
No	0.33	-0.28, 0.94	0.3
Gender
Male	—	—
Female	1.5	0.94, 2.1	<0.001
Ethnicity
NonHispanic	—	—
Hispanic	-0.11	-2.2, 2.0	>0.9
Soceioeconomic Status	3.4	1.7, 5.1	<0.001
Household Composition	1.4	-0.01, 2.8	0.052
Minority Status and Language	0.40	-0.66, 1.5	0.5
Housing and Transportation	0.14	-1.1, 1.3	0.8
1 CI = Confidence Interval

## Model performance
model_performance(phq2)

# Indices of model performance

AIC       |       BIC |    R2 | R2 (adj.) |  RMSE | Sigma
---------------------------------------------------------
14441.005 | 14548.919 | 0.089 |     0.081 | 6.746 | 6.774

performance::check_model(phq2, panel = TRUE)

## Margins
cplot(phq2, "RPL_THEME1", what = "prediction", main = "Predicted PHQ9 Given THEME1")

cplot(phq2, "RPL_THEME2", what = "prediction", main = "Predicted PHQ9 Given THEME2")

cplot(phq2, "RPL_THEME3", what = "prediction", main = "Predicted PHQ9 Given THEME3")

cplot(phq2, "RPL_THEME4", what = "prediction", main = "Predicted PHQ9 Given THEME4")

PHQ9 by SVI Quartile

phq3 <- lm(phq_score ~ mstat_5 + lang_3 + age_yrs + race_5 + relig_affil
               + gender + ethnic_3 + RPL_4,
              data = mh_vax_co_sub_mm)
summary(phq3)

Call:
lm(formula = phq_score ~ mstat_5 + lang_3 + age_yrs + race_5 + 
    relig_affil + gender + ethnic_3 + RPL_4, data = mh_vax_co_sub_mm)

Residuals:
    Min      1Q  Median      3Q     Max 
-12.339  -5.258  -1.702   4.495  21.083 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)               7.766282   0.644306  12.054  < 2e-16 ***
mstat_5Unknown            0.810313   0.453053   1.789  0.07382 .  
mstat_5Unmarried          1.160919   0.388606   2.987  0.00285 ** 
mstat_5DivorcedSeparated  4.585023   0.669272   6.851 9.53e-12 ***
mstat_5Widow              1.754147   0.767431   2.286  0.02237 *  
lang_3Other               0.247023   1.500234   0.165  0.86923    
age_yrs                  -0.044317   0.008419  -5.264 1.55e-07 ***
race_5Black              -0.570699   0.610981  -0.934  0.35037    
race_5Other              -0.118773   0.881508  -0.135  0.89283    
race_5Asian              -1.884071   0.790407  -2.384  0.01723 *  
race_5Native              1.084526   1.999640   0.542  0.58763    
relig_affilNo             0.331700   0.309731   1.071  0.28432    
genderFemale              1.540919   0.300247   5.132 3.12e-07 ***
ethnic_3Hispanic         -0.173384   1.075912  -0.161  0.87199    
RPL_4Second               1.307218   0.339261   3.853  0.00012 ***
RPL_4Third                1.221747   0.428423   2.852  0.00439 ** 
RPL_4Fourth               3.756220   0.657335   5.714 1.25e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.776 on 2172 degrees of freedom
  (13382 observations deleted due to missingness)
Multiple R-squared:  0.08563,   Adjusted R-squared:  0.07889 
F-statistic: 12.71 on 16 and 2172 DF,  p-value: < 2.2e-16

broom::glance(phq3)

broom::tidy(phq3)

tbl_regression(phq3, label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "English Speaking", relig_affil ~ "Any Religious Affiliation", RPL_4 ~ "SVI Quartile", mstat_5 ~ "Marital Status"))

Characteristic	Beta	95% CI1	p-value
Marital Status
Married	—	—
Unknown	0.81	-0.08, 1.7	0.074
Unmarried	1.2	0.40, 1.9	0.003
DivorcedSeparated	4.6	3.3, 5.9	<0.001
Widow	1.8	0.25, 3.3	0.022
English Speaking
English	—	—
Other	0.25	-2.7, 3.2	0.9
Age	-0.04	-0.06, -0.03	<0.001
Race
White	—	—
Black	-0.57	-1.8, 0.63	0.4
Other	-0.12	-1.8, 1.6	0.9
Asian	-1.9	-3.4, -0.33	0.017
Native	1.1	-2.8, 5.0	0.6
Any Religious Affiliation
Yes	—	—
No	0.33	-0.28, 0.94	0.3
Gender
Male	—	—
Female	1.5	0.95, 2.1	<0.001
Ethnicity
NonHispanic	—	—
Hispanic	-0.17	-2.3, 1.9	0.9
SVI Quartile
First	—	—
Second	1.3	0.64, 2.0	<0.001
Third	1.2	0.38, 2.1	0.004
Fourth	3.8	2.5, 5.0	<0.001
1 CI = Confidence Interval

## Model performance
model_performance(phq3)

# Indices of model performance

AIC       |       BIC |    R2 | R2 (adj.) |  RMSE | Sigma
---------------------------------------------------------
14607.730 | 14710.171 | 0.086 |     0.079 | 6.749 | 6.776

performance::check_model(phq3, panel = TRUE)

## Margins
cplot(phq3, "RPL_4", what = "prediction", main = "Predicted PHQ9 Given SVI Quartile")

Depression model including PHQ

phq4 <- glm(depression_2 ~ mstat_5 + lang_3 + age_yrs + race_5 + relig_affil
               + gender + ethnic_3 + phq_score + RPL_THEMES,
              family = "binomial",
              data = mh_vax_co_sub_mm)
summary(phq4)

Call:
glm(formula = depression_2 ~ mstat_5 + lang_3 + age_yrs + race_5 + 
    relig_affil + gender + ethnic_3 + phq_score + RPL_THEMES, 
    family = "binomial", data = mh_vax_co_sub_mm)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.6204  -0.7580  -0.5150   0.8024   2.6245  

Coefficients:
                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)              -2.043277   0.243010  -8.408  < 2e-16 ***
mstat_5Unknown           -0.091251   0.163241  -0.559  0.57616    
mstat_5Unmarried          0.064883   0.138343   0.469  0.63907    
mstat_5DivorcedSeparated  0.222710   0.247903   0.898  0.36899    
mstat_5Widow             -0.326932   0.272607  -1.199  0.23042    
lang_3Other               0.605330   0.568774   1.064  0.28721    
age_yrs                   0.001931   0.002996   0.645  0.51921    
race_5Black              -0.364582   0.215549  -1.691  0.09076 .  
race_5Other              -0.438687   0.335453  -1.308  0.19096    
race_5Asian              -1.449504   0.348870  -4.155 3.25e-05 ***
race_5Native              0.193995   0.717037   0.271  0.78674    
relig_affilNo            -0.102208   0.111344  -0.918  0.35865    
genderFemale              0.284843   0.108454   2.626  0.00863 ** 
ethnic_3Hispanic         -0.641765   0.422532  -1.519  0.12880    
phq_score                 0.200864   0.009491  21.164  < 2e-16 ***
RPL_THEMES                0.359290   0.230273   1.560  0.11869    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 3013.8  on 2188  degrees of freedom
Residual deviance: 2242.6  on 2173  degrees of freedom
  (13382 observations deleted due to missingness)
AIC: 2274.6

Number of Fisher Scoring iterations: 4

broom::glance(phq4)

broom::tidy(phq4, exponentiate = TRUE)

model_performance(phq4)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
2274.575 | 2365.634 |     0.318 | 0.411 | 1.016 |    0.512 |      -Inf |       6.657e-04 | 0.662

tbl_regression(phq4, label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "English Speaking", relig_affil ~ "Any Religious Affiliation", RPL_THEMES ~ "Total SVI", mstat_5 ~ "Marital Status", phq_score ~ "PHQ9 Score"), exponentiate = TRUE)  

Characteristic	OR1	95% CI1	p-value
Marital Status
Married	—	—
Unknown	0.91	0.66, 1.26	0.6
Unmarried	1.07	0.81, 1.40	0.6
DivorcedSeparated	1.25	0.77, 2.04	0.4
Widow	0.72	0.42, 1.22	0.2
English Speaking
English	—	—
Other	1.83	0.58, 5.46	0.3
Age	1.00	1.00, 1.01	0.5
Race
White	—	—
Black	0.69	0.45, 1.06	0.091
Other	0.64	0.33, 1.23	0.2
Asian	0.23	0.12, 0.45	<0.001
Native	1.21	0.28, 4.90	0.8
Any Religious Affiliation
Yes	—	—
No	0.90	0.73, 1.12	0.4
Gender
Male	—	—
Female	1.33	1.08, 1.65	0.009
Ethnicity
NonHispanic	—	—
Hispanic	0.53	0.22, 1.18	0.13
PHQ9 Score	1.22	1.20, 1.25	<0.001
Total SVI	1.43	0.91, 2.25	0.12
1 OR = Odds Ratio, CI = Confidence Interval

## Performance 
model_performance(phq4)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
2274.575 | 2365.634 |     0.318 | 0.411 | 1.016 |    0.512 |      -Inf |       6.657e-04 | 0.662

performance::check_model(phq4, panel = TRUE)

## Margins 
cplot(phq4, "phq_score", what = "prediction", main = "Percent Likelihood of Depression Given PHQ9")

PHQ9 logistic regression (none-mild, moderate-severe)

phq5 <- glm(phq_2 ~age_yrs + gender + race_5 + ethnic_3 + lang_3 + mstat_5 + relig_affil + RPL_THEMES,
              family = "binomial",
              data = mh_vax_co_sub_mm)
summary(phq5)

Call:
glm(formula = phq_2 ~ age_yrs + gender + race_5 + ethnic_3 + 
    lang_3 + mstat_5 + relig_affil + RPL_THEMES, family = "binomial", 
    data = mh_vax_co_sub_mm)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6882  -1.0252  -0.7853   1.1814   1.8728  

Coefficients:
                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)              -0.094071   0.216594  -0.434  0.66405    
age_yrs                  -0.016547   0.002811  -5.887 3.92e-09 ***
genderFemale              0.305373   0.101716   3.002  0.00268 ** 
race_5Black              -0.328907   0.195855  -1.679  0.09309 .  
race_5Other              -0.422853   0.295996  -1.429  0.15313    
race_5Asian              -0.442551   0.275703  -1.605  0.10846    
race_5Native             -0.522829   0.656620  -0.796  0.42589    
ethnic_3Hispanic          0.229178   0.356853   0.642  0.52073    
lang_3Other               0.152941   0.512041   0.299  0.76518    
mstat_5Unknown            0.208070   0.150212   1.385  0.16600    
mstat_5Unmarried          0.215968   0.126776   1.704  0.08847 .  
mstat_5DivorcedSeparated  0.905666   0.211025   4.292 1.77e-05 ***
mstat_5Widow              0.476371   0.250399   1.902  0.05711 .  
relig_affilNo             0.055918   0.102965   0.543  0.58707    
RPL_THEMES                0.917973   0.206723   4.441 8.97e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2574.9  on 1888  degrees of freedom
Residual deviance: 2455.7  on 1874  degrees of freedom
  (13682 observations deleted due to missingness)
AIC: 2485.7

Number of Fisher Scoring iterations: 4

broom::glance(phq5)

broom::tidy(phq5, exponentiate = TRUE)

model_performance(phq5)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
2485.742 | 2568.899 |     0.062 | 0.479 | 1.145 |    0.650 |      -Inf |       5.326e-04 | 0.542

tbl_regression(phq5, label = list(age_yrs ~ "Age", gender~ "Gender", race_5 ~ "Race", ethnic_3 ~ "Ethnicity", lang_3 ~ "Preferred Language", relig_affil ~ "Any Religious Affiliation", RPL_THEMES ~ "Total SVI", mstat_5 ~ "Marital Status"), exponentiate = TRUE)  

Characteristic	OR1	95% CI1	p-value
Age	0.98	0.98, 0.99	<0.001
Gender
Male	—	—
Female	1.36	1.11, 1.66	0.003
Race
White	—	—
Black	0.72	0.49, 1.05	0.093
Other	0.66	0.36, 1.16	0.2
Asian	0.64	0.37, 1.09	0.11
Native	0.59	0.15, 2.07	0.4
Ethnicity
NonHispanic	—	—
Hispanic	1.26	0.62, 2.54	0.5
Preferred Language
English	—	—
Other	1.17	0.41, 3.12	0.8
Marital Status
Married	—	—
Unknown	1.23	0.92, 1.65	0.2
Unmarried	1.24	0.97, 1.59	0.088
DivorcedSeparated	2.47	1.64, 3.76	<0.001
Widow	1.61	0.98, 2.62	0.057
Any Religious Affiliation
Yes	—	—
No	1.06	0.86, 1.29	0.6
Total SVI	2.50	1.67, 3.76	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

## Performance 
model_performance(phq5)

# Indices of model performance

AIC      |      BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
------------------------------------------------------------------------------------------------
2485.742 | 2568.899 |     0.062 | 0.479 | 1.145 |    0.650 |      -Inf |       5.326e-04 | 0.542

performance::check_model(phq5, panel = TRUE)