library(car)
library(stargazer)
## 
## Please cite as:
##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.1. 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)
load("C:/Users/Colton/Google Drive/PhD/Stats 2/brfss16_mmsa.Rdata")
newnames<-names(brfss16m)
head(newnames, n=10)
##  [1] "DISPCODE" "HHADULT"  "GENHLTH"  "PHYSHLTH" "MENTHLTH" "POORHLTH"
##  [7] "HLTHPLN1" "PERSDOC2" "MEDCOST"  "CHECKUP1"
bestnames<-tolower(gsub(pattern = "_", replacement = "", x = newnames))
names(brfss16m)<-bestnames
head(bestnames, n=10)
##  [1] "dispcode" "hhadult"  "genhlth"  "physhlth" "menthlth" "poorhlth"
##  [7] "hlthpln1" "persdoc2" "medcost"  "checkup1"
#Educational Level
brfss16m$educ <- car::recode(brfss16m$educag, 
                           recodes = 
                             "1='0 Less than HS';
                              2='1 HS';
                              3='2 Some Col'; 
                              4='3 Col Grad'; 
                              else=NA", as.factor.result = T)

brfss16m$educ <- relevel(brfss16m$educ, ref = '3 Col Grad')

#General Health
brfss16m$badhealth <- car::recode(brfss16m$genhlth, 
                                recodes = 
                                  "1:2='2 Good'; 
                                  3='1 Average'; 
                                  4:5='0 Poor'; 
                                  else=NA", as.factor.result = T)

brfss16m$badhealth <- relevel(brfss16m$badhealth, ref = '2 Good')

#Mental Health 
brfss16m$badmh <- car::recode(brfss16m$ment14d, 
                            recodes = 
                              "1='0 Good MH'; 
                              2:3='1 Bad MH'; 
                              else=NA")


brfss16m$ordmh <- car::recode(brfss16m$ment14d, 
                              recodes = 
                              "1 ='0 Good MH';
                              2 ='1 Some MH Issues'; 
                              3 ='2 Severe MH Issues'; 
                              else = NA", as.factor.result = T)

brfss16m$othmh <- car::recode(brfss16m$ment14d, 
                              recodes = 
                                "1 ='0 Good MH';
                              2 ='1 Some MH Issues'; 
                              3 ='2 Severe MH Issues'; 
                              else = NA", as.factor.result = F)

#Race/Ethnicity
brfss16m$race_eth <- recode(brfss16m$racegr3, 
                          recodes=
                          "1='nh white'; 
                          2='nh black'; 
                          3='nh other';
                          4='nh multirace'; 
                          5='hispanic'; 
                          else=NA", as.factor.result = T) 

brfss16m$race_eth <- relevel(brfss16m$race_eth, ref = "nh white")

#Age
brfss16m$age <- car::recode(brfss16m$ageg,
                            recodes = 
                              "1 = '0 18-24';
                              2 = '1 25-34';
                              3 = '2 35-44';
                              4 = '3 45-54';
                              5 = '4 55-64';
                              6 = '5 65 & Older';
                              else = NA", as.factor.result = T)

brfss16m$age <- relevel(brfss16m$age, ref = "0 18-24")

#Sex
brfss16m$rsex <- car::recode(brfss16m$sex,
                             recodes = 
                               "1 = '0 Male';
                                2 = '1 Female';
                                else = NA", as.factor.result = T)


brfss16m$rsex <- relevel(brfss16m$rsex, ref = "0 Male")

#Sleep Quality
brfss16m$sleep <- car::recode(brfss16m$sleptim1,
                              recodes = 
                                "1:4 = '0 Little Sleep';
                                5:7 = '1 Some Sleep';
                                8 = '2 Good Sleep';
                                9:12 = '3 Extra Sleep';
                                else = NA", as.factor.result = T)

brfss16m$sleep <- relevel(brfss16m$sleep, ref = "0 Little Sleep")

#Selected out
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
subdata<-brfss16m%>% 
  select(ordmh, othmh, badhealth, educ, sleep, age, rsex, race_eth, mmsawt, ststr) %>% 
  filter( complete.cases(.))

#Survey Design
options(survey.lonely.psu = "adjust")
mhdes<-svydesign(ids=~1, strata=~ststr, weights=~mmsawt, data = subdata)

Question 1 - I will look at the ordinal variable “ment14d”. This is an ordinal variable that assess how often in the past 30 days that individuals have dealth with poor mental health. I will recode this in a new variable with zero being equal to good mental health; one representing some mental health issues; two indicating severe mental health issues; else with be calculated as missing.

Question 2 - What factor operates as a stronger indicater of mental health, self-reported general health or level of educational attainment? Other factors like quality of sleep, respondents sex, age, and race/ethnicity wil all be controlled for.

mod1 <- svyolr(ordmh ~ badhealth + educ + sleep + age + rsex + race_eth, mhdes)
summary(mod1)
## Call:
## svyolr(ordmh ~ badhealth + educ + sleep + age + rsex + race_eth, 
##     mhdes)
## 
## Coefficients:
##                            Value Std. Error     t value
## badhealth0 Poor       1.50858015 0.02816939  53.5538846
## badhealth1 Average    0.46909525 0.02055384  22.8227507
## educ0 Less than HS   -0.01854953 0.03837574  -0.4833659
## educ1 HS             -0.11506182 0.02371407  -4.8520494
## educ2 Some Col        0.04976724 0.02128617   2.3380081
## sleep1 Some Sleep    -0.85200051 0.04884136 -17.4442413
## sleep2 Good Sleep    -1.21410096 0.05082587 -23.8874601
## sleep3 Extra Sleep   -0.87314647 0.05986675 -14.5848320
## age1 25-34           -0.36824108 0.03314643 -11.1095253
## age2 35-44           -0.67163893 0.03403841 -19.7317938
## age3 45-54           -0.78949180 0.03367708 -23.4430016
## age4 55-64           -1.08370034 0.03456693 -31.3507807
## age5 65 & Older      -1.69824660 0.03491945 -48.6332544
## rsex1 Female          0.48222692 0.01826552  26.4009413
## race_ethhispanic     -0.61176433 0.03088122 -19.8102369
## race_ethnh black     -0.32473795 0.02934300 -11.0669656
## race_ethnh multirace  0.24695305 0.06664162   3.7056877
## race_ethnh other     -0.47263420 0.04386113 -10.7756950
## 
## Intercepts:
##                                     Value    Std. Error t value 
## 0 Good MH|1 Some MH Issues           -0.5884   0.0579   -10.1578
## 1 Some MH Issues|2 Severe MH Issues   1.0247   0.0576    17.8043
exp(coef(mod1))
##                     badhealth0 Poor                  badhealth1 Average 
##                           4.5203081                           1.5985472 
##                  educ0 Less than HS                            educ1 HS 
##                           0.9816215                           0.8913110 
##                      educ2 Some Col                   sleep1 Some Sleep 
##                           1.0510264                           0.4265607 
##                   sleep2 Good Sleep                  sleep3 Extra Sleep 
##                           0.2969769                           0.4176354 
##                          age1 25-34                          age2 35-44 
##                           0.6919503                           0.5108706 
##                          age3 45-54                          age4 55-64 
##                           0.4540755                           0.3383412 
##                     age5 65 & Older                        rsex1 Female 
##                           0.1830041                           1.6196773 
##                    race_ethhispanic                    race_ethnh black 
##                           0.5423931                           0.7227167 
##                race_ethnh multirace                    race_ethnh other 
##                           1.2801190                           0.6233581 
##          0 Good MH|1 Some MH Issues 1 Some MH Issues|2 Severe MH Issues 
##                           0.5551939                           2.7862651
mod1$deviance+2*length(mod1$coefficients)
## [1] 376008.9

Proportional Odds Assumption

mod2 <- svyglm(I(othmh>0) ~ badhealth + educ + sleep + age + rsex + race_eth, mhdes, family = "binomial")
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
## Warning: glm.fit: algorithm did not converge
mod2$deviance+2*length(mod2$coefficients)
## [1] 38
mod3 <- svyglm(I(othmh>1) ~ badhealth + educ + sleep + age + rsex + race_eth, mhdes, family = "binomial")
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
mod3$deviance+2*length(mod3$coefficients)
## [1] 279119.1
plot(coef(mod2)[-1], ylim = c(-2, 2), type="l", xaxt= "n",
       xlab = "", ylab = "Beta", main = c("Comparison for Betas for", "Proportional Odds Assumption"))
lines(coef(mod3)[-1], col=1, lty=2)
axis(side = 1, at=1:18, labels = F)
text(x = 1:18, y = -2.50, srt = 90, pos = 1, xpd = T, cex = .8,
     labels = c( "Poor Health", "Average Health", "Less than HS", "HS", "Some College", 
                 "Some Sleep", "Extra Sleep", "Good Sleep", "Age 25-34", "Age 35-44", "Age 45-54",
                 "Age 55-64", "Age 65 & Older", "Female", "Hispanic", "Black", "Multirace", "Other"))
legend("topright",col=c(1,1),lty=c(1,2), legend=c
       (">0", ">1"))

exp(coef(mod2)) 
##          (Intercept)      badhealth0 Poor   badhealth1 Average 
##         4.362282e+11         9.652075e-01         9.840770e-01 
##   educ0 Less than HS             educ1 HS       educ2 Some Col 
##         1.362971e+00         1.133299e+00         1.158143e+00 
##    sleep1 Some Sleep    sleep2 Good Sleep   sleep3 Extra Sleep 
##         1.033814e+00         1.024119e+00         1.021186e+00 
##           age1 25-34           age2 35-44           age3 45-54 
##         9.452035e-01         8.966405e-01         8.289281e-01 
##           age4 55-64      age5 65 & Older         rsex1 Female 
##         7.650621e-01         7.128808e-01         9.622590e-01 
##     race_ethhispanic     race_ethnh black race_ethnh multirace 
##         1.254483e+00         1.151268e+00         9.554734e-01 
##     race_ethnh other 
##         1.394303e+00
exp(coef(mod3)) 
##          (Intercept)      badhealth0 Poor   badhealth1 Average 
##            1.6833420            3.7604130            1.5557656 
##   educ0 Less than HS             educ1 HS       educ2 Some Col 
##            0.9068982            0.8314859            0.9992219 
##    sleep1 Some Sleep    sleep2 Good Sleep   sleep3 Extra Sleep 
##            0.5246132            0.3655541            0.4934456 
##           age1 25-34           age2 35-44           age3 45-54 
##            0.6466520            0.4756014            0.4153658 
##           age4 55-64      age5 65 & Older         rsex1 Female 
##            0.3051856            0.1726553            1.6398144 
##     race_ethhispanic     race_ethnh black race_ethnh multirace 
##            0.5580529            0.7302345            1.2293104 
##     race_ethnh other 
##            0.6211178

I find that the assumptions of proportional odds model are not met. This is confirmed by the plotted graph of the models indicating that they are drastically different. Additionally, the odds ratios of the variables within the models vary by a great degree.

stargazer(mod3, apply.coef = exp, t.auto = F, p.auto = F, report = "vctsp*", type = "html", style = "demography", 
          covariate.labels = c("Poor Health", "Average Health", "Less than HS", "HS", "Some College", 
                              "Some Sleep", "Extra Sleep", "Good Sleep", "Age 25-34", "Age 35-44", "Age 45-54",
                               "Age 55-64", "Age 65 Older", "Female", "Hispanic", "Black", "Multirace", "Other"))
I(othmh > 1)
Poor Health 3.760
t = 46.717
(0.028)
p = 0.000***
Average Health 1.556
t = 20.475
(0.022)
p = 0.000***
Less than HS 0.907
t = -2.556
(0.038)
p = 0.011*
HS 0.831
t = -7.548
(0.024)
p = 0.000***
Some College 0.999
t = -0.034
(0.023)
p = 0.973
Some Sleep 0.525
t = -13.474
(0.048)
p = 0.000***
Extra Sleep 0.366
t = -20.176
(0.050)
p = 0.000***
Good Sleep 0.493
t = -12.145
(0.058)
p = 0.000***
Age 25-34 0.647
t = -11.862
(0.037)
p = 0.000***
Age 35-44 0.476
t = -19.905
(0.037)
p = 0.000***
Age 45-54 0.415
t = -23.726
(0.037)
p = 0.000***
Age 55-64 0.305
t = -32.181
(0.037)
p = 0.000***
Age 65 Older 0.173
t = -47.488
(0.037)
p = 0.000***
Female 1.640
t = 26.380
(0.019)
p = 0.000***
Hispanic 0.558
t = -18.697
(0.031)
p = 0.000***
Black 0.730
t = -10.433
(0.030)
p = 0.000***
Multirace 1.229
t = 3.028
(0.068)
p = 0.003**
Other 0.621
t = -10.461
(0.046)
p = 0.000***
Constant 1.683
t = 8.953
(0.058)
p = 0.000***
N 236,236
Log Likelihood -130,403.700
AIC 260,845.300
p < .05; p < .01; p < .001

Mod3 present the best fitting model. I found mod3 to present the best fitted model due to it having the lowest AIC score. In this model I assess the effect of multiple covariates impact on mental health. The findings reveal that those who report poor self-rated general health have 3.8 higher odds of reporting poor mental health in the past 30 days compared to those with good health. Interestingly, individuals that indicate attaining either less than high school or a high school education are at lower odds of reporting mental health when compared to college graduates. Regarding my initial research question, it seems that self-rated health tends to carry a stronger impact on mental health, rather than educational attainment. Furthermore, those who get some to extra hours of sleep have lower odds of reporting poor mental health compared to those that get little to no sleep. Age seems to share an inverse relationship with poor mental health, in that, as age increases odds of reporting poor mental health decrease. Females are at higher odds of reporting poor mental health compared to males. Regarding race, there is an interesting finding. Multirace has higher odds of reporting poor mental health than any other race, while blacks and hispanics, historical minorities, have lower odds of reporting poor mental health. ```