Research question: Do growing up with a problem drinker or alcoholic or growing up with anyone who was depressed, mentally ill, or suicidal predict the mental health days by gender?

library(car)
## Loading required package: carData
library(stargazer)
## 
## Please cite as:
##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.2. https://CRAN.R-project.org/package=stargazer
library(haven)
library(survey)
## Loading required package: grid
## Loading required package: Matrix
## Loading required package: survival
## 
## Attaching package: 'survey'
## The following object is masked from 'package:graphics':
## 
##     dotchart
library(questionr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.3     v purrr   0.3.4
## v tibble  3.0.5     v stringr 1.4.0
## v tidyr   1.1.2     v forcats 0.5.0
## v readr   1.4.0
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x tidyr::expand() masks Matrix::expand()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## x tidyr::pack()   masks Matrix::pack()
## x dplyr::recode() masks car::recode()
## x purrr::some()   masks car::some()
## x tidyr::unpack() masks Matrix::unpack()
BRFSS19<-read_xpt("C:/Users/adolp/Desktop/Statistics 2/BRFSS/LLCP2019.XPT")
nams<-names(BRFSS19)
head(BRFSS19, n=10)
## # A tibble: 10 x 342
##    `_STATE` FMONTH IDATE IMONTH IDAY  IYEAR DISPCODE SEQNO `_PSU` CTELENM1
##       <dbl>  <dbl> <chr> <chr>  <chr> <chr>    <dbl> <chr>  <dbl>    <dbl>
##  1        1      1 0118~ 01     18    2019      1100 2019~ 2.02e9        1
##  2        1      1 0113~ 01     13    2019      1100 2019~ 2.02e9        1
##  3        1      1 0118~ 01     18    2019      1100 2019~ 2.02e9        1
##  4        1      1 0118~ 01     18    2019      1200 2019~ 2.02e9        1
##  5        1      1 0104~ 01     04    2019      1100 2019~ 2.02e9        1
##  6        1      1 0118~ 01     18    2019      1200 2019~ 2.02e9        1
##  7        1      1 0104~ 01     04    2019      1100 2019~ 2.02e9        1
##  8        1      1 0123~ 01     23    2019      1100 2019~ 2.02e9        1
##  9        1      1 0124~ 01     24    2019      1100 2019~ 2.02e9        1
## 10        1      1 0113~ 01     13    2019      1100 2019~ 2.02e9        1
## # ... with 332 more variables: PVTRESD1 <dbl>, COLGHOUS <dbl>, STATERE1 <dbl>,
## #   CELPHONE <dbl>, LADULT1 <dbl>, COLGSEX <dbl>, NUMADULT <dbl>,
## #   LANDSEX <dbl>, NUMMEN <dbl>, NUMWOMEN <dbl>, RESPSLCT <dbl>,
## #   SAFETIME <dbl>, CTELNUM1 <dbl>, CELLFON5 <dbl>, CADULT1 <dbl>,
## #   CELLSEX <dbl>, PVTRESD3 <dbl>, CCLGHOUS <dbl>, CSTATE1 <dbl>,
## #   LANDLINE <dbl>, HHADULT <dbl>, SEXVAR <dbl>, GENHLTH <dbl>, PHYSHLTH <dbl>,
## #   MENTHLTH <dbl>, POORHLTH <dbl>, HLTHPLN1 <dbl>, PERSDOC2 <dbl>,
## #   MEDCOST <dbl>, CHECKUP1 <dbl>, BPHIGH4 <dbl>, BPMEDS <dbl>, CHOLCHK2 <dbl>,
## #   TOLDHI2 <dbl>, CHOLMED2 <dbl>, CVDINFR4 <dbl>, CVDCRHD4 <dbl>,
## #   CVDSTRK3 <dbl>, ASTHMA3 <dbl>, ASTHNOW <dbl>, CHCSCNCR <dbl>,
## #   CHCOCNCR <dbl>, CHCCOPD2 <dbl>, ADDEPEV3 <dbl>, CHCKDNY2 <dbl>,
## #   DIABETE4 <dbl>, DIABAGE3 <dbl>, HAVARTH4 <dbl>, ARTHEXER <dbl>,
## #   ARTHEDU <dbl>, LMTJOIN3 <dbl>, ARTHDIS2 <dbl>, JOINPAI2 <dbl>,
## #   MARITAL <dbl>, EDUCA <dbl>, RENTHOM1 <dbl>, NUMHHOL3 <dbl>, NUMPHON3 <dbl>,
## #   CPDEMO1B <dbl>, VETERAN3 <dbl>, EMPLOY1 <dbl>, CHILDREN <dbl>,
## #   INCOME2 <dbl>, WEIGHT2 <dbl>, HEIGHT3 <dbl>, PREGNANT <dbl>, DEAF <dbl>,
## #   BLIND <dbl>, DECIDE <dbl>, DIFFWALK <dbl>, DIFFDRES <dbl>, DIFFALON <dbl>,
## #   SMOKE100 <dbl>, SMOKDAY2 <dbl>, STOPSMK2 <dbl>, LASTSMK2 <dbl>,
## #   USENOW3 <dbl>, ALCDAY5 <dbl>, AVEDRNK3 <dbl>, DRNK3GE5 <dbl>,
## #   MAXDRNKS <dbl>, EXERANY2 <dbl>, EXRACT11 <dbl>, EXEROFT1 <dbl>,
## #   EXERHMM1 <dbl>, EXRACT21 <dbl>, EXEROFT2 <dbl>, EXERHMM2 <dbl>,
## #   STRENGTH <dbl>, FRUIT2 <dbl>, FRUITJU2 <dbl>, FVGREEN1 <dbl>,
## #   FRENCHF1 <dbl>, POTATOE1 <dbl>, VEGETAB2 <dbl>, FLUSHOT7 <dbl>,
## #   FLSHTMY3 <dbl>, TETANUS1 <dbl>, PNEUVAC4 <dbl>, HIVTST7 <dbl>, ...
newnames<-tolower(gsub(pattern = "_",replacement =  "",x =  nams))
names(BRFSS19)<-newnames
BRFSS19$Alc<-Recode(BRFSS19$acedrink, recodes ="7:9=NA; 1=1;2=0")
BRFSS19$Dep<-Recode(BRFSS19$acedeprs, recodes ="7:9=NA; 1=1;2=0") 
BRFSS19$Gen<-Recode(BRFSS19$birthsex, recodes = "7:9=NA; 1=1;2=0")
BRFSS19$Mhlth<-Recode(BRFSS19$ment14d, recodes ="7:9=NA; 1=0;2=1;3=1") 
sub<-BRFSS19 %>%
  select(Alc,Dep,Gen,Mhlth,llcpwt, ststr) %>%
  filter( complete.cases( . ))
options(survey.lonely.psu = "adjust")
des<-svydesign(ids=~1, strata=~ststr, weights=~llcpwt, data = sub)
cat<-svyby(formula = ~Mhlth,
           by = ~Alc,
           design = des,
           FUN = svymean,
           na.rm=T)
svychisq(~Mhlth+Alc,
         design = des)
## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(~Mhlth + Alc, design = des)
## F = 112.41, ndf = 1, ddf = 5336, p-value < 2.2e-16
cat%>%
  ggplot()+
  geom_point(aes(x=Alc,y=Mhlth))+
  geom_errorbar(aes(x=Alc, ymin = Mhlth-1.96*se, 
                    ymax= Mhlth+1.96*se),
                width=.25)+
  labs(title = "Percent % of US Adults' Mental Health Status by Growing Up With A Problem Drinker Or Alcoholic", 
       caption = "Source: CDC BRFSS - Annual Data, 2019 \n Calculations by Adolph J. Delgado, M.Ed.",
       x = "Growing Up With A Problem Drinker Or Alcoholic",
       y = "%  Mental Health Status")+
  theme_minimal()

cat2<-svyby(formula = ~Mhlth,
           by = ~Dep,
           design = des,
           FUN = svymean,
           na.rm=T)
svychisq(~Mhlth+Dep,
         design = des)
## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(~Mhlth + Dep, design = des)
## F = 193.07, ndf = 1, ddf = 5336, p-value < 2.2e-16
cat2%>%
  ggplot()+
  geom_point(aes(x=Dep,y=Mhlth))+
  geom_errorbar(aes(x=Dep, ymin = Mhlth-1.96*se, 
                    ymax= Mhlth+1.96*se),
                width=.25)+
  labs(title = "Percent % of US Adults' Mental Health Status by Growing Up With Anyone Who Was Depressed, Mentally Ill, Or Suicidal", 
       caption = "Source: CDC BRFSS - Annual Data, 2019 \n Calculations by Adolph J. Delgado, M.Ed.",
       x = "Growing Up With anyone who was depressed, mentally ill, or suicidal",
       y = "%  Mental Health Status")+
  theme_minimal()

cat3<-svyby(formula = ~Mhlth,
           by = ~Gen,
           design = des,
           FUN = svymean,
           na.rm=T)
svychisq(~Mhlth+Gen,
         design = des)
## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(~Mhlth + Gen, design = des)
## F = 44.794, ndf = 1, ddf = 5336, p-value = 2.414e-11
cat3%>%
  ggplot()+
  geom_point(aes(x=Gen,y=Mhlth))+
  geom_errorbar(aes(x=Gen, ymin = Mhlth-1.96*se, 
                    ymax= Mhlth+1.96*se),
                width=.25)+
  labs(title = "Percent % of US Adults' Mental Health Status by Gender", 
       caption = "Source: CDC BRFSS - Annual Data, 2019 \n Calculations by Adolph J. Delgado, M.Ed.",
       x = "Gender",
       y = "%  Mental Health Status")+
  theme_minimal()

fit.logit<-svyglm(Mhlth ~ Alc + Dep + Gen,
                  design = des,
                  family = binomial)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
library(broom)
fit.logit%>%
  tidy()%>%
  knitr::kable(digits = 3)
term estimate std.error statistic p.value
(Intercept) -0.631 0.059 -10.616 0
Alc 0.555 0.089 6.209 0
Dep 1.034 0.097 10.631 0
Gen -0.424 0.077 -5.528 0 
fit.logit%>%
  tidy()%>%
  mutate(OR = exp(estimate))%>%
  knitr::kable(digits = 3)
term estimate std.error statistic p.value OR
(Intercept) -0.631 0.059 -10.616 0 0.532
Alc 0.555 0.089 6.209 0 1.741
Dep 1.034 0.097 10.631 0 2.812
Gen -0.424 0.077 -5.528 0 0.654 
fit.logit%>%
  tidy()%>%
  mutate(OR = exp(estimate),
         LowerOR_Ci = exp(estimate - 1.96*std.error),
         UpperOR_Ci = exp(estimate + 1.96*std.error))%>%
  knitr::kable(digits = 3)
term estimate std.error statistic p.value OR LowerOR_Ci UpperOR_Ci
(Intercept) -0.631 0.059 -10.616 0 0.532 0.474 0.598
Alc 0.555 0.089 6.209 0 1.741 1.462 2.075
Dep 1.034 0.097 10.631 0 2.812 2.324 3.402
Gen -0.424 0.077 -5.528 0 0.654 0.563 0.760 
library(sjPlot)
## Registered S3 methods overwritten by 'lme4':
##   method                          from
##   cooks.distance.influence.merMod car 
##   influence.merMod                car 
##   dfbeta.influence.merMod         car 
##   dfbetas.influence.merMod        car
## Learn more about sjPlot with 'browseVignettes("sjPlot")'.
plot_model(fit.logit,
           title = "Odds ratios for Mental Health")

fit.probit<-svyglm(Mhlth~Alc+Dep+Gen,
                   design=des,
                   family=binomial(link= "probit"))
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
myexp<-function(x){exp(x)}
stargazer(fit.logit, fit.probit,
          type = "html",
          style="demography",
          covariate.labels=c("Alc", "Dep", "Gen", "Mhlth"),
          ci=T,
          column.sep.width = "10pt")
Mhlth
survey-weighted survey-weighted
logistic probit
Model 1 Model 2
Alc 0.555*** 0.341***
(0.380, 0.730) (0.233, 0.449)
Dep 1.034*** 0.641***
(0.843, 1.224) (0.524, 0.759)
Gen -0.424*** -0.258***
(-0.575, -0.274) (-0.349, -0.166)
Mhlth -0.631*** -0.392***
(-0.747, -0.514) (-0.463, -0.320)
N 5,360 5,360
Log Likelihood -3,322.921 -3,323.067
AIC 6,653.841 6,654.133
p < .05; p < .01; p < .001 
library(emmeans)
rg<-ref_grid(fit.logit)

marg_logit<-emmeans(object = rg,
                    specs = c( "Alc", "Dep", "Gen"),
                    type="response",data= sub )

knitr::kable(marg_logit,  digits = 4)
Alc Dep Gen prob SE df asymp.LCL asymp.UCL
0 0 0 0.3473 0.0135 Inf 0.3214 0.3742
1 0 0 0.4810 0.0223 Inf 0.4375 0.5248
0 1 0 0.5994 0.0239 Inf 0.5517 0.6453
1 1 0 0.7227 0.0199 Inf 0.6820 0.7600
0 0 1 0.2583 0.0118 Inf 0.2359 0.2820
1 0 1 0.3775 0.0218 Inf 0.3357 0.4211
0 1 1 0.4947 0.0261 Inf 0.4438 0.5456
1 1 1 0.6303 0.0246 Inf 0.5808 0.6772 
library(emmeans)
rg<-ref_grid(fit.probit)

marg_logit<-emmeans(object = rg,
                    specs = c( "Alc", "Dep", "Gen"),
                    type="response",data= sub )

knitr::kable(marg_logit,  digits = 4)
Alc Dep Gen prob SE df asymp.LCL asymp.UCL
0 0 0 0.3477 0.0135 Inf 0.3217 0.3745
1 0 0 0.4799 0.0221 Inf 0.4368 0.5232
0 1 0 0.5986 0.0238 Inf 0.5514 0.6445
1 1 0 0.7227 0.0203 Inf 0.6817 0.7610
0 0 1 0.2582 0.0119 Inf 0.2354 0.2820
1 0 1 0.3791 0.0218 Inf 0.3372 0.4224
0 1 1 0.4969 0.0257 Inf 0.4467 0.5472
1 1 1 0.6306 0.0245 Inf 0.5817 0.6774 
comps<-as.data.frame(marg_logit)
comps<-comps[is.na(comps$prob)==F,]

comps
##   Alc Dep Gen      prob         SE  df asymp.LCL asymp.UCL
## 1   0   0   0 0.3477102 0.01348707 Inf 0.3216649 0.3744951
## 2   1   0   0 0.4799139 0.02209037 Inf 0.4368205 0.5232435
## 3   0   1   0 0.5986093 0.02379980 Inf 0.5513683 0.6444500
## 4   1   1   0 0.7227027 0.02025984 Inf 0.6816684 0.7609646
## 5   0   0   1 0.2581560 0.01187362 Inf 0.2354390 0.2819594
## 6   1   0   1 0.3790789 0.02179118 Inf 0.3371859 0.4224436
## 7   0   1   1 0.4968955 0.02570932 Inf 0.4466645 0.5471758
## 8   1   1   1 0.6305680 0.02448800 Inf 0.5816737 0.6774354

Research question: Do growing up with a problem drinker or alcoholic or growing up with anyone who was depressed, mentally ill, or suicidal predict mental health days as an adult by gender? For females, growing up with a problem drinker or alcoholic or growing up with anyone who was depressed increases the probability mental health days as an adult from 35% to 72%. For males, growing up with a problem drinker or alcoholic or growing up with anyone who was depressed increases the probability mental health days as an adult from 26% to 63%. Thus, based on the self-reported responses on the 2019 BRFSS, females’ mental health days as an adult are equally impacted (37% increase) by growing up with a problem drinker or alcoholic or growing up with anyone who was depressed as males.