Updated CARBayes
Fin Mooney
2023-03-13
Scaling the variables
EJScreen <- st_read("/Users/finmooney/Downloads/EJScreenFinal/EJScreen032022.shp")## Reading layer `EJScreen032022' from data source
## `/Users/finmooney/Downloads/EJScreenFinal/EJScreen032022.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 2190 features and 56 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -8212581 ymin: 4949772 xmax: -7998050 ymax: 5044227
## Projected CRS: WGS 84 / Pseudo-MercatorEJScreen$PTRAF_3 <- ifelse(EJScreen$PTRAF_3 == 1, 3, ifelse(EJScreen$PTRAF_3 == 3, 1, EJScreen$PTRAF_3))
EJScreen$PRMP_3 <- ifelse(EJScreen$PRMP_3 == 1, 3, ifelse(EJScreen$PRMP_3 == 3, 1, EJScreen$PRMP_3))
EJScreen$PTSDF_3 <- ifelse(EJScreen$PTSDF_3 == 1, 3, ifelse(EJScreen$PTSDF_3 == 3, 1, EJScreen$PTSDF_3))
EJScreen$UST_3 <- ifelse(EJScreen$UST_3 == 1, 3, ifelse(EJScreen$UST_3 == 3, 1, EJScreen$UST_3))
EJScreen_transform = EJScreen %>%
mutate(
PNPL_Log = log(PNPL),
MINORPCT_S = scale(MINORPCT), # 1 SD unit increase
LOWINCPCT_S = scale(LOWINCPCT),
LESSHSPCT_S = scale(LESSHSPCT),
LINGISOPCT_S = scale(LINGISOPCT),
UNDER5PCT_S = scale(UNDER5PCT),
OVER64PCT_S = scale(OVER64PCT),
PRE1960PCT_S = scale(PRE1960PCT),
PTRAF_F = factor(PTRAF_3),
PRMP_F = factor(PRMP_3),
PTSDF_F = factor(PTSDF_3),
PWDIS_F = factor(PWDIS_3),
UST_F = factor(UST_3),
OZONE_S = scale(OZONE),
PM25_S = scale(PM25),
PNPL_MEDIAN = median(PNPL),
PNPL_ABOVE_MEDIAN = ifelse(PNPL > PNPL_MEDIAN, 1, 0)
)Creating Adjacency
# Create neighborhood structure with queen method
queenNB <- poly2nb(EJScreen_transform, queen = TRUE) # Queen adjacency of irregular lattice
coords <- st_coordinates(EJScreen_transform) # Obtain the centroids of the features
ids <- row.names(sf::st_set_geometry(EJScreen_transform, NULL))
# Convert queen neighbor list to a binary regular matrix object
W_mat <- nb2mat(queenNB, zero.policy = TRUE, style = "B")
# Show Queen Adjacency
plot(queenNB, st_geometry(EJScreen_transform), pch=20, cex=0.5, col="gray"); title("Queen Adjacency")
CARLeroux without Rho = 1
leroux_dem <- S.CARleroux(formula = PNPL_LOGAR ~ MINORPCT_S + LOWINCPCT_S + LESSHSPCT_S +
LINGISOPCT_S + OVER64PCT_S,
data = EJScreen_transform,
family = "gaussian",
W = W_mat,
MALA = FALSE,
burnin = 50000,
n.sample = 300000,
thin = 10,
prior.mean.beta = rep(0, times = 6),
prior.var.beta = c(rep((100^2), times = 6)),
prior.nu2 = c(0.01, 0.01),
prior.tau2 = c(0.01, 0.01),
verbose = TRUE)## Setting up the model.## Warning in mat2listw(W): style is M (missing); style should be set to a valid
## value## Generating 25000 post burnin and thinned (if requested) samples.
##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|=== | 4%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 24%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 90%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 100%
## Summarising results.
## Finished in 149.1 seconds.# Create a data frame to store the results
result <- data.frame(variable = character(),
mean = double(),
ci_lower = double(),
ci_upper = double(),
stringsAsFactors = FALSE)
# Exponentiate coefficients for MINORPCT
exp.coef_MINORPCT <- exp(1 * leroux_dem$samples$beta[, 2] / sd(EJScreen_transform$MINORPCT))
# Compute mean and 95% confidence interval
exp.coef.mean <- mean(exp.coef_MINORPCT)
exp.coef.ci <- quantile(exp.coef_MINORPCT, c(0.025, 0.975))
# Add results to data frame
result <- rbind(result, data.frame(variable = "MINORPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
# Repeat for each variable
exp.coef_LOWINCPCT <- exp(1 * leroux_dem$samples$beta[, 3] / sd(EJScreen_transform$LOWINCPCT))
exp.coef.mean <- mean(exp.coef_LOWINCPCT)
exp.coef.ci <- quantile(exp.coef_LOWINCPCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "LOWINCPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_LESSHSPCT <- exp(1 * leroux_dem$samples$beta[, 4] / sd(EJScreen_transform$LESSHSPCT))
exp.coef.mean <- mean(exp.coef_LESSHSPCT)
exp.coef.ci <- quantile(exp.coef_LESSHSPCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "LESSHSPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_LINGISOPCT <- exp(1 * leroux_dem$samples$beta[, 5] / sd(EJScreen_transform$LINGISOPCT))
exp.coef.mean <- mean(exp.coef_LINGISOPCT)
exp.coef.ci <- quantile(exp.coef_LINGISOPCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "LINGISOPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_OVER64PCT <- exp(1 * leroux_dem$samples$beta[, 6] / sd(EJScreen_transform$OVER64PCT))
exp.coef.mean <- mean(exp.coef_OVER64PCT)
exp.coef.ci <- quantile(exp.coef_OVER64PCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "OVER64PCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
print(result)## variable mean ci_lower ci_upper
## 2.5% MINORPCT 1.0315913 0.9637841 1.103836
## 2.5%1 LOWINCPCT 1.0366937 0.9481906 1.133225
## 2.5%2 LESSHSPCT 0.9858920 0.8447280 1.151514
## 2.5%3 LINGISOPCT 1.0404753 0.8661624 1.238192
## 2.5%4 OVER64PCT 0.9885946 0.8895695 1.094294# Fit Gaussian CAR model
leroux <- S.CARleroux(formula = PNPL_LOGAR ~ MINORPCT_S + LOWINCPCT_S + LESSHSPCT_S +
LINGISOPCT_S + OVER64PCT_S
+ PRE1960PCT_S + PTRAF_F
+ PWDIS_F + PRMP_F + PTSDF_F + UST_F + OZONE_S + PM25_S,
data = EJScreen_transform,
family = "gaussian",
W = W_mat,
MALA = FALSE,
burnin = 50000,
n.sample = 300000,
thin = 10,
prior.mean.beta = rep(0, times = 19),
prior.var.beta = c(rep((100^2), times = 19)),
prior.nu2 = c(0.01, 0.01),
prior.tau2 = c(0.01, 0.01),
verbose = TRUE)## Setting up the model.## Warning in mat2listw(W): style is M (missing); style should be set to a valid
## value## Generating 25000 post burnin and thinned (if requested) samples.
##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|=== | 4%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 24%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 90%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 100%
## Summarising results.
## Finished in 184.5 seconds.leroux$summary.results## Mean 2.5% 97.5% n.sample % accept n.effective Geweke.diag
## (Intercept) -1.6507 -1.7430 -1.5588 25000 100.0 396.0 -2.5
## MINORPCT_S 0.0023 -0.0156 0.0201 25000 100.0 1437.7 0.5
## LOWINCPCT_S 0.0054 -0.0060 0.0169 25000 100.0 2493.5 0.1
## LESSHSPCT_S -0.0033 -0.0176 0.0110 25000 100.0 2104.8 0.0
## LINGISOPCT_S -0.0004 -0.0123 0.0115 25000 100.0 2414.6 -1.6
## OVER64PCT_S -0.0013 -0.0121 0.0096 25000 100.0 2116.6 0.8
## PRE1960PCT_S -0.0042 -0.0205 0.0119 25000 100.0 1484.7 0.9
## PTRAF_F2 0.0050 -0.0228 0.0327 25000 100.0 2077.1 0.9
## PTRAF_F3 0.0062 -0.0270 0.0401 25000 100.0 1674.8 0.8
## PWDIS_F2 0.0896 0.0147 0.1622 25000 100.0 736.5 0.3
## PWDIS_F3 0.0558 -0.0175 0.1301 25000 100.0 423.2 0.0
## PRMP_F2 0.1769 0.1036 0.2502 25000 100.0 377.7 3.3
## PRMP_F3 0.3367 0.2431 0.4342 25000 100.0 282.0 2.8
## PTSDF_F2 0.0908 0.0477 0.1341 25000 100.0 750.4 1.3
## PTSDF_F3 0.1942 0.1354 0.2520 25000 100.0 541.6 0.5
## UST_F2 0.0061 -0.0194 0.0312 25000 100.0 2071.4 0.6
## UST_F3 0.0068 -0.0227 0.0357 25000 100.0 1850.9 0.5
## OZONE_S -0.1662 -0.2371 -0.0876 25000 100.0 66.5 -1.2
## PM25_S 0.5367 0.3485 0.6860 25000 100.0 19.6 0.9
## nu2 0.0011 0.0007 0.0016 25000 100.0 2374.5 -0.4
## tau2 0.3039 0.2860 0.3227 25000 100.0 25000.0 -2.4
## rho 0.9989 0.9967 0.9999 25000 41.8 22369.7 -0.9# Create a data frame to store the results
result <- data.frame(variable = character(),
mean = double(),
ci_lower = double(),
ci_upper = double(),
stringsAsFactors = FALSE)
# Exponentiate coefficients for MINORPCT
exp.coef_MINORPCT <- exp(1 * leroux$samples$beta[, 2] / sd(EJScreen_transform$MINORPCT))
# Compute mean and 95% confidence interval
exp.coef.mean <- mean(exp.coef_MINORPCT)
exp.coef.ci <- quantile(exp.coef_MINORPCT, c(0.025, 0.975))
# Add results to data frame
result <- rbind(result, data.frame(variable = "MINORPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
# Repeat for each variable
exp.coef_LOWINCPCT <- exp(1 * leroux$samples$beta[, 3] / sd(EJScreen_transform$LOWINCPCT))
exp.coef.mean <- mean(exp.coef_LOWINCPCT)
exp.coef.ci <- quantile(exp.coef_LOWINCPCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "LOWINCPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_LESSHSPCT <- exp(1 * leroux$samples$beta[, 4] / sd(EJScreen_transform$LESSHSPCT))
exp.coef.mean <- mean(exp.coef_LESSHSPCT)
exp.coef.ci <- quantile(exp.coef_LESSHSPCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "LESSHSPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_LINGISOPCT <- exp(1 * leroux$samples$beta[, 5] / sd(EJScreen_transform$LINGISOPCT))
exp.coef.mean <- mean(exp.coef_LINGISOPCT)
exp.coef.ci <- quantile(exp.coef_LINGISOPCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "LINGISOPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_OVER64PCT <- exp(1 * leroux$samples$beta[, 6] / sd(EJScreen_transform$OVER64PCT))
exp.coef.mean <- mean(exp.coef_OVER64PCT)
exp.coef.ci <- quantile(exp.coef_OVER64PCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "OVER64PCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PRE1960PCT <- exp(1 * leroux$samples$beta[, 7] / sd(EJScreen_transform$PRE1960PCT))
exp.coef.mean <- mean(exp.coef_PRE1960PCT)
exp.coef.ci <- quantile(exp.coef_PRE1960PCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PRE1960PCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PTRAF <- exp(1 * leroux$samples$beta[, 8] / sd(EJScreen_transform$PTRAF))
exp.coef.mean <- mean(exp.coef_PTRAF)
exp.coef.ci <- quantile(exp.coef_PTRAF, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PTRAF",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PTRAF2 <- exp(1 * leroux$samples$beta[, 9] / sd(EJScreen_transform$PTRAF))
exp.coef.mean <- mean(exp.coef_PTRAF2)
exp.coef.ci <- quantile(exp.coef_PTRAF2, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PTRAF2",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PWDIS <- exp(1 * leroux$samples$beta[, 10] / sd(EJScreen_transform$PWDIS))
exp.coef.mean <- mean(exp.coef_PWDIS)
exp.coef.ci <- quantile(exp.coef_PWDIS, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PWDIS",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PWDIS2 <- exp(1 * leroux$samples$beta[, 11] / sd(EJScreen_transform$PWDIS))
exp.coef.mean <- mean(exp.coef_PWDIS2)
exp.coef.ci <- quantile(exp.coef_PWDIS2, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PWDIS2",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PRMP <- exp(1 * leroux$samples$beta[, 12] / sd(EJScreen_transform$PRMP))
exp.coef.mean <- mean(exp.coef_PRMP)
exp.coef.ci <- quantile(exp.coef_PRMP, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PRMP",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PRMP2 <- exp(1 * leroux$samples$beta[, 13] / sd(EJScreen_transform$PRMP))
exp.coef.mean <- mean(exp.coef_PRMP2)
exp.coef.ci <- quantile(exp.coef_PRMP2, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PRMP2",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PTSDF <- exp(1 * leroux$samples$beta[, 14] / sd(EJScreen_transform$PTSDF))
exp.coef.mean <- mean(exp.coef_PTSDF)
exp.coef.ci <- quantile(exp.coef_PTSDF, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PTSDF",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PTSDF2 <- exp(1 * leroux$samples$beta[, 15] / sd(EJScreen_transform$PTSDF))
exp.coef.mean <- mean(exp.coef_PTSDF2)
exp.coef.ci <- quantile(exp.coef_PTSDF2, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PTSDF2",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_USTF <- exp(1 * leroux$samples$beta[, 16] / sd(EJScreen_transform$UST))
exp.coef.mean <- mean(exp.coef_USTF)
exp.coef.ci <- quantile(exp.coef_USTF, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "UST",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_USTF2 <- exp(1 * leroux$samples$beta[, 17] / sd(EJScreen_transform$UST))
exp.coef.mean <- mean(exp.coef_USTF2)
exp.coef.ci <- quantile(exp.coef_USTF2, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "UST",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PM25 <- exp(1 * leroux$samples$beta[, 18] / sd(EJScreen_transform$PM25))
exp.coef.mean <- mean(exp.coef_PM25)
exp.coef.ci <- quantile(exp.coef_PM25, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "Ozone",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_Ozone <- exp(1 * leroux$samples$beta[, 19] / sd(EJScreen_transform$OZONE))
exp.coef.mean <- mean(exp.coef_Ozone)
exp.coef.ci <- quantile(exp.coef_Ozone, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PM 2.5",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
print(result)## variable mean ci_lower ci_upper
## 2.5% MINORPCT 1.0089514 0.9451837 1.0754652
## 2.5%1 LOWINCPCT 1.0426086 0.9556793 1.1360550
## 2.5%2 LESSHSPCT 0.9691095 0.8308569 1.1228248
## 2.5%3 LINGISOPCT 0.9984320 0.8347226 1.1838013
## 2.5%4 OVER64PCT 0.9900260 0.8974896 1.0891523
## 2.5%5 PRE1960PCT 0.9860751 0.9321695 1.0416812
## 2.5%6 PTRAF 1.0000072 0.9999667 1.0000478
## 2.5%7 PTRAF2 1.0000090 0.9999606 1.0000584
## 2.5%8 PWDIS 1.6280474 1.0797091 2.3333934
## 2.5%9 PWDIS2 1.3646808 0.9125207 1.9728830
## 2.5%10 PRMP 2.4080213 1.6576169 3.3867960
## 2.5%11 PRMP2 5.3122943 3.2726468 8.3110542
## 2.5%12 PTSDF 1.0537806 1.0278617 1.0803238
## 2.5%13 PTSDF2 1.1185155 1.0810914 1.1562556
## 2.5%14 UST 1.0030912 0.9902767 1.0158915
## 2.5%15 UST 1.0034745 0.9886113 1.0181745
## 2.5%16 Ozone 0.7016788 0.6002532 0.8281719
## 2.5%17 PM 2.5 2.4443675 1.7751137 3.0945803rm(result)
residuals <- resid(leroux)
hist(residuals, main="Histogram of Residuals", xlab="Residuals")
print(leroux)##
## #################
## #### Model fitted
## #################
## Likelihood model - Gaussian (identity link function)
## Random effects model - Leroux CAR
## Regression equation - PNPL_LOGAR ~ MINORPCT_S + LOWINCPCT_S + LESSHSPCT_S + LINGISOPCT_S +
## OVER64PCT_S + PRE1960PCT_S + PTRAF_F + PWDIS_F + PRMP_F +
## PTSDF_F + UST_F + OZONE_S + PM25_S
## Number of missing observations - 0
##
## ############
## #### Results
## ############
## Posterior quantities and DIC
##
## Mean 2.5% 97.5% n.effective Geweke.diag
## (Intercept) -1.6507 -1.7430 -1.5588 396.0 -2.5
## MINORPCT_S 0.0023 -0.0156 0.0201 1437.7 0.5
## LOWINCPCT_S 0.0054 -0.0060 0.0169 2493.5 0.1
## LESSHSPCT_S -0.0033 -0.0176 0.0110 2104.8 0.0
## LINGISOPCT_S -0.0004 -0.0123 0.0115 2414.6 -1.6
## OVER64PCT_S -0.0013 -0.0121 0.0096 2116.6 0.8
## PRE1960PCT_S -0.0042 -0.0205 0.0119 1484.7 0.9
## PTRAF_F2 0.0050 -0.0228 0.0327 2077.1 0.9
## PTRAF_F3 0.0062 -0.0270 0.0401 1674.8 0.8
## PWDIS_F2 0.0896 0.0147 0.1622 736.5 0.3
## PWDIS_F3 0.0558 -0.0175 0.1301 423.2 0.0
## PRMP_F2 0.1769 0.1036 0.2502 377.7 3.3
## PRMP_F3 0.3367 0.2431 0.4342 282.0 2.8
## PTSDF_F2 0.0908 0.0477 0.1341 750.4 1.3
## PTSDF_F3 0.1942 0.1354 0.2520 541.6 0.5
## UST_F2 0.0061 -0.0194 0.0312 2071.4 0.6
## UST_F3 0.0068 -0.0227 0.0357 1850.9 0.5
## OZONE_S -0.1662 -0.2371 -0.0876 66.5 -1.2
## PM25_S 0.5367 0.3485 0.6860 19.6 0.9
## nu2 0.0011 0.0007 0.0016 2374.5 -0.4
## tau2 0.3039 0.2860 0.3227 25000.0 -2.4
## rho 0.9989 0.9967 0.9999 22369.7 -0.9
##
## DIC = -6785.826 p.d = 2088.8 LMPL = 2558.77Stratify by County
# Nassau
NassauEJ <- subset(EJScreen_transform, CNTY_NAME == "Nassau County")
# Create neighborhood structure with queen method
queenNB2 <- poly2nb(NassauEJ, queen = TRUE) # Queen adjacency of irregular lattice
coords2 <- st_coordinates(NassauEJ) # Obtain the centroids of the features
ids2 <- row.names(sf::st_set_geometry(NassauEJ, NULL))
# Convert queen neighbor list to a binary regular matrix object
W_mat2 <- nb2mat(queenNB2, zero.policy = TRUE, style = "B")
# Show Queen Adjacency
plot(queenNB2, st_geometry(NassauEJ), pch=20, cex=0.5, col="gray"); title("Queen Adjacency for Nassau County")
# Fit Gaussian CAR Model for Nassau County
leroux_nassau <- S.CARleroux(formula = PNPL_LOGAR ~ MINORPCT_S + LOWINCPCT_S + LESSHSPCT_S +
LINGISOPCT_S + OVER64PCT_S
+ PRE1960PCT_S + PTRAF_F + PWDIS_F + PRMP_F + PTSDF_F + UST_F + OZONE_S + PM25_S,
data = NassauEJ,
family = "gaussian",
W = W_mat2,
MALA = FALSE,
burnin = 50000,
n.sample = 300000,
thin = 10,
prior.mean.beta = rep(0, times = 19),
prior.var.beta = c(rep((100^2), times = 19)),
prior.nu2 = c(0.01, 0.01),
prior.tau2 = c(0.01, 0.01),
verbose = TRUE)## Setting up the model.## Warning in mat2listw(W): style is M (missing); style should be set to a valid
## value## Generating 25000 post burnin and thinned (if requested) samples.
##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|=== | 4%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 24%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 90%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 100%
## Summarising results.
## Finished in 97.8 seconds.print(leroux_nassau)##
## #################
## #### Model fitted
## #################
## Likelihood model - Gaussian (identity link function)
## Random effects model - Leroux CAR
## Regression equation - PNPL_LOGAR ~ MINORPCT_S + LOWINCPCT_S + LESSHSPCT_S + LINGISOPCT_S +
## OVER64PCT_S + PRE1960PCT_S + PTRAF_F + PWDIS_F + PRMP_F +
## PTSDF_F + UST_F + OZONE_S + PM25_S
## Number of missing observations - 0
##
## ############
## #### Results
## ############
## Posterior quantities and DIC
##
## Mean 2.5% 97.5% n.effective Geweke.diag
## (Intercept) -1.5205 -1.7285 -1.3021 344.4 -0.7
## MINORPCT_S -0.0098 -0.0347 0.0149 1616.7 0.3
## LOWINCPCT_S 0.0084 -0.0087 0.0252 3264.6 -0.4
## LESSHSPCT_S -0.0038 -0.0237 0.0161 2946.9 1.8
## LINGISOPCT_S 0.0000 -0.0143 0.0143 3383.7 1.6
## OVER64PCT_S 0.0038 -0.0139 0.0208 3582.8 1.0
## PRE1960PCT_S 0.0132 -0.0091 0.0357 2690.2 1.5
## PTRAF_F2 0.0195 -0.0230 0.0632 2342.8 -0.2
## PTRAF_F3 0.0236 -0.0258 0.0715 1847.1 -0.1
## PWDIS_F2 0.0268 -0.0705 0.1265 780.0 -0.3
## PWDIS_F3 0.0730 -0.0166 0.1650 673.8 -0.4
## PRMP_F2 0.0886 -0.0678 0.2394 790.4 1.5
## PRMP_F3 0.2598 0.0846 0.4249 760.9 1.9
## PTSDF_F2 0.0707 0.0049 0.1346 1186.2 0.5
## PTSDF_F3 0.1268 0.0428 0.2105 725.0 1.1
## UST_F2 0.0193 -0.0201 0.0593 2893.8 0.6
## UST_F3 0.0136 -0.0261 0.0531 2063.1 1.2
## OZONE_S -0.0837 -0.2337 0.0468 55.6 -0.9
## PM25_S 0.3967 0.1995 0.5977 63.8 -0.4
## nu2 0.0014 0.0008 0.0023 3277.9 1.6
## tau2 0.2927 0.2694 0.3186 21935.6 -0.7
## rho 0.9984 0.9951 0.9999 15872.7 2.6
##
## DIC = -3195.131 p.d = 1062.267 LMPL = 1181.45nassau_residuals <- resid(leroux_nassau)
hist(nassau_residuals, main="Histogram of Nassau County Residuals", xlab="Residuals")
#exponentiate
# Create a data frame to store the results
result <- data.frame(variable = character(),
mean = double(),
ci_lower = double(),
ci_upper = double(),
stringsAsFactors = FALSE)
# Exponentiate coefficients for MINORPCT
exp.coef_MINORPCT <- exp(1 * leroux_nassau$samples$beta[, 2] / sd(NassauEJ$MINORPCT))
# Compute mean and 95% confidence interval
exp.coef.mean <- mean(exp.coef_MINORPCT)
exp.coef.ci <- quantile(exp.coef_MINORPCT, c(0.025, 0.975))
# Add results to data frame
result <- rbind(result, data.frame(variable = "MINORPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
# Repeat for each variable
exp.coef_LOWINCPCT <- exp(1 * leroux_nassau$samples$beta[, 3] / sd(EJScreen_transform$LOWINCPCT))
exp.coef.mean <- mean(exp.coef_LOWINCPCT)
exp.coef.ci <- quantile(exp.coef_LOWINCPCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "LOWINCPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_LESSHSPCT <- exp(1 * leroux_nassau$samples$beta[, 4] / sd(EJScreen_transform$LESSHSPCT))
exp.coef.mean <- mean(exp.coef_LESSHSPCT)
exp.coef.ci <- quantile(exp.coef_LESSHSPCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "LESSHSPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_LINGISOPCT <- exp(1 * leroux_nassau$samples$beta[, 5] / sd(EJScreen_transform$LINGISOPCT))
exp.coef.mean <- mean(exp.coef_LINGISOPCT)
exp.coef.ci <- quantile(exp.coef_LINGISOPCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "LINGISOPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_OVER64PCT <- exp(1 * leroux_nassau$samples$beta[, 6] / sd(EJScreen_transform$OVER64PCT))
exp.coef.mean <- mean(exp.coef_OVER64PCT)
exp.coef.ci <- quantile(exp.coef_OVER64PCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "OVER64PCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PRE1960PCT <- exp(1 * leroux_nassau$samples$beta[, 7] / sd(EJScreen_transform$PRE1960PCT))
exp.coef.mean <- mean(exp.coef_PRE1960PCT)
exp.coef.ci <- quantile(exp.coef_PRE1960PCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PRE1960PCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PTRAF <- exp(1 * leroux_nassau$samples$beta[, 8] / sd(EJScreen_transform$PTRAF))
exp.coef.mean <- mean(exp.coef_PTRAF)
exp.coef.ci <- quantile(exp.coef_PTRAF, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PTRAF",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PTRAF2 <- exp(1 * leroux_nassau$samples$beta[, 9] / sd(EJScreen_transform$PTRAF))
exp.coef.mean <- mean(exp.coef_PTRAF2)
exp.coef.ci <- quantile(exp.coef_PTRAF2, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PTRAF2",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PWDIS <- exp(1 * leroux_nassau$samples$beta[, 10] / sd(EJScreen_transform$PWDIS))
exp.coef.mean <- mean(exp.coef_PWDIS)
exp.coef.ci <- quantile(exp.coef_PWDIS, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PWDIS",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PWDIS2 <- exp(1 * leroux_nassau$samples$beta[, 11] / sd(EJScreen_transform$PWDIS))
exp.coef.mean <- mean(exp.coef_PWDIS2)
exp.coef.ci <- quantile(exp.coef_PWDIS2, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PWDIS2",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PRMP <- exp(1 * leroux_nassau$samples$beta[, 12] / sd(EJScreen_transform$PRMP))
exp.coef.mean <- mean(exp.coef_PRMP)
exp.coef.ci <- quantile(exp.coef_PRMP, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PRMP",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PRMP2 <- exp(1 * leroux_nassau$samples$beta[, 13] / sd(EJScreen_transform$PRMP))
exp.coef.mean <- mean(exp.coef_PRMP2)
exp.coef.ci <- quantile(exp.coef_PRMP2, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PRMP2",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PTSDF <- exp(1 * leroux_nassau$samples$beta[, 14] / sd(EJScreen_transform$PTSDF))
exp.coef.mean <- mean(exp.coef_PTSDF)
exp.coef.ci <- quantile(exp.coef_PTSDF, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PTSDF",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PTSDF2 <- exp(1 * leroux_nassau$samples$beta[, 15] / sd(EJScreen_transform$PTSDF))
exp.coef.mean <- mean(exp.coef_PTSDF2)
exp.coef.ci <- quantile(exp.coef_PTSDF2, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PTSDF2",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_USTF <- exp(1 * leroux_nassau$samples$beta[, 16] / sd(EJScreen_transform$UST))
exp.coef.mean <- mean(exp.coef_USTF)
exp.coef.ci <- quantile(exp.coef_USTF, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "UST",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_USTF2 <- exp(1 * leroux_nassau$samples$beta[, 17] / sd(EJScreen_transform$UST))
exp.coef.mean <- mean(exp.coef_USTF2)
exp.coef.ci <- quantile(exp.coef_USTF2, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "UST",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_Ozone <- exp(1 * leroux_nassau$samples$beta[, 18] / sd(EJScreen_transform$OZONE))
exp.coef.mean <- mean(exp.coef_Ozone)
exp.coef.ci <- quantile(exp.coef_Ozone, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "Ozone",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PM25 <- exp(1 * leroux_nassau$samples$beta[, 19] / sd(EJScreen_transform$PM25))
exp.coef.mean <- mean(exp.coef_PM25)
exp.coef.ci <- quantile(exp.coef_PM25, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PM 2.5",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
print(result)## variable mean ci_lower ci_upper
## 2.5% MINORPCT 0.9676339 0.8867494 1.053027
## 2.5%1 LOWINCPCT 1.0673939 0.9364139 1.208930
## 2.5%2 LESSHSPCT 0.9665078 0.7788830 1.184580
## 2.5%3 LINGISOPCT 1.0050328 0.8106335 1.233738
## 2.5%4 OVER64PCT 1.0373003 0.8835148 1.204565
## 2.5%5 PRE1960PCT 1.0472621 0.9692092 1.130188
## 2.5%6 PTRAF 1.0000285 0.9999665 1.000092
## 2.5%7 PTRAF2 1.0000344 0.9999624 1.000104
## 2.5%8 PWDIS 1.1905898 0.6918703 1.936643
## 2.5%9 PWDIS2 1.5078026 0.9169971 2.368575
## 2.5%10 PRMP 1.6570402 0.7183690 3.214165
## 2.5%11 PRMP2 3.8743832 1.5107540 7.942776
## 2.5%12 PTSDF 1.0417577 1.0028424 1.080646
## 2.5%13 PTSDF2 1.0760906 1.0249598 1.128936
## 2.5%14 UST 1.0098454 0.9899115 1.030352
## 2.5%15 UST 1.0069644 0.9869347 1.027140
## 2.5%16 Ozone 0.8771968 0.6805547 1.080087
## 2.5%17 PM 2.5 2.4014042 1.5363799 3.619540# Suffolk
SuffolkEJ <- subset(EJScreen_transform, CNTY_NAME == "Suffolk County")
# Create neighborhood structure with queen method
queenNB3 <- poly2nb(SuffolkEJ, queen = TRUE) # Queen adjacency of irregular lattice
coords3 <- st_coordinates(SuffolkEJ) # Obtain the centroids of the features
ids3 <- row.names(sf::st_set_geometry(SuffolkEJ, NULL))
# Convert queen neighbor list to a binary regular matrix object
W_mat3 <- nb2mat(queenNB3, zero.policy = TRUE, style = "B")
# Show Queen Adjacency
plot(queenNB3, st_geometry(SuffolkEJ), pch=20, cex=0.5, col="gray"); title("Queen Adjacency for Suffolk County")
# Fit Gaussian CAR Model for Suffolk County
leroux_suffolk <- S.CARleroux(formula = PNPL_LOGAR ~ MINORPCT_S + LOWINCPCT_S + LESSHSPCT_S +
LINGISOPCT_S + OVER64PCT_S
+ PRE1960PCT_S + PTRAF_F + PRMP_F + PTSDF_F + UST_F + OZONE_S + PM25_S,
data = SuffolkEJ,
family = "gaussian",
W = W_mat3,
MALA = FALSE,
burnin = 50000,
n.sample = 300000,
thin = 10,
prior.mean.beta = rep(0, times = 17),
prior.var.beta = c(rep((100^2), times = 17)),
prior.nu2 = c(0.01, 0.01),
prior.tau2 = c(0.01, 0.01),
verbose = TRUE)## Setting up the model.## Warning in mat2listw(W): style is M (missing); style should be set to a valid
## value## Generating 25000 post burnin and thinned (if requested) samples.
##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|=== | 4%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 24%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 90%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 100%
## Summarising results.
## Finished in 94.3 seconds.print(leroux_suffolk)##
## #################
## #### Model fitted
## #################
## Likelihood model - Gaussian (identity link function)
## Random effects model - Leroux CAR
## Regression equation - PNPL_LOGAR ~ MINORPCT_S + LOWINCPCT_S + LESSHSPCT_S + LINGISOPCT_S +
## OVER64PCT_S + PRE1960PCT_S + PTRAF_F + PRMP_F + PTSDF_F +
## UST_F + OZONE_S + PM25_S
## Number of missing observations - 0
##
## ############
## #### Results
## ############
## Posterior quantities and DIC
##
## Mean 2.5% 97.5% n.effective Geweke.diag
## (Intercept) -1.3158 -1.5045 -1.1362 168.0 2.7
## MINORPCT_S 0.0114 -0.0153 0.0381 2486.2 0.4
## LOWINCPCT_S 0.0016 -0.0144 0.0179 3231.8 -0.6
## LESSHSPCT_S -0.0028 -0.0232 0.0176 3327.6 0.6
## LINGISOPCT_S 0.0044 -0.0164 0.0253 4016.1 -0.1
## OVER64PCT_S -0.0039 -0.0186 0.0108 2853.3 1.2
## PRE1960PCT_S -0.0187 -0.0427 0.0060 2357.9 -0.2
## PTRAF_F2 -0.0113 -0.0483 0.0259 3291.8 0.5
## PTRAF_F3 -0.0146 -0.0638 0.0348 2623.0 -0.2
## PRMP_F2 0.1941 0.1123 0.2804 617.7 -0.4
## PRMP_F3 0.3224 0.1856 0.4639 363.5 0.0
## PTSDF_F2 0.0941 0.0386 0.1503 1748.8 -1.7
## PTSDF_F3 0.2380 0.1631 0.3135 1228.4 -1.4
## UST_F2 -0.0047 -0.0383 0.0301 3011.4 -1.4
## UST_F3 -0.0017 -0.0504 0.0479 2727.9 -1.8
## OZONE_S -0.3193 -0.4357 -0.1956 190.9 -2.6
## PM25_S 0.8612 0.6303 1.0663 109.7 2.5
## nu2 0.0016 0.0009 0.0026 3288.0 0.0
## tau2 0.3031 0.2775 0.3309 23268.4 -0.4
## rho 0.9968 0.9902 0.9998 9267.5 1.5
##
## DIC = -2890.41 p.d = 983.6684 LMPL = 1068.11suffolk_residuals <- resid(leroux_suffolk)
hist(suffolk_residuals, main="Histogram of Suffolk County Residuals", xlab="Residuals")
#exponentiate
# Create a data frame to store the results
result <- data.frame(variable = character(),
mean = double(),
ci_lower = double(),
ci_upper = double(),
stringsAsFactors = FALSE)
# Exponentiate coefficients for MINORPCT
exp.coef_MINORPCT <- exp(1 * leroux_suffolk$samples$beta[, 2] / sd(SuffolkEJ$MINORPCT))
# Compute mean and 95% confidence interval
exp.coef.mean <- mean(exp.coef_MINORPCT)
exp.coef.ci <- quantile(exp.coef_MINORPCT, c(0.025, 0.975))
# Add results to data frame
result <- rbind(result, data.frame(variable = "MINORPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
# Repeat for each variable
exp.coef_LOWINCPCT <- exp(1 * leroux_suffolk$samples$beta[, 3] / sd(SuffolkEJ$LOWINCPCT))
exp.coef.mean <- mean(exp.coef_LOWINCPCT)
exp.coef.ci <- quantile(exp.coef_LOWINCPCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "LOWINCPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_LESSHSPCT <- exp(1 * leroux_suffolk$samples$beta[, 4] / sd(SuffolkEJ$LESSHSPCT))
exp.coef.mean <- mean(exp.coef_LESSHSPCT)
exp.coef.ci <- quantile(exp.coef_LESSHSPCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "LESSHSPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_LINGISOPCT <- exp(1 * leroux_suffolk$samples$beta[, 5] / sd(SuffolkEJ$LINGISOPCT))
exp.coef.mean <- mean(exp.coef_LINGISOPCT)
exp.coef.ci <- quantile(exp.coef_LINGISOPCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "LINGISOPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_OVER64PCT <- exp(1 * leroux_suffolk$samples$beta[, 6] / sd(SuffolkEJ$OVER64PCT))
exp.coef.mean <- mean(exp.coef_OVER64PCT)
exp.coef.ci <- quantile(exp.coef_OVER64PCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "OVER64PCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PRE1960PCT <- exp(1 * leroux_suffolk$samples$beta[, 7] / sd(SuffolkEJ$PRE1960PCT))
exp.coef.mean <- mean(exp.coef_PRE1960PCT)
exp.coef.ci <- quantile(exp.coef_PRE1960PCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PRE1960PCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PTRAF <- exp(1 * leroux_suffolk$samples$beta[, 8] / sd(SuffolkEJ$PTRAF))
exp.coef.mean <- mean(exp.coef_PTRAF)
exp.coef.ci <- quantile(exp.coef_PTRAF, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PTRAF",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PTRAF2 <- exp(1 * leroux_suffolk$samples$beta[, 9] / sd(SuffolkEJ$PTRAF))
exp.coef.mean <- mean(exp.coef_PTRAF2)
exp.coef.ci <- quantile(exp.coef_PTRAF2, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PTRAF2",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PRMP <- exp(1 * leroux_suffolk$samples$beta[, 10] / sd(SuffolkEJ$PRMP))
exp.coef.mean <- mean(exp.coef_PRMP)
exp.coef.ci <- quantile(exp.coef_PRMP, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PRMP",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PRMP2 <- exp(1 * leroux_suffolk$samples$beta[, 11] / sd(SuffolkEJ$PRMP))
exp.coef.mean <- mean(exp.coef_PRMP2)
exp.coef.ci <- quantile(exp.coef_PRMP2, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PRMP2",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PTSDF <- exp(1 * leroux_suffolk$samples$beta[, 12] / sd(SuffolkEJ$PTSDF))
exp.coef.mean <- mean(exp.coef_PTSDF)
exp.coef.ci <- quantile(exp.coef_PTSDF, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PTSDF",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PTSDF2 <- exp(1 * leroux_suffolk$samples$beta[, 13] / sd(SuffolkEJ$PTSDF))
exp.coef.mean <- mean(exp.coef_PTSDF2)
exp.coef.ci <- quantile(exp.coef_PTSDF2, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PTSDF2",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_USTF <- exp(1 * leroux_suffolk$samples$beta[, 14] / sd(SuffolkEJ$UST))
exp.coef.mean <- mean(exp.coef_USTF)
exp.coef.ci <- quantile(exp.coef_USTF, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "UST",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_USTF2 <- exp(1 * leroux_suffolk$samples$beta[, 15] / sd(SuffolkEJ$UST))
exp.coef.mean <- mean(exp.coef_USTF2)
exp.coef.ci <- quantile(exp.coef_USTF2, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "UST2",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PM25 <- exp(1 * leroux_suffolk$samples$beta[, 16] / sd(SuffolkEJ$PM25))
exp.coef.mean <- mean(exp.coef_PM25)
exp.coef.ci <- quantile(exp.coef_PM25, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "Ozone",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_Ozone <- exp(1 * leroux_suffolk$samples$beta[, 17] / sd(SuffolkEJ$OZONE))
exp.coef.mean <- mean(exp.coef_Ozone)
exp.coef.ci <- quantile(exp.coef_Ozone, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PM 2.5",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
print(result)## variable mean ci_lower ci_upper
## 2.5% MINORPCT 1.0469780 0.9420616 1.1607454
## 2.5%1 LOWINCPCT 1.0142325 0.8975217 1.1439550
## 2.5%2 LESSHSPCT 0.9768257 0.7872075 1.1982314
## 2.5%3 LINGISOPCT 1.1102087 0.7291245 1.6241839
## 2.5%4 OVER64PCT 0.9726198 0.8699170 1.0841730
## 2.5%5 PRE1960PCT 0.9206743 0.8249529 1.0271872
## 2.5%6 PTRAF 0.9999756 0.9998958 1.0000560
## 2.5%7 PTRAF2 0.9999684 0.9998622 1.0000752
## 2.5%8 PRMP 5.5264633 2.5892136 10.7490106
## 2.5%9 PRMP2 18.4110096 4.8158518 50.8603921
## 2.5%10 PTSDF 1.0580299 1.0233273 1.0940438
## 2.5%11 PTSDF2 1.1533164 1.1024576 1.2062364
## 2.5%12 UST 0.9951109 0.9596135 1.0328851
## 2.5%13 UST2 0.9985922 0.9471590 1.0529178
## 2.5%14 Ozone 0.4366253 0.3172238 0.5971324
## 2.5%15 PM 2.5 3.4880400 2.4722875 4.6236983glm <- S.glm(formula = PNPL_LOGAR ~ MINORPCT_S + LOWINCPCT_S + LESSHSPCT_S +
LINGISOPCT_S + OVER64PCT_S
+ PRE1960PCT_S + PTRAF_F + PWDIS_F + PRMP_F + PTSDF_F + UST_F + OZONE_S + PM25_S,
data = EJScreen_transform,
family = "gaussian",
MALA = FALSE,
burnin = 50000,
n.sample = 300000,
thin = 10,
prior.mean.beta = rep(0, times = 19),
prior.var.beta = c(rep((100^2), times = 19)),
prior.nu2 = c(0.01, 0.01),
verbose = TRUE)## Setting up the model.
## Generating 25000 post burnin and thinned (if requested) samples.
##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|=== | 4%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 24%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 90%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 100%
## Summarising results.
## Finished in 89.3 seconds.print(glm)##
## #################
## #### Model fitted
## #################
## Likelihood model - Gaussian (identity link function)
## Random effects model - None
## Regression equation - PNPL_LOGAR ~ MINORPCT_S + LOWINCPCT_S + LESSHSPCT_S + LINGISOPCT_S +
## OVER64PCT_S + PRE1960PCT_S + PTRAF_F + PWDIS_F + PRMP_F +
## PTSDF_F + UST_F + OZONE_S + PM25_S
## Number of missing observations - 0
##
## ############
## #### Results
## ############
## Posterior quantities and DIC
##
## Mean 2.5% 97.5% n.effective Geweke.diag
## (Intercept) -2.4914 -2.6007 -2.3822 25819.6 -0.1
## MINORPCT_S -0.1227 -0.1610 -0.0843 25000.0 -1.0
## LOWINCPCT_S 0.0263 -0.0059 0.0586 25953.3 0.5
## LESSHSPCT_S 0.0198 -0.0176 0.0571 25000.0 0.3
## LINGISOPCT_S 0.0321 -0.0005 0.0647 25000.0 1.3
## OVER64PCT_S -0.0305 -0.0593 -0.0019 25000.0 0.5
## PRE1960PCT_S -0.0506 -0.0860 -0.0150 25637.5 0.4
## PTRAF_F2 -0.0855 -0.1580 -0.0131 25000.0 -0.4
## PTRAF_F3 -0.0873 -0.1639 -0.0101 25592.2 -0.3
## PWDIS_F2 0.3199 0.2058 0.4346 25000.0 -0.8
## PWDIS_F3 0.6696 0.5822 0.7575 25000.0 -0.6
## PRMP_F2 0.3348 0.2510 0.4189 25000.0 0.1
## PRMP_F3 0.8524 0.7619 0.9441 23909.5 0.3
## PTSDF_F2 0.1971 0.1239 0.2701 25000.0 1.1
## PTSDF_F3 0.5688 0.4925 0.6454 25000.0 0.4
## UST_F2 0.0449 -0.0209 0.1113 25000.0 1.3
## UST_F3 0.0299 -0.0401 0.1009 25816.3 0.5
## OZONE_S -0.0998 -0.1344 -0.0651 25000.0 0.0
## PM25_S 0.3746 0.3259 0.4243 26636.9 -0.7
## nu2 0.4149 0.3910 0.4404 25000.0 1.1
##
## DIC = 4307.52 p.d = 19.99882 LMPL = -2154.07# Create a data frame to store the results
result <- data.frame(variable = character(),
mean = double(),
ci_lower = double(),
ci_upper = double(),
stringsAsFactors = FALSE)
# Exponentiate coefficients for MINORPCT
exp.coef_MINORPCT <- exp(1 * glm$samples$beta[, 2] / sd(EJScreen_transform$MINORPCT))
# Compute mean and 95% confidence interval
exp.coef.mean <- mean(exp.coef_MINORPCT)
exp.coef.ci <- quantile(exp.coef_MINORPCT, c(0.025, 0.975))
# Add results to data frame
result <- rbind(result, data.frame(variable = "MINORPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
# Repeat for each variable
exp.coef_LOWINCPCT <- exp(1 * glm$samples$beta[, 3] / sd(EJScreen_transform$LOWINCPCT))
exp.coef.mean <- mean(exp.coef_LOWINCPCT)
exp.coef.ci <- quantile(exp.coef_LOWINCPCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "LOWINCPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_LESSHSPCT <- exp(1 * glm$samples$beta[, 4] / sd(EJScreen_transform$LESSHSPCT))
exp.coef.mean <- mean(exp.coef_LESSHSPCT)
exp.coef.ci <- quantile(exp.coef_LESSHSPCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "LESSHSPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_LINGISOPCT <- exp(1 * glm$samples$beta[, 5] / sd(EJScreen_transform$LINGISOPCT))
exp.coef.mean <- mean(exp.coef_LINGISOPCT)
exp.coef.ci <- quantile(exp.coef_LINGISOPCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "LINGISOPCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_OVER64PCT <- exp(1 * glm$samples$beta[, 6] / sd(EJScreen_transform$OVER64PCT))
exp.coef.mean <- mean(exp.coef_OVER64PCT)
exp.coef.ci <- quantile(exp.coef_OVER64PCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "OVER64PCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PRE1960PCT <- exp(1 * glm$samples$beta[, 7] / sd(EJScreen_transform$PRE1960PCT))
exp.coef.mean <- mean(exp.coef_PRE1960PCT)
exp.coef.ci <- quantile(exp.coef_PRE1960PCT, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PRE1960PCT",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PTRAF <- exp(1 * glm$samples$beta[, 8] / sd(EJScreen_transform$PTRAF))
exp.coef.mean <- mean(exp.coef_PTRAF)
exp.coef.ci <- quantile(exp.coef_PTRAF, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PTRAF",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PTRAF2 <- exp(1 * glm$samples$beta[, 9] / sd(EJScreen_transform$PTRAF))
exp.coef.mean <- mean(exp.coef_PTRAF2)
exp.coef.ci <- quantile(exp.coef_PTRAF2, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PTRAF2",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PWDIS <- exp(1 * glm$samples$beta[, 10] / sd(EJScreen_transform$PWDIS))
exp.coef.mean <- mean(exp.coef_PWDIS)
exp.coef.ci <- quantile(exp.coef_PWDIS, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PWDIS",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PWDIS2 <- exp(1 * glm$samples$beta[, 11] / sd(EJScreen_transform$PWDIS))
exp.coef.mean <- mean(exp.coef_PWDIS2)
exp.coef.ci <- quantile(exp.coef_PWDIS2, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PWDIS2",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PRMP <- exp(1 * glm$samples$beta[, 12] / sd(EJScreen_transform$PRMP))
exp.coef.mean <- mean(exp.coef_PRMP)
exp.coef.ci <- quantile(exp.coef_PRMP, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PRMP",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PRMP2 <- exp(1 * glm$samples$beta[, 13] / sd(EJScreen_transform$PRMP))
exp.coef.mean <- mean(exp.coef_PRMP2)
exp.coef.ci <- quantile(exp.coef_PRMP2, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PRMP2",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PTSDF <- exp(1 * glm$samples$beta[, 14] / sd(EJScreen_transform$PTSDF))
exp.coef.mean <- mean(exp.coef_PTSDF)
exp.coef.ci <- quantile(exp.coef_PTSDF, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PTSDF",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PTSDF2 <- exp(1 * glm$samples$beta[, 15] / sd(EJScreen_transform$PTSDF))
exp.coef.mean <- mean(exp.coef_PTSDF2)
exp.coef.ci <- quantile(exp.coef_PTSDF2, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PTSDF2",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_USTF <- exp(1 * glm$samples$beta[, 16] / sd(EJScreen_transform$UST))
exp.coef.mean <- mean(exp.coef_USTF)
exp.coef.ci <- quantile(exp.coef_USTF, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "UST",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_USTF2 <- exp(1 * glm$samples$beta[, 17] / sd(EJScreen_transform$UST))
exp.coef.mean <- mean(exp.coef_USTF2)
exp.coef.ci <- quantile(exp.coef_USTF2, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "UST2",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_PM25 <- exp(1 * glm$samples$beta[, 18] / sd(EJScreen_transform$PM25))
exp.coef.mean <- mean(exp.coef_PM25)
exp.coef.ci <- quantile(exp.coef_PM25, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "Ozone",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
exp.coef_Ozone <- exp(1 * glm$samples$beta[, 19] / sd(EJScreen_transform$OZONE))
exp.coef.mean <- mean(exp.coef_Ozone)
exp.coef.ci <- quantile(exp.coef_Ozone, c(0.025, 0.975))
result <- rbind(result, data.frame(variable = "PM 2.5",
mean = exp.coef.mean,
ci_lower = exp.coef.ci[1],
ci_upper = exp.coef.ci[2],
stringsAsFactors = FALSE))
print(result)## variable mean ci_lower ci_upper
## 2.5% MINORPCT 0.6432830 0.5586957 0.7373228
## 2.5%1 LOWINCPCT 1.2285068 0.9566245 1.5550134
## 2.5%2 LESSHSPCT 1.2566861 0.8310571 1.8260186
## 2.5%3 LINGISOPCT 1.6503294 0.9929766 2.5857446
## 2.5%4 OVER64PCT 0.7682387 0.5891937 0.9830644
## 2.5%5 PRE1960PCT 0.8421577 0.7444880 0.9499431
## 2.5%6 PTRAF 0.9998753 0.9997696 0.9999808
## 2.5%7 PTRAF2 0.9998726 0.9997609 0.9999852
## 2.5%8 PWDIS 5.5709854 2.9307695 9.6852193
## 2.5%9 PWDIS2 33.9726527 20.9435888 52.3226927
## 2.5%10 PRMP 5.2308168 3.4000906 7.7103548
## 2.5%11 PRMP2 65.5096538 41.0675104 99.8886332
## 2.5%12 PTSDF 1.1204880 1.0740024 1.1683137
## 2.5%13 PTSDF2 1.3880207 1.3280277 1.4502809
## 2.5%14 UST 1.0230829 0.9894936 1.0577964
## 2.5%15 UST2 1.0153672 0.9799857 1.0522455
## 2.5%16 Ozone 0.8072610 0.7487588 0.8692968
## 2.5%17 PM 2.5 1.8547217 1.7102062 2.0110361