newdata<-read.csv("Data/sampleforglmmTMB.csv")
################################ SHORT TERM ##############################################
library(glmmTMB)
library(emmeans)
## Warning: package 'emmeans' was built under R version 3.6.3
m <- glmmTMB(
Deaths ~ meanUVAKJ100
# air pollution
+ scale(mean_pm25)
# other weather
+ scale(mean_winter_temp)
+ scale(UVvitdyearmean)
# demographics
+ scale(pct_blk) + scale(older_pecent) + scale(popdensity) + scale(hispanic)
# Deprivation
+ scale(PCA1)
# health care factors
+ scale(statepositivetestsperpop)
# urban rural
+ factor(RUCC_2013_CAT)
# population
+ offset(log(population)) + (1 | state),
ziformula = ~ meanUVAKJ100
+ offset(log(population)),
family = nbinom2,
data = newdata
)
## show estimates
m
## Formula:
## Deaths ~ meanUVAKJ100 + scale(mean_pm25) + scale(mean_winter_temp) +
## scale(UVvitdyearmean) + scale(pct_blk) + scale(older_pecent) +
## scale(popdensity) + scale(hispanic) + scale(PCA1) + scale(statepositivetestsperpop) +
## factor(RUCC_2013_CAT) + offset(log(population)) + (1 | state)
## Zero inflation: ~meanUVAKJ100 + offset(log(population))
## Data: newdata
## AIC BIC logLik df.resid
## 7746.928 7843.472 -3857.464 3068
## Random-effects (co)variances:
##
## Conditional model:
## Groups Name Std.Dev.
## state (Intercept) 0.6951
##
## Number of obs: 3084 / Conditional model: state, 49
##
## Overdispersion parameter for nbinom2 family (): 0.839
##
## Fixed Effects:
##
## Conditional model:
## (Intercept) meanUVAKJ100
## -10.52086 -0.16097
## scale(mean_pm25) scale(mean_winter_temp)
## 0.14183 -0.15514
## scale(UVvitdyearmean) scale(pct_blk)
## 0.33808 0.39985
## scale(older_pecent) scale(popdensity)
## 0.08541 0.08320
## scale(hispanic) scale(PCA1)
## 0.11972 -0.25780
## scale(statepositivetestsperpop) factor(RUCC_2013_CAT)1
## 0.39972 0.21842
##
## Zero-inflation model:
## (Intercept) meanUVAKJ100
## -8.682 -1.787
# create marginal effects
t<-emmeans(m, c("meanUVAKJ100"), at = list(meanUVAKJ100 = c(4,5,6,7,8,9)), type = "re", offset=c(log(1000000)))
t
## meanUVAKJ100 response SE df lower.CL upper.CL
## 4 15.80 4.07 3068 9.52 26.2
## 5 13.45 2.35 3068 9.54 19.0
## 6 11.45 1.40 3068 9.01 14.5
## 7 9.75 1.35 3068 7.42 12.8
## 8 8.30 1.73 3068 5.51 12.5
## 9 7.06 2.10 3068 3.94 12.6
##
## Results are averaged over the levels of: RUCC_2013_CAT
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
#plot(t, by = NULL, horizontal = FALSE, colors = "darkgreen")
######################################## LONG TERM ##############################################
m <- glmmTMB(
Deaths ~ UVAyearmean
# air pollution
+ scale(mean_pm25)
# other weather
+ scale(mean_winter_temp)
+ scale(UVvitdyearmean)
# demographics
+ scale(pct_blk) + scale(older_pecent) + scale(popdensity)+ scale(hispanic)
# Deprivation
+ scale(PCA1)
# health care factors
+ scale(statepositivetestsperpop)
# urban rural
+ factor(RUCC_2013_CAT)
# population
+ offset(log(population)) + (1 | state),
ziformula = ~ UVAyearmean
+ offset(log(population)),
family = nbinom2,
data = newdata
)
## show estimates
m
## Formula:
## Deaths ~ UVAyearmean + scale(mean_pm25) + scale(mean_winter_temp) +
## scale(UVvitdyearmean) + scale(pct_blk) + scale(older_pecent) +
## scale(popdensity) + scale(hispanic) + scale(PCA1) + scale(statepositivetestsperpop) +
## factor(RUCC_2013_CAT) + offset(log(population)) + (1 | state)
## Zero inflation: ~UVAyearmean + offset(log(population))
## Data: newdata
## AIC BIC logLik df.resid
## 7756.021 7852.591 -3862.011 3073
## Random-effects (co)variances:
##
## Conditional model:
## Groups Name Std.Dev.
## state (Intercept) 0.7062
##
## Number of obs: 3089 / Conditional model: state, 49
##
## Overdispersion parameter for nbinom2 family (): 0.84
##
## Fixed Effects:
##
## Conditional model:
## (Intercept) UVAyearmean
## -10.32926 -0.44186
## scale(mean_pm25) scale(mean_winter_temp)
## 0.16966 -0.07455
## scale(UVvitdyearmean) scale(pct_blk)
## 0.29147 0.39629
## scale(older_pecent) scale(popdensity)
## 0.08595 0.08097
## scale(hispanic) scale(PCA1)
## 0.11027 -0.25350
## scale(statepositivetestsperpop) factor(RUCC_2013_CAT)1
## 0.40750 0.22277
##
## Zero-inflation model:
## (Intercept) UVAyearmean
## -10.451 -3.025
#####################
#### eemeans
t<-emmeans(m, c("UVAyearmean"),at = list(UVAyearmean = c(2.3,2.6,2.9, 3.2)), type = "re", offset=c(log(1000000)))
t
## UVAyearmean response SE df lower.CL upper.CL
## 2.3 13.22 5.00 3073 6.29 27.8
## 2.6 11.58 1.87 3073 8.44 15.9
## 2.9 10.14 1.90 3073 7.02 14.6
## 3.2 8.88 3.67 3073 3.95 20.0
##
## Results are averaged over the levels of: RUCC_2013_CAT
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
############################################# interaction
# 1. . Include a product (interaction) term in the model if you want to investigate whether the effect of one independent variable is contingent on the level of a second independent variable
m <- glmmTMB(
Deaths ~ meanUVAKJ100*UVAyearmean
# air pollution
+ scale(mean_pm25)
# other weather
+ scale(mean_winter_temp)
+ scale(UVvitdyearmean)
# demographics
+ scale(pct_blk) + scale(older_pecent) + scale(popdensity) + scale(hispanic)
# Deprivation
+ scale(PCA1)
# health care factors
+ scale(statepositivetestsperpop)
# urban rural
+ factor(RUCC_2013_CAT)
# population
+ offset(log(population)) + (1 | state),
ziformula = ~ meanUVAKJ100*UVAyearmean
+ offset(log(population)),
family = nbinom2,
data = newdata
)
## show estimates
m
## Formula:
## Deaths ~ meanUVAKJ100 * UVAyearmean + scale(mean_pm25) + scale(mean_winter_temp) +
## scale(UVvitdyearmean) + scale(pct_blk) + scale(older_pecent) +
## scale(popdensity) + scale(hispanic) + scale(PCA1) + scale(statepositivetestsperpop) +
## factor(RUCC_2013_CAT) + offset(log(population)) + (1 | state)
## Zero inflation: ~meanUVAKJ100 * UVAyearmean + offset(log(population))
## Data: newdata
## AIC BIC logLik df.resid
## 7747.582 7868.262 -3853.791 3064
## Random-effects (co)variances:
##
## Conditional model:
## Groups Name Std.Dev.
## state (Intercept) 0.6784
##
## Number of obs: 3084 / Conditional model: state, 49
##
## Overdispersion parameter for nbinom2 family (): 0.843
##
## Fixed Effects:
##
## Conditional model:
## (Intercept) meanUVAKJ100
## -4.62111 -0.82633
## UVAyearmean scale(mean_pm25)
## -2.19857 0.16334
## scale(mean_winter_temp) scale(UVvitdyearmean)
## -0.15288 0.51022
## scale(pct_blk) scale(older_pecent)
## 0.40222 0.07566
## scale(popdensity) scale(hispanic)
## 0.08077 0.08923
## scale(PCA1) scale(statepositivetestsperpop)
## -0.25594 0.37694
## factor(RUCC_2013_CAT)1 meanUVAKJ100:UVAyearmean
## 0.20557 0.24477
##
## Zero-inflation model:
## (Intercept) meanUVAKJ100 UVAyearmean
## -394.36 74.08 165.45
## meanUVAKJ100:UVAyearmean
## -32.56
# 2. Ignore the coefficient of the product term: It does not necessarily provide accurate information about the significance, magnitude, or even the direction of the underlying interaction effect on the predictions.
#####################################
# 3. Plot the predictions to determine the nature of the underlying interaction effect on the metric of interest.
library(ggplot2)
p<-emmip(m,UVAyearmean~meanUVAKJ100, at = list(meanUVAKJ100 =c(5,6,7,8), UVAyearmean = c(2.3,2.6,2.9, 3.2)), type="re", CIs=T, offset=c(log(1000000))) + theme_classic(base_size = 15)
p
p + facet_grid(rows = vars(UVAyearmean))
#####################################
# 4.Determine the size and significance of the effects of interest using marginal effects, not regression coefficients.
df<-emmip(m,UVAyearmean~meanUVAKJ100, at = list(meanUVAKJ100 =c(5,6,7,8), UVAyearmean = c(2.3,2.6,2.9, 3.2)), type="re", plotit=F, CIs=T, offset=c(log(1000000)))
df
### There is a difference, almost treble the predicted deaths, the confidence intervals lap, so not a strong difference.
# 5. Use tests of second differences (whether two marginal effects are equal) to determine whether an interaction effect is significant for specific values of interest of your independent variables.
test(emtrends(m, ~ UVAyearmean, at=list(UVAyearmean = c(2.3)), var="meanUVAKJ100"))
test(emtrends(m, ~ UVAyearmean, at=list(UVAyearmean = c(2.6)), var="meanUVAKJ100"))
test(emtrends(m, ~ UVAyearmean, at=list(UVAyearmean = c(2.9)), var="meanUVAKJ100"))
test(emtrends(m, ~ UVAyearmean, at=list(UVAyearmean = c(3.2)), var="meanUVAKJ100"))
# 6. Absent substantive or theoretically interesting values of the focal independent or control variables to test the interaction effect, use average marginal effects to summarize whether there is an interaction effect present on average in the data.
# trend p-value
test(emtrends(m, pairwise ~ UVAyearmean, var="meanUVAKJ100"))
## $emtrends
## UVAyearmean meanUVAKJ100.trend SE df t.ratio p.value
## 2.72 -0.161 0.103 3064 -1.557 0.1196
##
## Results are averaged over the levels of: RUCC_2013_CAT
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## (nothing) nonEst NA NA NA NA
##
## Results are averaged over the levels of: RUCC_2013_CAT