---

title: "SNF Statistical Analysis"

author: "Amira Burns"

date: "`r Sys.Date()`"

output:

html_document:

df_print: paged

word_document:

keep_md: yes

pdf_document:

self_contained: no

---

```{r setup, include=FALSE}

knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE)

# -- Load Libraries

library(readxl) # read in Excel files

library(car) # Type III Sums of Squares Anova

library(emmeans) # multiple comparisons

library(tidyverse) # data processing, visualization

library(lubridate) # process dates

library(ggsci) # colors

library(rstatix)

library(gridExtra) # make and arrange tables

library(gtable)

library(kableExtra)

library(nnet) # multinomial model fitting

library(lme4) # multilevel model fitting

```

## Introduction

Rehospitalization of skilled nursing facility (SNF) residents remains a pertinent area of interest. Current financial penalties for rehospitalization of SNF patients are based in part on the studies by Ouslander et al., 2011 and Mor et al., 2010, demonstrating that many SNF hospitalizations were avoidable. With increasing age and complex illness severity and increasing use of SNFs for subacute rehabilitation/short-term stay, readmission metrics and financial penalties based on previous data may be due for reevaluation.

### Overview of Analysis Plan

- Load and format the data

- Exploratory tables and visuals, demographic disaggregation.

- Check underlying assumptions regarding the data to determine what analysis tools to use.

- Analysis of variance to determine demographic association with scores.

- Explore post-hoc group contrasts for Anovas.

- Define and fit general linear model (GLM) to assess readmission risk.

+ To determine current viability of calculated scores, use factors that contribute to those calculated scores and evaluate consistency between contribution to scores and contribution to GLM.

### Load and review the data

```{r load}

snf_data <- read_xlsx("/data/SNF_CODED_CALCULATED_SCORES_04042022.xlsx",

sheet = "CODED DATA AND CALCULATIONS")

```

**Select and Format Variables**

Convert dates, factors to correct their correct format. Remove duplicated fields. Combine two variables for liver disease into single ordered variable; combine two variables for diabetes into single ordered variable. Clean values for `Race` and `InsurancePlanName`. DOB fields are not formatted correctly, year keeps defaulting to 20- instead of 19-; use patient age instead. Remove years 2010 - 2013 due to concerns about missing data repopulated as 0's instead of NAs.

```{r format_vars}

snf_data <- snf_data %>%

transmute(SternClientID = `Stern Client Id`,

MRN = MRN,

Year = WorkSheet,

CurrentAdmit = as.Date(`Current Admit`, format = "%m/%d/%y"),

DischargeDate = as.Date(`Discharge Date`, format = "%m/%d/%y"),

Gender = factor(Gender),

Race = Race,

Ethnicity = factor(Ethnicity),

InsurancePlanName = (InsurancePlanName),

ED6moPrior = factor(`LACE Num_of_ED_visits_prior 6months`),

PatientAge = PatientAge,

AdmissionHospital = factor(AdmissionHospital),

CCIAgeScore = ordered(`CCI AGE SCORE`),

LACELOSDaysScore = ordered(`LACE LOSDays SCORE`),

LOSDays = LOSDays,

LOS5Days = factor(`LOS>=5days`),

HOSPITALLOS5Days = factor(`HOSPITAL LOS>=5days`),

VisitType = factor(VisitType),

AdmitDtm = as.Date(AdmitDtm, format = "%m/%d/%y"),

DischargeDtm = as.Date(DischargeDtm, format = "%m/%d/%y"),

NumAdmitPastYear = `NumberOfAdmissionInPastYear`,

PriorAdmin = ordered(`HOSPITAL # OF admission prior yr`),

Cancer = factor(`HOSPITAL SCORE D/C ONC SERVICE`),

AcuteAdmin = factor(`LACE ACUTE ADMISSION`),

AdminType = factor(`HOSPITAL SCORE INDEX ADMISSION TYPE`),

Hemoglobin = `Hemoglobin at discharge`,

HemoglobinLevel = factor(`HOSPITAL SCORE Hemoglobin at discharge`),

Sodium = `Sodium, Serum at discharge`,

SodiumLevel = factor(`HOSPITAL SCORE Sodium`),

LACECCI = `LACE CCI SCORE`,

MyoInfarc = factor(`Myocardial infarction`),

CongHrtFail = factor(`Congestive heart failure`),

PeriphVascDis = factor(`Peripheral vascular disease`),

CerebroVascDis = factor(`Cerebrovascular disease`),

Dementia = factor(Dementia),

ChronPulmoDis = factor(`Chronic pulmonary disease`),

RheuConTisDis = factor(`Rheumatic disease=connective tissue disease`),

PepUlcDis = factor(`Peptic ulcer disease`),

MildLiverDis = factor(`MILD liver disease`), # combine into one liver disease var

ModSevLiverDis = factor(`MODERATE or SEVERE liver disease`),

DiabetesWOChrCom = factor(`Diabetes WITHOUT chronic complication`), # combine into one diabetes var

DiabetesWChrCom = factor(`Diabetes WITH chronic complication`),

HemiParaPlegia = factor(`Hemiplegia or paraplegia`),

RenalDisCKD = factor(`Renal disease = mod to severe CKD`),

Malignancy = factor(`Any Malignancy`),

MetaSolidTumor = factor(`Metastatic solid tumour`),

AidsHIV = factor(`AIDS/HIV`),

CCI = `CALCULATED CCI SCORE`,

LACE = `CALCULATED LACE SCORE`,

HOSPITAL = `CALCULATED HOSPITAL SCORE`,

Readmit30 = `READMISSION <30DAYS`) %>%

mutate(StayDurationDays = as.numeric((DischargeDate - CurrentAdmit)+1),

DtmDuration = DischargeDtm - AdmitDtm)

```

```{r standardize}

# Exclude records before 2014 due to reporting changes

snf_data <- snf_data %>% filter(Year > 2013)

# Create new column for whether or not record is a readmission

snf_data <- snf_data %>%

mutate(Readmit30 = ifelse(Readmit30 == "N/A", "O", Readmit30)) %>%

rename(Readmit = Readmit30)

# Standardize values of Race, Insurance, Liver Disease, Diabetes

snf_data <- snf_data %>%

mutate(Race = as.factor(ifelse(Race == "White", "Caucasian/White",

ifelse((Race == "Declined" | Race == "Not Specified" | Race == "Unknown"), "Unavailable/Unknown",

ifelse(Race == "Multiracial", "Other/Multiracial", Race)))),

InsurancePlanCleaned = ifelse(grepl("MEDICARE", InsurancePlanName), "MCARE_MCAID",

ifelse(grepl("MCARE", InsurancePlanName), "MCARE_MCAID",

ifelse(grepl("MEDICAID", InsurancePlanName), "MCARE_MCAID",

ifelse(grepl("MCAID", InsurancePlanName), "MCARE_MCAID",

"NO_MCARE_MCAID")))),

LiverDis = ordered(ifelse(ModSevLiverDis == 3, 3,

ifelse(MildLiverDis == 1, 1, 0))),

Diabetes = ordered(ifelse(DiabetesWChrCom == 2, 2,

ifelse(DiabetesWOChrCom == 1, 1, 0))),

Readmit = as.factor(Readmit))

# relevel Readmit factor

snf_data$Readmit <- relevel(snf_data$Readmit, ref = "O")

```

**Check for Missing Data**

```{r, nas}

# Summarize missing values by variable

nas <- map(snf_data, ~sum(is.na(.)))

# nas %>% as_tibble(nas)

# Look at records with missing values

NAlook <- function(x) rowSums(x) > 0

NAs <- snf_data %>% filter(NAlook( across( .cols = everything(), .fns = ~ is.na(.x))))

snf_data <- snf_data %>% na.omit()

```

33 records have at least one NA value, a total of 0.24% of total records. The three variables with missing values contribute to determining whether a patient was readmitted within 30 days (`DischargeDate`) or Hospital score (`HemoglobinLevel`, `SodiumLevel`). These records were removed from the subsequent analysis.

***

## Exploratory Data Analysis

### Aggregated Scores Over Time

**Score Histograms**

```{r hist1}

score_piv <- snf_data %>%

dplyr::select(SternClientID, Year, CCI, LACE, HOSPITAL, Readmit) %>%

pivot_longer(!c(SternClientID, Year, Readmit), names_to = "score", values_to = "val")

ggplot(score_piv, aes(x = val, fill = score)) +

geom_histogram(bins = 20, alpha = 0.5) +

scale_fill_manual(values = c("red","cornflowerblue","goldenrod")) +

theme_minimal() +

labs(x = "score", title = "Score Histograms")

```

```{r hist2}

ggplot(score_piv, aes(x = val, fill = score)) +

geom_histogram(bins = 18, alpha = 0.5) +

facet_wrap(vars(as.factor(Year))) +

scale_fill_manual(values = c("red","cornflowerblue","goldenrod")) +

theme_minimal() +

labs(x = "score", title = "Score Histograms by Year")

```

No apparent changes in distribution for CCI and HOSPITAL scores over time; it looks like LACE score becomes more left-skewed over time. Distributions of the scores are a little skewed, but the sample size is large enough that we can expect assumptions of normality to hold. We can look at each score in separate plots to verify patterns.

```{r histCCI}

CCI_piv <- score_piv %>% filter(score == "CCI")

ggplot(CCI_piv, aes(x = val)) +

geom_histogram(bins = 18, alpha = 0.7, fill = "red") +

theme_minimal() +

labs(x = "CCI", fill = "CCI", title = "CCI Score Histogram")

ggplot(CCI_piv, aes(x = val)) +

geom_histogram(bins = 18, alpha = 0.7, fill = "red") +

theme_minimal() + facet_wrap(vars(Year)) +

labs(x = "CCI", fill = "CCI", title = "CCI Score by Year Histogram")

```

```{r histLACE}

LACE_piv <- score_piv %>% filter(score == "LACE")

ggplot(LACE_piv, aes(x = val)) +

geom_histogram(bins = 18, alpha = 0.7, fill = "cornflowerblue") +

theme_minimal() +

labs(x = "LACE", fill = "LACE", title = "LACE Score Histogram")

ggplot(LACE_piv, aes(x = val)) +

geom_histogram(bins = 18, fill = "cornflowerblue", alpha = 0.7) +

theme_minimal() + facet_wrap(vars(Year)) +

labs(x = "LACE", fill = "LACE", title = "LACE Score by Year Histogram")

```

```{r histHOSP}

HOSP_piv <- score_piv %>% filter(score == "HOSPITAL")

ggplot(HOSP_piv, aes(x = val)) +

geom_histogram(bins = 13, alpha = 0.7, fill = "goldenrod") +

theme_minimal() +

labs(x = "HOSP", fill = "HOSP", title = "HOSP Score Histogram")

ggplot(HOSP_piv, aes(x = val)) +

geom_histogram(bins = 13, fill = "goldenrod", alpha = 0.7) +

theme_minimal() + facet_wrap(vars(Year)) +

labs(x = "HOSP", fill = "HOSP", title = "HOSP Score by Year Histogram")

```

**Density plots**

Distribution of calculated scores by year. These plots tell us how score patterns differ year to year, accounting for the number of records each year. Distribution for each score looks relatively similar year over year, with a few notes.

```{r densplot}

# Density Plot CCI by Year

ggplot(snf_data, aes(x = CCI, after_stat(density), colour = factor(Year))) +

geom_freqpoly(bins = 17) + scale_color_d3() + theme_minimal() +

labs(title = "Density Plot : CCI Score by Year",

colour = "Year")

# Density Plot LACE by Year

ggplot(snf_data, aes(x = LACE, after_stat(density), colour = factor(Year))) +

geom_freqpoly(bins = 19) + scale_color_d3() + theme_minimal() +

labs(title = "Density Plot : LACE Score by Year",

colour = "Year")

# Density Plot Hospital by Year

ggplot(snf_data, aes(x = HOSPITAL, after_stat(density), colour = factor(Year))) +

geom_freqpoly(bins = 17) + scale_color_d3() + theme_minimal() +

labs(title = "Density Plot : HOSPITAL Score by Year",

colour = "Year")

```

We can also review QQ-plots of each score's sampling distribution against a theoretical normal distribution to assess any violations of normality assumptions.

```{r qq_scores}

qqC <- ggplot(snf_data, aes(sample = CCI)) +

stat_qq() + stat_qq_line() +

theme_minimal() + ggtitle("CCI QQNorm Plot")

qqL <- ggplot(snf_data, aes(sample = LACE)) +

stat_qq() + stat_qq_line() +

theme_minimal() + ggtitle("LACE QQNorm Plot")

qqH <- ggplot(snf_data, aes(sample = HOSPITAL)) +

stat_qq() + stat_qq_line() +

theme_minimal() + ggtitle("Hospital QQNorm Plot")

grid.arrange(qqC, qqL, qqH, ncol = 3)

```

Deviations from strict adherance to normality are to be expected after looking at the histogram plots. CCI looks a little wonky in the upper tail; comparing results from compatible parametric and non-parametric tests will make results more reliable. LACE and HOSPITAL qq-plots look "normal enough" at this sample size.

**Summary Statistics for Score**

```{r score_sums}

m_scores <- snf_data %>% group_by(Year) %>%

summarise(across(CCI:HOSPITAL, list(n = ~n(), mean = mean, median = median, sd = sd))) %>% dplyr::select(!c(LACE_n, HOSPITAL_n)) %>%

rename(N = CCI_n)

m_scores %>%

kbl(caption = "Num Records, Mean, Median Scores by Year", digits = 2,

col.names = c("Year", "N", rep(c("Mean", "Median", "SD"), 3))) %>%

kable_styling("striped", full_width = F) %>%

add_header_above(c(" " = 2, "CCI" = 3, "LACE" = 3, "HOSPITAL" = 3))

colors <- c("CCI" = "red", "LACE" = "cornflowerblue", "HOSPITAL" = "goldenrod")

linetypes <- c("mean" = 1, "median" = 22)

ggplot(m_scores, aes(x = Year,)) +

geom_line(aes(y = CCI_mean, color = "CCI", linetype = "mean"), lwd = 1.2) +

geom_line(aes(y = CCI_median, color = "CCI", linetype = "median"), lwd = 1.2) +

geom_line(aes(y = LACE_mean, color = "LACE", linetype = "mean"), lwd = 1.2) +

geom_line(aes(y = LACE_median, color = "LACE", linetype = "median"), lwd = 1.2) +

geom_line(aes(y = HOSPITAL_mean, color = "HOSPITAL", linetype = "mean"), lwd = 1.2) +

geom_line(aes(y = HOSPITAL_median, color = "HOSPITAL", linetype = "median"), lwd = 1.2) +

labs(title = "Mean, Median Scores by Year",

y = "Summarised Score", color = "Score") +

scale_color_manual(values = colors) +

scale_linetype_manual(name = "Statistic", values = linetypes) +

theme_minimal() + theme(text = element_text(size = 24)) +

scale_x_continuous(breaks = seq(2014, 2019, 1)) +

scale_y_continuous(limits = c(2, 12), breaks = seq(0, 12, 2))

```

Median scores are not significantly different from the mean scores, indicating symmetrical structure to all three scores. All three scores increase slightly over the years but move with the same pattern, indicating that these evaluation metrics have been assigned consistently over time. To confirm, variables can be scaled and tested for significant differences in means.

```{r scores_over_time}

CCI_lm <- lm(CCI ~ scale(Year), data = snf_data)

summary(CCI_lm)

CCI_ci <- confint(CCI_lm)[2,]

LACE_lm <- lm(LACE ~ scale(Year), data = snf_data)

summary(LACE_lm)

LACE_ci <- confint(LACE_lm)[2,]

HOS_lm <- lm(HOSPITAL ~ scale(Year), data = snf_data)

summary(HOS_lm)

HOS_ci <- confint(HOS_lm)[2,]

ests <- c(summary(CCI_lm)$coef[2,1], summary(LACE_lm)$coef[2,1], summary(HOS_lm)$coef[2,1])

low <- c(CCI_ci[1], LACE_ci[1], HOS_ci[1])

high <- c(CCI_ci[2], LACE_ci[2], HOS_ci[2])

lm_tab <- rbind(ests, low, high)

colnames(lm_tab) <- c("CCI", "LACE", "HOSPITAL")

rownames(lm_tab) <- c("Est Annual Change", "Low 95%", "High 95%")

kbl(lm_tab, caption = "Estimated Annual Increase in Calculated Scores",

digits = 3) %>%

kable_classic(font_size = 12, html_font = "Cambria") %>%

kable_styling(full_width = F)

```

```{r lmplot}

# Create fitted lines from lms

x <- seq(0, 6, length.out = 200)

fitline <- function(fit){

a <- fit$coefficients[1]

b <- fit$coefficients[2]

y <- a + b*x

return(y)

}

fitlow <- function(fit){

a <- confint(fit)[1,1]

b <- confint(fit)[2,1]

y <- a + b*x

return(y)

}

fithigh <- function(fit){

a <- confint(fit)[1,2]

b <- confint(fit)[2,2]

y <- a + b*x

return(y)

}

CCI_y <- fitline(CCI_lm)

CCI_low <- fitlow(CCI_lm)

CCI_high <- fithigh(CCI_lm)

LACE_y <- fitline(LACE_lm)

LACE_low <- fitlow(LACE_lm)

LACE_high <- fithigh(LACE_lm)

HOS_y <- fitline(HOS_lm)

HOS_low <- fitlow(HOS_lm)

HOS_high <- fithigh(HOS_lm)

lm_df <- data.frame(x, CCI_y, CCI_low, CCI_high,

LACE_y, LACE_low, LACE_high,

HOS_y, HOS_low, HOS_high)

colors <- c("CCI" = "red", "LACE" = "cornflowerblue", "HOSPITAL" = "goldenrod")

linetypes <- c("empirical mean" = 22, "estimated mean" = 1, "CI" = 6)

ggplot(snf_data, aes(x = Year)) +

geom_smooth(mapping = aes(y = CCI, color = "CCI"), method = lm, se=TRUE) +

geom_smooth(mapping = aes(y = LACE, color = "LACE"), method = lm, se = TRUE) +

geom_smooth(mapping = aes(y = HOSPITAL, color = "HOSPITAL"), method = lm, se = TRUE) +

labs(title = "Calculated Scores by Year", x = "Year",

y = "Score", color = "Score") +

scale_color_manual(values = colors) +

scale_linetype_manual(name = "Statistic", values = linetypes) +

theme_minimal() + theme(text = element_text(size = 24)) +

scale_x_continuous(breaks = seq(2014, 2019, 1)) +

scale_y_continuous(limits = c(2, 12), breaks = seq(2, 12, 2))

```

Linear models assess changes in score over time, averaged across all other predictors. Positive estimates associated with small p-values indicate that all scores increase over time. However, the low adjusted $R^2$ referring to proportion of total variability explained by Year lead to the conclusion that other factors are more strongly associated with score changes over time. Additionally, the relatively small estimates of annual score increase may imply limited clinically significant increases in scores over time.

***

### Demographics

**Gender**

Histograms and density plots of patient age by gender show no obvious changes in patient age or gender distributions over time. The data is also sufficiently symmetric.

```{r meangender1}

snf_data %>% group_by(Gender) %>%

summarise(across(PatientAge, list(n = ~n(), mean = mean, median = median, sd = sd))) %>%

rename(N = PatientAge_n, Mean = PatientAge_mean, Median = PatientAge_median,

SD = PatientAge_sd) %>%

kbl(caption = "Patient Age by Gender", digits = 2) %>%

kable_classic(font_size = 12, html_font = "Cambria") %>%

kable_styling(latex_options = c("striped", "HOLD_position"), full_width=F)

```

```{r genderhist1}

age <- snf_data %>%

select(SternClientID, Year, Gender, PatientAge) %>%

filter(Gender != "Unspecified")

ggplot(age, aes(x = PatientAge, fill = Gender)) +

geom_histogram(bins = 30, alpha = 0.7) +

scale_fill_manual(values = c("red","cornflowerblue")) +

theme_minimal() +

labs(title = "Patient Age by Gender")

```

```{r genderhist2}

ggplot(age, aes(x = PatientAge, fill = Gender)) +

geom_histogram(bins = 30, alpha = 0.7) +

facet_wrap(vars(as.factor(Year))) +

scale_fill_manual(values = c("red","cornflowerblue")) +

theme_minimal() +

labs(title = "Patient Age by Gender, by Year")

```

```{r genderhist3}

age_mf <- filter(age, Gender != "Unspecified")

ggplot(age_mf, aes(x = PatientAge, after_stat(density), colour = as.factor(Year))) +

geom_freqpoly(bins = 20) +

facet_wrap(vars(Gender)) +

scale_color_d3() + theme_minimal() +

labs(title = "Patient Age by Year Density Plot, by Gender",

color = "Year")

```

There are consistently more female patients, and age is slightly higher in female patients. Histograms show similar age distributions across genders and little changes over time.

```{r meangender2}

ageyr <- snf_data %>% group_by(Gender, Year) %>%

summarise(across(PatientAge, list(n = ~n(), mean = mean, median = median, sd = sd))) %>%

rename(N = PatientAge_n, Mean = PatientAge_mean, Median = PatientAge_median,

SD = PatientAge_sd)

ageyr[,2:6] %>%

kbl(caption = "Patient Age by Gender, Year", digits = 1) %>%

kable_classic(full_width = F, html_font = "Cambria") %>%

kable_styling(latex_options = c("striped", "HOLD_position")) %>%

pack_rows("Female", 1, 6) %>%

pack_rows("Male", 7, 12) %>%

pack_rows("Unspecified", 13, 14)

```

```{r meangender3}

ageyr <- ageyr %>% filter(Gender != "Unspecified")

linetypes <- c("Mean" = 1, "Median" = 2)

ggplot(ageyr, aes(x = Year)) +

geom_line(aes(y = Mean, color = Gender), lwd = 1.2) +

geom_line(aes(y = Median, color = Gender), lty = 2, lwd = 1.2) +

labs(title = "Mean, Median Age by Gender, Year",

y = "Age", color = "Gender") + theme_minimal() +

scale_color_manual(values = c("red","cornflowerblue")) +

scale_linetype_manual(name = "Statistic", values = linetypes) +

scale_x_continuous(breaks = seq(2014, 2019, 1))

# update legend with lty

```

Graph of mean and median age over time shows a slight increase for both genders; however, median age over time across gender diverges, with median Male age decreasing after 2013. This indicates a possible interaction between Gender and Age that must be checked in the model-building phase.

**Race and Ethnicity**

```{r race}

# Table of counts, frequencies

snf_data %>% count(Race) %>%

mutate(Freq = round(n/sum(n), 3)) %>%

kbl(caption = "Patient Counts and Proportion of Total by Race",

col.names = c("Race", "N", "Frequency")) %>%

kable_classic(full_width = F, html_font = "Cambria") %>%

kable_styling(latex_options = c("striped", "HOLD_position"))

```

```{r avg_race_score}

race_scores <- snf_data %>%

group_by(Race) %>%

summarise(across(CCI:HOSPITAL, list(mean = mean, median = median, sd = sd)))

race_scores %>%

kbl(caption = "Score Summaries by Race", digits = 1,

col.names = c("Race", rep(c("Mean", "Median", "SD"), 3))) %>%

kable_classic(full_width = F, html_font = "Cambria") %>%

kable_styling(latex_options = c("striped", "HOLD_position", "scale_down")) %>%

add_header_above(c(" " = 1, "CCI" = 3, "LACE" = 3, "HOSPITAL" = 3))

```

Patterns across score by race are also consistent.

```{r ethnicity}

# snf_data %>% count(Ethnicity) %>% mutate(freq = round(n/sum(n), 3))

#snf_data %>% group_by(Race, Ethnicity) %>% summarise(n = n()) %>% ungroup() %>%

# mutate(freq = round(n/sum(n), 3))

eth_scores <- snf_data %>% group_by(Ethnicity) %>%

summarise(across(CCI:HOSPITAL, list(mean = mean, median = median, sd = sd)))

eth_scores %>%

kbl(caption = "Score Summary by Ethnicity", digits = 1,

col.names = c("Ethnicity", rep(c("Mean", "Median", "SD"), 3))) %>%

kable_classic(html_font = "Cambria", full_width = F) %>%

kable_styling(latex_options = c("striped", "HOLD_position")) %>%

add_header_above(c(" " = 1, "CCI" = 3, "LACE" = 3, "HOSPITAL" = 3))

```

```{r raceeth}

# Incorporate Hispanic ethnicity into Race variable

snf_data <- snf_data %>%

mutate(Race_Eth = factor(ifelse(Ethnicity == "Hispanic or Latino", "Hispanic/Latino", as.character(Race))))

snf_data %>% group_by(Race_Eth, Year) %>%

summarise(n = n()) %>%

ungroup() %>%

mutate(freq = round(n/sum(n), 3)) %>%

pivot_wider(names_from = Race_Eth, values_from = c(n, freq), names_vary = "slowest") %>%

kbl(caption = "Counts and Proportions of Patients by Race/Ethnicity and Year",

col.names = c("Year", rep(c("N", "Freq"), 7))) %>%

kable_classic(html_font = "Cambria", full_width = F) %>%

kable_styling(latex_options = c("striped", "HOLD_position")) %>%

add_header_above(c(" " = 1, "African American/Black" = 2, "Asian" = 2,

"Caucasian/White" = 2, "Hispanic/Latino" = 2, "Native American/Alaskan" = 2,

"Other/Multiracial" = 2, "Unavailable/Unknown" = 2))

```

```{r race_scores}

race_score_piv <- snf_data %>%

select(SternClientID, Year, Race_Eth, CCI, LACE, HOSPITAL) %>%

pivot_longer(!c(SternClientID, Year, Race_Eth), names_to = "score", values_to = "val")

ggplot(race_score_piv, aes(x = val, fill = Race_Eth)) +

geom_histogram(bins = 15, alpha = 0.7) +

facet_wrap(vars(score)) + scale_fill_d3() +

theme_minimal() + labs(x = "score",

title = "Score Histogram by Race/Ethnicity",

fill = "Race/Ethnicity")

race_CCI_piv <- race_score_piv %>% filter(score == "CCI")

race_LACE_piv <- race_score_piv %>% filter(score == "LACE")

race_HOSP_piv <- race_score_piv %>% filter(score == "HOSPITAL")

ggplot(race_CCI_piv, aes(x = val, fill = Race_Eth)) +

geom_histogram(bins = 12, alpha = 0.7) + scale_fill_d3() +

theme_minimal() + labs(x = "score", title = "CCI by Race/Ethnicity",

fill = "Race/Ethnicity")

ggplot(race_LACE_piv, aes(x = val, fill = Race_Eth)) +

geom_histogram(bins = 12, alpha = 0.7) + scale_fill_d3() +

theme_minimal() + labs(x = "score", title = "LACE by Race/Ethnicity",

fill = "Race/Ethnicity")

ggplot(race_HOSP_piv, aes(x = val, fill = Race_Eth)) +

geom_histogram(bins = 12, alpha = 0.7) + scale_fill_d3() +

theme_minimal() + labs(x = "score", title = "HOSPITAL by Race/Ethnicity",

fill = "Race/Ethnicity")

```

```{r racescores2}

ggplot(race_score_piv, aes(x = val, after_stat(density), color = Race_Eth)) +

geom_freqpoly(bins = 10) + facet_wrap(vars(score)) +

theme_minimal() + scale_color_d3() +

labs(x = "score", title = "Density Plot : Score by Race")

```

```{r racescore_time}

raceyr <- snf_data %>% group_by(Race_Eth, Year) %>%

summarise(across(CCI:HOSPITAL, list(n = ~n(), mean = mean, median = median, sd = sd))) %>%

select(!c(LACE_n, HOSPITAL_n)) %>%

rename(Race_Ethnicity = Race_Eth, n = CCI_n,

meanCCI = CCI_mean, medCCI = CCI_median, sdCCI = CCI_sd,

meanLACE = LACE_mean, medLACE = LACE_median, sdLACE = LACE_sd,

meanHOSP = HOSPITAL_mean, medHOSP = HOSPITAL_median, sdHOSP = HOSPITAL_sd)

raceyr[,2:12] %>%

kbl(caption = "Score Summary by Race/Ethnicity, Year", digits = 1,

col.names = c("Year", "N", rep(c("Mean", "Median", "SD"), 3))) %>%

kable_styling(latex_options = c("striped", "HOLD_position", "scale_down")) %>%

kable_classic(full_width = F, html_font = "Cambria") %>%

add_header_above(c(" " = 2, "CCI" = 3, "LACE" = 3, "HOSP" = 3)) %>%

pack_rows("African American / Black", 1 , 6) %>%

pack_rows("Asian", 7, 12) %>%

pack_rows("Caucasian/White", 13, 18) %>%

pack_rows("Hispanic/Latino", 19, 24) %>%

pack_rows("Native American / Alaskan", 25, 30) %>%

pack_rows("Other / Multiracial", 31, 36) %>%

pack_rows("Unavailable / Unknown", 37, 42)

```

```{r raceyr}

# pivot table for faceting

ryr <- raceyr %>%

select(-c(n, medCCI, sdCCI, medLACE, sdLACE, medHOSP, sdHOSP)) %>%

rename(CCI = meanCCI, LACE = meanLACE, HOSPITAL = meanHOSP) %>%

pivot_longer(!c(Race_Ethnicity, Year),

names_to = "Score", values_to = "Value")

ggplot(ryr, aes(x = Year, y = Value, color = Race_Ethnicity)) +

geom_line(size = 1.2) +

facet_wrap(~Score) +

labs(title = "Empirical Mean Score by Race/Ethnicity, Year",

y = "Score", color = "Race/Ethnicity") + theme_minimal() +

scale_color_d3() +

theme(axis.text.x = element_text(angle = 45))

```

Plots showing both mean and median scores over time for all races were too busy to read clearly. Plots of means over time show a greater increase in all scores for Native Amer/Alaskan individuals, but the sample size for this group is in the single digits for each year observed.

**Insurance**

```{r ins}

ins <- snf_data %>% group_by(InsurancePlanCleaned, Year) %>%

summarise(across(CCI:HOSPITAL, list(n = ~n(), mean = mean,

median = median, sd = sd)))

ggplot(snf_data, aes(x = InsurancePlanCleaned, y = CCI,

fill = InsurancePlanCleaned)) +

geom_boxplot() + facet_wrap(vars(Year)) +

scale_fill_d3() + theme_minimal() +

scale_x_discrete(labels = c(rep(" ", 2))) +

labs(title = "Boxplot of CCI by Insurance Type, by Year",

fill = "Insurance")

ggplot(snf_data, aes(x = InsurancePlanCleaned, y = LACE,

fill = InsurancePlanCleaned)) +

geom_boxplot() + facet_wrap(vars(Year)) +

scale_fill_d3() + theme_minimal() +

scale_x_discrete(labels = c(rep(" ", 2))) +

labs(title = "Boxplot of LACE by Insurance Type, by Year",

fill = "Insurance")

ggplot(snf_data, aes(x = InsurancePlanCleaned, y = HOSPITAL,

fill = InsurancePlanCleaned)) +

geom_boxplot() + facet_wrap(vars(Year)) +

scale_fill_d3() + theme_minimal() +

scale_x_discrete(labels = c(rep(" ", 2))) +

labs(title = "Boxplot of HOSPITAL by Insurance Type, by Year",

fill = "Insurance")

```

***

# Analysis

Analysis of Variance is conducted for each of the three calculated scores (CCI, LACE, and HOSPITAL) for association with the following variables:

* **Year**: 2014 - 2019

* **Patient Age**: Patient age at time of hospitalization, ranging from 18 to 106. This variable is condensed into an ordered categorical variable and used in calculating CCI and LACE scores, so it is excluded from those analyses. Since association with PatientAge is inherent to these calculated scores, including the term would unnecessarily complicate interaction analysis without adding additional information.

* **Gender**: Male or Female. Three observations with `Gender` = "Unspecified" are excluded from analysis; the small group size produces erroneous errors when testing for interactions.

* **Race/Ethnicity**: Factor with six levels: "African American/Black", "Asian", Caucasian/White", "Hispanic/Latino", "Other/Multiracial", and "Unavailable/Unknown". 27 observations with `Race\Eth` = "Native American/Alaskan" are combined with "Other/Multiracial" since the small group size produces erroneous errors when testing for interactions.

* **Insurance**: Factor with two levels: "Medicare/Medicaid" and "Not Medicare/Medicaid". 235 unique values for Insurance Plan were combed for words "Medicare", "MCare", "Medicaid", "MCaid" and categorized as "Medicare/Medicaid"; all other observations are categorized as "Not Medicare/Medicaid".

All p-values within each ANOVA are adjusted using the Benjamini, Hochberg (1995) correction for False Discovery Rate. Since each score is calculated differently and is used as a different type of evaluation metric, p-value adjustment is not performed between ANOVAs.

Current post-hoc analysis examines differences of means for main effects and any identified interactions.

```{r data_prep}

snf_data_a <- snf_data %>%

filter(Gender != "Unspecified") %>%

mutate(Race_Eth = fct_recode(Race_Eth, "Other/Multiracial" = "Native Amer/Alaskan"),

Insurance = as.factor(InsurancePlanCleaned))

snf_data_a$Gender <- droplevels(snf_data_a$Gender)

```

## ANOVA

### CCI

After correcting for multiple comparisons, `Year` and `Race/Ethnicity`, and `Insurance` are strongly associated with changes in CCI scores. No interactions between demographic factors or time are indicated. Gender is not associated with meaningful differences in CCI scores.

```{r cci_aov}

CCI_mod <- lm(CCI ~ Year * Gender * Race_Eth * Insurance,

data = snf_data_a)

CCI_a <- Anova(CCI_mod, type = "II")

bhp <- round(p.adjust(CCI_a$`Pr(>F)`, method = "BH"), digits = 3)[1:15]

CCI_aov <- anova_summary(CCI_a, detailed = TRUE) %>%

select(c(1, 2, 4, 6)) %>%

mutate(`p-value` = bhp) %>%

mutate_if(is.numeric, round, 3)

CCI_aov %>%

kbl(caption = "ANOVA Table : CCI by Year, Demographics") %>%

kable_classic(full_width = F, html_font = "Cambria") %>%

kable_styling(latex_options = "HOLD_position")

```

### LACE

All main effects appear associated with LACE scores. LACE also incorporates age into calculation so `PatientAge` is excluded from the analysis. Among possible interactions, `Year`:`Insurance` is identified for further exploration.

```{r lace_aov}

LACE_mod <- lm(LACE ~ Year * PatientAge * Gender * Race_Eth * Insurance,

data = snf_data_a)

LACE_a <- Anova(LACE_mod, type = "II")

bhp <- round(p.adjust(LACE_a$`Pr(>F)`, method = "BH"), digits = 3)[1:31]

LACE_aov <- anova_summary(LACE_a, detailed = TRUE) %>%

select(c(1, 2, 4, 6)) %>%

mutate(p_adj = bhp) %>%

mutate_if(is.numeric, round, 3)

LACE_aov %>%

kbl(caption = "ANOVA Table : LACE by Year, Demographics") %>%

kable_classic(full_width = F, html_font = "Cambria") %>%

kable_styling(latex_options = "HOLD_position")

```

### HOSPITAL

We explored formatting PatientAge as a categorical variable with the same groups used to calculate LACE and CCI scores, but there were too many group combinations so PatientAge remains a continuous covatiate.

All main effects of demographic variables examined show association with changes in HOSPITAL scores. Additional two-way interactions to explore include `Year`:`Race/Ethnicity`, `PatientAge`:`Insurance`, `Year`:`PatientAge`, `Race/Ethnicity`:`Insurance`, along with a three-way interaction between `PatientAge`:`Race/Ethnicity`:`Insurance`.

```{r hosp_aov}

HOSP_mod <- lm(HOSPITAL ~ Year * PatientAge * Gender * Race_Eth * Insurance,

data = snf_data_a)

HOSP_a <- Anova(HOSP_mod, type = "II")

bhp <- round(p.adjust(HOSP_a$`Pr(>F)`, method = "BH"), digits = 3)[1:31]

HOSP_aov <- anova_summary(HOSP_a, detailed = TRUE) %>%

select(c(1, 2, 4, 6)) %>%

mutate(p_adj = bhp) %>%

mutate_if(is.numeric, round, 3)

HOSP_aov %>%

kbl(caption = "ANOVA Table : HOSPITAL by Year, Demographics") %>%

kable_classic(full_width = F, html_font = "Cambria") %>%

kable_styling(latex_options = "HOLD_position")

```

## GROUP CONTRASTS

```{r CCI_groupcontr}

# Include all main effects

CCI_cont <- lm(CCI ~ Year + Gender + Race_Eth + Insurance, data = snf_data_a)

```

```{r lace_grpcontr}

LACE_cont <- lm(LACE ~ Year + Gender + Race_Eth + Insurance + Race_Eth:Insurance +

Year:Insurance + Year:Race_Eth:Insurance, data = snf_data_a)

```

```{r hosp_grcontr}

HOSP_cont <- lm(HOSPITAL ~ Year + PatientAge + Gender + Race_Eth + Insurance +

Year:PatientAge + Year:Race_Eth + PatientAge:Insurance + Race_Eth:Insurance, data = snf_data_a)

```

### Gender

#### CCI

```{r CCI_gen}

# Marginal means and CIs

CCI_gen_emm <- emmeans(CCI_cont, "Gender")

# Test for significance

CCI_gen_pair <- pairs(CCI_gen_emm, adjust = "tukey")

# est marginal means and CIs

CCI_gen_cont <- confint(CCI_gen_emm, calc = c(n = ~.wgt.))

```

#### LACE

```{r lace_gen}

# Marginal means and CIs

LACE_gen_emm <- emmeans(LACE_cont, "Gender")

# Test for significance

LACE_gen_pair <- pairs(LACE_gen_emm, adjust = "tukey")

# est marginal means and CIs

LACE_gen_cont <- confint(LACE_gen_emm, calc = c(n = ~.wgt.))

```

#### HOSPITAL

```{r hosp_gen}

# Marginal means and CIs

HOSP_gen_emm <- emmeans(HOSP_cont, "Gender")

# Test of significance

HOSP_gen_pair <- pairs(HOSP_gen_emm, adjust = "tukey")

# Est marginal means and CIs

plot(HOSP_gen_emm, comparisons = TRUE) +

ggtitle("HOSPITAL : Estimated Marginal Means by Gender, with direction and CI")

HOS_gen_cont <- confint(HOSP_gen_emm, calc = c(n = ~.wgt.))

```

### Race/Ethnicity

#### CCI

```{r CCI_re}

# Marginal means and CIs

CCI_reth_emm <- emmeans(CCI_cont, "Race_Eth")

# pairwise comparison test of significance

CCI_reth_pair <- pairs(CCI_reth_emm)

# plots

plot(CCI_reth_emm, comparisons = TRUE) +

ggtitle("CCI : Est Marginal Means by Race/Ethnicity, with direction and CI")

emmip(CCI_cont, Race_Eth ~ Year, cov.reduce = range) +

theme_minimal() +

labs(title = "Annual Trend of CCI Score by Race/Ethnicity",

color = "Race/Ethnicity")

```

#### LACE

```{r lace_reth}

# est marginal means and CIs

LACE_reth_emm <- emmeans(LACE_cont, "Race_Eth")

# pairwise comparisons

LACE_reth_pair <- pairs(LACE_reth_emm)

# plots

plot(LACE_reth_emm, comparisons = TRUE) +

ggtitle("LACE : Est Marginal Means by Race/Ethnicity, with direction and CI")

```

#### HOSPITAL

```{r hosp_reth}

# emmeans and CIs

HOSP_reth_emm <- emmeans(HOSP_cont, "Race_Eth")

# pairwise comparisons

HOSP_reth_pair <- pairs(HOSP_reth_emm)

# plots

plot(HOSP_reth_emm, comparisons = TRUE) +

ggtitle("HOSPITAL : Estimated Marginal Means by Race/Ethnicity, with direction and CI")

```

### Insurance

#### CCI

```{r CCI_ins}

# emmeans and CI

CCI_ins_emm <- emmeans(CCI_cont, "Insurance")

# pairwise comparisons

CCI_ins_pair <- pairs(CCI_ins_emm)

# plots

plot(CCI_ins_emm, comparisons = TRUE) +

ggtitle("CCI : Est Marginal Means by Insurance, with direction and CI")

```

#### LACE

```{r lace_ins}

# emmeans and CIs

LACE_ins_emm <- emmeans(LACE_cont, "Insurance")

# pairwise comparisons

LACE_ins_pair <- pairs(LACE_ins_emm)

# plots

plot(LACE_ins_emm, comparisons = TRUE) +

ggtitle("LACE : Est Marginal Means by Insurance, with direction and CI")

```

#### HOSPITAL

```{r hosp_ins}

# emmeans and CIs

HOSP_ins_emm <- emmeans(HOSP_cont, "Insurance")

# pairwise comparisons

HOSP_ins_pair <- pairs(HOSP_ins_emm)

# plots

plot(HOSP_ins_emm, comparisons = TRUE) +

ggtitle("HOSPITAL : Estimated Marginal Means by Insurance, with direction and CI")

```

#### Main Effects Table

```{r maintab}

# Gender

CCI_cdf <- data.frame(CCI_gen_emm)

LACE_cdf <- data.frame(LACE_gen_emm)

HOS_cdf <- data.frame(HOSP_gen_emm)

fem <- c(round(CCI_cdf[1,2], digits = 2),

paste0("(", round(CCI_cdf[1,5], digits = 2), ", ", round(CCI_cdf[1,6], digits = 2), ")"),

round(LACE_cdf[1,2], digits = 2),

paste0("(", round(LACE_cdf[1,5], digits = 2), " , ", round(LACE_cdf[1,6], digits = 2), ")"),

round(HOS_cdf[1,2], digits = 2),

paste0("(", round(HOS_cdf[1,5], digits = 2), " , ", round(HOS_cdf[1,6], digits = 2), ")"))

male <- c(round(CCI_cdf[2,2], digits = 2),

paste0("(", round(CCI_cdf[2,5], digits = 2), ", ", round(CCI_cdf[2,6], digits = 2), ")"),

round(LACE_cdf[2,2], digits = 2),

paste0("(", round(LACE_cdf[2,5], digits = 2), " , ", round(LACE_cdf[2,6], digits = 2), ")"),

round(HOS_cdf[2,2], digits = 2),

paste0("(", round(HOS_cdf[2,5], digits = 2), " , ", round(HOS_cdf[2,6], digits = 2), ")"))

# Race/Ethnicity

CCI_r_cdf <- data.frame(CCI_reth_emm)

LACE_r_cdf <- data.frame(LACE_reth_emm)

HOS_r_cdf <- data.frame(HOSP_reth_emm)

bl <- c(round(CCI_r_cdf[1,2], digits = 2),

paste0("(", round(CCI_cdf[1,5], digits = 2), ", ", round(CCI_r_cdf[1,6], digits = 2), ")"),

round(LACE_r_cdf[1,2], digits = 2),

paste0("(", round(LACE_r_cdf[1,5], digits = 2), " , ", round(LACE_r_cdf[1,6], digits = 2), ")"),

round(HOS_r_cdf[1,2], digits = 2),

paste0("(", round(HOS_r_cdf[1,5], digits = 2), " , ", round(HOS_r_cdf[1,6], digits = 2), ")"))

as <- c(round(CCI_r_cdf[2,2], digits = 2),

paste0("(", round(CCI_r_cdf[2,5], digits = 2), ", ", round(CCI_r_cdf[2,6], digits = 2), ")"),

round(LACE_r_cdf[2,2], digits = 2),

paste0("(", round(LACE_r_cdf[2,5], digits = 2), " , ", round(LACE_r_cdf[2,6], digits = 2), ")"),

round(HOS_r_cdf[2,2], digits = 2),

paste0("(", round(HOS_r_cdf[2,5], digits = 2), " , ", round(HOS_r_cdf[2,6], digits = 2), ")"))

cau <- c(round(CCI_r_cdf[3,2], digits = 2),

paste0("(", round(CCI_r_cdf[3,5], digits = 2), ", ", round(CCI_r_cdf[3,6], digits = 2), ")"),

round(LACE_r_cdf[3,2], digits = 2),

paste0("(", round(LACE_r_cdf[3,5], digits = 2), " , ", round(LACE_r_cdf[3,6], digits = 2), ")"),

round(HOS_r_cdf[3,2], digits = 2),

paste0("(", round(HOS_r_cdf[3,5], digits = 2), " , ", round(HOS_r_cdf[3,6], digits = 2), ")"))

his <- c(round(CCI_r_cdf[4,2], digits = 2),

paste0("(", round(CCI_r_cdf[4,5], digits = 2), ", ", round(CCI_r_cdf[4,6], digits = 2), ")"),

round(LACE_r_cdf[4,2], digits = 2),

paste0("(", round(LACE_r_cdf[4,5], digits = 2), " , ", round(LACE_r_cdf[4,6], digits = 2), ")"),

round(HOS_r_cdf[4,2], digits = 2),

paste0("(", round(HOS_r_cdf[4,5], digits = 2), " , ", round(HOS_r_cdf[4,6], digits = 2), ")"))

oth <- c(round(CCI_r_cdf[5,2], digits = 2),

paste0("(", round(CCI_r_cdf[5,5], digits = 2), ", ", round(CCI_r_cdf[5,6], digits = 2), ")"),

round(LACE_r_cdf[5,2], digits = 2),

paste0("(", round(LACE_r_cdf[5,5], digits = 2), " , ", round(LACE_r_cdf[5,6], digits = 2), ")"),

round(HOS_r_cdf[5,2], digits = 2),

paste0("(", round(HOS_r_cdf[5,5], digits = 2), " , ", round(HOS_r_cdf[5,6], digits = 2), ")"))

unk <- c(round(CCI_r_cdf[6,2], digits = 2),

paste0("(", round(CCI_r_cdf[6,5], digits = 2), ", ", round(CCI_r_cdf[6,6], digits = 2), ")"),

round(LACE_r_cdf[6,2], digits = 2),

paste0("(", round(LACE_r_cdf[6,5], digits = 2), " , ", round(LACE_r_cdf[6,6], digits = 2), ")"),

round(HOS_r_cdf[6,2], digits = 2),

paste0("(", round(HOS_r_cdf[6,5], digits = 2), " , ", round(HOS_r_cdf[6,6], digits = 2), ")"))

# Insurance

CCI_i_cdf <- data.frame(CCI_ins_emm)

LACE_i_cdf <- data.frame(LACE_ins_emm)

HOS_i_cdf <- data.frame(HOSP_ins_emm)

care <- c(round(CCI_i_cdf[1,2], digits = 2),

paste0("(", round(CCI_i_cdf[1,5], digits = 2), ", ", round(CCI_i_cdf[1,6], digits = 2), ")"),

round(LACE_i_cdf[1,2], digits = 2),

paste0("(", round(LACE_i_cdf[1,5], digits = 2), " , ", round(LACE_i_cdf[1,6], digits = 2), ")"),

round(HOS_i_cdf[1,2], digits = 2),

paste0("(", round(HOS_i_cdf[1,5], digits = 2), " , ", round(HOS_i_cdf[1,6], digits = 2), ")"))

nocare <- c(round(CCI_i_cdf[2,2], digits = 2),

paste0("(", round(CCI_i_cdf[2,5], digits = 2), ", ", round(CCI_i_cdf[2,6], digits = 2), ")"),

round(LACE_i_cdf[2,2], digits = 2),

paste0("(", round(LACE_i_cdf[2,5], digits = 2), " , ", round(LACE_i_cdf[2,6], digits = 2), ")"),

round(HOS_i_cdf[2,2], digits = 2),

paste0("(", round(HOS_i_cdf[2,5], digits = 2), " , ", round(HOS_i_cdf[2,6], digits = 3), ")"))

# make dataframe