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
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. ```