library(haven)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(janitor)
##
## Attaching package: 'janitor'
##
## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.test
library(sjPlot)
## #refugeeswelcome
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## some
library(ggplot2)
library(patchwork)
library(ggpubr)
library(gtsummary)
library(haven)
library(stats)
library(survey)
## Loading required package: grid
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## Loading required package: survival
##
## Attaching package: 'survey'
##
## The following object is masked from 'package:graphics':
##
## dotchart
library(marginaleffects)
library(readr)
usa_00043 <- read_csv("usa_00043.csv")
## Rows: 16545837 Columns: 28
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (28): YEAR, MULTYEAR, SAMPLE, SERIAL, CBSERIAL, HHWT, CLUSTER, STRATA, G...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
dat1 <- usa_00043
colnames(dat1) <- tolower(colnames(dat1))
dat1 <- na.omit(dat1)
dat1$educd <- as.numeric(dat1$educd)
dat1$degree <- Recode(
dat1$educd,
recodes = "000:061='No High School Degree'; 062:064='High School Degree'; 065:100='Some College'; 100:116='BA/BS Degree'"
)
dat1$degree <- factor(
dat1$degree,
levels = c('No High School Degree', 'High School Degree', 'Some College', 'BA/BS Degree')
)
degreetable <- tabyl(dat1$degree)
dat1$empstat <- as.numeric(dat1$empstat)
dat1$employmentstatus <- recode(dat1$empstat, "1='Employed'; 2='Unemployed'; 3='Not in labor force'; else=NA")
dat1$employmentstatus <- factor(dat1$employmentstatus, levels = c("Employed", "Unemployed", "Not in labor force"), exclude = NULL)
dat1$insured <- Recode(dat1$hcovany, recodes = "2=1;else=0")
tabyl(dat1$insured)
## dat1$insured n percent
## 0 1953809 0.1180846
## 1 14592028 0.8819154
blkamerica <- dat1 %>%
filter((raced == 200 | raced == 210) & (bpld == 99900 | bpld < 99900))
blkamerica$group <- "usablack"
blkafrican <- dat1 %>%
filter((raced == 200 | raced == 210) & bpld >= 60000 & bpld < 69999 & bpld != 99900)
blkafrican$group <- "africanblack"
blkcaribbean <- dat1 %>%
filter((raced == 200 | raced == 210) &
(bpld >= 21000 & bpld < 39999 |
bpld >= 9000 & bpld < 10000 |
bpld >= 200 & bpld < 300)
)
blkcaribbean$group <- "blackcaribbean"
blkcombined <- rbind(blkamerica, blkafrican, blkcaribbean)
tabyl(blkcombined$group)
## blkcombined$group n percent
## africanblack 83657 0.04511008
## blackcaribbean 114209 0.06158453
## usablack 1656642 0.89330539
ageprobs<- blkcombined %>%
group_by(age, group) %>%
summarize(p=mean(insured),n=n())
## `summarise()` has grouped output by 'age'. You can override using the `.groups`
## argument.
degreeprobs<- blkcombined %>%
group_by(degree, group) %>%
summarize(p=mean(insured),n=n())
## `summarise()` has grouped output by 'degree'. You can override using the
## `.groups` argument.
empstprobs<- blkcombined %>%
group_by(employmentstatus, group) %>%
summarize(p=mean(insured),n=n())
## `summarise()` has grouped output by 'employmentstatus'. You can override using
## the `.groups` argument.
summary_table <- blkcombined %>%
select(insured, degree, ,employmentstatus, age, group) %>%
tbl_summary(
by = group,
statistic = list(insured ~ "{p}% ({n})",
degree ~ "{p}% ({n})",
age ~ "{mean} ({sd})")
) %>%
modify_header(label ~ "**Variable**") %>%
modify_spanning_header(all_stat_cols() ~ "**Group**") %>%
modify_caption("**Statistics for Heath Insured, Employment Status, Educational Attainment , and Age of Native Blacks, And Nonnative Black African, Carribean/Latin American **")
summary_table
**Statistics for Heath Insured, Employment Status, Educational Attainment , and Age of Native Blacks, And Nonnative Black African, Carribean/Latin American **
| Variable |
Group
|
| africanblack, N = 83,657 |
blackcaribbean, N = 114,209 |
usablack, N = 1,656,642 |
| insured |
83% (69,506) |
84% (95,469) |
83% (1,380,076) |
| degree |
|
|
|
| Â Â Â Â No High School Degree |
11% (9,043) |
15% (16,911) |
14% (224,848) |
| Â Â Â Â High School Degree |
19% (16,070) |
29% (33,490) |
31% (518,815) |
| Â Â Â Â Some College |
28% (23,282) |
30% (34,139) |
33% (543,575) |
| Â Â Â Â BA/BS Degree |
42% (35,262) |
26% (29,669) |
22% (369,404) |
| employmentstatus |
|
|
|
| Â Â Â Â Employed |
64,922 (78%) |
85,873 (75%) |
1,025,561 (62%) |
| Â Â Â Â Unemployed |
4,345 (5.2%) |
6,244 (5.5%) |
95,568 (5.8%) |
| Â Â Â Â Not in labor force |
14,390 (17%) |
22,092 (19%) |
535,513 (32%) |
| age |
43 (11) |
47 (11) |
46 (12) |
model1 <- glm(insured ~ group, data = blkcombined, family = binomial)
summary(model1)
##
## Call:
## glm(formula = insured ~ group, family = binomial, data = blkcombined)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.591628 0.009222 172.588 < 2e-16 ***
## groupblackcaribbean 0.036514 0.012202 2.992 0.00277 **
## groupusablack 0.015817 0.009455 1.673 0.09434 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1672337 on 1854507 degrees of freedom
## Residual deviance: 1672327 on 1854505 degrees of freedom
## AIC: 1672333
##
## Number of Fisher Scoring iterations: 3
coefficients <- coef(model1)
mod1logodds <- exp(coefficients)
mod1logodds
## (Intercept) groupblackcaribbean groupusablack
## 4.911738 1.037188 1.015942
model2 <- glm(insured ~ age + degree + employmentstatus, data = blkcombined, family = binomial)
summary(model2)
##
## Call:
## glm(formula = insured ~ age + degree + employmentstatus, family = binomial,
## data = blkcombined)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.1532242 0.0093088 -16.46 <2e-16 ***
## age 0.0374505 0.0001762 212.50 <2e-16 ***
## degreeHigh School Degree 0.2227031 0.0057881 38.48 <2e-16 ***
## degreeSome College 0.6192147 0.0061396 100.86 <2e-16 ***
## degreeBA/BS Degree 1.2163981 0.0076884 158.21 <2e-16 ***
## employmentstatusUnemployed -1.1152404 0.0073965 -150.78 <2e-16 ***
## employmentstatusNot in labor force -0.7976759 0.0045826 -174.07 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1672337 on 1854507 degrees of freedom
## Residual deviance: 1529256 on 1854501 degrees of freedom
## AIC: 1529270
##
## Number of Fisher Scoring iterations: 5
coefficients <- coef(model2)
mod2logodds <- exp(coefficients)
mod2logodds
## (Intercept) age
## 0.8579373 1.0381606
## degreeHigh School Degree degreeSome College
## 1.2494496 1.8574687
## degreeBA/BS Degree employmentstatusUnemployed
## 3.3750094 0.3278365
## employmentstatusNot in labor force
## 0.4503745
ageprobs$num <-ageprobs$p*(1-ageprobs$p)
ageprobs$sep <- sqrt(ageprobs$num/degreeprobs$n)
## Warning in ageprobs$num/degreeprobs$n: longer object length is not a multiple
## of shorter object length
ageprobs$me <- 2*ageprobs$sep
ggplot(data =ageprobs, aes(x = age , y = p, ymin=p-me, ymax=p+me, fill=group,color = group, group = group)) +
geom_line() + ylim(0,1) + geom_ribbon(alpha=0.3,aes(color=NULL)) +
labs(title="Health Insured By Age",
subtitle="Native Black/Non-Native Black African/Carribean/Latin American Respondents",
caption="Source: 2017 & 2022 ACS 5-Year File") +
xlab(label="Age") +
ylab(label="Insured Rate")
degreeprobs$num <- degreeprobs$p*(1-degreeprobs$p)
degreeprobs$sep <- sqrt(degreeprobs$num/degreeprobs$n)
degreeprobs$me <- 2*degreeprobs$sep
ggplot(data =degreeprobs, aes(x = degree , y = p, ymin=p-me, ymax=p+me, fill=group,color = group, group = group)) +
geom_line() + ylim(0,1) + geom_ribbon(alpha=0.3,aes(color=NULL)) +
labs(title="Health Insured Rates By Education Level",
subtitle="Native Black/Non-Native Black African/Carribean/Latin American Respondents",
caption="Source: 2017 & 2022 ACS 5-Year File") +
xlab(label="Education Attainment") +
ylab(label="Insured Rate")
empstprobs$num <- empstprobs$p*(1-empstprobs$p)
empstprobs$sep <- sqrt(empstprobs$num/empstprobs$n)
empstprobs$me <- 2*empstprobs$sep
ggplot(data =empstprobs, aes(x = employmentstatus , y = p, ymin=p-me, ymax=p+me, fill=group,color = group, group = group)) +
geom_line() + ylim(0,1) + geom_ribbon(alpha=0.3,aes(color=NULL)) +
labs(title="Health Insured Rates by Employment Status",
subtitle="Native Black/Non-Native African/Carribean Respondents Ages 25 to 65",
caption="Source: 2017 & 2022 ACS 5-Year File") +
xlab(label="Education Attainment") +
ylab(label="Insured Rate")
p2 <- plot_predictions(model2, condition = c("employmentstatus", "degree")) +
geom_line() +
ylim(0, 1) +
geom_ribbon(alpha = 0.4, aes(color = NULL)) +
labs(
title = "Health Insured Rates by Employment Status and Educational Attainment",
subtitle = "Native Black/Non-Native African/Caribbean/Latin American Respondents Ages 25 to 65",
caption = "Source: 2017 & 2022 ACS 5-Year File",
xlab = "Age",
ylab = "Insured Rate"
)
p2
