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)
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$employed <- Recode(dat1$empstat, recodes="1='Employed'; 2='Unemployed'; 3='Not in labor force'; else=NA")
tabyl(dat1$employed)
##       dat1$employed        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 < 70000 & bpld != 99900)
blkafrican$group <- "africanblack"
blkcaribbean <- dat1 %>%
  filter(raced == 200 & 
         (bpld >= 25000 & bpld < 26000 |  
          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.051779531
##     blackcaribbean    2275 0.001458894
##           usablack 1476380 0.946761575
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.
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") 

blkcombined$degree <- factor(blkcombined$degree)
summary_table <- blkcombined %>%
  select(insured, degree, 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 Insured, Degree, and Age by of African and African Americans**")
summary_table
Statistics for Insured, Degree, and Age by of African and African Americans
Variable Group
africanblack, N = 80,7451 blackcaribbean, N = 2,2751 usablack, N = 1,476,3801
insured 84% (67,544) 78% (1,782) 87% (1,281,962)
degree


    No High School Degree 11% (8,525) 15% (333) 11% (165,019)
    High School Degree 19% (15,466) 30% (690) 30% (443,455)
    Some College 28% (22,500) 27% (621) 34% (505,471)
    BA/BS Degree 42% (34,254) 28% (631) 25% (362,435)
age 43 (11) 45 (11) 46 (12)
1 % (n); Mean (SD) 
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.245186   0.252472  -0.971   0.3315    
## age                       0.022172   0.004682   4.735 2.19e-06 ***
## degreeHigh School Degree  0.298027   0.148983   2.000   0.0455 *  
## degreeSome College        0.755754   0.159804   4.729 2.25e-06 ***
## degreeBA/BS Degree        1.001284   0.165043   6.067 1.30e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2378.3  on 2274  degrees of freedom
## Residual deviance: 2311.6  on 2270  degrees of freedom
## AIC: 2321.6
## 
## Number of Fisher Scoring iterations: 4

#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.347504   0.051770  -6.712 1.91e-11 ***
## 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: 1225392  on 1559399  degrees of freedom
## Residual deviance: 1224640  on 1559397  degrees of freedom
## AIC: 1224646
## 
## 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.4478217  0.0143653 -31.174  < 2e-16 ***
## groupblackcaribbean      -0.2663153  0.0531119  -5.014 5.32e-07 ***
## groupusablack             0.3310423  0.0101088  32.748  < 2e-16 ***
## age                       0.0316356  0.0002049 154.358  < 2e-16 ***
## degreeHigh School Degree  0.2908994  0.0071927  40.443  < 2e-16 ***
## degreeSome College        0.7245341  0.0073643  98.385  < 2e-16 ***
## degreeBA/BS Degree        1.3364915  0.0086336 154.800  < 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: 1225392  on 1559399  degrees of freedom
## Residual deviance: 1171446  on 1559393  degrees of freedom
## AIC: 1171460
## 
## 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.5804563  0.2569285  -2.259
## 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.0014803  0.0047715   0.310
## groupusablack:age                             0.0115897  0.0009412  12.314
## groupblackcaribbean:degreeHigh School Degree  0.2007380  0.1525477   1.316
## groupusablack:degreeHigh School Degree        0.2075933  0.0336060   6.177
## groupblackcaribbean:degreeSome College        0.3963718  0.1628997   2.433
## groupusablack:degreeSome College              0.3873479  0.0325046  11.917
## groupblackcaribbean:degreeBA/BS Degree        0.2270357  0.1679347   1.352
## groupusablack:degreeBA/BS Degree              0.6173315  0.0323312  19.094
##                                              Pr(>|z|)    
## (Intercept)                                  1.96e-12 ***
## groupblackcaribbean                            0.0239 *  
## groupusablack                                 < 2e-16 ***
## age                                           < 2e-16 ***
## degreeHigh School Degree                       0.0030 ** 
## degreeSome College                            < 2e-16 ***
## degreeBA/BS Degree                            < 2e-16 ***
## groupblackcaribbean:age                        0.7564    
## groupusablack:age                             < 2e-16 ***
## groupblackcaribbean:degreeHigh School Degree   0.1882    
## groupusablack:degreeHigh School Degree       6.52e-10 ***
## groupblackcaribbean:degreeSome College         0.0150 *  
## groupusablack:degreeSome College              < 2e-16 ***
## groupblackcaribbean:degreeBA/BS Degree         0.1764    
## 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: 1225392  on 1559399  degrees of freedom
## Residual deviance: 1170869  on 1559385  degrees of freedom
## AIC: 1170899
## 
## 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 Pr(>|z|)    
## (Intercept)                  1.463e+00  2.532e-02  57.787  < 2e-16 ***
## groupblackcaribbean         -5.921e-01  1.381e-01  -4.288 1.80e-05 ***
## groupusablack               -8.707e-02  2.621e-02  -3.322 0.000894 ***
## empstat                     -3.024e-01  1.183e-02 -25.549  < 2e-16 ***
## ftotinc                      7.642e-06  1.793e-07  42.634  < 2e-16 ***
## groupblackcaribbean:empstat  2.615e-01  6.288e-02   4.159 3.19e-05 ***
## groupusablack:empstat        2.839e-01  1.217e-02  23.324  < 2e-16 ***
## groupblackcaribbean:ftotinc -3.968e-07  1.069e-06  -0.371 0.710444    
## groupusablack:ftotinc        4.854e-07  1.869e-07   2.598 0.009382 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1225392  on 1559399  degrees of freedom
## Residual deviance: 1184564  on 1559391  degrees of freedom
## AIC: 1184582
## 
## 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.

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