#Data Prep
nams<-names(brfss_17)
head(nams, n=10)
## [1] "dispcode" "statere1" "safetime" "hhadult" "genhlth" "physhlth"
## [7] "menthlth" "poorhlth" "hlthpln1" "persdoc2"
newnames<-tolower(gsub(pattern = "_",replacement = "",x = nams))
names(brfss_17)<-newnames
##Get TX Cities
brfss_17$tx<-NA
brfss_17$tx[grep(pattern = "TX", brfss_17$mmsaname)]<-1
##Remove Non-Responses
brfss_17<-brfss_17%>%
filter(tx==1, is.na(mmsawt)==F, is.na(hlthpln1)==F)
##Recodes
brfss_17$hlthpln1<-Recode(brfss_17$hlthpln1, recodes ="7:9=NA; 1=1;2=0")
brfss_17$educ<-Recode(brfss_17$educa, recodes="1:3='0hs'; 4='2hsgrad'; 5='3somecol'; 6='4colgrad';9=NA", as.factor=T)
##Gender variable recode
brfss_17$gender<-as.factor(ifelse(brfss_17$sex==1, "Male", "Female"))
View(brfss_17)
library(tableone)
des<-svydesign(ids=~1, strata=~ststr, weights=~mmsawt, data = brfss_17)
options(survey.lonely.psu = "adjust")
st1<-svyCreateTableOne(vars = c("gender", "educ"), strata = "hlthpln1", test = T, data = des)
#st1<-print(st1, format="p")
print(st1, format="p")
## Stratified by hlthpln1
## 0 1 p test
## n 3579519.3 11755198.9
## gender = Male (%) 56.2 46.4 0.001
## educ (%) <0.001
## 0hs 36.0 10.5
## 2hsgrad 31.9 22.6
## 3somecol 24.7 32.9
## 4colgrad 7.5 34.0
#Homework #3
#Create a research question and perform logistic regression:
##Research question: How does gender and education impact health insurance status?
library(survey)
fit.logit<-svyglm(hlthpln1 ~ gender + educ,
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)
(Intercept) |
0.157 |
0.147 |
1.069 |
0.285 |
genderMale |
-0.403 |
0.124 |
-3.236 |
0.001 |
educ2hsgrad |
0.901 |
0.167 |
5.398 |
0.000 |
educ3somecol |
1.513 |
0.176 |
8.594 |
0.000 |
educ4colgrad |
2.763 |
0.186 |
14.824 |
0.000 |
#Present results from a model with sample weights and design effects, if your data allow for this.
library(broom)
fit.logit%>%
tidy()%>%
knitr::kable(digits = 3)
(Intercept) |
0.157 |
0.147 |
1.069 |
0.285 |
genderMale |
-0.403 |
0.124 |
-3.236 |
0.001 |
educ2hsgrad |
0.901 |
0.167 |
5.398 |
0.000 |
educ3somecol |
1.513 |
0.176 |
8.594 |
0.000 |
educ4colgrad |
2.763 |
0.186 |
14.824 |
0.000 |
fit.logit%>%
tidy()%>%
mutate(OR = exp(estimate))%>%
knitr::kable(digits = 3)
(Intercept) |
0.157 |
0.147 |
1.069 |
0.285 |
1.170 |
genderMale |
-0.403 |
0.124 |
-3.236 |
0.001 |
0.669 |
educ2hsgrad |
0.901 |
0.167 |
5.398 |
0.000 |
2.461 |
educ3somecol |
1.513 |
0.176 |
8.594 |
0.000 |
4.539 |
educ4colgrad |
2.763 |
0.186 |
14.824 |
0.000 |
15.854 |
#Present the results in tabular form, with Parameter estimates, odds ratios (if using the logit model) and confidence intervals for the odds ratios.
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)
(Intercept) |
0.157 |
0.147 |
1.069 |
0.285 |
1.170 |
0.877 |
1.560 |
genderMale |
-0.403 |
0.124 |
-3.236 |
0.001 |
0.669 |
0.524 |
0.853 |
educ2hsgrad |
0.901 |
0.167 |
5.398 |
0.000 |
2.461 |
1.775 |
3.414 |
educ3somecol |
1.513 |
0.176 |
8.594 |
0.000 |
4.539 |
3.215 |
6.410 |
educ4colgrad |
2.763 |
0.186 |
14.824 |
0.000 |
15.854 |
11.002 |
22.847 |
#Generate predicted probabilities for some “interesting” cases from your analysis, to highlight the effects from the model and your stated research question.
library(emmeans)
rg<-ref_grid(fit.logit, data = brfss_17)
marg_logit<-emmeans(object = rg,
specs = c( "gender", "educ"),
type="response" )
knitr::kable(marg_logit, digits = 4)
Female |
0hs |
0.5392 |
0.0365 |
Inf |
0.4673 |
0.6094 |
Male |
0hs |
0.4389 |
0.0348 |
Inf |
0.3722 |
0.5079 |
Female |
2hsgrad |
0.7423 |
0.0238 |
Inf |
0.6929 |
0.7861 |
Male |
2hsgrad |
0.6582 |
0.0268 |
Inf |
0.6039 |
0.7086 |
Female |
3somecol |
0.8416 |
0.0181 |
Inf |
0.8028 |
0.8739 |
Male |
3somecol |
0.7803 |
0.0226 |
Inf |
0.7327 |
0.8214 |
Female |
4colgrad |
0.9488 |
0.0073 |
Inf |
0.9325 |
0.9614 |
Male |
4colgrad |
0.9254 |
0.0099 |
Inf |
0.9035 |
0.9426 |