library(tidyverse)
library(dplyr)
library(tidycensus)
library(tmap)
library(stringr)
library(mediation)
library(hdmed)causalinferenceproject
1. Load Libraries
2. Load Fertility Data
#load fertility data here
#2010 data is weird
#t2010new <- read_csv("t2010new")
#t2010new <- t2010new %>% select(FIPScode,n)
#2011 data is weird
#t2011new <- read_csv("t2011new")
#t2011new <- t2011new %>% select(FIPScode,n)
#2012
t2012new <- read_csv("t2012new")
t2012new <- t2012new %>% dplyr::select(FIPScode, n)
#2013
t2013new <- read_csv("t2013new")
t2013new <- t2013new %>% dplyr::select(FIPScode, n)3. Load Mortality Data
#08 three year (07,08,09)
ml083year <- read.delim("ml083year.txt")
ml083year <- ml083year %>% dplyr::select(County.Code, Age.Adjusted.Rate) %>% dplyr::rename(threeyear08 = Age.Adjusted.Rate, FIPScode = County.Code)
ml083year$FIPScode <- str_pad(ml083year$FIPScode,5,side = c("left"),pad=0)
ml083year$threeyear08 <- ml083year$threeyear08 %>% as.numeric()
#09 three year (08,09,10)
ml093year <- read.delim("ml093year.txt")
ml093year <- ml093year %>% dplyr::select(County.Code, Age.Adjusted.Rate) %>% dplyr::rename(threeyear09 = Age.Adjusted.Rate,FIPScode=County.Code)
ml093year$FIPScode <- str_pad(ml093year$FIPScode,5,side = c("left"),pad=0)
ml093year$threeyear09 <- ml093year$threeyear09 %>% as.numeric()
#10 three year (09, 10, 11)
ml103year <- read.delim("ml103year.txt")
ml103year <- ml103year %>% dplyr::select(County.Code, Age.Adjusted.Rate) %>% dplyr::rename(threeyear10 = Age.Adjusted.Rate,FIPScode=County.Code)
ml103year$FIPScode <- str_pad(ml103year$FIPScode,5,side = c("left"),pad=0)
ml103year$threeyear10 <- ml103year$threeyear10 %>% as.numeric()
#11 three year (10,11,12)
ml113year <- read.delim("ml113year.txt")
ml113year <- ml113year %>% dplyr::select(County.Code, Age.Adjusted.Rate) %>% dplyr::rename(threeyear11 = Age.Adjusted.Rate,FIPScode=County.Code)
ml113year$FIPScode <- str_pad(ml113year$FIPScode,5,side = c("left"),pad=0)
ml113year$threeyear11 <- ml113year$threeyear11 %>% as.numeric()Load SVI Data
SVI_2010_US_county <- read_csv("SVI_2010_US_county.csv")
SVI_2010_US_county <- SVI_2010_US_county %>% dplyr::select(FIPS, E_P_POV, E_P_UNEMP, E_P_PCI, S_PL_THEME1, S_PL_THEME2, S_PL_THEME3, S_PL_THEME4, S_PL_THEMES, PL_SNGPRNT, PL_AGE17, PL_AGE65, P_MINORITY, PL_MINORITY, R_PL_THEME1) %>% dplyr::rename(propbelowpov = E_P_POV, propunemployed = E_P_UNEMP, percapinc = E_P_PCI, sumofsestheme = S_PL_THEME1, sumofhouseholdcomptheme = S_PL_THEME2, sumofminoritytheme = S_PL_THEME3, sumofhousingtransporttheme = S_PL_THEME4, sumofallthemes = S_PL_THEMES, FIPScode = FIPS, singleparenthh = PL_SNGPRNT, under17 = PL_AGE17, over65 = PL_AGE65, minority=P_MINORITY, pl_minority=PL_MINORITY, percentileSES=R_PL_THEME1)4. get ACS women data
#2012
nevermarried12 <- get_acs(
geography = "county",
variables = c('15-17'= "B12002_097", '18-19'="B12002_098", '20-24'="B12002_099", '25-29'= "B12002_100", '30-34'="B12002_101", '35-39'= "B12002_102", '40-44'="B12002_103", '45-49'="B12002_104"),
year = 2012)
totalwomen12 <- aggregate(nevermarried12$estimate, by=list(nevermarried12$GEOID), FUN=sum) %>% dplyr::rename(FIPScode = Group.1, numberofwomen = x)
t2012new <- left_join(t2012new, totalwomen12, by="FIPScode")
t2012new$nmfr <- (t2012new$n/t2012new$numberofwomen)*1000
#2013
nevermarried13 <- get_acs(
geography = "county",
variables = c('15-17'= "B12002_097", '18-19'="B12002_098", '20-24'="B12002_099", '25-29'= "B12002_100", '30-34'="B12002_101", '35-39'= "B12002_102", '40-44'="B12002_103", '45-49'="B12002_104"),
year = 2013)
totalwomen13 <- aggregate(nevermarried13$estimate, by=list(nevermarried13$GEOID), FUN=sum) %>% dplyr::rename(FIPScode = Group.1, numberofwomen = x)
t2013new <- left_join(t2013new, totalwomen13, by="FIPScode")
t2013new$nmfr <- (t2013new$n/t2013new$numberofwomen)*10005. get geometry shapefiles
geography <- get_acs(
geography = "county",
variables = c('15-17'= "B12002_097"),
year = 2012,
geometry = T) %>% dplyr::rename(FIPScode = GEOID)
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 41%
|
|============================= | 41%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%t2012new <- left_join(t2012new, geography, by="FIPScode")6. merge mortality and svi with fertility
#fertility 2012 join
fert2012 <- left_join(t2012new, ml083year, by="FIPScode") %>% left_join(., ml093year, by="FIPScode") %>% left_join(.,ml103year, by="FIPScode") %>% left_join(.,ml113year, by="FIPScode") %>% left_join(., SVI_2010_US_county, by="FIPScode")
fert2012 <- fert2012 %>% dplyr::select(FIPScode, nmfr, threeyear08, threeyear09, threeyear10, threeyear11, propbelowpov, propunemployed, percapinc, sumofsestheme, sumofhouseholdcomptheme, sumofminoritytheme, sumofhousingtransporttheme,sumofallthemes, minority, percentileSES)
fert2012 <- fert2012 %>% na.omit() %>% filter(sumofallthemes > -1)
#fertility 2013 join
fert2013 <- left_join(t2013new, ml083year, by="FIPScode") %>% left_join(., ml093year, by="FIPScode") %>% left_join(.,ml103year, by="FIPScode") %>% left_join(.,ml113year, by="FIPScode") %>% left_join(., SVI_2010_US_county, by="FIPScode")
fert2013 <- fert2013 %>% dplyr::select(FIPScode, nmfr, threeyear08, threeyear09, threeyear10, threeyear11, propbelowpov, propunemployed, percapinc, sumofsestheme, sumofhouseholdcomptheme, sumofminoritytheme, sumofhousingtransporttheme,sumofallthemes, singleparenthh, under17, over65, minority, percentileSES)
fert2013 <- fert2013 %>% na.omit() %>% filter(sumofallthemes > -1)7. bivariate plots for 2012 fertility
ggplot(data=fert2012, aes(x=threeyear08, y=nmfr)) +
geom_point() +
geom_smooth(method=lm)+
ggtitle("Midlife Mortality and Nonmarital Fertility - 2012 (Five Year Lag)")+
labs(caption = "2012 NM Ferility and 2008 Centered Mortality Rates \n Sources: ACS and CDC WONDER")+
xlab("Midlife Mortality Rate")+
ylab("Nonmarital Fertility rate")
ggplot(data=fert2012, aes(x=threeyear09, y=nmfr)) +
geom_point() +
geom_smooth(method=lm)+
ggtitle("Midlife Mortality and Nonmarital Fertility - 2012 (Four Year Lag)")+
labs(caption = "2012 NM Ferility and 2009 Centered Mortality Rates \n Sources: ACS and CDC WONDER")+
xlab("Midlife Mortality Rate")+
ylab("Nonmarital Fertility rate")
ggplot(data=fert2012, aes(x=threeyear10, y=nmfr)) +
geom_point() +
geom_smooth(method=lm)+
ggtitle("Midlife Mortality and Nonmarital Fertility - 2012 (Three Year Lag)")+
labs(caption = "2012 NM Ferility and 2010 Centered Mortality Rates \n Sources: ACS and CDC WONDER")+
xlab("Midlife Mortality Rate")+
ylab("Nonmarital Fertility rate")
ggplot(data=fert2012, aes(x=threeyear11, y=nmfr)) +
geom_point() +
geom_smooth(method=lm)+
ggtitle("Midlife Mortality and Nonmarital Fertility - 2012 (Two Year Lag)")+
labs(caption = "2012 NM Ferility and 2011 Centered Mortality Rates \n Sources: ACS and CDC WONDER")+
xlab("Midlife Mortality Rate")+
ylab("Nonmarital Fertility rate")
8. Bivariate plots for 2013
ggplot(data=fert2013, aes(x=threeyear08, y=nmfr)) +
geom_point() +
geom_smooth(method=lm)+
ggtitle("Midlife Mortality and Nonmarital Fertility - 2013 (Five Year Lag)")+
labs(caption = "2013 NM Ferility and 2008 Centered Mortality Rates \n Sources: ACS and CDC WONDER")+
xlab("Midlife Mortality Rate")+
ylab("Nonmarital Fertility rate")
ggplot(data=fert2013, aes(x=threeyear09, y=nmfr)) +
geom_point() +
geom_smooth(method=lm)+
ggtitle("Midlife Mortality and Nonmarital Fertility - 2013 (Four Year Lag)")+
labs(caption = "2013 NM Ferility and 2009 Centered Mortality Rates \n Sources: ACS and CDC WONDER")+
xlab("Midlife Mortality Rate")+
ylab("Nonmarital Fertility rate")
ggplot(data=fert2013, aes(x=threeyear10, y=nmfr)) +
geom_point() +
geom_smooth(method=lm)+
ggtitle("Midlife Mortality and Nonmarital Fertility - 2013 (Three Year Lag)")+
labs(caption = "2013 NM Ferility and 2010 Centered Mortality Rates \n Sources: ACS and CDC WONDER")+
xlab("Midlife Mortality Rate")+
ylab("Nonmarital Fertility rate")
ggplot(data=fert2013, aes(x=threeyear11, y=nmfr)) +
geom_point() +
geom_smooth(method=lm)+
ggtitle("Midlife Mortality and Nonmarital Fertility - 2013 (Two Year Lag)")+
labs(caption = "2013 NM Ferility and 2011 Centered Mortality Rates \n Sources: ACS and CDC WONDER")+
xlab("Midlife Mortality Rate")+
ylab("Nonmarital Fertility rate")
ggplot(data=fert2013, aes(x= sumofsestheme, y=nmfr)) +
geom_point() +
geom_smooth(method=lm)+
ggtitle("SVI SES Theme and Nonmarital Fertility - 2013")+
labs(caption = "2013 NM Ferility and 2010 SVI SES \n Sources: ACS and CDC WONDER")+
xlab("SVI SES")+
ylab("Nonmarital Fertility rate")
img <- ggplot(data=fert2013, aes(x= sumofhouseholdcomptheme, y=nmfr)) +
geom_point() +
geom_smooth(method=lm)+
ggtitle("SVI Household Composition Theme and Nonmarital Fertility - 2013")+
labs(caption = "2013 NM Ferility and 2010 SVI HH Comp \n Sources: ACS and CDC WONDER")+
xlab("SVI Household Composition")+
ylab("Nonmarital Fertility rate")
ggplot(data=fert2013, aes(x= fert2013$sumofminoritytheme, y=nmfr)) +
geom_point() +
geom_smooth(method=lm)+
ggtitle("SVI Minority Theme and Nonmarital Fertility - 2013 (Five Year Lag)")+
labs(caption = "2013 NM Ferility and 2010 SVI Minority \n Sources: ACS and CDC WONDER")+
xlab("SVI Minority")+
ylab("Nonmarital Fertility rate")
ggplot(data=fert2013, aes(x= fert2013$sumofhousingtransporttheme, y=nmfr)) +
geom_point() +
geom_smooth(method=lm)+
ggtitle("SVI Housing Theme and Nonmarital Fertility - 2013 (Five Year Lag)")+
labs(caption = "2013 NM Ferility and 2010 SVI Housing \n Sources: ACS and CDC WONDER")+
xlab("SVI Housing")+
ylab("Nonmarital Fertility rate")
ggplot(data=fert2013, aes(x= sumofallthemes, y=nmfr)) +
geom_point() +
geom_smooth(method=lm)+
ggtitle("SVI and Nonmarital Fertility - 2013")+
labs(caption = "2013 NM Ferility and 2010 SVI \n Sources: ACS and CDC WONDER")+
xlab("SVI")+
ylab("Nonmarital Fertility rate")
ggplot(data=fert2013, aes(x= fert2013$propbelowpov, y=nmfr)) +
geom_point() +
geom_smooth(method=lm)+
ggtitle("Proportion below Poverty and Nonmarital Fertility - 2013 (Five Year Lag)")+
labs(caption = "2013 NM Ferility and 2010 Poverty \n Sources: ACS and CDC WONDER")+
xlab("Proportion below Poverty Line")+
ylab("Nonmarital Fertility rate")
ggplot(data=fert2013, aes(x= fert2013$propunemployed, y=nmfr)) +
geom_point() +
geom_smooth(method=lm)+
ggtitle("Proportion Unemployed and Nonmarital Fertility - 2013 (Five Year Lag)")+
labs(caption = "2013 NM Ferility and 2010 Unemployment \n Sources: ACS and CDC WONDER")+
xlab("Proportion Unemployment")+
ylab("Nonmarital Fertility rate")
ggplot(data=fert2013, aes(x= fert2013$percapinc, y=nmfr)) +
geom_point() +
geom_smooth(method=lm)+
ggtitle("Per Capita Income and Nonmarital Fertility - 2013 (Five Year Lag)")+
labs(caption = "2013 NM Ferility and 2010 Per Capita Income \n Sources: ACS and CDC WONDER")+
xlab("Per Capita Income")+
ylab("Nonmarital Fertility rate")
ggplot(data=fert2013, aes(x= singleparenthh, y=nmfr)) +
geom_point() +
geom_smooth(method=lm)+
ggtitle("SVI Single Parent HH and Nonmarital Fertility - 2013")+
labs(caption = "2013 NM Ferility and 2010 SVI Single Parent HH \n Sources: ACS and CDC WONDER")+
xlab("SVI SES")+
ylab("Nonmarital Fertility rate")
ggplot(data=fert2013, aes(x= under17, y=nmfr)) +
geom_point() +
geom_smooth(method=lm)+
ggtitle("SVI Under 17 and Nonmarital Fertility - 2013")+
labs(caption = "2013 NM Ferility and 2010 SVI Under 17 \n Sources: ACS and CDC WONDER")+
xlab("SVI U17")+
ylab("Nonmarital Fertility rate")
ggplot(data=fert2013, aes(x= over65, y=nmfr)) +
geom_point() +
geom_smooth(method=lm)+
ggtitle("SVI Over 65 and Nonmarital Fertility - 2013")+
labs(caption = "2013 NM Ferility and 2010 SVI Over 65 \n Sources: ACS and CDC WONDER")+
xlab("SVI O 65")+
ylab("Nonmarital Fertility rate")
8. simple models
#causal mediation model for 2008
med.simple08 <- lm(threeyear08 ~ sumofallthemes, data=fert2013)
out.simple08 <- lm(nmfr~sumofallthemes + threeyear08, data=fert2013)
med.out.simple08 <- mediate(med.simple08, out.simple08, treat = "sumofallthemes", mediator = "threeyear08", boot = T, sims = 1000, boot.ci.type = "bca")Running nonparametric bootstrapsummary(med.out.simple08)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the BCa Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 2.733 2.322 3.18 <2e-16 ***
ADE 2.703 2.172 3.22 <2e-16 ***
Total Effect 5.435 5.052 5.86 <2e-16 ***
Prop. Mediated 0.503 0.424 0.59 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 2899
Simulations: 1000 summary(med.simple08)
Call:
lm(formula = threeyear08 ~ sumofallthemes, data = fert2013)
Residuals:
Min 1Q Median 3Q Max
-281.52 -58.84 -4.54 50.64 883.11
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 105.0340 6.8392 15.36 <2e-16 ***
sumofallthemes 40.2476 0.9293 43.31 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 92.13 on 2897 degrees of freedom
Multiple R-squared: 0.393, Adjusted R-squared: 0.3928
F-statistic: 1876 on 1 and 2897 DF, p-value: < 2.2e-16summary(out.simple08)
Call:
lm(formula = nmfr ~ sumofallthemes + threeyear08, data = fert2013)
Residuals:
Min 1Q Median 3Q Max
-59.727 -11.166 -1.010 8.381 132.123
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.724450 1.363991 9.329 <2e-16 ***
sumofallthemes 2.702661 0.228759 11.814 <2e-16 ***
threeyear08 0.067896 0.003563 19.055 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 17.67 on 2896 degrees of freedom
Multiple R-squared: 0.3087, Adjusted R-squared: 0.3082
F-statistic: 646.6 on 2 and 2896 DF, p-value: < 2.2e-16mod1 <- lm(nmfr~scale(sumofallthemes), data = fert2013)
summary(mod1)
Call:
lm(formula = nmfr ~ scale(sumofallthemes), data = fert2013)
Residuals:
Min 1Q Median 3Q Max
-54.050 -11.845 -0.519 9.632 132.780
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 58.5842 0.3481 168.32 <2e-16 ***
scale(sumofallthemes) 10.0093 0.3481 28.75 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 18.74 on 2897 degrees of freedom
Multiple R-squared: 0.222, Adjusted R-squared: 0.2217
F-statistic: 826.7 on 1 and 2897 DF, p-value: < 2.2e-16mod2 <- lm(threeyear08~scale(sumofallthemes), data = fert2013)
summary(mod2)
Call:
lm(formula = threeyear08 ~ scale(sumofallthemes), data = fert2013)
Residuals:
Min 1Q Median 3Q Max
-281.52 -58.84 -4.54 50.64 883.11
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 391.812 1.711 228.99 <2e-16 ***
scale(sumofallthemes) 74.118 1.711 43.31 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 92.13 on 2897 degrees of freedom
Multiple R-squared: 0.393, Adjusted R-squared: 0.3928
F-statistic: 1876 on 1 and 2897 DF, p-value: < 2.2e-16mod3 <- lm(nmfr~scale(sumofallthemes)+scale(threeyear08), data = fert2013)
summary(mod3)
Call:
lm(formula = nmfr ~ scale(sumofallthemes) + scale(threeyear08),
data = fert2013)
Residuals:
Min 1Q Median 3Q Max
-59.727 -11.166 -1.010 8.381 132.123
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 58.5842 0.3282 178.53 <2e-16 ***
scale(sumofallthemes) 4.9771 0.4213 11.81 <2e-16 ***
scale(threeyear08) 8.0273 0.4213 19.05 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 17.67 on 2896 degrees of freedom
Multiple R-squared: 0.3087, Adjusted R-squared: 0.3082
F-statistic: 646.6 on 2 and 2896 DF, p-value: < 2.2e-16#2009
med.simple09 <- lm(threeyear09 ~ sumofallthemes, data=fert2013)
out.simple09 <- lm(nmfr~sumofallthemes + threeyear09, data=fert2013)
med.out.simple09 <- mediate(med.simple09, out.simple09, treat = "sumofallthemes", mediator = "threeyear09", boot = T, sims = 1000, boot.ci.type = "bca")Running nonparametric bootstrapsummary(med.out.simple09)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the BCa Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 2.687 2.338 3.06 <2e-16 ***
ADE 2.749 2.242 3.27 <2e-16 ***
Total Effect 5.435 5.042 5.84 <2e-16 ***
Prop. Mediated 0.494 0.422 0.57 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 2899
Simulations: 1000 summary(med.simple09)
Call:
lm(formula = threeyear09 ~ sumofallthemes, data = fert2013)
Residuals:
Min 1Q Median 3Q Max
-294.83 -61.93 -5.90 50.80 689.68
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 112.5730 7.0105 16.06 <2e-16 ***
sumofallthemes 39.2974 0.9526 41.25 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 94.44 on 2897 degrees of freedom
Multiple R-squared: 0.3701, Adjusted R-squared: 0.3698
F-statistic: 1702 on 1 and 2897 DF, p-value: < 2.2e-16summary(out.simple09)
Call:
lm(formula = nmfr ~ sumofallthemes + threeyear09, data = fert2013)
Residuals:
Min 1Q Median 3Q Max
-53.335 -10.893 -0.956 8.227 131.653
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.159329 1.363148 8.92 <2e-16 ***
sumofallthemes 2.748588 0.223632 12.29 <2e-16 ***
threeyear09 0.068369 0.003462 19.75 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 17.6 on 2896 degrees of freedom
Multiple R-squared: 0.3144, Adjusted R-squared: 0.3139
F-statistic: 663.9 on 2 and 2896 DF, p-value: < 2.2e-16med.simple10 <- lm(threeyear10 ~ sumofallthemes, data=fert2013)
out.simple10 <- lm(nmfr~sumofallthemes + threeyear10, data=fert2013)
med.out.simple10 <- mediate(med.simple10, out.simple10, treat = "sumofallthemes", mediator = "threeyear10", boot = T, sims = 1000, boot.ci.type = "bca")Running nonparametric bootstrapsummary(med.out.simple10)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the BCa Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 2.718 2.365 3.13 <2e-16 ***
ADE 2.717 2.229 3.21 <2e-16 ***
Total Effect 5.435 5.033 5.84 <2e-16 ***
Prop. Mediated 0.500 0.433 0.58 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 2899
Simulations: 1000 summary(med.simple10)
Call:
lm(formula = threeyear10 ~ sumofallthemes, data = fert2013)
Residuals:
Min 1Q Median 3Q Max
-302.65 -62.15 -7.79 51.49 599.34
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 113.3371 7.2052 15.73 <2e-16 ***
sumofallthemes 39.3939 0.9791 40.24 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 97.06 on 2897 degrees of freedom
Multiple R-squared: 0.3585, Adjusted R-squared: 0.3583
F-statistic: 1619 on 1 and 2897 DF, p-value: < 2.2e-16summary(out.simple10)
Call:
lm(formula = nmfr ~ sumofallthemes + threeyear10, data = fert2013)
Residuals:
Min 1Q Median 3Q Max
-52.153 -10.797 -1.038 8.391 129.125
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.035591 1.353919 8.889 <2e-16 ***
sumofallthemes 2.717139 0.220474 12.324 <2e-16 ***
threeyear10 0.069000 0.003351 20.591 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 17.51 on 2896 degrees of freedom
Multiple R-squared: 0.3214, Adjusted R-squared: 0.3209
F-statistic: 685.7 on 2 and 2896 DF, p-value: < 2.2e-16med.simple11 <- lm(threeyear11 ~ sumofallthemes, data=fert2013)
out.simple11 <- lm(nmfr~sumofallthemes + threeyear11, data=fert2013)
med.out.simple11 <- mediate(med.simple11, out.simple11, treat = "sumofallthemes", mediator = "threeyear11", boot = T, sims = 5000, boot.ci.type = "bca")Running nonparametric bootstrapsummary(med.out.simple11)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the BCa Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 2.618 2.264 2.99 <2e-16 ***
ADE 2.818 2.311 3.29 <2e-16 ***
Total Effect 5.435 5.054 5.82 <2e-16 ***
Prop. Mediated 0.482 0.416 0.56 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 2899
Simulations: 5000 summary(med.simple11)
Call:
lm(formula = threeyear11 ~ sumofallthemes, data = fert2013)
Residuals:
Min 1Q Median 3Q Max
-314.31 -65.53 -9.20 53.89 656.93
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 120.5189 7.3493 16.40 <2e-16 ***
sumofallthemes 38.5243 0.9986 38.58 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 99 on 2897 degrees of freedom
Multiple R-squared: 0.3394, Adjusted R-squared: 0.3391
F-statistic: 1488 on 1 and 2897 DF, p-value: < 2.2e-16summary(out.simple11)
Call:
lm(formula = nmfr ~ sumofallthemes + threeyear11, data = fert2013)
Residuals:
Min 1Q Median 3Q Max
-57.308 -10.793 -0.970 8.332 126.752
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.666474 1.357644 8.593 <2e-16 ***
sumofallthemes 2.817558 0.217115 12.977 <2e-16 ***
threeyear11 0.067951 0.003283 20.697 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 17.49 on 2896 degrees of freedom
Multiple R-squared: 0.3223, Adjusted R-squared: 0.3218
F-statistic: 688.5 on 2 and 2896 DF, p-value: < 2.2e-16med.simple11SES <- lm(threeyear11 ~ sumofsestheme, data=fert2013)
out.simple11SES <- lm(nmfr~sumofsestheme + threeyear11, data=fert2013)
med.out.simple11SES <- mediate(med.simple11SES, out.simple11SES, treat = "sumofsestheme", mediator = "threeyear11", boot = T, sims = 5000, boot.ci.type = "bca")Running nonparametric bootstrapsummary(med.out.simple11SES)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the BCa Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 6.257 5.229 7.33 <2e-16 ***
ADE 4.367 3.202 5.62 <2e-16 ***
Total Effect 10.624 9.923 11.33 <2e-16 ***
Prop. Mediated 0.589 0.488 0.69 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 2899
Simulations: 5000 summary(med.simple11SES)
Call:
lm(formula = threeyear11 ~ sumofsestheme, data = fert2013)
Residuals:
Min 1Q Median 3Q Max
-347.80 -51.79 -4.49 44.88 613.07
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 205.529 3.656 56.22 <2e-16 ***
sumofsestheme 92.420 1.615 57.23 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 83.45 on 2897 degrees of freedom
Multiple R-squared: 0.5306, Adjusted R-squared: 0.5304
F-statistic: 3275 on 1 and 2897 DF, p-value: < 2.2e-16summary(out.simple11SES)
Call:
lm(formula = nmfr ~ sumofsestheme + threeyear11, data = fert2013)
Residuals:
Min 1Q Median 3Q Max
-57.192 -10.672 -1.127 8.756 128.756
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.888249 1.125430 20.337 <2e-16 ***
sumofsestheme 4.367328 0.501820 8.703 <2e-16 ***
threeyear11 0.067697 0.003955 17.116 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 17.77 on 2896 degrees of freedom
Multiple R-squared: 0.3011, Adjusted R-squared: 0.3006
F-statistic: 623.9 on 2 and 2896 DF, p-value: < 2.2e-169. modeling with controls
#upload Area Health Resource Files (2010, 2011 and 2012)
AHRF<-read_csv("AHRFhospop.csv")Rows: 3201 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): hospitalsper1000
dbl (4): FIPS, Year, Total Number Hospitals, Census Population
ℹ 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.#merge with dataframe for analysis
datacontrols<-merge(fert2013, AHRF, by.x = "FIPScode", by.y = "FIPS")
datacontrols$hospitals<-datacontrols$`Total Number Hospitals`
datacontrols$hospitalsper1000<-as.numeric(datacontrols$hospitalsper1000)
datacontrols$hospitalsper1000scaled<-scale(datacontrols$hospitalsper1000)med11SES <- lm(threeyear11 ~ percentileSES + minority + hospitalsper1000scaled, data=datacontrols)
out11SES <- lm(nmfr~percentileSES + threeyear11 + minority + hospitalsper1000scaled, data=datacontrols)
med.out.11SES <- mediate(med11SES, out11SES, treat = "percentileSES", mediator = "threeyear11", boot = T, sims = 5000, boot.ci.type = "bca")Running nonparametric bootstrapsummary(med.out.11SES)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the BCa Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 19.862 16.684 22.97 <2e-16 ***
ADE 11.247 7.435 15.24 <2e-16 ***
Total Effect 31.109 28.774 33.41 <2e-16 ***
Prop. Mediated 0.638 0.530 0.75 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 2622
Simulations: 5000 summary(med11SES)
Call:
lm(formula = threeyear11 ~ percentileSES + minority + hospitalsper1000scaled,
data = datacontrols)
Residuals:
Min 1Q Median 3Q Max
-247.94 -50.95 -4.07 44.06 481.87
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 239.885 3.338 71.857 < 2e-16 ***
percentileSES 308.507 6.118 50.427 < 2e-16 ***
minority -17.188 9.083 -1.892 0.0586 .
hospitalsper1000scaled 12.199 1.603 7.610 3.8e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 81.99 on 2618 degrees of freedom
Multiple R-squared: 0.53, Adjusted R-squared: 0.5295
F-statistic: 984.1 on 3 and 2618 DF, p-value: < 2.2e-16summary(out11SES)
Call:
lm(formula = nmfr ~ percentileSES + threeyear11 + minority +
hospitalsper1000scaled, data = datacontrols)
Residuals:
Min 1Q Median 3Q Max
-59.463 -10.141 -0.293 8.777 116.071
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 24.284946 1.193890 20.341 < 2e-16 ***
percentileSES 11.246655 1.781818 6.312 3.23e-10 ***
threeyear11 0.064382 0.004054 15.881 < 2e-16 ***
minority 17.802014 1.885359 9.442 < 2e-16 ***
hospitalsper1000scaled 4.480203 0.336182 13.327 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 17.01 on 2617 degrees of freedom
Multiple R-squared: 0.3593, Adjusted R-squared: 0.3583
F-statistic: 366.8 on 4 and 2617 DF, p-value: < 2.2e-16#sensitivity analysis
s.out<-medsens(med.out.11SES)
summary(s.out)
Mediation Sensitivity Analysis for Average Causal Mediation Effect
Sensitivity Region
Rho ACME 95% CI Lower 95% CI Upper R^2_M*R^2_Y* R^2_M~R^2_Y~
[1,] 0.3 -0.2594 -2.7089 2.1901 0.09 0.0271
Rho at which ACME = 0: 0.3
R^2_M*R^2_Y* at which ACME = 0: 0.09
R^2_M~R^2_Y~ at which ACME = 0: 0.0271 plot(s.out, "rho")
plot(s.out, "R2")
#outcome model
outcome<-lm(nmfr~percentileSES + minority + hospitalsper1000scaled, data=datacontrols)
summary(outcome)
Call:
lm(formula = nmfr ~ percentileSES + minority + hospitalsper1000scaled,
data = datacontrols)
Residuals:
Min 1Q Median 3Q Max
-72.264 -10.455 0.078 9.436 121.933
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 39.7293 0.7250 54.802 <2e-16 ***
percentileSES 31.1091 1.3286 23.416 <2e-16 ***
minority 16.6954 1.9724 8.465 <2e-16 ***
hospitalsper1000scaled 5.2656 0.3481 15.126 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 17.81 on 2618 degrees of freedom
Multiple R-squared: 0.2975, Adjusted R-squared: 0.2967
F-statistic: 369.6 on 3 and 2618 DF, p-value: < 2.2e-16#descriptive table
library(gtsummary)
Attaching package: 'gtsummary'The following object is masked from 'package:MASS':
selecttable1 <-
datacontrols |>
tbl_summary(
include = c(nmfr, percentileSES, threeyear11, minority, hospitalsper1000),
by = NULL,
label = list(
nmfr = "Nonmarital fertility rate (2013) (per 1,000)",
percentileSES = "Socioeconomic SVI theme (2010) (sum score percentile)",
threeyear11 = "Midlife mortality rate (2010-2012) (per 100,00)",
minority = "Proportion of county population is a minority racial/ethnic group (2010)",
hospitalsper1000 = "Hospitals per county (2010) (per 1,000)"
),
statistic = list(
all_continuous() ~ "{mean} ({min}, {max})"
)
) %>%
modify_header(label = "**Variable**") %>%
modify_caption("**Sample distribution of dependent and independent variables at the county-level**")
table1| Variable | N = 2,6221 |
|---|---|
| Nonmarital fertility rate (2013) (per 1,000) | 59 (2, 195) |
| Socioeconomic SVI theme (2010) (sum score percentile) | 0.51 (0.00, 1.00) |
| Midlife mortality rate (2010-2012) (per 100,00) | 392 (151, 989) |
| Proportion of county population is a minority racial/ethnic group (2010) | 0.21 (0.01, 0.97) |
| Hospitals per county (2010) (per 1,000) | 0.05 (0.00, 0.60) |
| 1 Mean (Min, Max) | |