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)
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
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<-filter(dat1, ftotinc>=0)
dat1 <- dat1 %>%
filter(ftotinc != 9999998 & ftotinc != 9999999)
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, recodes="1='Employed'; 2='Unemployed'; 3='Not in labor force'; else=NA")
dat1$empstat <- factor(dat1$empstat, levels = c(1, 2, 3), labels = c("Employed", "Unemployed", "Not in labor force"))
tabyl(dat1$employmentstatus)
## dat1$employmentstatus n percent
## Employed 11813954 0.7413898
## Not in labor force 3567047 0.2238516
## Unemployed 553874 0.0347586
dat1$insured <- Recode(dat1$hcovany, recodes = "2=1;else=0")
tabyl(dat1$insured)
## dat1$insured n percent
## 0 1715794 0.1076754
## 1 14219081 0.8923246
blkamerica <- dat1 %>%
filter(raced == 200 & (bpld == 99900 | bpld < 99900))
blkamerica$group <- "usablack"
blkafrican <- dat1 %>%
filter(raced == 200 & bpld >= 60000 & bpld < 69999 & bpld != 99900)
blkafrican$group <- "africanblack"
blkcaribbean <- dat1 %>%
filter(raced == 200 &
(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 80745 0.04843516
## blackcaribbean 109949 0.06595328
## usablack 1476380 0.88561156
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.
degree_summary <- blkamerica %>%
group_by(degree, insured) %>%
summarize(count = n()) %>%
group_by(degree) %>%
mutate(percent = count / sum(count) * 100)
## `summarise()` has grouped output by 'degree'. You can override using the
## `.groups` argument.
ggplot(degree_summary, aes(x = degree, y = percent, fill = factor(insured))) +
geom_bar(position = "dodge", stat = "identity") +
labs(title = "Health Insurance of Black Americans by Degree",
x = "Degree",
y = "Percentage",
fill = "Health Insurance") +
scale_fill_manual(values = c("blue", "red"), labels = c("Insured", "Uninsured")) +
theme_minimal()
degree_summary <- blkafrican %>%
group_by(degree, insured) %>%
summarize(count = n()) %>%
group_by(degree) %>%
mutate(percent = count / sum(count) * 100)
## `summarise()` has grouped output by 'degree'. You can override using the
## `.groups` argument.
ggplot(degree_summary, aes(x = degree, y = percent, fill = factor(insured))) +
geom_bar(position = "dodge", stat = "identity") +
labs(title = "Health Insurance of Black Africans by Degree",
x = "Degree",
y = "Percentage",
fill = "Health Insurance") +
scale_fill_manual(values = c("blue", "red"), labels = c("Insured", "Uninsured")) +
theme_minimal()
degree_summary <- blkcaribbean %>%
group_by(degree, insured) %>%
summarize(count = n()) %>%
group_by(degree) %>%
mutate(percent = count / sum(count) * 100)
## `summarise()` has grouped output by 'degree'. You can override using the
## `.groups` argument.
ggplot(degree_summary, aes(x = degree, y = percent, fill = factor(insured))) +
geom_bar(position = "dodge", stat = "identity") +
labs(title = "Health Insurance of Black Carribeans/South/Central Americans by Degree",
x = "Degree",
y = "Percentage",
fill = "Health Insurance") +
scale_fill_manual(values = c("blue", "red"), labels = c("Insured", "Uninsured")) +
theme_minimal()
#These models were built by dividing those identifying as the race black that were born in the United States from non-native blacks born in Africa. Black Carribean is the combined total of those identifying as black that were born in the Carribean, as well as South and Central America. There is a limitiation because my sample size for Black Carribean is so low.
model1 <- glm(insured ~ age + degree, data = blkamerica, family = binomial)
summary(model1)
##
## Call:
## glm(formula = insured ~ age + degree, family = binomial, data = blkamerica)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.1647912 0.0113747 -14.49 <2e-16 ***
## age 0.0322817 0.0002105 153.38 <2e-16 ***
## degreeHigh School Degree 0.3048819 0.0073790 41.32 <2e-16 ***
## degreeSome College 0.7467304 0.0075793 98.52 <2e-16 ***
## degreeBA/BS Degree 1.3915801 0.0090882 153.12 <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: 1150332 on 1476379 degrees of freedom
## Residual deviance: 1098130 on 1476375 degrees of freedom
## AIC: 1098140
##
## Number of Fisher Scoring iterations: 5
#Model 1 logistic regression analysis shows that there is a positive statistically significant when accounting for increaseing in age and educational attainment and having health insurance coverage for native blacks.
model2 <- glm(insured ~ age + degree, data = blkafrican, family = binomial)
summary(model2)
##
## Call:
## glm(formula = insured ~ age + degree, family = binomial, data = blkafrican)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.3352707 0.0476435 7.037 1.96e-12 ***
## age 0.0206921 0.0009174 22.556 < 2e-16 ***
## degreeHigh School Degree 0.0972886 0.0327858 2.967 0.003 **
## degreeSome College 0.3593825 0.0316086 11.370 < 2e-16 ***
## degreeBA/BS Degree 0.7742486 0.0310274 24.954 < 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: 71930 on 80744 degrees of freedom
## Residual deviance: 70427 on 80740 degrees of freedom
## AIC: 70437
##
## Number of Fisher Scoring iterations: 4
#Model 2 logistic regression analysis shows that there is a positive statistically significant when accounting for increaseing in age and educational attainment and having health insurance coverage for non-natice Africa blacks but some degree has a 36% higher chance than those with just high school.
model3 <- glm(insured ~ age + degree, data = blkcaribbean, family = binomial)
summary(model3)
##
## Call:
## glm(formula = insured ~ age + degree, family = binomial, data = blkcaribbean)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3671837 0.0421544 -8.71 <2e-16 ***
## age 0.0299463 0.0007712 38.83 <2e-16 ***
## degreeHigh School Degree 0.3970372 0.0237947 16.69 <2e-16 ***
## degreeSome College 0.8567849 0.0251054 34.13 <2e-16 ***
## degreeBA/BS Degree 1.4359595 0.0288202 49.83 <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: 94754 on 109948 degrees of freedom
## Residual deviance: 90570 on 109944 degrees of freedom
## AIC: 90580
##
## Number of Fisher Scoring iterations: 5
#Model 3 logistic regression analysis shows that there is a positive statistically significant when accounting for increaseing in age and educational attainment and having health insurance coverage for non-natice Carribean blacks but some degree has a 77% higher chance than those with just high school.
model4 <- glm(insured ~ group, data = blkcombined, family = binomial)
summary(model4)
##
## Call:
## glm(formula = insured ~ group, family = binomial, data = blkcombined)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.632487 0.009516 171.549 < 2e-16 ***
## groupblackcaribbean 0.065169 0.012652 5.151 2.6e-07 ***
## groupusablack 0.253650 0.009822 25.823 < 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: 1318039 on 1667073 degrees of freedom
## Residual deviance: 1317015 on 1667071 degrees of freedom
## AIC: 1317021
##
## Number of Fisher Scoring iterations: 4
#Model 4 is the combined Black native and non=native health insurance and Black Caribbean decreases the chancces of being insured compared to the African Black group and usa black group. These is no significant statistically measures between the three groups.
model5 <- glm(insured ~ group + age + degree , data = blkcombined, family = binomial)
summary(model5)
##
## Call:
## glm(formula = insured ~ group + age + degree, family = binomial,
## data = blkcombined)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.4517736 0.0141039 -32.032 < 2e-16 ***
## groupblackcaribbean 0.0967518 0.0129786 7.455 9.01e-14 ***
## groupusablack 0.3312829 0.0101060 32.781 < 2e-16 ***
## age 0.0315234 0.0001982 159.047 < 2e-16 ***
## degreeHigh School Degree 0.3003563 0.0068879 43.606 < 2e-16 ***
## degreeSome College 0.7356751 0.0070654 104.124 < 2e-16 ***
## degreeBA/BS Degree 1.3463650 0.0082769 162.665 < 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: 1318039 on 1667073 degrees of freedom
## Residual deviance: 1259729 on 1667067 degrees of freedom
## AIC: 1259743
##
## Number of Fisher Scoring iterations: 5
model6 <- glm(insured ~ group * (age + degree) , data = blkcombined, family = binomial)
summary(model6)
##
## Call:
## glm(formula = insured ~ group * (age + degree), family = binomial,
## data = blkcombined)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) 0.3352707 0.0476438 7.037
## groupblackcaribbean -0.7024544 0.0636154 -11.042
## groupusablack -0.5000619 0.0489828 -10.209
## age 0.0206921 0.0009174 22.556
## degreeHigh School Degree 0.0972886 0.0327859 2.967
## degreeSome College 0.3593825 0.0316086 11.370
## degreeBA/BS Degree 0.7742486 0.0310276 24.954
## groupblackcaribbean:age 0.0092543 0.0011985 7.722
## groupusablack:age 0.0115897 0.0009412 12.314
## groupblackcaribbean:degreeHigh School Degree 0.2997485 0.0405105 7.399
## groupusablack:degreeHigh School Degree 0.2075933 0.0336060 6.177
## groupblackcaribbean:degreeSome College 0.4974024 0.0403657 12.322
## groupusablack:degreeSome College 0.3873479 0.0325046 11.917
## groupblackcaribbean:degreeBA/BS Degree 0.6617109 0.0423475 15.626
## groupusablack:degreeBA/BS Degree 0.6173315 0.0323312 19.094
## Pr(>|z|)
## (Intercept) 1.96e-12 ***
## groupblackcaribbean < 2e-16 ***
## groupusablack < 2e-16 ***
## age < 2e-16 ***
## degreeHigh School Degree 0.003 **
## degreeSome College < 2e-16 ***
## degreeBA/BS Degree < 2e-16 ***
## groupblackcaribbean:age 1.15e-14 ***
## groupusablack:age < 2e-16 ***
## groupblackcaribbean:degreeHigh School Degree 1.37e-13 ***
## groupusablack:degreeHigh School Degree 6.52e-10 ***
## groupblackcaribbean:degreeSome College < 2e-16 ***
## groupusablack:degreeSome College < 2e-16 ***
## groupblackcaribbean:degreeBA/BS Degree < 2e-16 ***
## groupusablack:degreeBA/BS Degree < 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: 1318039 on 1667073 degrees of freedom
## Residual deviance: 1259127 on 1667059 degrees of freedom
## AIC: 1259157
##
## Number of Fisher Scoring iterations: 5
#Model 6 is similar to the previous models
model7 <- glm(insured ~ group * (empstat + ftotinc) , data = blkcombined, family = binomial)
summary(model7)
##
## Call:
## glm(formula = insured ~ group * (empstat + ftotinc), family = binomial,
## data = blkcombined)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) 1.201e+00 1.760e-02 68.237
## groupblackcaribbean 4.185e-02 2.388e-02 1.752
## groupusablack 2.728e-01 1.824e-02 14.954
## empstatUnemployed -8.779e-01 3.602e-02 -24.371
## empstatNot in labor force -5.466e-01 2.431e-02 -22.485
## ftotinc 7.485e-06 1.788e-07 41.854
## groupblackcaribbean:empstatUnemployed -1.867e-01 4.673e-02 -3.997
## groupusablack:empstatUnemployed -2.491e-01 3.689e-02 -6.751
## groupblackcaribbean:empstatNot in labor force 1.177e-01 3.221e-02 3.654
## groupusablack:empstatNot in labor force 5.842e-01 2.504e-02 23.333
## groupblackcaribbean:ftotinc 1.399e-07 2.395e-07 0.584
## groupusablack:ftotinc 9.030e-08 1.863e-07 0.485
## Pr(>|z|)
## (Intercept) < 2e-16 ***
## groupblackcaribbean 0.079696 .
## groupusablack < 2e-16 ***
## empstatUnemployed < 2e-16 ***
## empstatNot in labor force < 2e-16 ***
## ftotinc < 2e-16 ***
## groupblackcaribbean:empstatUnemployed 6.42e-05 ***
## groupusablack:empstatUnemployed 1.47e-11 ***
## groupblackcaribbean:empstatNot in labor force 0.000258 ***
## groupusablack:empstatNot in labor force < 2e-16 ***
## groupblackcaribbean:ftotinc 0.559098
## groupusablack:ftotinc 0.627919
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1318039 on 1667073 degrees of freedom
## Residual deviance: 1251655 on 1667062 degrees of freedom
## AIC: 1251679
##
## Number of Fisher Scoring iterations: 5
#Model 7, being in the Black Caribbean group decreases the chances of being insured compared to the African Black group, while being in the USA Black group increases the chances.
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 Age and Education Level",
subtitle="African American/African Respondents Ages 25 to 65",
caption="Source: 2017 & 2022 ACS 5-Year File") +
xlab(label="Age at Survey") +
ylab(label="Insured Rate")
p1<-plot_predictions(model5, condition = c("age","group"))
p1
#The smaller sample size for Carribean blacks has effected the results
but is consitant with models showing an increase in age increases
insurance chances in all groups.
p2<-plot_predictions(model6, condition = c("degree","group"))
p2
#USA blacks show much higher increase in coverage with educationtional
attainment