Spatial GLM
Jairemy Edwards
10/18/2021
library(haven)
library(foreign)
library(readr)
library(dplyr)
library(ggplot2)
library(broom)
library(car)
library(MASS)
library(lmtest)
library(zoo)
library(nortest)
library(plotrix)
library(scales)
library(tableone)
library(Weighted.Desc.Stat)
library(mitools)
library(survey)
library(VGAM)
library(stargazer)
library(sandwich)
library(pastecs)
library(muhaz)
library(ggpubr)
library(survminer)
library(eha)
library(reshape2)
library(data.table)
library(magrittr)
library(tidyverse)
library(sjmisc)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(weights)
library(GGally)
library(tigris)
library(RColorBrewer)
library(patchwork)
library(tidycensus)
library(censusapi)
library(spdep)
library(classInt)
library(sf)
library(spatialreg)Data extraction, define outcome and predictors (AHRF 2018-2019)
arf<-read_sas("C://Users//Jaire//OneDrive//Desktop//Spatial Demography//Data//ahrf2019.sas7bdat")# Outcome/offset - suicides (2015-2017)
# Denominator - total deaths (2015-2017)
# Predictors:
# number of mental health facilities
# number of long-term psychiatric hospitals
# number of short-term psychiatric hospitals
# number of VA hospital admissionsset.seed(123)
drop_na(arf)
arf$cofips=arf$f00002
arf$coname=arf$f00010
arf$state = arf$f00008
arf$deaths1517=arf$f1194115
arf$deaths0810=arf$f1194108
arf$suicide1517=arf$f1316615
arf$suicide0810=arf$f1316608
arf$mentfac10=arf$f1322110
arf$psychlt=arf$f1067817
arf$psychst=arf$f1318517
arf$vaadm=arf$f0892017
scale(arf$mentfac10)
scale(arf$psychst)
scale(arf$psychlt)
scale(arf$vaadm)options(tigris_class="sf")
usco<-counties(cb=T, year=2010)
usco$cofips<-substr(usco$GEO_ID, 10, 15)
sts<-states(cb = T, year=2010)
sts<-st_boundary(sts)
arf.g<-geo_join(usco, arf, by_sp="cofips", by_df="cofips",how="left" )## Warning: We recommend using the dplyr::*_join() family of functions instead.## Warning: `group_by_()` was deprecated in dplyr 0.7.0.
## Please use `group_by()` instead.
## See vignette('programming') for more help
Poisson model and interpretation
arf.use<-filter(arf.g, is.na(suicide1517)==F)
fit_pois<- glm(suicide1517 ~ offset(log(deaths1517))+mentfac10+psychst+psychlt+vaadm,
family=poisson,
data=arf.use)
summary(fit_pois)##
## Call:
## glm(formula = suicide1517 ~ offset(log(deaths1517)) + mentfac10 +
## psychst + psychlt + vaadm, family = poisson, data = arf.use)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -8.4179 -0.8057 0.2048 1.0687 10.5588
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.056e+00 6.374e-03 -636.424 < 2e-16 ***
## mentfac10 -2.989e-03 7.747e-04 -3.859 0.000114 ***
## psychst 2.414e-02 2.943e-03 8.202 2.36e-16 ***
## psychlt -1.075e-01 1.251e-02 -8.597 < 2e-16 ***
## vaadm -1.382e-05 1.363e-06 -10.142 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 3602.6 on 928 degrees of freedom
## Residual deviance: 3374.5 on 924 degrees of freedom
## AIC: 8154.2
##
## Number of Fisher Scoring iterations: 4Check for overdispersion
scale<-sqrt(fit_pois$deviance/fit_pois$df.residual)
scale## [1] 1.911023exp(coef(fit_pois))## (Intercept) mentfac10 psychst psychlt vaadm
## 0.01731308 0.99701520 1.02443047 0.89804966 0.99998618The model shows that an increase in the number of mental health facilities in an area is associated with .003% decrease in the mean number of suicides (p<0.05). An increase in the number of short term psychiatric facilities in an area is associated with an 2.4% increase in the mean number of suicides (p<0.05). An increase in the number of long term psychiatric facilities in an area is associated with an 11% decrease in the mean number of suicides (p<0.05). An increase in the number of admissions to VA hospitals in an area is associated with an .00002% decrease in the mean number of suicides (p<0.05). Over-dispersion is present.
Test for residual autocorrelation
nbs<-knearneigh(coordinates(as(arf.use, "Spatial")), k = 4, longlat = T)
nbs<-knn2nb(nbs, sym = T)
us.wt4<-nb2listw(nbs, style="W", zero.policy = TRUE)
length(nbs)## [1] 929can.be.simmed(us.wt4)## [1] TRUElength(us.wt4)## [1] 3lm.morantest(fit_pois, listw = us.wt4, zero.policy = TRUE)##
## Global Moran I for regression residuals
##
## data:
## model: glm(formula = suicide1517 ~ offset(log(deaths1517)) + mentfac10
## + psychst + psychlt + vaadm, family = poisson, data = arf.use)
## weights: us.wt4
##
## Moran I statistic standard deviate = 25.882, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 5.472563e-01 -3.278262e-05 4.471312e-04There is moderate positive autocorrelation among residuals.
Spatial filtering analysis of Eigenvectors
# nb model
fit_nb<-glm.nb(suicide1517 ~ offset(log(deaths1517))+mentfac10+psychlt+psychst+vaadm,
data=arf.use)
summary(fit_nb)##
## Call:
## glm.nb(formula = suicide1517 ~ offset(log(deaths1517)) + mentfac10 +
## psychlt + psychst + vaadm, data = arf.use, init.theta = 17.214375,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.9671 -0.6635 -0.0325 0.5086 3.8584
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.994e+00 1.196e-02 -333.951 < 2e-16 ***
## mentfac10 -5.466e-03 3.031e-03 -1.804 0.0713 .
## psychlt -1.145e-01 2.825e-02 -4.053 5.06e-05 ***
## psychst 9.630e-03 1.016e-02 0.948 0.3433
## vaadm -1.033e-05 4.049e-06 -2.552 0.0107 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(17.2144) family taken to be 1)
##
## Null deviance: 903.26 on 928 degrees of freedom
## Residual deviance: 871.72 on 924 degrees of freedom
## AIC: 6594.8
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 17.21
## Std. Err.: 1.19
##
## 2 x log-likelihood: -6582.77fit_nb$theta## [1] 17.21438me.fit<-ME(suicide1517 ~ offset(log(deaths1517))+mentfac10+psychlt+psychst+vaadm,
data=arf.use,
family=negative.binomial(17.21438),
listw = us.wt4,
verbose=T,
alpha=.05)## eV[,2], I: 0.3085447 ZI: NA, pr(ZI): 0.01
## eV[,34], I: 0.2824097 ZI: NA, pr(ZI): 0.01
## eV[,13], I: 0.2555259 ZI: NA, pr(ZI): 0.01
## eV[,18], I: 0.2313331 ZI: NA, pr(ZI): 0.01
## eV[,9], I: 0.2164284 ZI: NA, pr(ZI): 0.01
## eV[,10], I: 0.1848292 ZI: NA, pr(ZI): 0.01
## eV[,3], I: 0.1619137 ZI: NA, pr(ZI): 0.01
## eV[,90], I: 0.1418652 ZI: NA, pr(ZI): 0.01
## eV[,93], I: 0.1290831 ZI: NA, pr(ZI): 0.01
## eV[,55], I: 0.1161711 ZI: NA, pr(ZI): 0.01
## eV[,63], I: 0.1035144 ZI: NA, pr(ZI): 0.01
## eV[,24], I: 0.09463813 ZI: NA, pr(ZI): 0.01
## eV[,807], I: 0.08673058 ZI: NA, pr(ZI): 0.01
## eV[,1], I: 0.07940913 ZI: NA, pr(ZI): 0.01
## eV[,25], I: 0.07290248 ZI: NA, pr(ZI): 0.01
## eV[,65], I: 0.06647661 ZI: NA, pr(ZI): 0.01
## eV[,426], I: 0.06065521 ZI: NA, pr(ZI): 0.01
## eV[,8], I: 0.0543628 ZI: NA, pr(ZI): 0.03
## eV[,140], I: 0.04758679 ZI: NA, pr(ZI): 0.02
## eV[,618], I: 0.04095998 ZI: NA, pr(ZI): 0.07me.fit## Eigenvector ZI pr(ZI)
## 0 NA NA 0.01
## 1 2 NA 0.01
## 2 34 NA 0.01
## 3 13 NA 0.01
## 4 18 NA 0.01
## 5 9 NA 0.01
## 6 10 NA 0.01
## 7 3 NA 0.01
## 8 90 NA 0.01
## 9 93 NA 0.01
## 10 55 NA 0.01
## 11 63 NA 0.01
## 12 24 NA 0.01
## 13 807 NA 0.01
## 14 1 NA 0.01
## 15 25 NA 0.01
## 16 65 NA 0.01
## 17 426 NA 0.01
## 18 8 NA 0.03
## 19 140 NA 0.02
## 20 618 NA 0.07fits<-data.frame(fitted(me.fit))
arf.use$me1<-fits$vec1
arf.use$me13<-fits$vec13
Cleaned negative binomial:
clean.nb<-glm.nb(suicide1517 ~ offset(log(deaths1517)) + mentfac10 +psychlt+psychst+vaadm+fitted(me.fit), arf.use)
summary(clean.nb)##
## Call:
## glm.nb(formula = suicide1517 ~ offset(log(deaths1517)) + mentfac10 +
## psychlt + psychst + vaadm + fitted(me.fit), data = arf.use,
## init.theta = 63.76636108, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.3386 -0.5681 0.0034 0.6201 3.3779
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.024e+00 8.764e-03 -459.118 < 2e-16 ***
## mentfac10 -4.441e-03 1.773e-03 -2.505 0.012243 *
## psychlt -4.241e-02 1.941e-02 -2.185 0.028911 *
## psychst 1.346e-02 6.431e-03 2.093 0.036343 *
## vaadm -1.098e-05 2.593e-06 -4.234 2.30e-05 ***
## fitted(me.fit)vec2 4.865e+00 2.186e-01 22.256 < 2e-16 ***
## fitted(me.fit)vec34 5.403e-01 1.990e-01 2.715 0.006632 **
## fitted(me.fit)vec13 -1.677e+00 2.316e-01 -7.244 4.37e-13 ***
## fitted(me.fit)vec18 -1.198e+00 2.228e-01 -5.376 7.62e-08 ***
## fitted(me.fit)vec9 7.850e-01 2.121e-01 3.701 0.000215 ***
## fitted(me.fit)vec10 -2.679e+00 2.165e-01 -12.375 < 2e-16 ***
## fitted(me.fit)vec3 1.431e+00 2.338e-01 6.122 9.22e-10 ***
## fitted(me.fit)vec90 6.179e-01 1.977e-01 3.125 0.001780 **
## fitted(me.fit)vec93 2.360e-01 2.065e-01 1.143 0.253182
## fitted(me.fit)vec55 -9.074e-01 2.280e-01 -3.980 6.89e-05 ***
## fitted(me.fit)vec63 -9.369e-01 2.192e-01 -4.275 1.92e-05 ***
## fitted(me.fit)vec24 -5.864e-01 2.163e-01 -2.711 0.006712 **
## fitted(me.fit)vec807 -2.346e-01 2.087e-01 -1.124 0.260847
## fitted(me.fit)vec1 7.695e-01 2.122e-01 3.627 0.000287 ***
## fitted(me.fit)vec25 -6.729e-01 2.163e-01 -3.112 0.001861 **
## fitted(me.fit)vec65 5.276e-01 2.237e-01 2.359 0.018346 *
## fitted(me.fit)vec426 2.674e-01 2.222e-01 1.203 0.228945
## fitted(me.fit)vec8 -7.748e-01 2.248e-01 -3.447 0.000568 ***
## fitted(me.fit)vec140 -2.949e-01 2.021e-01 -1.459 0.144621
## fitted(me.fit)vec618 1.565e-01 1.970e-01 0.794 0.427081
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(63.7664) family taken to be 1)
##
## Null deviance: 1793.31 on 928 degrees of freedom
## Residual deviance: 824.78 on 904 degrees of freedom
## AIC: 6037.9
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 63.77
## Std. Err.: 7.53
##
## 2 x log-likelihood: -5985.867lm.morantest(fit_nb, listw=us.wt4)##
## Global Moran I for regression residuals
##
## data:
## model: glm.nb(formula = suicide1517 ~ offset(log(deaths1517)) +
## mentfac10 + psychlt + psychst + vaadm, data = arf.use, init.theta =
## 17.214375, link = log)
## weights: us.wt4
##
## Moran I statistic standard deviate = 27.351, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.5778601597 -0.0001166645 0.0004465395lm.morantest(clean.nb, listw=us.wt4)##
## Global Moran I for regression residuals
##
## data:
## model: glm.nb(formula = suicide1517 ~ offset(log(deaths1517)) +
## mentfac10 + psychlt + psychst + vaadm + fitted(me.fit), data = arf.use,
## init.theta = 63.76636108, link = log)
## weights: us.wt4
##
## Moran I statistic standard deviate = 9.274, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.1988180487 -0.0008419133 0.0004634955The filtering approach filtered out a large amount of spatial autocorrelation in the negative binomial model.