Assignment 3

1. Define a binary outcome variable of your choosing and define how you recode the original variable.

Proposed topic: Advserse health conditions associated with opportunity youths

Outcome variabels: Health conditions such Health status, depression, anxiety, smoking etc. The purpose of this assignment I would focus on self rated health status

2 State a research question about what factors you believe will affect your outcome variable.

Qst: Is Opportunity youth status associated with Fair/Poor Health status?

Define at least 2 predictor variables, based on your research question

– Opportunity youth status – Gender – Urban-rural status – Racial ethnicity

fit1<-lm(badhealth~opportunity_youth_cat+urban_rural,
         data=data)

fit2<-lm(badhealth~opportunity_youth_cat+urban_rural,
         data=data,
         weights = sampweight)

fit3<-svyglm(badhealth~opportunity_youth_cat+urban_rural,
             design = des, 
             family=gaussian)

stargazer(fit1, fit2, fit3,
          style="demography", type="html",
          column.labels = c("OLS", "Weights Only", "Survey Design"),
          title = "Regression models for Self rated health using survey data - IPUMS Health survey", 
          covariate.labels=c("Opportunity youth", "Not opportunity youth", "urban", "rural"),
          
          keep.stat="n",
          model.names=F, 
          align=T,
          ci=T)

Regression models for Self rated health using survey data - IPUMS Health survey

	badhealth
	OLS	Weights Only	Survey Design
	Model 1	Model 2	Model 3
Opportunity youth	-0.077***	-0.060***	-0.060***
	(-0.101, -0.054)	(-0.081, -0.039)	(-0.093, -0.026)
Not opportunity youth	-0.026*	-0.030**	-0.030
	(-0.047, -0.004)	(-0.052, -0.008)	(-0.060, 0.001)
urban	0.146***	0.129***	0.129***
	(0.117, 0.175)	(0.102, 0.156)	(0.089, 0.169)
N	3,771	3,771	3,771
p < .05; p < .01; p < .001

library(ggplot2)
library(dplyr)
coefs<-data.frame(coefs=c(coef(fit1)[-1], coef(fit3)[-1]),
                  mod=c(rep("Non Survey Model", 8),rep("Survey Model", 8)),
                  effect=rep(names(coef(fit1)[-1]), 2))

coefs%>%
  ggplot()+
  geom_point(aes( x=effect,
                  y=coefs,
                  group=effect,
                  color=effect,
                  shape=mod),
             position=position_jitterdodge(jitter.width = 1),
             size=2)+
  ylab("Regression Coefficient")+
  xlab("Beta")+
  geom_abline(intercept = 0, slope=0)+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  ggtitle(label = "Comparison of Survey and Non-Survey Regression effects")

Analysis

sub<-data %>%
  select(badhealth, opportunity_youth_cat_num, opportunity_youth_cat, race_eth,educ,whitemajority, otherminority,urban_rural, male,sampweight,male, strata) %>%
  filter( complete.cases( . ))


#cat<-sample(1:nrow(sub), size = 1000, replace = FALSE)

#sub<-sub[cat, ]

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

cat<-svyby(formula = ~badhealth,
           by = ~opportunity_youth_cat,
           design = des,
           FUN = svymean,
           na.rm=T)

svychisq(~badhealth+opportunity_youth_cat,
         design = des)

## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(~badhealth + opportunity_youth_cat, design = des)
## F = 25.042, ndf = 1, ddf = 3710, p-value = 5.87e-07

cat%>%
  ggplot()+
  geom_point(aes(x=opportunity_youth_cat,y=badhealth))+
  geom_errorbar(aes(x=opportunity_youth_cat, ymin = badhealth-1.96*se, 
                    ymax= badhealth+1.96*se),
                width=.25)+
   labs(title = "Percent % of Young Adults Ages 16-24 with Fair/Poor Health by Opportunity youth category", 
        caption = "Source: IPUMS Health Survey Data, 2019-2020 \n Calculations by Joseph Jaiyeola",
       x = "Opportunity youth category",
       y = "%  Fair/Poor Health")+
  theme_minimal()

dog<-svyby(formula = ~badhealth,
           by = ~male, 
           design = des, 
           FUN = svymean,
           na.rm=T)

svychisq(~badhealth+male,
         design = des)

## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(~badhealth + male, design = des)
## F = 3.4527, ndf = 1, ddf = 3710, p-value = 0.06323

dog%>%
  ggplot()+
  geom_point(aes(x=male,y=badhealth))+
  geom_errorbar(aes(x=male, ymin = badhealth-1.96*se, 
                    ymax= badhealth+1.96*se),
                width=.25)+
   labs(title = "Percent % of US  with Fair/Poor Health by Gender", 
        caption = "Source: IPUMS Health Survey Data, 2019-2020 \n Calculations by Joseph Jaiyeola",
       x = "Gender",
       y = "%  Fair/Poor Health")+
  theme_minimal()

catdog<-svyby(formula = ~badhealth,
              by = ~male+opportunity_youth_cat,
              design = des,
              FUN = svymean,
              na.rm=T)

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

catdog%>%
  ggplot()+
  #geom_point(aes(x=educ, y = badhealth, color=race_eth, group=race_eth), position="dodge")+ 
  geom_errorbar(aes(x=opportunity_youth_cat,y = badhealth,
                    ymin = badhealth-1.96*se, 
                   ymax= badhealth+1.96*se,
                   color=male,
                   group=male),
                width=.25,
                position="dodge")+
  #facet_wrap(~ race_eth, nrow = 3)+
  labs(title = "Percent % of US  with Fair/Poor Health by Gender and Oppor Youth Cat", 
        caption = "Source: IPUMS Health Survey Data, 2019-2020 \n Calculations by Joseph Jaiyeola",
       x = "Opportunity Youth Category",
       y = "%  Fair/Poor Health")+
  theme_minimal()

#Logit model
fit.logit<-svyglm(badhealth ~ opportunity_youth_cat + urban_rural + male  + educ, 
                  design = des,
                  family = binomial)

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

summary(fit.logit)

## 
## Call:
## svyglm(formula = badhealth ~ opportunity_youth_cat + urban_rural + 
##     male + educ, design = des, family = binomial)
## 
## Survey design:
## svydesign(ids = ~1, strata = ~strata, weights = ~sampweight, 
##     data = sub)
## 
## Coefficients:
##                                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                 -1.2283     0.3312  -3.709 0.000211
## opportunity_youth_catNot opportunity youth  -0.9007     0.2078  -4.334  1.5e-05
## urban_ruralurban                            -0.4207     0.2384  -1.765 0.077698
## maleMale                                    -0.3444     0.1763  -1.953 0.050848
## educ1somehs                                  0.1592     0.4232   0.376 0.706758
## educ2hsgrad                                 -0.5323     0.2878  -1.849 0.064469
## educ3somecol                                -0.4732     0.2911  -1.626 0.104073
## educ4colgrad                                -1.0339     0.3843  -2.690 0.007175
##                                               
## (Intercept)                                ***
## opportunity_youth_catNot opportunity youth ***
## urban_ruralurban                           .  
## maleMale                                   .  
## educ1somehs                                   
## educ2hsgrad                                .  
## educ3somecol                                  
## educ4colgrad                               ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1.001159)
## 
## Number of Fisher Scoring iterations: 6

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

Characteristic	OR1	95% CI1	p-value
opportunity_youth_cat
Opportunity youth	—	—
Not opportunity youth	0.41	0.27, 0.61	<0.001
urban_rural
rural	—	—
urban	0.66	0.41, 1.05	0.078
male
Female	—	—
Male	0.71	0.50, 1.00	0.051
educ
0Prim	—	—
1somehs	1.17	0.51, 2.69	0.7
2hsgrad	0.59	0.33, 1.03	0.064
3somecol	0.62	0.35, 1.10	0.10
4colgrad	0.36	0.17, 0.76	0.007
1 OR = Odds Ratio, CI = Confidence Interval

library(sjPlot)
plot_model(fit.logit,
           axis.lim = c(.1, 10), #you may need to modify these
           title = "Odds ratios for Poor Self Rated Health")

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

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

knitr::kable(marg_logit,  digits = 4)

opportunity_youth_cat	urban_rural	male	educ	prob	SE	df	asymp.LCL	asymp.UCL
Opportunity youth	rural	Female	0Prim	0.2265	0.0580	Inf	0.1327	0.3591
Not opportunity youth	rural	Female	0Prim	0.1063	0.0311	Inf	0.0590	0.1842
Opportunity youth	urban	Female	0Prim	0.1612	0.0431	Inf	0.0933	0.2642
Not opportunity youth	urban	Female	0Prim	0.0724	0.0193	Inf	0.0426	0.1206
Opportunity youth	rural	Male	0Prim	0.1718	0.0449	Inf	0.1005	0.2781
Not opportunity youth	rural	Male	0Prim	0.0777	0.0211	Inf	0.0452	0.1306
Opportunity youth	urban	Male	0Prim	0.1199	0.0335	Inf	0.0681	0.2024
Not opportunity youth	urban	Male	0Prim	0.0524	0.0133	Inf	0.0317	0.0855
Opportunity youth	rural	Female	1somehs	0.2556	0.0797	Inf	0.1312	0.4383
Not opportunity youth	rural	Female	1somehs	0.1224	0.0440	Inf	0.0588	0.2376
Opportunity youth	urban	Female	1somehs	0.1840	0.0591	Inf	0.0944	0.3278
Not opportunity youth	urban	Female	1somehs	0.0839	0.0278	Inf	0.0431	0.1570
Opportunity youth	rural	Male	1somehs	0.1957	0.0653	Inf	0.0973	0.3544
Not opportunity youth	rural	Male	1somehs	0.0900	0.0322	Inf	0.0437	0.1762
Opportunity youth	urban	Male	1somehs	0.1377	0.0476	Inf	0.0678	0.2596
Not opportunity youth	urban	Male	1somehs	0.0609	0.0204	Inf	0.0313	0.1154
Opportunity youth	rural	Female	2hsgrad	0.1467	0.0341	Inf	0.0916	0.2268
Not opportunity youth	rural	Female	2hsgrad	0.0653	0.0170	Inf	0.0389	0.1076
Opportunity youth	urban	Female	2hsgrad	0.1014	0.0198	Inf	0.0687	0.1474
Not opportunity youth	urban	Female	2hsgrad	0.0439	0.0078	Inf	0.0309	0.0618
Opportunity youth	rural	Male	2hsgrad	0.1086	0.0266	Inf	0.0664	0.1726
Not opportunity youth	rural	Male	2hsgrad	0.0472	0.0118	Inf	0.0288	0.0764
Opportunity youth	urban	Male	2hsgrad	0.0741	0.0164	Inf	0.0476	0.1135
Not opportunity youth	urban	Male	2hsgrad	0.0315	0.0057	Inf	0.0221	0.0447
Opportunity youth	rural	Female	3somecol	0.1543	0.0390	Inf	0.0921	0.2469
Not opportunity youth	rural	Female	3somecol	0.0690	0.0178	Inf	0.0413	0.1131
Opportunity youth	urban	Female	3somecol	0.1070	0.0216	Inf	0.0714	0.1571
Not opportunity youth	urban	Female	3somecol	0.0464	0.0065	Inf	0.0352	0.0610
Opportunity youth	rural	Male	3somecol	0.1145	0.0324	Inf	0.0646	0.1949
Not opportunity youth	rural	Male	3somecol	0.0499	0.0134	Inf	0.0293	0.0837
Opportunity youth	urban	Male	3somecol	0.0782	0.0195	Inf	0.0476	0.1260
Not opportunity youth	urban	Male	3somecol	0.0333	0.0059	Inf	0.0235	0.0471
Opportunity youth	rural	Female	4colgrad	0.0943	0.0328	Inf	0.0468	0.1809
Not opportunity youth	rural	Female	4colgrad	0.0406	0.0145	Inf	0.0200	0.0807
Opportunity youth	urban	Female	4colgrad	0.0640	0.0200	Inf	0.0344	0.1161
Not opportunity youth	urban	Female	4colgrad	0.0270	0.0077	Inf	0.0154	0.0470
Opportunity youth	rural	Male	4colgrad	0.0687	0.0255	Inf	0.0327	0.1386
Not opportunity youth	rural	Male	4colgrad	0.0291	0.0106	Inf	0.0142	0.0587
Opportunity youth	urban	Male	4colgrad	0.0462	0.0160	Inf	0.0233	0.0897
Not opportunity youth	urban	Male	4colgrad	0.0193	0.0059	Inf	0.0106	0.0349

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

comps<-as.data.frame(marg_logit)

comps[comps$opportunity_youth_cat=="Opportunity youth" & comps$educ == "4colgrad" & comps$urban_rural == "urban" , ]

comps[comps$opportunity_youth_cat=="Not opportunity youth" & comps$educ == "4colgrad" & comps$urban_rural == "urban" , ]

According to the model, holding staying in a urban environment and being a college grad constant, opportunity youths have higher probability of reporting fair/poor self health rating compared to connected youth. And this is this for both male and female