---
title: "Is green land cover associated with less health care spending? Promising findings from county-level Medicare spending in the continental United States"
author: "Matthew Browning"
date: "February 22, 2018"
output: html_document
---
#Set up
```{r setup}
  #setup R Markdown
library(knitr)
knitr::opts_chunk$set(echo = TRUE)
  #load packges
library(readr) #for loading CSV file
library(maps) #for maps and map data
library(plyr) #for merging datasets and renaming variables
library(sjPlot) #for APA formatted tables
library(lme4) # GLMMs and VIFs
library(nlme) # I believe this is needed for GLMM too
library(Rmisc) #create quartiles
library(ggplot2) #plots
library(gridExtra) #side by side plots
library(gtable) #dummy plots
library(grid) #overlap dummy plot
library(cowplot) #multiple ggplots combined
library(readxl) #read excel file
library(nlme) #generalized linear least square models

  #load packages needed for GIS and spatial error models
library(maptools)
library(rgdal)
library(spdep)
library(ctv)
library(lmtest)

  #make results replicable
set.seed(3221)
  #set number of significant digits
options(digits=3)  
  #do not print scientific notation
options(scipen=999)

  #set working directory
setwd("~/Google Drive File Stream/My Drive/Projects/Tree$/Tree$_RA/Pilot_Study")
```
*****
```{r load data}
  #load master datafile
dat <- read_csv("Green_land_cover_Medicare_expenditures_data.csv", na = "NA")

  #transform dwellings per square mile
dat$HouseDensityLog <- log(dat$HouseDensity)

  #create education-income index
    #scale income 
dat$Income_scaled <- scale(dat$Income, center = T, scale = T)
    #scale education
dat$Ed_scaled <- scale(dat$PctCollegeDegree, center = T, scale = T)
    #create index
dat$EdIncome <- (dat$Income_scaled + dat$Ed_scaled) / 2

  #calculate FFS per total population ratio
dat$UserConcentration <- 100 * (dat$PopMedicare / dat$PopTotal)
```
******
#Descriptive Statistics
```{r descriptive statistics}
  #subset and rearrange columns
dat_desc_pre <- dat[c(6:9, 18:23, 25:29, 37, 40, 41)]
dat_desc <- dat_desc_pre[c(1,2,3,4,16,17,18,5,6,7,8,9,10,11,12,13,14,15)]
  #create and save table
sjt.df(dat_desc, altr.row.col = T, title = "Descriptive statistics")
```
******
#Bivariate correlations
```{r bivariate correlation}
  #subset and rearrange columns
dat_corr_pre <- dat[c(6:9, 18:21, 22, 23, 25:29, 40, 41)]
dat_corr <- dat_corr_pre[c(1,2,3,4,16,17,5,6,7,8,9,10,11,12,13,14,15)]
  #create and save table with fewer number of significant digits and no leading zeros
sjt.corr(dat_corr, triangle = "lower", title = "Bivariate correlations", digits = 2)
```
******
#Plots of greenness and expenditures
``` {r bivariate plotes of green land cover and expenditures}
#forest
dat_q_forest <- within(dat, quartile <- as.integer(cut(Forest, quantile(Forest, probs=0:4/4), include.lowest=TRUE)))
dat_q_forest_sum <- summarySE(dat_q_forest,measurevar="Expenditures",groupvars=c("quartile"))
Forest_plot <- ggplot(dat_q_forest_sum, aes(x=quartile, y=Expenditures)) +
  geom_errorbar(aes(ymin=Expenditures - (2*se),
                    ymax=Expenditures + (2*se),
                    width=.05)) + geom_line() +
  scale_y_continuous(limits = c(8725, 10000), labels = scales::dollar_format("$")) +
  xlab("Forest") +
  ylab("") +
  theme(axis.title.x=element_text(face="bold"),)

#shrubs
dat_q_shrub <- within(dat, quartile <- as.integer(cut(Shrub, quantile(Shrub, probs=0:4/4), include.lowest=TRUE)))
dat_q_shrub_sum <- summarySE(dat_q_shrub, measurevar="Expenditures", groupvars=c("quartile"))
Shrub_plot <- ggplot(dat_q_shrub_sum, aes(x=quartile, y=Expenditures)) +
  geom_errorbar(aes(ymin=Expenditures - (2*se),
                    ymax=Expenditures + (2*se),
                    width=.05)) + geom_line() +
  scale_y_continuous(limits = c(8725, 10000)) +
  xlab("Shrub") +
  ylab("") +
  theme(axis.text.y=element_blank(),) +
  theme(axis.title.x=element_text(face="bold"),)

#grass
dat_q_grass <- within(dat, quartile <- as.integer(cut(Grass, quantile(Grass, probs=0:4/4), include.lowest=TRUE)))
dat_q_grass_sum <- summarySE(dat_q_grass, measurevar="Expenditures", groupvars=c("quartile"))
Grass_plot <- ggplot(dat_q_grass_sum, aes(x=quartile, y=Expenditures)) +
  geom_errorbar(aes(ymin=Expenditures - (2*se),
                    ymax=Expenditures + (2*se),
                    width=.05)) + geom_line() +
  scale_y_continuous(limits = c(8725, 10000)) +
  xlab("Grass") +
  ylab("") +
  theme(axis.text.y=element_blank(),) +
  theme(axis.title.x=element_text(face="bold"),)

#urban vegetation
dat_q_urbanveg <- within(dat, quartile <- as.integer(cut(UrbanVeg, quantile(UrbanVeg, probs=0:4/4), include.lowest=TRUE)))
dat_q_urbanveg_sum <- summarySE(dat_q_urbanveg,measurevar="Expenditures",groupvars=c("quartile"))
Urbanveg_plot <- ggplot(dat_q_urbanveg_sum, aes(x=quartile, y=Expenditures)) +
  geom_errorbar(aes(ymin=Expenditures - (2*se),
                    ymax=Expenditures + (2*se),
                    width=.05)) + geom_line() +
  scale_y_continuous(limits = c(8725, 10000), labels = scales::dollar_format("$")) +
  xlab("Urban Veg") +
  ylab("") +
  theme(axis.text.y=element_blank(),) +
  theme(axis.title.x=element_text(face="bold"),)

#agriculture
dat_q_agriculture <- within(dat, quartile <- as.integer(cut(Agriculture, quantile(Agriculture, probs=0:4/4), include.lowest=TRUE)))
dat_q_agriculture_sum <- summarySE(dat_q_agriculture,measurevar="Expenditures",groupvars=c("quartile"))
Agriculture_plot <- ggplot(dat_q_agriculture_sum, aes(x=quartile, y=Expenditures)) +
  geom_errorbar(aes(ymin=Expenditures - (2*se),
                    ymax=Expenditures + (2*se),
                    width=.05)) + geom_line() +
  scale_y_continuous(limits = c(8725, 10000), labels = scales::dollar_format("$")) +
  xlab("Agriculture") +
  ylab("") +
  theme(axis.text.y=element_blank(),) +
  theme(axis.title.x=element_text(face="bold"),)

#plots
NLCD_plot <- plot_grid(Forest_plot, Shrub_plot, Grass_plot, Urbanveg_plot, Agriculture_plot, labels= "AUTO", ncol = 5)
save_plot("Bivariate plots of green land cover and expenditures.jpg", NLCD_plot, ncol = 5, base_aspect_ratio = 1.3)
```
******
#Checking statistical differences between bivariate plot means
```{r ANOVA for quartiles}
  #Forest
#set quartile value as a factor, not a numeric
dat_q_forest$quartile <- as.factor(dat_q_forest$quartile)
#run ANOVA
ANOVA_forest <- aov(Expenditures ~ quartile, data=dat_q_forest)
#print results
summary(ANOVA_forest)
TukeyHSD(ANOVA_forest)

  #Urban vegetation
#set quartile value as a factor, not a numeric
dat_q_urbanveg$quartile <- as.factor(dat_q_urbanveg$quartile)
#run ANOVA
ANOVA_urbanveg <- aov(Expenditures ~ quartile, data=dat_q_urbanveg)
#print results
summary(ANOVA_urbanveg)
TukeyHSD(ANOVA_urbanveg)

  #Agriculture
#set quartile value as a factor, not a numeric
dat_q_agriculture$quartile <- as.factor(dat_q_agriculture$quartile)
#run ANOVA
ANOVA_agriculture <- aov(Expenditures ~ quartile, data=dat_q_agriculture)
#print results
summary(ANOVA_agriculture)
TukeyHSD(ANOVA_agriculture)
```
******
#Run GLMMs
```{r run generalized linear mixed model}
#Run GLMM
model_GLMM_state <- lme4::lmer(Expenditures~Female+Age+White+HouseDensityLog+EdIncome+UserConcentration+
                                 Doctors+Hospitals+HospitalBeds+
                                 Particulates+Inactive+
                                 Forest+Shrub+Grass+UrbanVeg+
                                 (1 | State),
                               data=dat,
                               REML=F)
sjt.lmer(model_GLMM_state, show.aic=T, p.kr=F, emph.p=T, p.numeric=F, show.std = T)

  #Check multicollinearity
fit <- model_GLMM_state
vif.lme <- function (fit) {
  v <- vcov(fit)
  nam <- names(fixef(fit))
  ns <- sum(1 *
  (nam == "Intercept" | nam == "(Intercept)"))
  if (ns > 0) {
  v <- v[-(1:ns), -(1:ns), drop = FALSE]
  nam <- nam[-(1:ns)] }
  d <- diag(v)^0.5
  v <- diag(solve(v/(d %o% d)))
  names(v) <- nam
  v }
  #display VIFs
vif.lme(fit)

```
*****
#Quartiles of Education-Income
``` {r Models with education-income index quartiles}
#create new income dataframe with assigned quartiles
dat_edincome_q <- within(dat, quartile <- as.integer(cut(EdIncome, quantile(EdIncome, probs=0:4/4), include.lowest=TRUE)))

#separate dataframes by educ income quartile
dat_edincome_q1 <- dat[(dat_edincome_q$quartile == 1),]
dat_edincome_q2 <- dat[(dat_edincome_q$quartile == 2),]
dat_edincome_q3 <- dat[(dat_edincome_q$quartile == 3),]
dat_edincome_q4 <- dat[(dat_edincome_q$quartile == 4),]

##Education Income quartile 1
GLMM_edincome_q1 <- lme4::lmer(Expenditures~
                               Female+Age+White+HouseDensityLog+EdIncome+UserConcentration+
                               Doctors+Hospitals+HospitalBeds+
                               Particulates+Inactive+
                               Forest+Shrub+Grass+UrbanVeg+
                               (1 | State),
                               REML=F,
                             data=dat_edincome_q1)

##Education Income quartile 2
GLMM_edincome_q2 <- lme4::lmer(Expenditures~
                               Female+Age+White+HouseDensityLog+EdIncome+UserConcentration+
                               Doctors+Hospitals+HospitalBeds+
                               Particulates+Inactive+
                               Forest+Shrub+Grass+UrbanVeg+
                               (1 | State),
                               REML=F,
                             data=dat_edincome_q2)

##Education Income quartile 3
GLMM_edincome_q3 <- lme4::lmer(Expenditures~
                               Female+Age+White+HouseDensityLog+EdIncome+UserConcentration+
                               Doctors+Hospitals+HospitalBeds+
                               Particulates+Inactive+
                               Forest+Shrub+Grass+UrbanVeg+
                               (1 | State),
                               REML=F,
                             data=dat_edincome_q3)

##Education Income quartile 4
GLMM_edincome_q4 <- lme4::lmer(Expenditures~
                               Female+Age+White+HouseDensityLog+EdIncome+UserConcentration+
                               Doctors+Hospitals+HospitalBeds+
                               Particulates+Inactive+
                               Forest+Shrub+Grass+UrbanVeg+
                               (1 | State),
                               REML=F,
                             data=dat_edincome_q4)

sjt.lmer(GLMM_edincome_q1, GLMM_edincome_q2, GLMM_edincome_q3, GLMM_edincome_q4, show.aic=T, p.kr=F, emph.p=T, p.numeric=F, show.std = T)
```
*****
#Quartiles of Density
``` {r Models with density quartiles}
#create new income dataframe with assigned quartiles
dat_HouseDensity_q <- within(dat, quartile <- as.integer(cut(HouseDensityLog, quantile(HouseDensityLog, probs=0:4/4), include.lowest=TRUE)))

#separate dataframes by educ income quartile
dat_HouseDensity_q1 <- dat[(dat_HouseDensity_q$quartile == 1),]
dat_HouseDensity_q2 <- dat[(dat_HouseDensity_q$quartile == 2),]
dat_HouseDensity_q3 <- dat[(dat_HouseDensity_q$quartile == 3),]
dat_HouseDensity_q4 <- dat[(dat_HouseDensity_q$quartile == 4),]

##Density quartile 1
GLMM_HouseDensity_q1 <- lme4::lmer(Expenditures~
                               Female+Age+White+HouseDensityLog+HouseDensity+UserConcentration+
                               Doctors+Hospitals+HospitalBeds+
                               Particulates+Inactive+
                               Forest+Shrub+Grass+UrbanVeg+
                               (1 | State),
                               REML=F,
                             data=dat_HouseDensity_q1)

##Density quartile 2
GLMM_HouseDensity_q2 <- lme4::lmer(Expenditures~
                               Female+Age+White+HouseDensityLog+HouseDensity+UserConcentration+
                               Doctors+Hospitals+HospitalBeds+
                               Particulates+Inactive+
                               Forest+Shrub+Grass+UrbanVeg+
                               (1 | State),
                               REML=F,
                             data=dat_HouseDensity_q2)

##Density quartile 3
GLMM_HouseDensity_q3 <- lme4::lmer(Expenditures~
                               Female+Age+White+HouseDensityLog+HouseDensity+UserConcentration+
                               Doctors+Hospitals+HospitalBeds+
                               Particulates+Inactive+
                               Forest+Shrub+Grass+UrbanVeg+
                               (1 | State),
                               REML=F,
                             data=dat_HouseDensity_q3)

##Density quartile 4
GLMM_HouseDensity_q4 <- lme4::lmer(Expenditures~
                               Female+Age+White+HouseDensityLog+HouseDensity+UserConcentration+
                               Doctors+Hospitals+HospitalBeds+
                               Particulates+Inactive+
                               Forest+Shrub+Grass+UrbanVeg+
                               (1 | State),
                               REML=F,
                             data=dat_HouseDensity_q4)

sjt.lmer(GLMM_HouseDensity_q1, GLMM_HouseDensity_q2, GLMM_HouseDensity_q3, GLMM_HouseDensity_q4, show.aic=T, p.kr=F, emph.p=T, p.numeric=F, show.std = T)
```
******
#Other spatial models
``` {r compare GLMM to other spatial models}

 ###### prepare data ###### 

  #load shapefile
dat_shp_pre<-readOGR(dsn="Green_land_cover_Medicare_expenditures_shapefile.shp")

  #remove non-numeric variables
dat_shp <- dat_shp_pre[c(7:30)]

  #change variables that loaded as factors and should be numeric
dat_shp$Expend <- as.numeric(dat_shp$Expend)
dat_shp$Female <- as.numeric(dat_shp$Female)
dat_shp$Age <- as.numeric(dat_shp$Age)
dat_shp$White <- as.numeric(dat_shp$White)
dat_shp$Income <- as.numeric(dat_shp$Income)
dat_shp$PctColl <- as.numeric(dat_shp$PctColl)
dat_shp$PopTot <- as.numeric(dat_shp$PopTot)
dat_shp$PopMed <- as.numeric(dat_shp$PopMed)
dat_shp$PopDens <- as.numeric(dat_shp$PopDens)
dat_shp$Dwellings <- as.numeric(dat_shp$Dwellings)
dat_shp$HouseDens <- as.numeric(dat_shp$HouseDens)
dat_shp$Doctors <- as.numeric(dat_shp$Doctors)
dat_shp$Hospitals <- as.numeric(dat_shp$Hospitals)
dat_shp$HospitalB <- as.numeric(dat_shp$HospitalB)
dat_shp$Partic <- as.numeric(dat_shp$Partic)
dat_shp$Inactive <- as.numeric(dat_shp$Inactive)
dat_shp$Grass <- as.numeric(dat_shp$Grass)
dat_shp$Shrub <- as.numeric(dat_shp$Shrub)
dat_shp$Forest <- as.numeric(dat_shp$Forest)
dat_shp$Agricult <- as.numeric(dat_shp$Agricult)
dat_shp$UrbanVeg <- as.numeric(dat_shp$UrbanVeg)
dat_shp$Area <- as.numeric(dat_shp$Area)
dat_shp$Ozone <- as.numeric(dat_shp$Ozone)
dat_shp$HealthyFd <- as.numeric(dat_shp$HealthyFd)

  #transform housing density
dat_shp$HouseDensL <- log(dat_shp$HouseDens)

  #create education-income index
    #scale income 
dat_shp$Income <- as.numeric(dat_shp$Income)
dat_shp$Income_scaled <- scale(dat_shp$Income, center = T, scale = T)
    #scale education
dat_shp$Ed_scaled <- scale(dat_shp$PctColl, center = T, scale = T)
    #create index
dat_shp$EdIncome <- (dat_shp$Income_scaled + dat_shp$Ed_scaled) / 2

  #calculate FFS per total population ratio
dat_shp$UserConc <- 100 * (dat_shp$PopMed / dat_shp$PopTot)

  #assign coordinates
coords<-coordinates(dat_shp)


 ###### Prepare neighbors lists for shapefile ###### 

  #create neighbors file based on 'k' number of nearest neighbors
dat_shp_kn1 <- knn2nb(knearneigh(coords, k = 1))
dat_shp_kn2 <- knn2nb(knearneigh(coords, k = 2))
dat_shp_kn4 <- knn2nb(knearneigh(coords, k = 4))

  #plot links between counties
plot(dat_shp)
plot(dat_shp_kn2, coords, add = T)

  #determine appropraite threshold for distance-based neighbors
dist <- unlist(nbdists(dat_shp_kn1, coords))
summary(dist)

  #store maximum distance
max_k1 <- max(dist)

  #create distance-based neighbors lists
dat_shp_kd1 <- dnearneigh(coords, d1 = 0, d2 = 0.75 * max_k1)
dat_shp_kd2 <- dnearneigh(coords, d1 = 0, d2 = 1 * max_k1)
dat_shp_kd3 <- dnearneigh(coords, d1 = 0, d2 = 1.5 * max_k1)

  #assign Queen's case neighbors, the zero-policy ensures no counties have zero neighbors
dat_shp_kd2_w <- nb2listw(dat_shp_kd2, zero.policy = T)
dat_shp_kd2_w

  #alternatively, create neighbors file based on 'k' number of nearest neighbors
dat_shp_kn1 <- knn2nb(knearneigh(coords, k = 1))
dat_shp_kn2 <- knn2nb(knearneigh(coords, k = 2))
dat_shp_kn4 <- knn2nb(knearneigh(coords, k = 4))

 ###### calculate Moran's I and models results ###### 

  #GLMM with state as random effects variable
dat$model_GLMM_forest_state_residuals <- NA
dat$model_GLMM_state_residuals <- residuals(model_GLMM_state)
moran.test(dat$model_GLMM_state_residuals, dat_shp_kd2_w, randomisation=FALSE)

  #OLS
model_OLS<-lm(Expenditures~Female+Age+White+HouseDensityLog+EdIncome+UserConcentration+Doctors+Hospitals+HospitalBeds+Particulates+Inactive+Forest+Shrub+Grass+UrbanVeg,data=dat)
sjt.lm(model_OLS,show.aic=T, p.kr=F, emph.p=T, p.numeric=F)
dat$model_OLS_residuals <- NA
dat$model_OLS_residuals <- residuals(model_OLS)
moran.test(dat$model_OLS_residuals, dat_shp_kd2_w, randomisation=FALSE)

  #GLS with continuous autoccrelation structure
model_GLS<-gls(Expenditures~Female+Age+White+HouseDensityLog+EdIncome+UserConcentration+Doctors+Hospitals+HospitalBeds+Particulates+Inactive+Forest+Shrub+Grass+UrbanVeg,correlation = corCAR1(form = ~LON | LAT),data=dat)
summary(model_GLS)
coef(model_GLS)
dat$model_GLS_residuals <- NA
dat$model_GLS_residuals <- residuals(model_GLS)
moran.test(dat$model_GLS_residuals, dat_shp_kd2_w, randomisation=FALSE)
  
  #GLMM with ecoregion as random effects variable
model_GLMM_ecoregion <- lme4::lmer(Expenditures~Female+Age+White+HouseDensityLog+EdIncome+UserConcentration+Doctors+Hospitals+HospitalBeds+Particulates+Inactive+Forest+Shrub+Grass+UrbanVeg+
                                 (1 | Ecozone),
                               data=dat,
                               REML=F)
sjt.lmer(model_GLMM_ecoregion, show.aic=T, p.kr=F, emph.p=T, p.numeric=F, show.std = T)
dat$model_GLMM_ecoregion_residuals <- NA
dat$model_GLMM_ecoregion_residuals <- residuals(model_GLMM_ecoregion)
moran.test(dat$model_GLMM_ecoregion_residuals, dat_shp_kd2_w, randomisation=FALSE)

  #Spatial lag model
model_splag <- lagsarlm(Expend~Female+Age+White+HouseDensL+EdIncome+UserConc+Doctors+Hospitals+HospitalB+Partic+Inactive+Forest+Shrub+Grass+UrbanVeg,dat=dat_shp,listw=dat_shp_kd2_w)
summary(model_splag)
dat$model_splag_residuals <- NA
dat$model_splag_residuals <- residuals(model_splag)
moran.test(dat$model_splag_residuals, dat_shp_kd2_w, randomisation=FALSE)
```
******
```{run GLMM with other measures of Medicare expenditure}

  ###### Load data ###### 
  #need to retrieve data from github site, as it has been formatted to work in R: https://github.com/jjchern/ahrf#

#load data
load("/Volumes/GoogleDrive/My Drive/Projects/Tree$/Tree$_RA/Pilot_Study/Data/Variables Data/AHRF/ahrf_county.rda")
 
  ###### Actual Per Capita Costs from AHRF calculation ###### 
#retrieve 2011 "Actual Per Capita Mediare Cost" 
Expenditures_Actual_Per_Capita <- ahrf_county[, c("F00002", "F15297-11")]

#relabel variable and transform (AHRF includes cents but does not have decimal point)
Expenditures_Actual_Per_Capita$Expenditures_Actual_Per_Capita<-as.numeric(Expenditures_Actual_Per_Capita$`F15297-11`)
Expenditures_Actual_Per_Capita$Expenditures_Actual_Per_Capita <- Expenditures_Actual_Per_Capita$Expenditures_Actual_Per_Capita / 100
write.csv(Expenditures_Actual_Per_Capita, "AHRF_Expenditures_Actual_Per_Capita.csv")

#merge data
dat_Expenditures_Actual_Per_Capita_pre <- merge(dat,Expenditures_Actual_Per_Capita,by.x=c("FIPS"), by.y=c("F00002"))

#remove rows with missing data
dat_Expenditures_Actual_Per_Capita <- dat_Expenditures_Actual_Per_Capita_pre[complete.cases(dat_Expenditures_Actual_Per_Capita_pre),]

#Run GLMM with actual per capita
model_GLMM_state_actual_per_capita <- lme4::lmer(Expenditures_Actual_Per_Capita~Female+Age+White+HouseDensityLog+EdIncome+UserConcentration+Doctors+Hospitals+HospitalBeds+Particulates+Inactive+Forest+Shrub+Grass+UrbanVeg+(1 | State),data=dat_Expenditures_Actual_Per_Capita,REML=F)
summary(model_GLMM_state_actual_per_capita)

 ###### Actual Standardized Per Capita Costs from AHRF calculation ###### 

#retrieve 2011 "Standardized Per Capita Mediare Cost" 
Expenditures_Stand_Per_Capita <- ahrf_county[, c("F00002", "F15298-11")]

#relabel variable and transform (AHRF includes cents but does not have decimal point)
Expenditures_Stand_Per_Capita$Expenditures_Stand_Per_Capita<-as.numeric(Expenditures_Stand_Per_Capita$`F15298-11`)
Expenditures_Stand_Per_Capita$Expenditures_Stand_Per_Capita <- Expenditures_Stand_Per_Capita$Expenditures_Stand_Per_Capita / 100
write.csv(Expenditures_Stand_Per_Capita, "AHRF_Expenditures_Stand_Per_Capita.csv")

#merge data
dat_Expenditures_Stand_Per_Capita_pre <- merge(dat,Expenditures_Stand_Per_Capita,by.x=c("FIPS"), by.y=c("F00002"))

#remove rows with missing data
dat_Expenditures_Stand_Per_Capita <- dat_Expenditures_Stand_Per_Capita_pre[complete.cases(dat_Expenditures_Stand_Per_Capita_pre),]

#Run GLMM with actual per capita
model_GLMM_state_stand_per_capita <- lme4::lmer(Expenditures_Stand_Per_Capita~Female+Age+White+HouseDensityLog+EdIncome+UserConcentration+Doctors+Hospitals+HospitalBeds+Particulates+Inactive+Forest+Shrub+Grass+UrbanVeg+(1 | State),data=dat_Expenditures_Stand_Per_Capita,REML=F)
summary(model_GLMM_state_stand_per_capita)

  ######Compare models######

#Run original GLMM with standardized, risk-adjusted per capita expenditures measure
model_GLMM_state_stand_riskadj_per_capita <- lme4::lmer(Expenditures~Female+Age+White+HouseDensityLog+EdIncome+UserConcentration+
                                 Doctors+Hospitals+HospitalBeds+
                                 Particulates+Inactive+
                                 Forest+Shrub+Grass+UrbanVeg+
                                 (1 | State),
                               data=dat,
                               REML=F)

sjt.lmer(model_GLMM_state_stand_riskadj_per_capita, 
         model_GLMM_state_stand_per_capita, 
         model_GLMM_state_actual_per_capita, 
         depvar.labels = c("Standardized risk-adjusted","Standardized","Actual"),
         show.aic=T, p.kr=F, emph.p=T, p.numeric=F, show.std = T)

```
******
#Looking at differences in health expenditures by income-education index
```{r ANOVA for quartiles of income-education index}
#create quartile
dat_edincome_q <- within(dat, quartile <- as.integer(cut(EdIncome, quantile(EdIncome, probs=0:4/4), include.lowest=TRUE)))

#set quartile as factor
dat_edincome_q$quartile<-as.factor(dat_edincome_q$quartile)

#separate dataframes by educ income quartiles
dat_edincome_q1 <- dat[(dat_edincome_q$quartile == 1),]
dat_edincome_q2 <- dat[(dat_edincome_q$quartile == 2),]
dat_edincome_q3 <- dat[(dat_edincome_q$quartile == 3),]
dat_edincome_q4 <- dat[(dat_edincome_q$quartile == 4),]

#summary table
dat_q_edincome_sum <- summarySE(dat_edincome_q,measurevar="Expenditures",groupvars=c("quartile"))

#run ANOVA
ANOVA_edincome <- aov(Expenditures ~ quartile, data=dat_edincome_q)

#print results
summary(dat_q_edincome_sum)
summary(ANOVA_edincome)
TukeyHSD(ANOVA_edincome)
```
