#load brfss
library(car)
library(stargazer)
library(survey)
library(sjPlot)
library(ggplot2)
library(pander)
library(knitr)
library(lme4)
load("~/CSparks_StatsII/brfss16_mmsa.Rdata")
set.seed(12345)
#samps<-sample(1:nrow(brfss16m), size = 40000, replace=F)
#brfss16m<-brfss16m[samps,]
#The names in the data are very ugly, so I make them less ugly
nams<-names(brfss16m)
#we see some names are lower case, some are upper and some have a little _ in the first position. This is a nightmare.
newnames<-gsub(pattern = "_",replacement = "",x = nams)
names(brfss16m)<-tolower(newnames)
Recode variables
#sex
brfss16m$female<-ifelse(brfss16m$sex==2, 1, 0)
#BMI
brfss16m$bmi<-ifelse(is.na(brfss16m$bmi5)==T, NA, brfss16m$bmi5/100)
#Healthy days
brfss16m$healthdays<-recode(brfss16m$physhlth, recodes = "88=0; 77=NA; 99=NA")
#Healthy mental health days
brfss16m$healthmdays<-recode(brfss16m$menthlth, recodes = "88=0; 77=NA; 99=NA")
brfss16m$badhealth<-recode(brfss16m$genhlth, recodes="4:5=1; 1:3=0; else=NA")
#race/ethnicity
brfss16m$black<-recode(brfss16m$racegr3, recodes="2=1; 9=NA; else=0")
brfss16m$white<-recode(brfss16m$racegr3, recodes="1=1; 9=NA; else=0")
brfss16m$other<-recode(brfss16m$racegr3, recodes="3:4=1; 9=NA; else=0")
brfss16m$hispanic<-recode(brfss16m$racegr3, recodes="5=1; 9=NA; else=0")
brfss16m$race_eth<-recode(brfss16m$racegr3,
recodes="1='nhwhite'; 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 = "nhwhite")
#insurance
brfss16m$ins<-recode(brfss16m$hlthpln1, recodes ="7:9=NA; 1=1;2=0")
#income grouping
brfss16m$inc<-recode(brfss16m$incomg, recodes = "9= NA;1='1_lt15k'; 2='2_15-25k';3='3_25-35k';4='4_35-50k';5='5_50kplus'", as.factor.result = T)
brfss16m$inc<-as.ordered(brfss16m$inc)
#education level
brfss16m$educ<-recode(brfss16m$educa,
recodes="1:2='0Prim'; 3='1somehs'; 4='2hsgrad'; 5='3somecol'; 6='4colgrad';9=NA",
as.factor.result=T)
brfss16m$educ<-relevel(brfss16m$educ, ref='2hsgrad')
#employloyment
brfss16m$employ<-recode(brfss16m$employ1,
recodes="1:2='employloyed'; 2:6='nilf'; 7='retired'; 8='unable'; else=NA",
as.factor.result=T)
brfss16m$employ<-relevel(brfss16m$employ, ref='employloyed')
#marital status
brfss16m$marst<-recode(brfss16m$marital,
recodes="1='married'; 2='divorced'; 3='widowed'; 4='separated'; 5='nm';6='cohab'; else=NA",
as.factor.result=T)
brfss16m$marst<-relevel(brfss16m$marst, ref='married')
#Age cut into intervals
brfss16m$agec<-cut(brfss16m$age80, breaks=c(0,24,39,59,79,99))
#BMI, in the brfss16ma the bmi variable has 2 implied decimal places,
#so we must divide by 100 to get real bmi's
brfss16m$bmi<-brfss16m$bmi5/100
#smoking currently
brfss16m$smoke<-recode(brfss16m$smoker3,
recodes="1:2=1; 3:4=0; else=NA")
#brfss16m$smoke<-relevel(brfss16m$smoke, ref = "NeverSmoked")
brfss16m$obese<-ifelse(is.na(brfss16m$bmi)==T, NA,
ifelse(brfss16m$bmi>30,1,0))
First multilevel regression model with the MSA and determine length of MSA
head(data.frame(name=table(brfss16m$mmsaname),id=unique(brfss16m$mmsa)))
## name.Var1
## 1 Albany-Schenectady-Troy, NY, Metropolitan Statistical Area
## 2 Albuquerque, NM, Metropolitan Statistical Area
## 3 Allentown-Bethlehem-Easton, PA-NJ, Metropolitan Statistical Area
## 4 Anchorage, AK, Metropolitan Statistical Area
## 5 Atlanta-Sandy Springs-Roswell, GA, Metropolitan Statistical Area
## 6 Augusta-Richmond County, GA-SC, Metropolitan Statistical Area
## name.Freq id
## 1 2901 10580
## 2 1300 10740
## 3 675 10900
## 4 1028 11260
## 5 2444 12060
## 6 849 12260
length(table(brfss16m$mmsa))
## [1] 143
#MSAs
Higher level predictors
library(acs)
#Get 2010 ACS median household incomes for tracts in Texas
msaacs<-geo.make(msa="*")
#ACS tables B17 is poverty, B19 is Gini, B25001 is housing vacancy, B25035 is median year built
acsecon<-acs.fetch(key=mykey, endyear=2010, span=5, geography=msaacs, variable = c("B19083_001","B17001_001","B17001_002" , "B25002_001","B25002_003", "B25035_001"))
colnames(acsecon@estimate)
## [1] "B19083_001" "B17001_001" "B17001_002" "B25002_001" "B25002_003"
## [6] "B25035_001"
msaecon<-data.frame(gini=acsecon@estimate[, "B19083_001"],
ppoverty=acsecon@estimate[, "B17001_002"]/acsecon@estimate[, "B17001_001"],
giniz=scale(acsecon@estimate[, "B19083_001"]),
pvacant=acsecon@estimate[,"B25002_003"]/acsecon@estimate[, "B25002_001"],
ppovertyz=scale(acsecon@estimate[, "B17001_002"]/acsecon@estimate[, "B17001_001"]),
pvacantz=scale(acsecon@estimate[,"B25002_003"]/acsecon@estimate[, "B25002_001"]),
medhouse=acsecon@estimate[, "B25035_001" ],
medhousez=scale(acsecon@estimate[, "B25035_001" ]))
msaecon$ids<-paste(acsecon@geography$metropolitanstatisticalareamicropolitanstatisticalarea)
head(msaecon)
## gini ppoverty giniz
## Adjuntas, PR Micro Area 0.525 0.5925656 2.718785
## Aguadilla-Isabela-San Sebastián, PR Metro Area 0.529 0.5577311 2.847562
## Coamo, PR Micro Area 0.532 0.5633984 2.944145
## Fajardo, PR Metro Area 0.483 0.4309070 1.366623
## Guayama, PR Metro Area 0.518 0.4980518 2.493425
## Jayuya, PR Micro Area 0.502 0.5446137 1.978315
## pvacant ppovertyz
## Adjuntas, PR Micro Area 0.1996029 6.271934
## Aguadilla-Isabela-San Sebastián, PR Metro Area 0.2005956 5.762039
## Coamo, PR Micro Area 0.1500353 5.844994
## Fajardo, PR Metro Area 0.3143172 3.905632
## Guayama, PR Metro Area 0.1951184 4.888474
## Jayuya, PR Micro Area 0.1562109 5.570031
## pvacantz medhouse
## Adjuntas, PR Micro Area 0.9053869 1978
## Aguadilla-Isabela-San Sebastián, PR Metro Area 0.9199516 1979
## Coamo, PR Micro Area 0.1781186 1979
## Fajardo, PR Metro Area 2.5885036 1977
## Guayama, PR Metro Area 0.8395882 1978
## Jayuya, PR Micro Area 0.2687281 1978
## medhousez ids
## Adjuntas, PR Micro Area 0.4923898 10260
## Aguadilla-Isabela-San Sebastián, PR Metro Area 0.5982027 10380
## Coamo, PR Micro Area 0.5982027 17620
## Fajardo, PR Metro Area 0.3865770 21940
## Guayama, PR Metro Area 0.4923898 25020
## Jayuya, PR Micro Area 0.4923898 27580
summary(msaecon)
## gini ppoverty giniz pvacant
## Min. :0.3510 Min. :0.02425 Min. :-2.88303 Min. :0.04246
## 1st Qu.:0.4190 1st Qu.:0.12308 1st Qu.:-0.69381 1st Qu.:0.09454
## Median :0.4390 Median :0.15172 Median :-0.04993 Median :0.12076
## Mean :0.4406 Mean :0.16409 Mean : 0.00000 Mean :0.13790
## 3rd Qu.:0.4590 3rd Qu.:0.18746 3rd Qu.: 0.59396 3rd Qu.:0.15964
## Max. :0.5730 Max. :0.59257 Max. : 4.26411 Max. :0.63877
## ppovertyz pvacantz medhouse medhousez
## Min. :-2.0468 Min. :-1.4002 Min. :1940 Min. :-3.5285
## 1st Qu.:-0.6002 1st Qu.:-0.6361 1st Qu.:1968 1st Qu.:-0.5657
## Median :-0.1810 Median :-0.2514 Median :1975 Median : 0.1750
## Mean : 0.0000 Mean : 0.0000 Mean :1973 Mean : 0.0000
## 3rd Qu.: 0.3421 3rd Qu.: 0.3190 3rd Qu.:1980 3rd Qu.: 0.7040
## Max. : 6.2719 Max. : 7.3489 Max. :1998 Max. : 2.6086
## ids
## Length:955
## Class :character
## Mode :character
##
##
##
library(tigris)
## Warning: package 'tigris' was built under R version 3.4.4
## To enable
## caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.
##
## Attaching package: 'tigris'
## The following object is masked from 'package:graphics':
##
## plot
msa<-core_based_statistical_areas(cb=T)
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 6%
|
|==== | 7%
|
|===== | 7%
|
|===== | 8%
|
|====== | 9%
|
|====== | 10%
|
|======= | 10%
|
|======= | 11%
|
|======== | 12%
|
|======== | 13%
|
|========= | 13%
|
|========= | 14%
|
|========== | 15%
|
|========== | 16%
|
|=========== | 16%
|
|=========== | 17%
|
|=========== | 18%
|
|============ | 18%
|
|============ | 19%
|
|============= | 19%
|
|============= | 20%
|
|============= | 21%
|
|============== | 21%
|
|============== | 22%
|
|=============== | 23%
|
|=============== | 24%
|
|================ | 24%
|
|================ | 25%
|
|================= | 26%
|
|================= | 27%
|
|================== | 27%
|
|================== | 28%
|
|=================== | 29%
|
|=================== | 30%
|
|==================== | 30%
|
|==================== | 31%
|
|===================== | 32%
|
|===================== | 33%
|
|====================== | 33%
|
|====================== | 34%
|
|======================= | 35%
|
|======================= | 36%
|
|======================== | 36%
|
|======================== | 37%
|
|======================== | 38%
|
|========================= | 38%
|
|========================= | 39%
|
|========================== | 39%
|
|========================== | 40%
|
|========================== | 41%
|
|=========================== | 41%
|
|=========================== | 42%
|
|============================ | 43%
|
|============================ | 44%
|
|============================= | 44%
|
|============================= | 45%
|
|============================== | 46%
|
|============================== | 47%
|
|=============================== | 47%
|
|=============================== | 48%
|
|================================ | 49%
|
|================================ | 50%
|
|================================= | 50%
|
|================================= | 51%
|
|================================== | 52%
|
|================================== | 53%
|
|=================================== | 53%
|
|=================================== | 54%
|
|==================================== | 55%
|
|==================================== | 56%
|
|===================================== | 56%
|
|===================================== | 57%
|
|===================================== | 58%
|
|====================================== | 58%
|
|====================================== | 59%
|
|======================================= | 59%
|
|======================================= | 60%
|
|======================================= | 61%
|
|======================================== | 61%
|
|======================================== | 62%
|
|========================================= | 63%
|
|========================================= | 64%
|
|========================================== | 64%
|
|========================================== | 65%
|
|=========================================== | 66%
|
|=========================================== | 67%
|
|============================================ | 67%
|
|============================================ | 68%
|
|============================================= | 69%
|
|============================================= | 70%
|
|============================================== | 70%
|
|============================================== | 71%
|
|============================================== | 72%
|
|=============================================== | 72%
|
|=============================================== | 73%
|
|================================================ | 73%
|
|================================================ | 74%
|
|================================================= | 75%
|
|================================================= | 76%
|
|================================================== | 76%
|
|================================================== | 77%
|
|================================================== | 78%
|
|=================================================== | 78%
|
|=================================================== | 79%
|
|==================================================== | 79%
|
|==================================================== | 80%
|
|==================================================== | 81%
|
|===================================================== | 81%
|
|===================================================== | 82%
|
|====================================================== | 83%
|
|====================================================== | 84%
|
|======================================================= | 84%
|
|======================================================= | 85%
|
|======================================================== | 86%
|
|======================================================== | 87%
|
|========================================================= | 87%
|
|========================================================= | 88%
|
|========================================================== | 89%
|
|========================================================== | 90%
|
|=========================================================== | 90%
|
|=========================================================== | 91%
|
|=========================================================== | 92%
|
|============================================================ | 92%
|
|============================================================ | 93%
|
|============================================================= | 93%
|
|============================================================= | 94%
|
|============================================================== | 95%
|
|============================================================== | 96%
|
|=============================================================== | 96%
|
|=============================================================== | 97%
|
|=============================================================== | 98%
|
|================================================================ | 98%
|
|================================================================ | 99%
|
|=================================================================| 99%
|
|=================================================================| 100%
msa_ec<-geo_join(msa, msaecon, "CBSAFP", "ids", how="inner")
Merge the MSA data to the BRFSS data
Research Question: How do socioeconomic factors affect the odds of having poor mental health?
msa_ec$state<-str_sub(msa_ec$NAME, -2,-1)
msaecon<-merge(msaecon, msa_ec@data[, c(1:7,18)], by.x="ids", by.y="CBSAFP", all.x=T)
joindata<-merge(brfss16m, msaecon, by.x="mmsa",by.y="ids", all.x=T)
joindata$healthmdaysz<-as.numeric(scale(joindata$healthmdays, center=T, scale=T))
#joindata<-joindata[complete.cases(joindata[, c("bmiz", "race_eth", "agec", "educ", "gini")]),]
#and merge the data back to the kids data
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:acs':
##
## combine
## 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
joindata<-joindata%>%
select(healthmdaysz, marst, female, mmsa, agec, educ, race_eth, employ, healthmdays, badhealth,bmi,gini, ppoverty, ppovertyz,giniz,medhouse,medhousez, pvacant, pvacantz, mmsawt, mmsaname )%>%
filter(complete.cases(.))
head(joindata[, c("healthmdaysz", "marst", "female", "educ","ppoverty", "mmsa")])
## healthmdaysz marst female educ ppoverty mmsa
## 1 0.2095226 widowed 1 3somecol 0.1053253 10580
## 2 -0.4434798 married 1 2hsgrad 0.1053253 10580
## 3 -0.4434798 divorced 0 4colgrad 0.1053253 10580
## 4 -0.3128793 divorced 0 4colgrad 0.1053253 10580
## 5 -0.4434798 widowed 0 2hsgrad 0.1053253 10580
## 6 -0.4434798 widowed 1 3somecol 0.1053253 10580
meanhealthmdays<-mean(joindata$healthmdays, na.rm=T)
sdhealthmdays<-sd(joindata$healthmdays, na.rm=T)
MSA level variables
fit.mix2<-lmer(healthmdaysz~marst+female+educ+race_eth+ppovertyz+(1|mmsa), data=joindata)
arm::display(fit.mix2, detail=T)
## lmer(formula = healthmdaysz ~ marst + female + educ + race_eth +
## ppovertyz + (1 | mmsa), data = joindata)
## coef.est coef.se t value
## (Intercept) -0.10 0.01 -11.87
## marstcohab 0.27 0.01 19.76
## marstdivorced 0.24 0.01 33.25
## marstnm 0.26 0.01 38.44
## marstseparated 0.54 0.02 31.16
## marstwidowed 0.01 0.01 0.66
## female 0.15 0.00 30.83
## educ0Prim 0.15 0.02 8.52
## educ1somehs 0.22 0.01 17.64
## educ3somecol -0.01 0.01 -1.28
## educ4colgrad -0.16 0.01 -25.37
## race_ethhispanic -0.04 0.01 -4.16
## race_ethnh black -0.07 0.01 -7.39
## race_ethnh multirace 0.26 0.02 13.33
## race_ethnh other 0.01 0.01 0.88
## ppovertyz 0.03 0.01 3.45
##
## Error terms:
## Groups Name Std.Dev.
## mmsa (Intercept) 0.05
## Residual 0.99
## ---
## number of obs: 174058, groups: mmsa, 123
## AIC = 490341, DIC = 490063.3
## deviance = 490184.0
fit.mix4<-lmer(healthmdaysz~marst+female+educ+race_eth+ppovertyz*race_eth+(1|mmsa), joindata, REML=F)
arm::display(fit.mix4, detail=T)
## lmer(formula = healthmdaysz ~ marst + female + educ + race_eth +
## ppovertyz * race_eth + (1 | mmsa), data = joindata, REML = F)
## coef.est coef.se t value
## (Intercept) -0.09 0.01 -9.69
## marstcohab 0.27 0.01 19.77
## marstdivorced 0.24 0.01 33.25
## marstnm 0.26 0.01 38.46
## marstseparated 0.54 0.02 31.15
## marstwidowed 0.01 0.01 0.66
## female 0.15 0.00 30.84
## educ0Prim 0.15 0.02 8.51
## educ1somehs 0.22 0.01 17.63
## educ3somecol -0.01 0.01 -1.28
## educ4colgrad -0.16 0.01 -25.34
## race_ethhispanic -0.06 0.01 -5.06
## race_ethnh black -0.08 0.01 -5.75
## race_ethnh multirace 0.29 0.03 8.38
## race_ethnh other 0.01 0.02 0.47
## ppovertyz 0.05 0.01 4.50
## race_ethhispanic:ppovertyz -0.04 0.01 -2.86
## race_ethnh black:ppovertyz -0.02 0.02 -0.87
## race_ethnh multirace:ppovertyz 0.06 0.05 1.11
## race_ethnh other:ppovertyz 0.00 0.03 -0.05
##
## Error terms:
## Groups Name Std.Dev.
## mmsa (Intercept) 0.04
## Residual 0.99
## ---
## number of obs: 174058, groups: mmsa, 123
## AIC = 490218, DIC = 490174.2
## deviance = 490174.2
Hispanics and blacks in MSAs with higher than average poverty, have fewer mental health days, while the mental health days of other groups is affected by poverty (is this correct? Not sure I interpreted correctly)
Compare the models with a LRT
anova(fit.mix2, fit.mix4)
## refitting model(s) with ML (instead of REML)
## Data: joindata
## Models:
## fit.mix2: healthmdaysz ~ marst + female + educ + race_eth + ppovertyz +
## fit.mix2: (1 | mmsa)
## fit.mix4: healthmdaysz ~ marst + female + educ + race_eth + ppovertyz *
## fit.mix4: race_eth + (1 | mmsa)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit.mix2 18 490220 490401 -245092 490184
## fit.mix4 22 490218 490440 -245087 490174 9.868 4 0.04271 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Test if GLMM is justified
anova(glm(I(healthmdaysz>0)~ mmsa, joindata, family = binomial), test= "Chisq")
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: I(healthmdaysz > 0)
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 174057 177370
## mmsa 1 0.48749 174056 177370 0.485
Research question: How do sociodemographic factors like education and race and ethnicity affect the odds of someone experiencing poor mental health in U.S. cities?
Multilevel Logistic Model
joindata$mmsa_fac<-as.factor(joindata$mmsa)
fit.mix.bin<-glmer(I(healthmdaysz>0)~ female+educ+race_eth+(1|mmsa_fac), family=binomial, joindata,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
fit.mix.bin<-refit(fit.mix.bin)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00142437 (tol =
## 0.001, component 1)
arm::display(fit.mix.bin, detail=T)
## glmer(formula = I(healthmdaysz > 0) ~ female + educ + race_eth +
## (1 | mmsa_fac), data = joindata, family = binomial, control = glmerControl(optimizer = "bobyqa",
## optCtrl = list(maxfun = 2e+05)))
## coef.est coef.se z value Pr(>|z|)
## (Intercept) -1.53 0.02 -79.83 0.00
## female 0.44 0.01 35.49 0.00
## educ0Prim 0.16 0.04 3.93 0.00
## educ1somehs 0.38 0.03 13.77 0.00
## educ3somecol 0.03 0.02 2.21 0.03
## educ4colgrad -0.34 0.02 -22.31 0.00
## race_ethhispanic 0.09 0.02 3.80 0.00
## race_ethnh black 0.02 0.02 1.06 0.29
## race_ethnh multirace 0.62 0.04 14.72 0.00
## race_ethnh other 0.08 0.03 2.33 0.02
##
## Error terms:
## Groups Name Std.Dev.
## mmsa_fac (Intercept) 0.14
## Residual 1.00
## ---
## number of obs: 174058, groups: mmsa_fac, 123
## AIC = 174327, DIC = 173749
## deviance = 174027.1
#odds ratios
exp(fixef(fit.mix.bin)[-1])
## female educ0Prim educ1somehs
## 1.5469795 1.1787001 1.4618929
## educ3somecol educ4colgrad race_ethhispanic
## 1.0353841 0.7083097 1.0984831
## race_ethnh black race_ethnh multirace race_ethnh other
## 1.0235395 1.8518392 1.0841388
exp(confint(fit.mix.bin, method="Wald"))
## 2.5 % 97.5 %
## .sig01 NA NA
## (Intercept) 0.2089223 0.2252035
## female 1.5101450 1.5847124
## educ0Prim 1.0858029 1.2795452
## educ1somehs 1.3849979 1.5430571
## educ3somecol 1.0039277 1.0678261
## educ4colgrad 0.6871669 0.7301030
## race_ethhispanic 1.0465929 1.1529460
## race_ethnh black 0.9804323 1.0685420
## race_ethnh multirace 1.7059311 2.0102268
## race_ethnh other 1.0127997 1.1605028
fit.mix.bin2<-glmer(I(healthmdaysz>0)~ female+educ+race_eth+ ppovertyz+(1|mmsa_fac), family=binomial, joindata,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
fit.mix.bin2<-refit(fit.mix.bin2)
arm::display(fit.mix.bin2, detail=T)
## glmer(formula = I(healthmdaysz > 0) ~ female + educ + race_eth +
## ppovertyz + (1 | mmsa_fac), data = joindata, family = binomial,
## control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05)))
## coef.est coef.se z value Pr(>|z|)
## (Intercept) -1.52 0.02 -70.13 0.00
## female 0.44 0.01 35.46 0.00
## educ0Prim 0.16 0.04 3.94 0.00
## educ1somehs 0.38 0.03 13.76 0.00
## educ3somecol 0.03 0.02 2.21 0.03
## educ4colgrad -0.34 0.02 -22.31 0.00
## race_ethhispanic 0.09 0.02 3.74 0.00
## race_ethnh black 0.02 0.02 1.05 0.29
## race_ethnh multirace 0.62 0.04 14.71 0.00
## race_ethnh other 0.08 0.03 2.33 0.02
## ppovertyz 0.01 0.02 0.33 0.74
##
## Error terms:
## Groups Name Std.Dev.
## mmsa_fac (Intercept) 0.14
## Residual 1.00
## ---
## number of obs: 174058, groups: mmsa_fac, 123
## AIC = 174329, DIC = 173748.5
## deviance = 174026.8
#odds ratios
exp(fixef(fit.mix.bin2)[-1])
## female educ0Prim educ1somehs
## 1.5469098 1.1786720 1.4618351
## educ3somecol educ4colgrad race_ethhispanic
## 1.0353850 0.7083098 1.0973320
## race_ethnh black race_ethnh multirace race_ethnh other
## 1.0232986 1.8516580 1.0841354
## ppovertyz
## 1.0079409
exp(confint(fit.mix.bin, method="Wald"))
## 2.5 % 97.5 %
## .sig01 NA NA
## (Intercept) 0.2089223 0.2252035
## female 1.5101450 1.5847124
## educ0Prim 1.0858029 1.2795452
## educ1somehs 1.3849979 1.5430571
## educ3somecol 1.0039277 1.0678261
## educ4colgrad 0.6871669 0.7301030
## race_ethhispanic 1.0465929 1.1529460
## race_ethnh black 0.9804323 1.0685420
## race_ethnh multirace 1.7059311 2.0102268
## race_ethnh other 1.0127997 1.1605028
#LRT between models
anova(fit.mix.bin, fit.mix.bin2)
## Data: joindata
## Models:
## fit.mix.bin: I(healthmdaysz > 0) ~ female + educ + race_eth + (1 | mmsa_fac)
## fit.mix.bin2: I(healthmdaysz > 0) ~ female + educ + race_eth + ppovertyz +
## fit.mix.bin2: (1 | mmsa_fac)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit.mix.bin 11 174327 174438 -87153 174305
## fit.mix.bin2 12 174329 174450 -87153 174305 0.1098 1 0.7403
highmsa<-levels(joindata$mmsa_fac)[as.numeric(lapply(ranef(fit.mix.bin)$mmsa, which.max))]
lowmsa<-levels(joindata$mmsa_fac)[as.numeric(lapply(ranef(fit.mix.bin)$mmsa, which.min))]
anova(fit.mix.bin, fit.mix.bin2)
## Data: joindata
## Models:
## fit.mix.bin: I(healthmdaysz > 0) ~ female + educ + race_eth + (1 | mmsa_fac)
## fit.mix.bin2: I(healthmdaysz > 0) ~ female + educ + race_eth + ppovertyz +
## fit.mix.bin2: (1 | mmsa_fac)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit.mix.bin 11 174327 174438 -87153 174305
## fit.mix.bin2 12 174329 174450 -87153 174305 0.1098 1 0.7403
Interpretation: Individuals living in poverty, ethnic minorities and less than college degree affects the probablility of someone experiencing poor mental health, net of individual factors.