dag { bb=“0,0,1,1” “HIV test” [pos=“0.397,0.910”] “Health Insurance” [pos=“0.423,0.265”] “Heavy drinking” [pos=“0.704,0.764”] “Number of Sex_Partners” [pos=“0.115,0.692”] “Person Healthcare Provider” [pos=“0.570,0.612”] “Prep Awareness” [outcome,pos=“0.827,0.450”] “Race/ethnicity” [pos=“0.581,0.074”] “Sexual Identity” [exposure,pos=“0.101,0.452”] “sexually active” [pos=“0.332,0.692”] Age [pos=“0.399,0.068”] Birthsex [pos=“0.227,0.129”] condom [pos=“0.510,0.750”] education [pos=“0.706,0.216”] “HIV test” -> “Prep Awareness” “Health Insurance” -> “Prep Awareness” “Heavy drinking” -> “Prep Awareness” “Number of Sex_Partners” -> “sexually active” “Person Healthcare Provider” -> “Prep Awareness” “Race/ethnicity” -> “Prep Awareness” “Race/ethnicity” -> “Sexual Identity” “Sexual Identity” -> “Prep Awareness” “sexually active” -> “Prep Awareness” Age -> “Prep Awareness” Age -> “Sexual Identity” Birth-sex -> “Sexual Identity” condom -> “Prep Awareness” education -> “Prep Awareness” education -> “Sexual Identity” }
Dag 2 dag { bb=“0,0,1,1” “HIV test” [pos=“0.602,0.655”] “Health Insurance” [pos=“0.165,0.297”] “Heavy drinking” [pos=“0.220,0.770”] “Number of Sex_Partners” [pos=“0.213,0.635”] “Person Healthcare Provider” [exposure,pos=“0.145,0.449”] “Prep Awareness” [outcome,pos=“0.827,0.450”] “Race/ethnicity” [pos=“0.596,0.049”] “Sexual Identity” [pos=“0.710,0.771”] “sexually active” [pos=“0.468,0.883”] Age [pos=“0.399,0.068”] Birthsex [pos=“0.261,0.158”] condom [pos=“0.417,0.704”] education [pos=“0.492,0.214”] “HIV test” -> “Health Insurance” “HIV test” -> “Prep Awareness” “Health Insurance” -> “Prep Awareness” “Heavy drinking” -> “Number of Sex_Partners” “Heavy drinking” -> “Prep Awareness” “Heavy drinking” -> “sexually active” “Number of Sex_Partners” -> “Prep Awareness” “Person Healthcare Provider” -> “Prep Awareness” “Race/ethnicity” -> “Health Insurance” “Race/ethnicity” -> “Person Healthcare Provider” “Race/ethnicity” -> “Prep Awareness” “Sexual Identity” -> “Health Insurance” “Sexual Identity” -> “Person Healthcare Provider” “Sexual Identity” -> “Prep Awareness” “sexually active” -> “Prep Awareness” Age -> “Health Insurance” Age -> “Person Healthcare Provider” Age -> “Prep Awareness” Birthsex -> “Health Insurance” condom -> “Prep Awareness” education -> “Health Insurance” education -> “Person Healthcare Provider” }
# Install packages and load libraries for figure 1 and table1
pacman::p_load(tidyverse, readr, table1, DiagrammeR, rsvg, broom, odds.n.ends)# Loading the dataset
library(haven)
PrEP_Use <- read_sas("chs2020_public.sas7bdat")library(dplyr)
PrEP_Use_Sorted1 <- PrEP_Use %>%
dplyr::select(birthsex, agegroup6, newrace6, education, insuredgateway20,
everheardofprep, sexualid20, pcp20, hiv12months20, condom20,
bingenew, sexpartner, sexuallyactive20)# Selecting the required variables
library(dplyr)
PrEP_Use_Sorted <- PrEP_Use %>%
dplyr::select(birthsex, agegroup6, newrace6, education, insuredgateway20,
everheardofprep, sexualid20, pcp20, hiv12months20, condom20,
bingenew, sexpartner, sexuallyactive20) %>%
dplyr::mutate(
insuredgateway20 = case_when(
insuredgateway20 == 1 ~ 2,
insuredgateway20 == 2 ~ 1,
TRUE ~ insuredgateway20
),
everheardofprep = case_when(
everheardofprep == 1 ~ 2,
everheardofprep == 2 ~ 1,
TRUE ~ everheardofprep
),
pcp20 = case_when(
pcp20 == 1 ~ 2,
pcp20 == 2 ~ 1,
TRUE ~ pcp20))table(PrEP_Use_Sorted$birthsex)##
## 1 2
## 3850 4902table(PrEP_Use_Sorted$everheardofprep)##
## 1 2
## 2853 2360table(PrEP_Use_Sorted$newrace6)##
## 1 2 3 4 5 6
## 2763 1837 2457 1340 96 288head(PrEP_Use_Sorted)library(dplyr)
library(labelled)
PrEP_Use_Sorted <- PrEP_Use_Sorted %>%
mutate (birthsex = recode_factor(birthsex, `1` = "Male", `2` = "Female")) %>%
mutate(agegroup6 = recode_factor(agegroup6, `1` = "18-24", `2` = "25-29", `3` = "30-44", `4` = "45-64", `5` = "65-74", `6` = "75+")) %>%
mutate(newrace6 = recode_factor(newrace6, `1` = "White, non-Hispanic", `2` = "Black, non-Hispanic", `3` = "Hispanic", `4` = "Asian/PI, non-Hispanic", `5` = "N African/Mid Eastern, non-Hispanic", `6` = "Other, non-Hispanic")) %>%
mutate(education = recode_factor(education, `1` = "Less than high school", `2` = "High school graduate", `3` = "Some college/technical school", `4` = "College graduate")) %>%
mutate(insuredgateway20 = recode_factor(insuredgateway20, `1` = "No", `2` = "Yes")) %>%
mutate(everheardofprep = recode_factor(everheardofprep, `1` = "No", `2` = "Yes")) %>%
mutate(sexualid20 = recode_factor(sexualid20, `1` = "Gay/Lesbian", `2` = "Straight, that is not gay", `3` = "Bisexual", `4` = "Something else")) %>%
mutate(pcp20 = recode_factor(pcp20, `1` = "No", `2` = "Yes")) %>%
mutate(hiv12months20 = recode_factor(hiv12months20, `1` = "Yes", `2` = "No")) %>%
mutate(condom20 = recode_factor(condom20, `1` = "Yes", `2` = "No")) %>%
mutate(bingenew = recode_factor(bingenew, `1` = "Yes", `2` = "No")) %>%
mutate(sexpartner = recode_factor(sexpartner, `1` = "None or One", `2` = "None or One", `3` = "Two", `4` = "Three or more")) %>%
mutate(sexuallyactive20 = recode_factor(sexuallyactive20, `1` = "Yes", `2` = "No"))
# Labelled the variables for Table One
var_label(PrEP_Use_Sorted) <- list(
birthsex = "Birth Sex",
agegroup6 = "Age Group",
newrace6 = "Race/Ethnicity",
education = "Education Level",
insuredgateway20 = "Insured",
everheardofprep = "Ever Heard of PrEP",
sexualid20 = "Sexual Identity",
pcp20 = "Personal Healthcare Provider",
hiv12months20 = "HIV Test (12 Months)",
condom20 = "Condom Use",
bingenew = "Binge Drinking",
sexpartner = "Sex Partners",
sexuallyactive20 = "Sexually Active"
)
# View the labeled data frame
print(var_label(PrEP_Use_Sorted))## $birthsex
## [1] "Birth Sex"
##
## $agegroup6
## [1] "Age Group"
##
## $newrace6
## [1] "Race/Ethnicity"
##
## $education
## [1] "Education Level"
##
## $insuredgateway20
## [1] "Insured"
##
## $everheardofprep
## [1] "Ever Heard of PrEP"
##
## $sexualid20
## [1] "Sexual Identity"
##
## $pcp20
## [1] "Personal Healthcare Provider"
##
## $hiv12months20
## [1] "HIV Test (12 Months)"
##
## $condom20
## [1] "Condom Use"
##
## $bingenew
## [1] "Binge Drinking"
##
## $sexpartner
## [1] "Sex Partners"
##
## $sexuallyactive20
## [1] "Sexually Active"library(dplyr)
PrEPUse_Clean <- PrEP_Use_Sorted %>%
drop_na(everheardofprep)
PrEPUse_Clean2 <- PrEPUse_Clean %>%
drop_na(birthsex)
PrEPUse_Clean3 <- PrEPUse_Clean2 %>%
drop_na(agegroup6)
PrEPUse_Clean4 <- PrEPUse_Clean3 %>%
drop_na(condom20)
PrEPUse_Clean5 <- PrEPUse_Clean4%>%
drop_na(sexuallyactive20)
PrEPUse_Clean6 <- PrEPUse_Clean5%>%
drop_na(hiv12months20)
PrEPUse_Clean7 <- PrEPUse_Clean6%>%
drop_na(bingenew)
PrEPUse_Clean8 <- PrEPUse_Clean7%>%
drop_na(sexpartner)
PrEPUse_Clean9 <- PrEPUse_Clean8%>%
drop_na(sexualid20)
PrEPUse_Clean10 <- PrEPUse_Clean9%>%
drop_na(newrace6)
PrEPUse_Clean11 <- PrEPUse_Clean10%>%
drop_na(insuredgateway20)
PrEPUse_Clean12 <- PrEPUse_Clean11%>%
drop_na(pcp20)
PrEPUse_Clean13 <- PrEPUse_Clean12%>%
drop_na(education)head(PrEP_Use_Sorted)library(DiagrammeR)
grViz(diagram = "digraph flowchart {
node [fontname = Helvetica, shape = rectangle, fontsize=15]
node1 [label = '@@1'] # Starting number
node2 [label = '@@2'] # number after exclusion 1
node3 [label = '@@3'] # number after exclusion 2
node4 [label = '@@4'] # number after exclusion 3
node5 [label = '@@5'] # number after exclusion 4
node6 [label = '@@6'] # number after exclusion 5
node7 [label = '@@7'] # number after exclusion 6
node8 [label = '@@8'] # number after exclusion 7
node9 [label = '@@9'] # number after exclusion 8
node10 [label = '@@10'] # number after exclusion 9
node11 [label = '@@11'] # number after exclusion 10
node12 [label = '@@12'] # number after exclusion 11
node13 [label = '@@13'] # number after exclusion 12
node14 [label = '@@13'] # number after exclusion 14
node1 -> node2 -> node3 -> node4 -> node5 -> node6 -> node7 -> node8 -> node9 -> node10 -> node11 -> node12 -> node13 -> node14
}
[1]: 'Starting number of observations in PrEP_Use_Sorted n =8,781'
[2]: 'Excluding 3568 individuals with missing data on Ever heard of PrEP n=5213'
[3]: 'Excluding 0 individuals with missing data on Birthsex n=5213'
[4]: 'Excluding 2 individuals with missing data on Age group n=5211'
[5]: 'Excluding 136 individuals with missing data on Condom Use n=5075'
[6]: 'Excluding 0 individuals with missing data on Sexually active n=5075'
[7]: 'Excluding 64 individuals with missing data on HIV test (12 Months) n=5011'
[8]: 'Excluding 37 individuals with missing data on Binge dringking n=4974'
[9]: 'Excluding 0 individuals with missing data on Sex Partners n=4974'
[10]:'Excluding 109 individuals with missing data on Sexual Identity n=4865'
[11]:'Excluding 0 individuals with missing data on Race/Ethinicty n=4865'
[12]:'Excluding 12 individuals with missing data on Insured n=4853'
[13]:'Excluding 15 individuals with missing data on Personal Helathcare Provider n=4838'
[13]:'Excluding 11 individuals with missing data on Educational Level n=4827'
")# Created the table1 object with renamed variables
table1(~ birthsex + agegroup6 + newrace6 + education + insuredgateway20 + sexualid20 + pcp20 + hiv12months20 + condom20 + bingenew + sexpartner| everheardofprep, data = PrEPUse_Clean13)| No (N=2614) |
Yes (N=2213) |
Overall (N=4827) |
|
|---|---|---|---|
| Birth Sex | |||
| Male | 1370 (52.4%) | 1085 (49.0%) | 2455 (50.9%) |
| Female | 1244 (47.6%) | 1128 (51.0%) | 2372 (49.1%) |
| Age Group | |||
| 18-24 | 198 (7.6%) | 163 (7.4%) | 361 (7.5%) |
| 25-29 | 231 (8.8%) | 284 (12.8%) | 515 (10.7%) |
| 30-44 | 903 (34.5%) | 878 (39.7%) | 1781 (36.9%) |
| 45-64 | 951 (36.4%) | 702 (31.7%) | 1653 (34.2%) |
| 65-74 | 248 (9.5%) | 155 (7.0%) | 403 (8.3%) |
| 75+ | 83 (3.2%) | 31 (1.4%) | 114 (2.4%) |
| Race/Ethnicity | |||
| White, non-Hispanic | 675 (25.8%) | 917 (41.4%) | 1592 (33.0%) |
| Black, non-Hispanic | 427 (16.3%) | 555 (25.1%) | 982 (20.3%) |
| Hispanic | 971 (37.1%) | 477 (21.6%) | 1448 (30.0%) |
| Asian/PI, non-Hispanic | 441 (16.9%) | 143 (6.5%) | 584 (12.1%) |
| N African/Mid Eastern, non-Hispanic | 37 (1.4%) | 24 (1.1%) | 61 (1.3%) |
| Other, non-Hispanic | 63 (2.4%) | 97 (4.4%) | 160 (3.3%) |
| Education Level | |||
| Less than high school | 391 (15.0%) | 92 (4.2%) | 483 (10.0%) |
| High school graduate | 628 (24.0%) | 293 (13.2%) | 921 (19.1%) |
| Some college/technical school | 523 (20.0%) | 460 (20.8%) | 983 (20.4%) |
| College graduate | 1072 (41.0%) | 1368 (61.8%) | 2440 (50.5%) |
| Insured | |||
| No | 398 (15.2%) | 155 (7.0%) | 553 (11.5%) |
| Yes | 2216 (84.8%) | 2058 (93.0%) | 4274 (88.5%) |
| Sexual Identity | |||
| Gay/Lesbian | 14 (0.5%) | 151 (6.8%) | 165 (3.4%) |
| Straight, that is not gay | 2537 (97.1%) | 1900 (85.9%) | 4437 (91.9%) |
| Bisexual | 28 (1.1%) | 129 (5.8%) | 157 (3.3%) |
| Something else | 35 (1.3%) | 33 (1.5%) | 68 (1.4%) |
| Personal Healthcare Provider | |||
| No | 778 (29.8%) | 602 (27.2%) | 1380 (28.6%) |
| Yes | 1836 (70.2%) | 1611 (72.8%) | 3447 (71.4%) |
| HIV Test (12 Months) | |||
| Yes | 824 (31.5%) | 882 (39.9%) | 1706 (35.3%) |
| No | 1790 (68.5%) | 1331 (60.1%) | 3121 (64.7%) |
| Condom Use | |||
| Yes | 741 (28.3%) | 704 (31.8%) | 1445 (29.9%) |
| No | 1873 (71.7%) | 1509 (68.2%) | 3382 (70.1%) |
| Binge Drinking | |||
| Yes | 439 (16.8%) | 566 (25.6%) | 1005 (20.8%) |
| No | 2175 (83.2%) | 1647 (74.4%) | 3822 (79.2%) |
| Sex Partners | |||
| None or One | 2262 (86.5%) | 1661 (75.1%) | 3923 (81.3%) |
| Two | 158 (6.0%) | 193 (8.7%) | 351 (7.3%) |
| Three or more | 194 (7.4%) | 359 (16.2%) | 553 (11.5%) |
# Model 1: Adjusting for Age, Personal Healthcare Provider, Sexual-Identity, Race/ethnicity and Education
model1 <- glm(everheardofprep ~ pcp20 + insuredgateway20 + agegroup6 +birthsex+ newrace6 + sexualid20 + education,
data = PrEPUse_Clean13,
family = binomial)
model1_summary <- tidy(model1, conf.int = TRUE)
odds.n.ends(model1)## Waiting for profiling to be done...## $`Logistic regression model significance`
## Chi-squared d.f. p
## 936.169 19 <.001
##
## $`Contingency tables (model fit): frequency predicted`
## Number observed
## Number predicted 1 0 Sum
## 1 1369 743 2112
## 0 844 1871 2715
## Sum 2213 2614 4827
##
## $`Count R-squared (model fit): percent correctly predicted`
## [1] 67.12244
##
## $`Model sensitivity`
## [1] 0.6186173
##
## $`Model specificity`
## [1] 0.7157613
##
## $`Predictor odds ratios and 95% CI`
## OR 2.5 % 97.5 %
## (Intercept) 3.52521633 1.80903047 7.2398058
## pcp20Yes 1.08823832 0.93692766 1.2642124
## insuredgateway20Yes 1.82351973 1.45742860 2.2895235
## agegroup625-29 1.31562256 0.97543889 1.7768329
## agegroup630-44 1.04665645 0.81239001 1.3503223
## agegroup645-64 0.68111167 0.52562299 0.8833337
## agegroup665-74 0.47187691 0.34106699 0.6521640
## agegroup675+ 0.29033947 0.17420252 0.4749509
## birthsexFemale 1.09653895 0.96451439 1.2465712
## newrace6Black, non-Hispanic 1.33311273 1.11467015 1.5954743
## newrace6Hispanic 0.49189689 0.41280940 0.5857078
## newrace6Asian/PI, non-Hispanic 0.23428011 0.18516555 0.2949778
## newrace6N African/Mid Eastern, non-Hispanic 0.49739008 0.28484183 0.8527777
## newrace6Other, non-Hispanic 1.21477730 0.85720220 1.7309502
## sexualid20Straight, that is not gay 0.06965059 0.03754786 0.1194561
## sexualid20Bisexual 0.39976114 0.19054283 0.8071427
## sexualid20Something else 0.11636582 0.05199326 0.2495416
## educationHigh school graduate 1.64542508 1.23896514 2.2001777
## educationSome college/technical school 2.50837000 1.89678018 3.3415958
## educationCollege graduate 3.62950434 2.77736406 4.7827035# Model 2: Adjusting for Personal Healthcare Provider with an interaction Health insurance, Age, New race Sexual-Identity, Race/ethnicity and Education
model2 <- glm(everheardofprep ~ pcp20 * insuredgateway20 + agegroup6 + newrace6 + education +
sexualid20,
data = PrEPUse_Clean13,
family = binomial)
model2_summary <- tidy(model2, conf.int = TRUE)
odds.n.ends(model2)## Waiting for profiling to be done...## $`Logistic regression model significance`
## Chi-squared d.f. p
## 935.408 19 <.001
##
## $`Contingency tables (model fit): frequency predicted`
## Number observed
## Number predicted 1 0 Sum
## 1 1368 746 2114
## 0 845 1868 2713
## Sum 2213 2614 4827
##
## $`Count R-squared (model fit): percent correctly predicted`
## [1] 67.03957
##
## $`Model sensitivity`
## [1] 0.6181654
##
## $`Model specificity`
## [1] 0.7146136
##
## $`Predictor odds ratios and 95% CI`
## OR 2.5 % 97.5 %
## (Intercept) 3.15622964 1.58664775 6.5937526
## pcp20Yes 1.37325537 0.90128916 2.0876564
## insuredgateway20Yes 2.04637391 1.52475253 2.7643871
## agegroup625-29 1.32485497 0.98208567 1.7896854
## agegroup630-44 1.04997254 0.81491163 1.3546832
## agegroup645-64 0.68131639 0.52574833 0.8836446
## agegroup665-74 0.46961379 0.33956165 0.6487599
## agegroup675+ 0.28534511 0.17144585 0.4660420
## newrace6Black, non-Hispanic 1.34885858 1.12867924 1.6131493
## newrace6Hispanic 0.49940415 0.41934675 0.5943259
## newrace6Asian/PI, non-Hispanic 0.23482816 0.18558313 0.2956957
## newrace6N African/Mid Eastern, non-Hispanic 0.49351503 0.28280524 0.8455925
## newrace6Other, non-Hispanic 1.22458313 0.86390614 1.7454136
## educationHigh school graduate 1.64748463 1.24049308 2.2029508
## educationSome college/technical school 2.51078474 1.89866789 3.3446853
## educationCollege graduate 3.65619592 2.79866953 4.8164391
## sexualid20Straight, that is not gay 0.07288069 0.03938878 0.1245634
## sexualid20Bisexual 0.42458254 0.20309920 0.8538267
## sexualid20Something else 0.12395370 0.05550567 0.2651767
## pcp20Yes:insuredgateway20Yes 0.77697706 0.49744322 1.2162974model3 <- glm(everheardofprep ~ pcp20,
data = PrEPUse_Clean13 %>%filter (insuredgateway20=="Yes"),
family = binomial)
model3_summary <- tidy(model3, conf.int = TRUE)
odds.n.ends(model3)## Waiting for profiling to be done...## $`Logistic regression model significance`
## Chi-squared d.f. p
## 2.353 1.000 0.125
##
## $`Contingency tables (model fit): frequency predicted`
## Number observed
## Number predicted 1 0 Sum
## 1 513 508 1021
## 0 1545 1708 3253
## Sum 2058 2216 4274
##
## $`Count R-squared (model fit): percent correctly predicted`
## [1] 51.96537
##
## $`Model sensitivity`
## [1] 0.2492711
##
## $`Model specificity`
## [1] 0.7707581
##
## $`Predictor odds ratios and 95% CI`
## OR 2.5 % 97.5 %
## (Intercept) 1.0098425 0.8932302 1.141707
## pcp20Yes 0.8957503 0.7781778 1.031048model4 <- glm(everheardofprep ~ pcp20 +birthsex+ newrace6+ sexualid20 + education,
data = PrEPUse_Clean13 %>%filter (insuredgateway20=="No"),
family = binomial)
model4_summary <- tidy(model4, conf.int = TRUE)
odds.n.ends(model4)## Waiting for profiling to be done...## $`Logistic regression model significance`
## Chi-squared d.f. p
## 103.668 13 <.001
##
## $`Contingency tables (model fit): frequency predicted`
## Number observed
## Number predicted 1 0 Sum
## 1 57 40 97
## 0 98 358 456
## Sum 155 398 553
##
## $`Count R-squared (model fit): percent correctly predicted`
## [1] 75.04521
##
## $`Model sensitivity`
## [1] 0.3677419
##
## $`Model specificity`
## [1] 0.8994975
##
## $`Predictor odds ratios and 95% CI`
## OR 2.5 % 97.5 %
## (Intercept) 1.6318804 0.37402344 8.9558093
## pcp20Yes 1.1993653 0.77449879 1.8476233
## birthsexFemale 0.9236954 0.59894562 1.4182007
## newrace6Black, non-Hispanic 1.6212776 0.83627433 3.1684118
## newrace6Hispanic 0.4232120 0.23135335 0.7742368
## newrace6Asian/PI, non-Hispanic 0.2319726 0.08640276 0.5662972
## newrace6N African/Mid Eastern, non-Hispanic 1.1279677 0.25041356 4.8348262
## newrace6Other, non-Hispanic 1.0179726 0.25423861 3.9023257
## sexualid20Straight, that is not gay 0.1351886 0.02883756 0.4764161
## sexualid20Bisexual 0.7260345 0.11188654 4.0855752
## sexualid20Something else 0.1032829 0.01309385 0.6039255
## educationHigh school graduate 1.8030302 0.91115910 3.6962958
## educationSome college/technical school 3.2764849 1.59116643 6.9710762
## educationCollege graduate 4.2870998 2.18761012 8.7800654# Comparing models using AIC
AIC(model1, model2)# Likelihood Ratio Test
anova(model1, model2, test = "Chisq")# Load necessary libraries
library(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# Calculate Variance Inflation Factor (VIF)
vif(model2)## there are higher-order terms (interactions) in this model
## consider setting type = 'predictor'; see ?vif## GVIF Df GVIF^(1/(2*Df))
## pcp20 9.038462 1 3.006404
## insuredgateway20 1.930540 1 1.389439
## agegroup6 1.233489 5 1.021206
## newrace6 1.334791 5 1.029298
## education 1.250699 3 1.037987
## sexualid20 1.034260 3 1.005630
## pcp20:insuredgateway20 11.011191 1 3.318311# Load the necessary library
library(ResourceSelection)## Warning: package 'ResourceSelection' was built under R version 4.4.3## ResourceSelection 0.3-6 2023-06-27# Perform the Hosmer-Lemeshow test for both models
hl_test_model1 <- hoslem.test(model1$y, fitted(model1), g = 10)
hl_test_model2 <- hoslem.test(model2$y, fitted(model2), g = 10)
# Print the results of both tests
print("Hosmer-Lemeshow Test for Model 1:")## [1] "Hosmer-Lemeshow Test for Model 1:"print(hl_test_model1)##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: model1$y, fitted(model1)
## X-squared = 8.1073, df = 8, p-value = 0.4231print("Hosmer-Lemeshow Test for Model 2:")## [1] "Hosmer-Lemeshow Test for Model 2:"print(hl_test_model2)##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: model2$y, fitted(model2)
## X-squared = 7.7417, df = 8, p-value = 0.4591