Logit Model, HW 3

Research question: Does alcohol consumption vary with educational level and racial/ethnic background?

Predictor Variables: Educational level and racial/ethnic background

Loading libraries

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(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(haven)
library(ggplot2)

Loading data from BRFSS 2020

brfss20<- readRDS(url("https://github.com/coreysparks/DEM7283/blob/master/data/brfss20sm.rds?raw=true"))
### Fix variable names
names(brfss20) <- tolower(gsub(pattern = "_",replacement =  "",x =  names(brfss20)))

Recoding predictor and Outcome variables

#sex
brfss20$male<-as.factor(ifelse(brfss20$sex==1,
                                "Male",
                                "Female"))

#Age cut into intervals
brfss20$agec<-cut(brfss20$age80,
                   breaks=c(0,24,39,59,79,99))

brfss20$educ<-Recode(brfss20$educa,
                      recodes="1:2='0Prim'; 3='1somehs'; 4='2hsgrad'; 5='3somecol'; 6='4colgrad';9=NA",
                      as.factor=T)

#Racial/Ethnic Background 

brfss20$black<-Recode(brfss20$racegr3,
                       recodes="2=1; 9=NA; else=0")
brfss20$white<-Recode(brfss20$racegr3,
                       recodes="1=1; 9=NA; else=0")
brfss20$other<-Recode(brfss20$racegr3,
                       recodes="3:4=1; 9=NA; else=0")
brfss20$hispanic<-Recode(brfss20$racegr3,
                          recodes="5=1; 9=NA; else=0")

brfss20$race_eth<-Recode(brfss20$racegr3,
                          recodes="1='nhwhite'; 2='nh_black'; 3='nh_other';4='nh_multirace'; 5='hispanic'; else=NA",
                          as.factor = T)

#insurance
brfss20$ins<-Recode(brfss20$hlthpln1,
                     recodes ="7:9=NA; 1=1;2=0")

#income grouping
brfss20$inc<-ifelse(brfss20$incomg==9,
                     NA,
                     brfss20$incomg)
#employment
brfss20$employ<-Recode(brfss20$employ1,
                        recodes="1:2='Employed'; 2:6='nilf'; 7='retired'; 8='unable'; else=NA",
                        as.factor=T)

brfss20$employ<-relevel(brfss20$employ,
                         ref='Employed')

#marital status
brfss20$marst<-Recode(brfss20$marital,
                       recodes="1='married'; 2='divorced'; 3='widowed'; 4='separated'; 5='nm';6='cohab'; else=NA",
                       as.factor=T)

brfss20$marst<-relevel(brfss20$marst,
                        ref='married')
#BMI, in the brfss20a the bmi variable has 2 implied decimal places, so we must divide by 100 to get real bmi's

brfss20$bmi<-brfss20$bmi5/100

brfss20$obese<-ifelse(brfss20$bmi>=30,
                       1,
                       0)

Recoding outcome variable

I have recoded alcohol consumption variable in past 30 days (drnkany5) into binary variable. I put “1” for for those who consumed alcohol in the past 30 days and “0” for those who do not consumed alcohol in the past days. Hence, it has been recoded as-1: Yes and 0: No. For NA, I have recoded as-7 and 9:NA.

#Alcohol Intake within past 30 days
brfss20$alcohol<-Recode(brfss20$drnkany5, recodes="1=1; 2=0; else = NA")
table(brfss20$alcohol)

## 
##     0     1 
## 82998 98015

Analysis

Descriptive Analyses

library(dplyr)
df<-brfss20%>%
  
  select(drnkany5,mmsaname, bmi, agec,race_eth, marst, educ,white, black, hispanic, other, alcohol, ins, mmsawt, ststr) %>%
  filter( complete.cases(.))

#First we tell R our survey design
options(survey.lonely.psu = "adjust")
des<-svydesign(ids=~1, strata=~ststr, weights=~mmsawt, data =df )

First, we examine the % of US adults with alcohol consumption by education level, and do a survey-corrected chi-square test for independence.

Data<-svyby(formula = ~alcohol,
           by = ~educ,
           design = des,
           FUN = svymean,
           na.rm=T)

svychisq(~alcohol+educ,
         design = des)

## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(~alcohol + educ, design = des)
## F = 229.28, ndf = 3.4414e+00, ddf = 5.5440e+05, p-value < 2.2e-16

Plot of estimates with standard errors

Data%>%
  ggplot()+
  geom_point(aes(x=educ,y=alcohol))+
  geom_errorbar(aes(x=educ, ymin = alcohol-1.96*se, 
                    ymax= alcohol+1.96*se),
                width=.25)+
   labs(title = "Percent % of US Adults with Alcohol Consumption", 
        caption = "Source: CDC BRFSS - SMART Data, 2020 \n Calculations by Jyoti Nepal, MSW",
       x = "Education",
       y = "%  Alcohol Consumption in Past 30 Days")+
  theme_minimal()

The plot shows that with increasing education, the percentage of alcohol consumption increases.

Calculate race*alcohol cross tabulation, and plot it

Data<-svyby(formula = ~alcohol,
           by = ~race_eth, 
           design = des, 
           FUN = svymean,
           na.rm=T)

svychisq(~alcohol+race_eth,
         design = des)

## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(~alcohol + race_eth, design = des)
## F = 101.28, ndf = 3.4718e+00, ddf = 5.5930e+05, p-value < 2.2e-16

Data%>%
  ggplot()+
  geom_point(aes(x=race_eth,y=alcohol))+
  geom_errorbar(aes(x=race_eth, ymin = alcohol-1.96*se, 
                    ymax= alcohol+1.96*se),
                width=.25)+
   labs(title = "Percent % of US Adults with Alcohol Comsumption by Race/Ethnicity", 
        caption = "Source: CDC BRFSS - SMART Data, 2020 \n Calculations by Jyoti Nepal, MSW",
       x = "Race/Ethnicity",
       y = "%  Alcohol Consumption")+
  theme_minimal()

With respect to racial and ethnic population and their chance of drinking alcohol, the non-Hispanic Whites have the higher percent of alcohol consumption with the non-Hispanic multirace having second highest consumption. The non-Hispanic other have the least percent of drinking alcohol compared to other racial and ethnic groups.

Calculate race by education by alcohol consumtion cross tabulation, and plot it

Data<-svyby(formula = ~alcohol,
              by = ~race_eth+educ,
              design = des,
              FUN = svymean,
              na.rm=T) 

#this plot is a little more complicated, but facet_wrap() plots separate plots for groups

Data%>%
  ggplot()+
  #geom_point(aes(x=educ, y = alcohol, color=race_eth, group=race_eth), position="dodge")+ 
  geom_errorbar(aes(x=educ,y = alcohol,
                    ymin = alcohol-1.96*se, 
                   ymax= alcohol+1.96*se,
                   color=race_eth,
                   group=race_eth),
                width=.25,
                position="dodge")+
  #facet_wrap(~ race_eth, nrow = 3)+
  labs(title = "Percent % of US Adults with Alcohol Consumption by Race/Ethnicity and Education", 
        caption = "Source: CDC BRFSS - SMART Data, 2017 \n Calculations by Jyoti Nepal, MSW",
       x = "Education",
       y = "Alcohol")+
  theme_minimal()

Generally, higher education degree correlates to higher alcohol consumption, within the population with same education, nh white appears to have higher chance of consuming alcohol compared to other racial groups. Similarly, for a given education category, nh others had lowest percentage of alcohol consumption. It should be noted however that the overlapping error bars means that the data have high degree of variability.

Logit Model

#Logit model
fit.logit<-svyglm(alcohol~race_eth+educ,
                  design= des,
                  family=binomial)

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

summary(fit.logit)

## 
## Call:
## svyglm(formula = alcohol ~ race_eth + educ, design = des, family = binomial)
## 
## Survey design:
## svydesign(ids = ~1, strata = ~ststr, weights = ~mmsawt, data = df)
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -0.68038    0.09647  -7.053 1.76e-12 ***
## race_ethnh_black     -0.10046    0.04895  -2.052   0.0402 *  
## race_ethnh_multirace  0.08192    0.08145   1.006   0.3145    
## race_ethnh_other     -0.52640    0.06402  -8.222  < 2e-16 ***
## race_ethnhwhite       0.23944    0.03890   6.155 7.52e-10 ***
## educ1somehs           0.12067    0.11008   1.096   0.2730    
## educ2hsgrad           0.46781    0.09852   4.748 2.05e-06 ***
## educ3somecol          0.79344    0.09851   8.054 8.05e-16 ***
## educ4colgrad          1.30787    0.09720  13.455  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 0.9996027)
## 
## Number of Fisher Scoring iterations: 4

Here, we only see coefficients, the direction of relationship, and the significance. Now, we need to calculate odds ratios and confidence intervals.

Get odd ratios and confidence intervals for the estimates

library(gtsummary)
fit.logit%>%
  tbl_regression(exponentiate=TRUE )

Characteristic	OR1	95% CI1	p-value
race_eth
hispanic	—	—
nh_black	0.90	0.82, 1.00	0.040
nh_multirace	1.09	0.93, 1.27	0.3
nh_other	0.59	0.52, 0.67	<0.001
nhwhite	1.27	1.18, 1.37	<0.001
educ
0Prim	—	—
1somehs	1.13	0.91, 1.40	0.3
2hsgrad	1.60	1.32, 1.94	<0.001
3somecol	2.21	1.82, 2.68	<0.001
4colgrad	3.70	3.06, 4.47	<0.001
1 OR = Odds Ratio, CI = Confidence Interval

library(sjPlot)

## #refugeeswelcome

plot_model(fit.logit,
           axis.lim = c(.1, 10),
           title = "Odds Ratios for Alcohol Consumption")

Non-Hispanic Whites drink alcohol the most out of all of the races included in the study and interestingly, those who are college graduates drink the most.

knitr::kable(data.frame(OR = exp(coef(fit.logit)), ci = exp(confint(fit.logit)))) 

	OR	ci.2.5..	ci.97.5..
(Intercept)	0.5064253	0.4191804	0.6118285
race_ethnh_black	0.9044207	0.8216750	0.9954991
race_ethnh_multirace	1.0853721	0.9252210	1.2732445
race_ethnh_other	0.5907297	0.5210683	0.6697042
race_ethnhwhite	1.2705358	1.1772647	1.3711965
educ1somehs	1.1282508	0.9092992	1.3999242
educ2hsgrad	1.5964985	1.3161517	1.9365606
educ3somecol	2.2109971	1.8227750	2.6819044
educ4colgrad	3.6982840	3.0567651	4.4744375

Fitted values

#get a series of predicted probabilities for different "types" of people for each model
#ref_grid will generate all possible combinations of predictors from a model

library(emmeans)
rg<-ref_grid(fit.logit)

marg_logit<-emmeans(object = rg,
              specs = c( "race_eth",  "educ"),
              type="response" )

knitr::kable(marg_logit,  digits = 4)

race_eth	educ	prob	SE	df	asymp.LCL	asymp.UCL
hispanic	0Prim	0.3362	0.0215	Inf	0.2954	0.3796
nh_black	0Prim	0.3141	0.0217	Inf	0.2733	0.3581
nh_multirace	0Prim	0.3547	0.0274	Inf	0.3029	0.4101
nh_other	0Prim	0.2303	0.0193	Inf	0.1947	0.2702
nhwhite	0Prim	0.3915	0.0230	Inf	0.3475	0.4374
hispanic	1somehs	0.3636	0.0145	Inf	0.3357	0.3925
nh_black	1somehs	0.3407	0.0138	Inf	0.3142	0.3682
nh_multirace	1somehs	0.3828	0.0211	Inf	0.3424	0.4248
nh_other	1somehs	0.2524	0.0142	Inf	0.2255	0.2812
nhwhite	1somehs	0.4206	0.0136	Inf	0.3943	0.4474
hispanic	2hsgrad	0.4471	0.0102	Inf	0.4272	0.4671
nh_black	2hsgrad	0.4224	0.0091	Inf	0.4046	0.4404
nh_multirace	2hsgrad	0.4674	0.0187	Inf	0.4309	0.5042
nh_other	2hsgrad	0.3232	0.0122	Inf	0.2997	0.3477
nhwhite	2hsgrad	0.5067	0.0057	Inf	0.4956	0.5178
hispanic	3somecol	0.5282	0.0105	Inf	0.5077	0.5487
nh_black	3somecol	0.5032	0.0091	Inf	0.4853	0.5210
nh_multirace	3somecol	0.5486	0.0184	Inf	0.5123	0.5844
nh_other	3somecol	0.3981	0.0136	Inf	0.3717	0.4251
nhwhite	3somecol	0.5872	0.0052	Inf	0.5770	0.5974
hispanic	4colgrad	0.6519	0.0088	Inf	0.6345	0.6690
nh_black	4colgrad	0.6288	0.0082	Inf	0.6125	0.6448
nh_multirace	4colgrad	0.6703	0.0164	Inf	0.6374	0.7016
nh_other	4colgrad	0.5253	0.0133	Inf	0.4992	0.5512
nhwhite	4colgrad	0.7041	0.0035	Inf	0.6971	0.7110

From the above analyses, we see the relationship between alcohol consumption and each predictor variable is significant except for non-Hispanic (NH) multirace and some high school graduates. Furthermore, all of the odd ratios fall within 95% confidence intervals. NH Blacks are 10% less likely to drink alcohol compared to Hispanics whereas NH Whites are 27% more likely to drink alcohol compared to Hispanics. Moreover, college graduates are 270% more likely to drink alcohol compared to adults with primary education.

Generate predicted probabilities for some “interesting” cases from your analysis, to highlight the effects from the model and your stated research question.

## Case from college graduates Hispanics and NH Whites:

comps<-as.data.frame(marg_logit)

comps[comps$race_eth=="hispanic" & comps$educ == "4colgrad" , ]

comps[comps$race_eth=="nhwhite" & comps$educ == "4colgrad" ,  ]

According to the model for Hispanic and NH White college graduates, whites have greater probability (about 70%) of drinking alcohol compared to Hispanic graduates (about 65%). Hence, according to my research question, there are influence of races and educational background on probability of drinking alcohol.