Commit cb4e6a44 authored by Kara McCormack's avatar Kara McCormack

made table with exponentiated fixed effects from model

parent cf942e2c
......@@ -39,12 +39,13 @@ gr <- inla.read.graph(filename = "nyc_map.adj")
m0 <- ma_pos_7 ~ male_prop +
pop_non_hispanic_black_prop +
hispanic_prop +
age_over_65_prop
age_over_65_prop + mgmt_prof_prop
m1 <- ma_pos_7 ~ male_prop +
pop_non_hispanic_black_prop +
hispanic_prop +
age_over_65_prop +
mgmt_prof_prop +
f(idarea, model = "bym", graph = gr) +
f(idtime, model = "ar1")
......@@ -52,6 +53,7 @@ m2 <- ma_pos_7 ~ male_prop +
pop_non_hispanic_black_prop +
hispanic_prop +
age_over_65_prop +
mgmt_prof_prop +
f(com_class, model = "iid") +
f(acs_class, model = "iid") +
f(idarea, model = "bym", graph = gr) +
......@@ -61,6 +63,7 @@ m3 <- ma_pos_7 ~ male_prop +
pop_non_hispanic_black_prop +
hispanic_prop +
age_over_65_prop +
mgmt_prof_prop +
male_prop*pop_non_hispanic_black_prop +
male_prop*hispanic_prop +
f(com_class, model = "iid") +
......@@ -72,6 +75,7 @@ m4 <- ma_pos_7 ~ male_prop +
pop_non_hispanic_black_prop +
hispanic_prop +
age_over_65_prop +
mgmt_prof_prop +
male_prop*hispanic_prop +
f(com_class, model = "iid") +
f(acs_class, model = "iid") +
......@@ -82,22 +86,24 @@ m5 <- ma_pos_7 ~ male_prop +
pop_non_hispanic_black_prop +
hispanic_prop +
age_over_65_prop +
mgmt_prof_prop +
as.factor(com_class) +
as.factor(acs_class) +
f(idarea, model = "bym", graph = gr) +
f(idtime, model = "ar1")
list_models <- list(m0, m1, m2, m3, m4)
model_names <- c("m0", "m1", "m2", "m3", "m4")
list_models <- list(m0, m1, m2, m3, m4, m5)
model_names <- c("m0", "m1", "m2", "m3", "m4", "m5")
model_descriptions <- c("only fixed effects",
"fixed + S + T",
"fixed + S + T + latent classes",
"fixed + interactions + S + T + LC",
"m3 without male*black")
"m3 without male*black",
"no interax, random vars as both fixed+random")
## adding in proportion of management and professional.
m6 <-
```
```{r eval = F}
......@@ -124,13 +130,6 @@ saveRDS(res_results, file = "res_results.rds")
saveRDS(dic_results, file = "dic_results.rds")
saveRDS(waic_results, file = "waic_results.rds")
#additional model
res_m5 <- inla(m5,
family = "poisson",
data = d3,
E = expected_daily_ma,
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE))
```
Load the results so don't have to run it every time.
```{r}
......@@ -139,19 +138,23 @@ dic_results <- readRDS("dic_results.rds")
waic_results <- readRDS("waic_results.rds")
```
$$y_{ij} = \beta_0 + $$
Table to show all information.
```{r}
table_dat <- data.frame(model_name = model_names,
model_description = model_descriptions,
DIC = dic_results,
WAIC = waic_results)
# put results in table, bold the best model
table_dat %>%
arrange(DIC)%>%
kable %>%
kable_styling("striped", full_width = F) %>%
row_spec(3, bold = T, color = "black", background = "khaki1")
row_spec(1, bold = T, color = "black", background = "khaki1") %>%
row_spec(2, bold = T, color = "coral2", background = "khaki1")
```
Model 2 (the model with only fixed effects, spatial and temporal effects, and latent classes) looks to be our best bet here in terms of DIC and WAIC.
Model 5 (all fixed effects, random variables as both fixed and random) looks to be our best bet here in terms of DIC and WAIC. However, the parametrization doesn't make sense. Model 2 then is our final model.
Note: WAIC is an extension of the Akaike Information Criterion (AIC) that is more fully Bayesian than the Deviance Information Criterion (DIC). Compared to AIC and DIC, WAIC has the desirable property of averaging over the posterior distribution rather than conditioning on a point estimate. This is especially relevant in a predictive context, as WAIC is evaluating the predictions that are actually being used for new data in a Bayesian context.
......@@ -163,11 +166,47 @@ res_results[[3]]$summary.fixed[,c(1, 3, 5)] %>%
fixed_results <- res_results[[3]]$summary.fixed[,c(1, 3, 5)]
rownames(fixed_results) <- c("Intercept", "Male", "Non-Hispanic Black",
"Hispanic", "Age Over 65")
"Hispanic", "Age Over 65", "Mgmt/Prof")
colnames(fixed_results) <- c("Mean", ".025 Quant", ".975 Quant")
fixed_results %>%
kable() %>%
kable_styling("striped", full_width = F)
```
```{r}
mod6 <- res_results[[6]]
mod3 <- res_results[[3]]
names(mod3$marginals.fixed)
marg.male<- inla.tmarginal(exp, mod3$marginals.fixed[["male_prop"]])
inla.est_male <- inla.zmarginal(marg.male)
marg.nh_black<- inla.tmarginal(exp, mod3$marginals.fixed[["pop_non_hispanic_black_prop"]])
inla.est_nh_black <- inla.zmarginal(marg.nh_black)
marg.hispanic<- inla.tmarginal(exp, mod3$marginals.fixed[["hispanic_prop"]])
inla.est_hispanic <- inla.zmarginal(marg.hispanic)
marg.age65<- inla.tmarginal(exp, mod3$marginals.fixed[["age_over_65_prop"]])
inla.est_age65 <- inla.zmarginal(marg.age65)
marg.mgmt_prof<- inla.tmarginal(exp, mod3$marginals.fixed[["mgmt_prof_prop"]])
inla.est_mgmt_prof <- inla.zmarginal(marg.mgmt_prof)
mod_final_fixed <- rbind(inla.est_male,
inla.est_nh_black,
inla.est_hispanic,
inla.est_age65,
inla.est_mgmt_prof)
```
```{r}
mod_final_fixed %>%
as.data.frame() %>%
select(mean, sd, quant0.025, quant0.975) %>%
kable() %>%
kable_styling()
```
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment