PROLOG

PROJECT: HDS Final Project

PURPOSE: Analyze relationship between social determinants of health, obgyn access, and maternal mortality rate

DIR: C:/Users/cadye/OneDrive/Desktop/R Class Folders/Final Project

DATA: obgyn_access.csv

AUTHOR: Catherine Bailey

CREATED: 4/24/24

LAST EDITED: 5/2/24

NOTES

PROLOG

#Load tidyverse package
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.5.1     ✔ 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(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(knitr)
# read obgyn access data
data <- read.csv("obgyn_access_data.csv")
# Summarize the data

summary(data)
##  State.Territory    M.D..Obstetrics.and.Gynecology..County.Level.File.
##  Length:51          Min.   :  56.0                                    
##  Class :character   1st Qu.: 214.8                                    
##  Mode  :character   Median : 582.0                                    
##                     Mean   : 799.0                                    
##                     3rd Qu.:1025.5                                    
##                     Max.   :4928.0                                    
##                     NA's   :1                                         
##  Population.All..County.Level.File. Rate..per.100000.population.
##  Min.   :  582328                   Min.   : 6.86               
##  1st Qu.: 1854573                   1st Qu.:10.07               
##  Median : 4561284                   Median :11.15               
##  Mean   : 6575426                   Mean   :11.78               
##  3rd Qu.: 7625559                   3rd Qu.:12.91               
##  Max.   :39368078                   Max.   :19.20               
##  NA's   :1                          NA's   :1                   
##  MMR..per.100000.live.births. Uninsured....women.age.15.44.
##  Min.   :10.10                Min.   : 3.000               
##  1st Qu.:18.40                1st Qu.: 7.912               
##  Median :25.40                Median :10.012               
##  Mean   :25.23                Mean   :11.015               
##  3rd Qu.:30.20                3rd Qu.:13.506               
##  Max.   :43.50                Max.   :24.550               
##  NA's   :6                    NA's   :7                    
##  Median.Household.Income High.School.Completion.. X..of.C.sections
##  Min.   :45928           Min.   :83.30            Min.   :22.77   
##  1st Qu.:57343           1st Qu.:87.05            1st Qu.:28.07   
##  Median :62028           Median :89.85            Median :30.85   
##  Mean   :64738           Mean   :89.20            Mean   :30.72   
##  3rd Qu.:74823           3rd Qu.:91.33            3rd Qu.:33.71   
##  Max.   :86644           Max.   :93.60            Max.   :38.38   
##  NA's   :7               NA's   :7                NA's   :7       
##  Late.No.Prenatal.Care..   black_mmr     
##  Min.   : 3.450          Min.   :  3.60  
##  1st Qu.: 4.844          1st Qu.: 30.85  
##  Median : 5.975          Median : 54.90  
##  Mean   : 6.253          Mean   : 54.96  
##  3rd Qu.: 7.088          3rd Qu.: 74.15  
##  Max.   :11.750          Max.   :122.90  
##  NA's   :7               NA's   :1

Since there are some NA’s present and the column labels are lengthy, the data will be recoded

# recode the data
data.clean <- data %>%
  rename(state = State.Territory) %>%
  rename(total_obgyn = M.D..Obstetrics.and.Gynecology..County.Level.File.) %>%
  rename(population = Population.All..County.Level.File.) %>%
  rename(obgyn_per_100K = Rate..per.100000.population.) %>%
  rename(mmr = MMR..per.100000.live.births.) %>%
  rename(unins = Uninsured....women.age.15.44.) %>%
  rename(income = Median.Household.Income) %>%
  rename(education = High.School.Completion..) %>%
  rename(c.sect = X..of.C.sections) %>%
  rename(prenatal = Late.No.Prenatal.Care..) %>%
  drop_na()
# Repeat data summary
summary(data.clean)
##     state            total_obgyn       population       obgyn_per_100K 
##  Length:44          Min.   :  77.0   Min.   :  731158   Min.   : 6.86  
##  Class :character   1st Qu.: 277.8   1st Qu.: 3014588   1st Qu.:10.03  
##  Mode  :character   Median : 673.5   Median : 5437691   Median :11.15  
##                     Mean   : 890.4   Mean   : 7336484   Mean   :11.60  
##                     3rd Qu.:1151.8   3rd Qu.: 8663515   3rd Qu.:12.53  
##                     Max.   :4928.0   Max.   :39368078   Max.   :19.20  
##       mmr            unins            income        education    
##  Min.   :10.10   Min.   : 3.000   Min.   :45928   Min.   :83.30  
##  1st Qu.:18.38   1st Qu.: 7.912   1st Qu.:57343   1st Qu.:87.05  
##  Median :25.55   Median :10.012   Median :62028   Median :89.85  
##  Mean   :25.27   Mean   :11.015   Mean   :64738   Mean   :89.20  
##  3rd Qu.:30.23   3rd Qu.:13.506   3rd Qu.:74823   3rd Qu.:91.33  
##  Max.   :43.50   Max.   :24.550   Max.   :86644   Max.   :93.60  
##      c.sect         prenatal        black_mmr     
##  Min.   :22.77   Min.   : 3.450   Min.   :  3.60  
##  1st Qu.:28.07   1st Qu.: 4.844   1st Qu.: 33.77  
##  Median :30.85   Median : 5.975   Median : 56.50  
##  Mean   :30.72   Mean   : 6.253   Mean   : 58.08  
##  3rd Qu.:33.71   3rd Qu.: 7.088   3rd Qu.: 74.92  
##  Max.   :38.38   Max.   :11.750   Max.   :122.90

The data is much more polished and managable now

# Plot a histogram of obgyn access level
ggplot(data = data.clean, aes(x = obgyn_per_100K)) +
  geom_histogram(bins = 15)

Interpretation: Number of OBGYNs per 100K population is not normally distributed

# Plot histogram of maternal mortality rate

data.clean %>%
  ggplot(aes(x = mmr)) +
  geom_histogram(bins = 15)

Interpretation: Maternal mortality rate is not normally distributed across states

# Create a scatter plot of OBGYNs per 100K population and maternal mortality rate with a linear fit line

data.clean %>%
  ggplot(aes(x = obgyn_per_100K, y = mmr)) +
  geom_point(aes(size = "state"), color = "blue") +
  geom_smooth(aes(linetype = "Linear fit line"), method = "lm", se = FALSE, color = "gray60") +
  theme_minimal() +
  labs(x = "Number of OBGYNs per 100K Population", y = "Maternal Mortality Rate") + 
  ggtitle("Maternal mortality rate vs. OBGYNs per 100K")
## Warning: Using size for a discrete variable is not advised.
## `geom_smooth()` using formula = 'y ~ x'

Interpretation: The is a weakly negative correlation between # of OBGYNs per 100K population and MMR.

# Create a linear model comparing MMR to OBGYN access

prd.by.ob <- lm(formula = mmr ~ obgyn_per_100K, data = data.clean, na.action = na.exclude)
# Get linear model summary
summary(prd.by.ob)
## 
## Call:
## lm(formula = mmr ~ obgyn_per_100K, data = data.clean, na.action = na.exclude)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.7271  -6.1656   0.2045   3.9359  16.9358 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     34.6681     5.9911   5.787 8.05e-07 ***
## obgyn_per_100K  -0.8098     0.5044  -1.605    0.116    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.497 on 42 degrees of freedom
## Multiple R-squared:  0.05781,    Adjusted R-squared:  0.03538 
## F-statistic: 2.577 on 1 and 42 DF,  p-value: 0.1159

The y intercept is 34.6681. This estimates that a state with 0 obgyns would have 34.6681 pregnancy related deaths per 100K births.

The slope is -0.8098. This means if the number of obgyns per 100K population in a state goes up by 1, there would be 0.8098 fewer pregnancy related deaths per 100K births.

# Get the confidence intervals
confint(object = prd.by.ob)
##                    2.5 %     97.5 %
## (Intercept)    22.577608 46.7585049
## obgyn_per_100K -1.827799  0.2082004

Interpretation: A simple linear regression analysis found that a state’s number of OBGYNs is not a statistically significant predictor for number of pregnancy related deaths (b = -0.8098, p > 0.05). The slope value is likely between -1.828 and 0.208. The model was not statistically better from the baseline (F(1, 42) = 2.577, p > 0.05). These results do not suggest a relationship between number obgyns and number of pregnancy related deaths in a state.

Checking assumptions

1: Independent Observations * Number of pregancy related deaths are not reliant on any other pregnancy related death. * This assumption is met

2: Continuous outcome * Since NAs were eliminated, # of OBGYNs and # of pregnancy related deaths can be any value greater than or equal to 0 * This assumption is met

3: Linearity

# Plot a scatter plot of obgyn access vs mmr with a linear fit line and loess curve
data.clean %>%
  ggplot(aes(x = obgyn_per_100K, y = mmr)) +
  geom_point() +
  geom_smooth(aes(color = "Linear fit line"), method = "lm", se = FALSE) +
  geom_smooth(aes(color = "Loess curve"), se = FALSE) +
  theme_minimal() +
  labs(y = "Maternal Mortality Rate", x = "OBGYNs per 100K population") +
  scale_color_manual(values = c("gray60", "deeppink"), name = "") +
  ggtitle("Maternal mortality rate vs. OBGYNs per 100K population")
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Interpretation: The Loess curve does not match the linear fit line at all. The assumption is failed.

4: Homeoscadacity

# Conduct a bp test for homeoscadacity
const.var.test <- lmtest::bptest(formula = prd.by.ob)

const.var.test
## 
##  studentized Breusch-Pagan test
## 
## data:  prd.by.ob
## BP = 0.4322, df = 1, p-value = 0.5109

The Breush-Pagan test has a high p-value. The null hypothesis is accepted, and the assumption of homoescedasticity is met.

Assumption 5: Independence of residuals

# Perforn a Durbin-Watson test for independence of residuals
lmtest::dwtest(formula = prd.by.ob)
## 
##  Durbin-Watson test
## 
## data:  prd.by.ob
## DW = 1.8048, p-value = 0.2533
## alternative hypothesis: true autocorrelation is greater than 0

The D-W statistic is near 2 and the p-value is > 0.05, so the null hypothesis is accepted. This assumption is met.

Assumption 6: Normality of residuals

# Plot a histogram of the residuals
data.frame(prd.by.ob$residuals) %>%
  ggplot(aes(x = prd.by.ob.residuals)) +
  labs(x = "Residuals") +
  geom_histogram(bins = 15) +
  ggtitle("Residuals of MMR vs. OBGYNs per 100K")

The residuals are not normally distributed. This assumption is not met.

Conclusion: A linear model is not an appropriate way to describe this data

Since a linear fit model is not an appropriate choice for analysis, covariance will be tried instead.

Check the strength of the Pearson’s Product Moment Correlation Coefficient

# get the Pearson's product moment correlation coefficient
data.clean %>% summarize(cov.ob.prd = cor(x = obgyn_per_100K,
                                          y = mmr,
                                          use = "complete"))
##   cov.ob.prd
## 1 -0.2404432

Interpretation: the strength of the r is weakly negative

Check assumptions (see above for tests)

1: Observations are independent: met 2: Both variables are continuous: met 3: Both variables are normally distributed: not met 4: Linearity: not met 5: Homoscedasticity: met

Because one of the the main issue is the nonlinaerity of the relationship, I will first try transforming the data

# apply log and arcsin transformations to the data

data.clean <- data.clean %>%
  mutate(log.ob = log(x = (obgyn_per_100K/100)/(1-obgyn_per_100K/100))) %>%
  mutate(log.prd = log(x = (mmr/100)/(1-mmr/100))) %>%
  mutate(arcsin.ob = asin(x = sqrt(obgyn_per_100K/100))) %>%
  mutate(arcsin.prd = asin(x = sqrt(mmr/100)))

Recheck normality

# Check normality of log transformed obgyn access using a box plot
data.clean %>%
  ggplot(aes(x = log.ob)) +
  geom_histogram(binwidth = 0.05)

Interpretation: a log transformation does not make the obgyn access data more normally distributed

# Check normality of log transformed mmr with a histogram

data.clean %>%
  ggplot(aes(x = log.prd)) +
  geom_histogram(bins = 15)

Interpretation: Log transformation was unsuccessful in normalizing mmr data

# Check normality of arcsin transformed obgyn access data with a histogram
data.clean %>%
  ggplot(aes(x = arcsin.ob)) +
  geom_histogram(bins = 15)

Interpretation: Arsin transformation was unsuccessful in normalizing obgyn access data

# Check the normality of arcsin transformed mmr data with a histogram

data.clean %>%
  ggplot(aes(x = arcsin.prd)) +
  geom_histogram(bins = 15)

Interpretation: Arsin transformation was unsuccessful in normalizing mmr data

Next, folded power transformation will be attempted

# Get lambda tranformation value for obgyn access
p.ob <- rcompanion::transformTukey(x = data.clean$obgyn_per_100K,
                                   plotit = FALSE,
                                   quiet = TRUE,
                                   returnLambda = TRUE)

p.ob
## lambda 
##  -0.25
# Get lambda tranformation value for mmr
p.prd <- rcompanion::transformTukey(x = data.clean$mmr,
                                   plotit = FALSE,
                                   quiet = TRUE,
                                   returnLambda = TRUE)

p.prd
## lambda 
##  0.275
# Apply folded power transformation to data

data.clean <- data.clean %>%
  mutate(folded.ob = (obgyn_per_100K/100)^(1/p.ob)-(1-obgyn_per_100K/100)^(1/p.ob)) %>%
  mutate(folded.prd = (mmr/100)^(1/p.prd)-(1-mmr/100)^(1/p.prd))
# Check normality of folded power transformed obgyn access data for normality
data.clean %>%
  ggplot(aes(x = folded.ob)) +
  geom_histogram(bins = 15)

Interpretation: Folded power transformation was unsuccessful at normalizing obgyn access data

# Check normality of folded power transformed mmr data using a histogram
data.clean %>%
  ggplot(aes(x = folded.prd)) +
  geom_histogram(bins = 15)

Interpretation: Folded power transformation was unsuccessful at normalizing mmr data

Conclusion: Transforming the data was not successful. It may be better to transform the data into categories.

Recoding data as categorical

OBGYN access data will be recoded as categorical, high OBGYN access, mid OBGYN access, and low OBGYN access. I will do this by getting the range of the data, and dividing by 3. This will provide the low, mid, and high cutoffs.

# Define the cutoff value for 
cutoff <- (max(data.clean$obgyn_per_100K) - min(data.clean$obgyn_per_100K))/3;cutoff
## [1] 4.113333
# Define the minimum values for each obgyn access category
low.ob <- min(data.clean$obgyn_per_100K);low.ob
## [1] 6.86
mid.ob <- low.ob + cutoff;mid.ob
## [1] 10.97333
high.ob <- mid.ob + cutoff;high.ob
## [1] 15.08667
# Add a column to the dataset for obgyn access category
data.clean <- data.clean %>%
  mutate(ob.category = 0)

# Population the newly created ob.category variable based on OBGYNs per 100K
for(i in 1:nrow(data.clean)){
  if(data.clean$obgyn_per_100K[i] < mid.ob){
    data.clean$ob.category[i] <-  "low"}
  else if (data.clean$obgyn_per_100K[i] >= mid.ob & data.clean$obgyn_per_100K[i] < high.ob){
    data.clean$ob.category[i] <- "mid"
  } else{
    data.clean$ob.category[i] <- "high"
  }
  }

The appropriate parametric test for this data is ANOVA. First, assumptions will be tested.

1: Independence of observations: met (see above) 2: Continuous variables and independent groups: met 3: Normal distribution in each group

# Perform Shapiro Wilk test for normality among groups
data.clean %>%
  group_by(ob.category) %>%
  summarize(shapiro.pval = shapiro.test(x = mmr)$p.value)
## # A tibble: 3 × 2
##   ob.category shapiro.pval
##   <chr>              <dbl>
## 1 high              0.0367
## 2 low               0.563 
## 3 mid               0.569

For high obgyn access states, mmr is not normally distributed (Shapiro-Wilk test, p < 0.05). This assumption is failed.

4: Equal variances for each group

# Perform a Levene test for homogeny of variance

leveneTest(mmr ~ ob.category, data.clean)
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.2444 0.7843
##       41

A Levene test confirms variance is equal between groups (p >0.05).

Since the data fails the normality assumption, we will use the Kruskal-Wallis rank sum test.

# Run KW test comparing mmr to obgyn access category

kw.ob.by.prd <- kruskal.test(formula = mmr ~ ob.category, data = data.clean);kw.ob.by.prd
## 
##  Kruskal-Wallis rank sum test
## 
## data:  mmr by ob.category
## Kruskal-Wallis chi-squared = 2.8224, df = 2, p-value = 0.2438

Based on a Kruskal-Wallis test, we accept the null hypothesis that there is no difference in maternal mortality between states with low, mid, and high obgyn availability (X2 = 2.8224, p > 0.05).

# order the obgyn access categories as factors
data.clean$ob.category <- factor(data.clean$ob.category, levels = c("low", "mid", "high"))
# Create a box plot of mmr in states with low, mid, and high OBGYN access
data.clean %>% 
  ggplot(aes(x = ob.category, y = mmr)) +
  geom_boxplot(aes(fill = ob.category)) +
  theme_minimal() +
  scale_fill_brewer(palette = "Spectral", guide = "none") +
  scale_color_brewer(palette = "Spectral", guide = "none") +
  labs(x = "Level of obgyn access", y = "Maternal mortality rate") +
  ggtitle("Maternal mortality rate v. OBGYN access level")

Next, social determininants of health will be compared to OBGYN access level. Assumptions for ANOVA will be considered first, then an appropriate test will be determined.

How does percent of women age 15-44 in a state effect obgyn access?

1: Independece of observations: met

2: Continuous variable and independent groups: Met; percent uninsured can range from 0-100 continuously. OBGYN access levels are independent of other states.

3: Assumption of normality between groups

# Perform a Shapiro-Wilk test 
data.clean %>%
  group_by(ob.category) %>%
  summarize(shapiro.pval = shapiro.test(x = unins)$p.value)
## # A tibble: 3 × 2
##   ob.category shapiro.pval
##   <fct>              <dbl>
## 1 low                0.758
## 2 mid                0.128
## 3 high               0.668

All categories have p > 0.05 for the Shapiro-Wilk test. The null hypothesis that data is normally distributed between groups is met.

4: Homogeny of variance

# Perform a Levene test for homogeny of variance
leveneTest(unins ~ ob.category, data.clean)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  1.3371 0.2738
##       41

With p > 0.05, the Levene Test confirms that the variance between groups is homogenous. The assumption is met.

Based on this data, I will proceed with a Kruskall-Wallis test again

# Create a boxplot of percent of women ages 15-44 uninsured grouped by obgyn access level
data.clean %>% 
  ggplot(aes(x = ob.category, y = unins)) +
  geom_boxplot(aes(fill = ob.category)) +
  scale_fill_brewer(palette = "Spectral", guide = "none") +
  scale_color_brewer(palette = "Spectral", guide = "none") +
  theme_minimal() +
  labs(x = "Level of obgyn access", y = "% Women 15-44 uninsured") +
  ggtitle("% Women age 15-44 Uninsured vs. OBGYN access level")

# 
kw.ob.by.unins <- kruskal.test(formula =  unins ~ ob.category, data = data.clean);kw.ob.by.unins
## 
##  Kruskal-Wallis rank sum test
## 
## data:  unins by ob.category
## Kruskal-Wallis chi-squared = 7.3577, df = 2, p-value = 0.02525

Based on a Kruskal-Wallis test, the null hypothesis is rejected in favor of the alternate hypothesis that a relationship exists between percent of the population that is uninsured and access to obgyn care (X2 = 7.3577, p < 0.05).

# Run a dunn post hoc test on the relationship between the % population of uninsured women and obgyn access level

dunn.test::dunn.test(x = data.clean$unins, g = data.clean$ob.category, method = "bonferroni")
##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 7.3577, df = 2, p-value = 0.03
## 
## 
##                            Comparison of x by group                            
##                                  (Bonferroni)                                  
## Col Mean-|
## Row Mean |       high        low
## ---------+----------------------
##      low |  -2.709367
##          |    0.0101*
##          |
##      mid |  -2.077258   0.969535
##          |     0.0567     0.4984
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2

Based on the Dunn test, there is a significant difference in percent of the population uninsured in states with high vs. low obgyn access (p < 0.05).

Effect size

etaSquared.unins.ob <- (7.3577-3+1)/(44+3);etaSquared.unins.ob
## [1] 0.1139936

The effect size is medium

A Kruskal-Wallis test found a statistically significant difference in percent of the population uninsured across states with different levels of obgyn access (X2 = 7.3577, p < 0.05). Based on a Dunn’s post hoc test, states with low obgyn access had a significantly higher percent of the population uninsured than states with high obgyn access (p < 0.05). There was a medium effect size for the relationship between obgyn access level and ranked values of population uninsured (n2 = 0.114).

Does OBGYN access differ by median income

Check assumptions

1: Independent observations: met; one household’s income does not affect another’s 2: Continuous variable and independent groups: met; median household income can be any positive value. OBGYN access level is not dependent on access level of other states.

3: Normality between groups

# Perform a Shapiro-Wilk test 
data.clean %>%
  group_by(ob.category) %>%
  summarize(shapiro.pval = shapiro.test(x = income)$p.value)
## # A tibble: 3 × 2
##   ob.category shapiro.pval
##   <fct>              <dbl>
## 1 low               0.341 
## 2 mid               0.0527
## 3 high              0.258

All groups have p > 0.05 for a shapiro-wilk test, so the null hypothesis that data is normally distributed between groups is met.

4: Homogeny of variance

leveneTest(income ~ ob.category, data.clean)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.6546  0.525
##       41

The Levene’s test p-value is > 0.05, so the null hypothesis that variance is homogenous is met.

ANOVA is the appropriate analysis for this data.

H0: There is no relationship between median household income and obgyn access level in a state. HA: A relationship exists between median household income and obgyn access level in a state.

# Perform ANOVA

aov.inc.ob <- aov(income ~ ob.category, data = data.clean)

summary(aov.inc.ob)
##             Df    Sum Sq   Mean Sq F value  Pr(>F)   
## ob.category  2 1.282e+09 641163427   6.583 0.00332 **
## Residuals   41 3.993e+09  97393894                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Based on an ANOVA test, the null hypothesis is rejected in favor of the alternate hypothesis, that a relationship exists between median household income and obgyn access level in a state (F = 6.583, p < 0.01).

# Make a box plot showing distribution of income grouped by obgyn access
data.clean %>% 
  ggplot(aes(x = ob.category, y = income)) +
  geom_boxplot(aes(fill = ob.category)) +
  scale_fill_brewer(palette = "Spectral", guide = "none") +
  scale_color_brewer(palette = "Spectral", guide = "none") +
  theme_minimal() +
  labs(x = "Level of obgyn access", y = "Median Income") +
  ggtitle("Median income vs. level of OBGYN access")

As median household income increases, so does obgyn access.

A Bonferroni pairwise test will be run as a post hoc test

# Run Bonferroni pairwise post hoc test. 

pairwise.t.test(x = data.clean$income,
                g = data.clean$ob.category,
                p.adj = "bonf")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  data.clean$income and data.clean$ob.category 
## 
##      low   mid  
## mid  0.025 -    
## high 0.011 0.589
## 
## P value adjustment method: bonferroni

Based on the Bonferroni pairwise comparison, states with low obgyn access have lower median household income than states with high obgyn access (p < 0.05) and states with mid obgyn access (p < 0.05).

Effect size

# Calculate omega squared effect size

omega.inc.ob <- (6.583-1)/(6.583 + ((44-3+1)/3-1));omega.inc.ob
## [1] 0.2850942

Based on the omega squared value, the effect size is large (w2 = 0.285).

A one-way ANOVA found a statistically significant difference in percent of the population uninsured across states with different levels of obgyn access (F = 6.583, p < 0.01). Based on a Bonferroni post hoc test, states with low obgyn access have significantly lower median income than states with high obgyn access (p < 0.05) and states with mid obgyn access (p < 0.05). There was a large effect size for the relationship between obgyn access level and ranked values of median income (w2 = .285).

Is there are realationship between percent of deliveries done by c-section and obgyn access

1: Independent observations: met; one person’s delivery does not impact another’s 2: Continuous variable and independent groups: met; percent of deliveries by c-section can be any value 0-100. OBGYN access level is not dependent on access level of other states.

3: Normality between groups

# Perform a Shapiro-Wilk test 
data.clean %>%
  group_by(ob.category) %>%
  summarize(shapiro.pval = shapiro.test(x = c.sect)$p.value)
## # A tibble: 3 × 2
##   ob.category shapiro.pval
##   <fct>              <dbl>
## 1 low                0.598
## 2 mid                0.267
## 3 high               0.182

All groups have p > 0.05 for a shapiro-wilk test, so the null hypothesis that data is normally distributed between groups is met.

4: Homogeny of variance

# Perform a Levene test for homogeneity of variance
leveneTest(c.sect ~ ob.category, data.clean)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  2.3739 0.1058
##       41

The Levene’s test p-value is > 0.05, so the null hypothesis that variance is homogenous is met.

ANOVA is the appropriate analysis for this data.

H0: There is no relationship between percent of births caried out by c-section and obgyn access level in a state. HA: A relationship exists between percent of births caried out by c-section and obgyn access level in a state.

# Perform ANOVA
aov.csect.ob <- aov(c.sect ~ ob.category, data = data.clean)

summary(aov.csect.ob)
##             Df Sum Sq Mean Sq F value Pr(>F)
## ob.category  2   37.3   18.65   1.345  0.272
## Residuals   41  568.4   13.86

Based on a one way ANOVA, the null hypothesis is accepted, that no relationship between percent of births caried out by c-section and obgyn access level in a state.

Is there a relationship between percent of population who has completed high school and OBGYN access?

1: Independent observations: met; one person’s education status is not dependent on another’s 2: Continuous variable and independent groups: met; percent of the population that completed high school can be any value 0-100. OBGYN access level is not dependent on access level of other states.

3: Normality between groups

# Perform a Shapiro-Wilk test 
data.clean %>%
  group_by(ob.category) %>%
  summarize(shapiro.pval = shapiro.test(x = education)$p.value)
## # A tibble: 3 × 2
##   ob.category shapiro.pval
##   <fct>              <dbl>
## 1 low               0.0793
## 2 mid               0.0683
## 3 high              0.504

All groups have p > 0.05 for a shapiro-wilk test, so the null hypothesis that data is normally distributed between groups is met.

4: Homogeny of variance

# Perform a Levene test for homogeneity of variance
leveneTest(education ~ ob.category, data.clean)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.6183 0.5438
##       41

The Levene’s test p-value is > 0.05, so the null hypothesis that variance is homogenous is met.

ANOVA is the appropriate analysis for this data.

H0: There is no relationship between percent of the population with a high school education and obgyn access level in a state. HA: A relationship exists between percent of the population with a high school education and obgyn access level in a state.

# Perform ANOVA
aov.edu.ob <- aov(education ~ ob.category, data = data.clean)

summary(aov.edu.ob)
##             Df Sum Sq Mean Sq F value Pr(>F)
## ob.category  2   0.84   0.421   0.055  0.946
## Residuals   41 311.96   7.609

Based on a one way ANOVA, the null hypothesis that no relationship exists between percent of the population with a high school education and obgyn access level in a state is accepted (F = 0.055, p > 0.05).

Is there a relationship between OBGYN access and percent of pregnancies with late/no prenatal care

1: Independent observations: met; one person’s access to or choice to obtain prenatal care is not dependent on another person 2: Continuous variable and independent groups: met; percent of pregnancies with late/no prenatal carel can be any value 0-100. OBGYN access level is not dependent on access level of other states.

3: Normality between groups

# Perform a Shapiro-Wilk test 
data.clean %>%
  group_by(ob.category) %>%
  summarize(shapiro.pval = shapiro.test(x = prenatal)$p.value)
## # A tibble: 3 × 2
##   ob.category shapiro.pval
##   <fct>              <dbl>
## 1 low               0.0914
## 2 mid               0.802 
## 3 high              0.611

All groups have p > 0.05 for a shapiro-wilk test, so the null hypothesis that data is normally distributed between groups is met.

4: Homogeny of variance

# Perform a Levene test for homogeneity of variance
leveneTest(prenatal~ ob.category, data.clean)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.6028 0.5521
##       41

The Levene’s test p-value is > 0.05, so the null hypothesis that variance is homogenous is met.

ANOVA is the appropriate analysis for this data.

H0: There is no relationship between percent of pregnancies with late/no prenatal care and obgyn access level in a state. HA: A relationship exists between percent of pregnancies with late/no prenatal care education and obgyn access level in a state.

# Perform ANOVA
aov.prenat.ob <- aov(prenatal ~ ob.category, data = data.clean)

summary(aov.prenat.ob)
##             Df Sum Sq Mean Sq F value Pr(>F)
## ob.category  2   3.34   1.668   0.399  0.673
## Residuals   41 171.30   4.178

Based on one-way analysis of variance, the null hypothesis is accepted, that there is no relationship between percent of pregnancies with late/no prenatal care and obgyn access level in a state (F = 0.399, p > 0.05).

Does OBGYN access impact Black maternal mortality rates?

Check assumptions for ANOVA

1: Independent observations: met; no case of mortality is dependent on another 2: Continuous variable and independent groups: met; MMR can be any positive value. OBGYN access level is not dependent on access level of other states.

3: Normality between groups

# Perform a Shapiro-Wilk test 
data.clean %>%
  group_by(ob.category) %>%
  summarize(shapiro.pval = shapiro.test(x = black_mmr)$p.value)
## # A tibble: 3 × 2
##   ob.category shapiro.pval
##   <fct>              <dbl>
## 1 low               0.906 
## 2 mid               0.0918
## 3 high              0.890

All groups have p > 0.05 for a shapiro-wilk test, so the null hypothesis that data is normally distributed between groups is met.

4: Homogeny of variance

# Perform a Levene test for homogeneity of variance
leveneTest(black_mmr ~ ob.category, data.clean)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2   0.091 0.9132
##       41

The Levene’s test p-value is > 0.05, so the null hypothesis that variance is homogenous is met.

ANOVA is the appropriate analysis for this data.

H0: There is no relationship between Black maternal mortality and obgyn access level in a state. HA: A relationship exists between Black maternal mortality and obgyn access level in a state.

# Perform ANOVA
aov.bmmr.ob <- aov(black_mmr ~ ob.category, data = data.clean)

summary(aov.bmmr.ob)
##             Df Sum Sq Mean Sq F value Pr(>F)
## ob.category  2    552   275.8   0.322  0.727
## Residuals   41  35166   857.7

Based on one-way analysis of variance, the null hypothesis is accepted, that there is no relationship between Black maternal mortality and obgyn access level in a state (F = 0.322, p > 0.05).