Multilevel model 2 - age-specific mortality rates
# I have three data sets that contains the number of deaths - homicides and homicides committed with firearms too - occurred by municipalitiy (2457); year (26); sex (2); age group (19); and population
# I am going to use these three data sets to make another three but with information for the Border States. On these new data sets I'm going to use the model desing that we discussed to estimate age and sex specific mortality rates for each municipality. Hope it works.
# Data set for the Border States - all deaths, homicides and homicides with firearms
border.all<-rbind(bc.all,son.all,chi.all,coa.all,nl.all,tam.all)
border.hom<-rbind(bc.hom,son.hom,chi.hom,coa.hom,nl.hom,tam.hom)
border.fa<-rbind(bc.fa,son.fa,chi.fa,coa.fa,nl.fa,tam.fa)
save(border.all, file = "border.alldeaths.RData")
save(border.hom, file = "border.homicides.RData")
save(border.fa, file = "border.firearms.RData")
The first models uses death rate as the outcome variable with a random intercept.
# All-causes of death in the Border States
border<-glmer(cbind(Freq,pop2)~factor(age.n) + factor(sex) + factor(time) + (1|clave), data = border.all, family = binomial)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00234038 (tol =
## 0.001, component 1)
# Homicide rates in the Border States
borderhom<-glmer(cbind(Freq,pop2)~factor(age.n) + factor(sex) + factor(time) + (1|clave), data = border.hom, family = binomial)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.0117324 (tol =
## 0.001, component 1)
# Homicides committed with firearms rates in the Border States
borderfa<-glmer(cbind(Freq,pop2)~factor(age.n) + factor(sex) + factor(time) + (1|clave), data = border.fa, family = binomial)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.0537465 (tol =
## 0.001, component 1)
# Considering all municipalities in the Border States (276) and all causes of death, we have a variance of 0.1325 with an Intercept of -4.4184.
summary(border)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(Freq, pop2) ~ factor(age.n) + factor(sex) + factor(time) +
## (1 | clave)
## Data: border.all
##
## AIC BIC logLik deviance df.resid
## 178281.9 178503.6 -89116.0 178231.9 52326
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.663 -0.651 -0.287 0.364 47.244
##
## Random effects:
## Groups Name Variance Std.Dev.
## clave (Intercept) 0.1325 0.364
## Number of obs: 52351, groups: clave, 276
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.418421 0.023818 -185.505 < 2e-16 ***
## factor(age.n)2 -3.183841 0.016166 -196.941 < 2e-16 ***
## factor(age.n)3 -4.045693 0.021461 -188.511 < 2e-16 ***
## factor(age.n)4 -3.878884 0.020204 -191.985 < 2e-16 ***
## factor(age.n)5 -2.816486 0.012984 -216.915 < 2e-16 ***
## factor(age.n)6 -2.416984 0.011273 -214.400 < 2e-16 ***
## factor(age.n)7 -2.221597 0.010933 -203.207 < 2e-16 ***
## factor(age.n)8 -2.071285 0.010638 -194.708 < 2e-16 ***
## factor(age.n)9 -1.874007 0.010254 -182.754 < 2e-16 ***
## factor(age.n)10 -1.598539 0.009911 -161.289 < 2e-16 ***
## factor(age.n)11 -1.220503 0.009453 -129.117 < 2e-16 ***
## factor(age.n)12 -0.840982 0.008973 -93.720 < 2e-16 ***
## factor(age.n)13 -0.358359 0.008539 -41.965 < 2e-16 ***
## factor(age.n)14 0.021368 0.008245 2.592 0.00955 **
## factor(age.n)15 0.446357 0.008062 55.368 < 2e-16 ***
## factor(age.n)16 0.854639 0.007946 107.556 < 2e-16 ***
## factor(age.n)17 1.288866 0.007967 161.766 < 2e-16 ***
## factor(age.n)18 1.753319 0.008124 215.822 < 2e-16 ***
## factor(age.n)19 2.388583 0.007685 310.798 < 2e-16 ***
## factor(sex)2 -0.457160 0.003054 -149.677 < 2e-16 ***
## factor(time)2000 -0.039220 0.005336 -7.350 1.98e-13 ***
## factor(time)2005 -0.042063 0.005163 -8.147 3.74e-16 ***
## factor(time)2010 0.020996 0.004935 4.254 2.10e-05 ***
## factor(time)2015 -0.118449 0.004909 -24.131 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 24 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
## convergence code: 0
## Model failed to converge with max|grad| = 0.00234038 (tol = 0.001, component 1)
# Considering all municipalities in the Border States (272) and homicide as cause of death, we have a variance of 1.298 with an Intercept of -11.4935.
summary(borderhom)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(Freq, pop2) ~ factor(age.n) + factor(sex) + factor(time) +
## (1 | clave)
## Data: border.hom
##
## AIC BIC logLik deviance df.resid
## 31661.2 31882.4 -15805.6 31611.2 51567
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.576 -0.189 -0.076 -0.028 119.356
##
## Random effects:
## Groups Name Variance Std.Dev.
## clave (Intercept) 1.298 1.139
## Number of obs: 51592, groups: clave, 272
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.49356 0.19595 -58.656 < 2e-16 ***
## factor(age.n)2 -0.44790 0.20895 -2.144 0.0321 *
## factor(age.n)3 -1.14542 0.22472 -5.097 3.45e-07 ***
## factor(age.n)4 0.04685 0.19459 0.241 0.8097
## factor(age.n)5 2.47361 0.18060 13.696 < 2e-16 ***
## factor(age.n)6 3.06768 0.18000 17.043 < 2e-16 ***
## factor(age.n)7 3.23552 0.17994 17.981 < 2e-16 ***
## factor(age.n)8 3.19969 0.18000 17.776 < 2e-16 ***
## factor(age.n)9 3.03993 0.18019 16.870 < 2e-16 ***
## factor(age.n)10 2.83938 0.18066 15.717 < 2e-16 ***
## factor(age.n)11 2.65316 0.18141 14.625 < 2e-16 ***
## factor(age.n)12 2.37456 0.18281 12.989 < 2e-16 ***
## factor(age.n)13 2.24215 0.18471 12.139 < 2e-16 ***
## factor(age.n)14 1.99809 0.18791 10.633 < 2e-16 ***
## factor(age.n)15 2.00473 0.19076 10.509 < 2e-16 ***
## factor(age.n)16 1.63650 0.20049 8.162 3.28e-16 ***
## factor(age.n)17 1.70413 0.20950 8.134 4.14e-16 ***
## factor(age.n)18 1.41308 0.24421 5.786 7.19e-09 ***
## factor(age.n)19 1.92356 0.23258 8.270 < 2e-16 ***
## factor(sex)2 -2.21959 0.02387 -92.978 < 2e-16 ***
## factor(time)2000 0.16685 0.03694 4.516 6.29e-06 ***
## factor(time)2005 0.22409 0.03596 6.232 4.59e-10 ***
## factor(time)2010 1.87301 0.02935 63.818 < 2e-16 ***
## factor(time)2015 0.80247 0.03203 25.055 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 24 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
## convergence code: 0
## Model failed to converge with max|grad| = 0.0117324 (tol = 0.001, component 1)
# Considering all municipalities in the Border States (270) and homicide committed with firearms as cause of death, we have a variance of 1.633 with an Intercept of -14.0009.
summary(borderfa)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(Freq, pop2) ~ factor(age.n) + factor(sex) + factor(time) +
## (1 | clave)
## Data: border.fa
##
## AIC BIC logLik deviance df.resid
## 24081.8 24302.9 -12015.9 24031.8 51187
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.495 -0.144 -0.052 -0.017 115.038
##
## Random effects:
## Groups Name Variance Std.Dev.
## clave (Intercept) 1.633 1.278
## Number of obs: 51212, groups: clave, 270
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -14.00096 0.79659 -17.576 < 2e-16 ***
## factor(age.n)2 1.16710 0.81382 1.434 0.15155
## factor(age.n)3 0.93546 0.81438 1.149 0.25068
## factor(age.n)4 2.39944 0.79716 3.010 0.00261 **
## factor(age.n)5 4.90962 0.79205 6.199 5.70e-10 ***
## factor(age.n)6 5.54920 0.79185 7.008 2.42e-12 ***
## factor(age.n)7 5.71552 0.79183 7.218 5.27e-13 ***
## factor(age.n)8 5.66101 0.79185 7.149 8.73e-13 ***
## factor(age.n)9 5.49299 0.79191 6.936 4.02e-12 ***
## factor(age.n)10 5.22726 0.79209 6.599 4.13e-11 ***
## factor(age.n)11 4.96860 0.79239 6.270 3.60e-10 ***
## factor(age.n)12 4.62351 0.79307 5.830 5.55e-09 ***
## factor(age.n)13 4.40076 0.79393 5.543 2.97e-08 ***
## factor(age.n)14 4.06525 0.79580 5.108 3.25e-07 ***
## factor(age.n)15 3.98847 0.79752 5.001 5.70e-07 ***
## factor(age.n)16 3.50227 0.80373 4.358 1.32e-05 ***
## factor(age.n)17 3.47986 0.81109 4.290 1.78e-05 ***
## factor(age.n)18 2.77020 0.86082 3.218 0.00129 **
## factor(age.n)19 3.50011 0.83551 4.189 2.80e-05 ***
## factor(sex)2 -2.45268 0.03100 -79.124 < 2e-16 ***
## factor(time)2000 -0.29659 0.04267 -6.952 3.61e-12 ***
## factor(time)2005 -0.25221 0.04146 -6.083 1.18e-09 ***
## factor(time)2010 1.71379 0.03054 56.119 < 2e-16 ***
## factor(time)2015 0.35250 0.03542 9.953 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 24 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
## convergence code: 0
## Model failed to converge with max|grad| = 0.0537465 (tol = 0.001, component 1)
These results suggest that homicides with firearms rates vary more in the municipalities of the Border States than homicide rates and all-causes mortality rates.
The next models consider both random intercepts and random slopes. The random slope for the age structure took a long time, so I am just showing the random slope for sex.
# All-causes of death rates in the Border States
border2<-glmer(cbind(Freq,pop2)~factor(age.n) + factor(sex) + factor(time) + (1+factor(sex)|clave), data = border.all, family = binomial)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00854424 (tol =
## 0.001, component 1)
# Homicide rates in the Border States
borderhom2<-glmer(cbind(Freq,pop2)~factor(age.n) + factor(sex) + factor(time) + (1+factor(sex)|clave), data = border.hom, family = binomial)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00605949 (tol =
## 0.001, component 1)
# Homicide committed with firearms rates in the Border States
borderfa2<-glmer(cbind(Freq,pop2)~factor(age.n) + factor(sex) + factor(time) + (1+factor(sex)|clave), data = border.fa, family = binomial)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.0651098 (tol =
## 0.001, component 1)
summary(border2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(Freq, pop2) ~ factor(age.n) + factor(sex) + factor(time) +
## (1 + factor(sex) | clave)
## Data: border.all
##
## AIC BIC logLik deviance df.resid
## 177241.2 177480.6 -88593.6 177187.2 52324
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -10.442 -0.642 -0.283 0.380 45.418
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## clave (Intercept) 0.14014 0.3744
## factor(sex)2 0.03275 0.1810 -0.25
## Number of obs: 52351, groups: clave, 276
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.402708 0.024635 -178.715 < 2e-16 ***
## factor(age.n)2 -3.183976 0.016165 -196.966 < 2e-16 ***
## factor(age.n)3 -4.045806 0.021460 -188.528 < 2e-16 ***
## factor(age.n)4 -3.878905 0.020202 -192.004 < 2e-16 ***
## factor(age.n)5 -2.816510 0.012983 -216.931 < 2e-16 ***
## factor(age.n)6 -2.417041 0.011273 -214.413 < 2e-16 ***
## factor(age.n)7 -2.221771 0.010932 -203.236 < 2e-16 ***
## factor(age.n)8 -2.071839 0.010637 -194.770 < 2e-16 ***
## factor(age.n)9 -1.874667 0.010254 -182.823 < 2e-16 ***
## factor(age.n)10 -1.599111 0.009911 -161.354 < 2e-16 ***
## factor(age.n)11 -1.220903 0.009452 -129.165 < 2e-16 ***
## factor(age.n)12 -0.841198 0.008973 -93.746 < 2e-16 ***
## factor(age.n)13 -0.358220 0.008539 -41.950 < 2e-16 ***
## factor(age.n)14 0.021692 0.008244 2.631 0.00851 **
## factor(age.n)15 0.446579 0.008062 55.394 < 2e-16 ***
## factor(age.n)16 0.854470 0.007946 107.530 < 2e-16 ***
## factor(age.n)17 1.288606 0.007968 161.722 < 2e-16 ***
## factor(age.n)18 1.752965 0.008125 215.749 < 2e-16 ***
## factor(age.n)19 2.387992 0.007687 310.653 < 2e-16 ***
## factor(sex)2 -0.508730 0.014734 -34.528 < 2e-16 ***
## factor(time)2000 -0.039324 0.005336 -7.370 1.70e-13 ***
## factor(time)2005 -0.042065 0.005163 -8.148 3.71e-16 ***
## factor(time)2010 0.020913 0.004935 4.238 2.26e-05 ***
## factor(time)2015 -0.118565 0.004908 -24.155 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 24 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
## convergence code: 0
## Model failed to converge with max|grad| = 0.00854424 (tol = 0.001, component 1)
summary(borderhom2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(Freq, pop2) ~ factor(age.n) + factor(sex) + factor(time) +
## (1 + factor(sex) | clave)
## Data: border.hom
##
## AIC BIC logLik deviance df.resid
## 31513.8 31752.8 -15729.9 31459.8 51565
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.496 -0.185 -0.075 -0.028 91.656
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## clave (Intercept) 1.3916 1.1797
## factor(sex)2 0.2592 0.5091 -0.83
## Number of obs: 51592, groups: clave, 272
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.50604 0.19575 -58.781 < 2e-16 ***
## factor(age.n)2 -0.45047 0.20761 -2.170 0.03 *
## factor(age.n)3 -1.14838 0.22344 -5.140 2.75e-07 ***
## factor(age.n)4 0.04401 0.19328 0.228 0.82
## factor(age.n)5 2.47072 0.17924 13.784 < 2e-16 ***
## factor(age.n)6 3.06483 0.17863 17.157 < 2e-16 ***
## factor(age.n)7 3.23274 0.17857 18.103 < 2e-16 ***
## factor(age.n)8 3.19675 0.17863 17.896 < 2e-16 ***
## factor(age.n)9 3.03707 0.17883 16.983 < 2e-16 ***
## factor(age.n)10 2.83632 0.17930 15.819 < 2e-16 ***
## factor(age.n)11 2.64996 0.18006 14.717 < 2e-16 ***
## factor(age.n)12 2.37117 0.18146 13.067 < 2e-16 ***
## factor(age.n)13 2.23876 0.18337 12.209 < 2e-16 ***
## factor(age.n)14 1.99413 0.18657 10.688 < 2e-16 ***
## factor(age.n)15 2.00051 0.18945 10.560 < 2e-16 ***
## factor(age.n)16 1.63177 0.19925 8.190 2.62e-16 ***
## factor(age.n)17 1.69913 0.20824 8.159 3.37e-16 ***
## factor(age.n)18 1.40820 0.24302 5.795 6.85e-09 ***
## factor(age.n)19 1.91820 0.23142 8.289 < 2e-16 ***
## factor(sex)2 -2.18731 0.05476 -39.941 < 2e-16 ***
## factor(time)2000 0.16690 0.03692 4.520 6.18e-06 ***
## factor(time)2005 0.22428 0.03594 6.241 4.35e-10 ***
## factor(time)2010 1.87321 0.02933 63.859 < 2e-16 ***
## factor(time)2015 0.80282 0.03201 25.079 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 24 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
## convergence code: 0
## Model failed to converge with max|grad| = 0.00605949 (tol = 0.001, component 1)
summary(borderfa2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(Freq, pop2) ~ factor(age.n) + factor(sex) + factor(time) +
## (1 + factor(sex) | clave)
## Data: border.fa
##
## AIC BIC logLik deviance df.resid
## 24000.8 24239.6 -11973.4 23946.8 51185
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.365 -0.139 -0.050 -0.016 133.041
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## clave (Intercept) 1.6976 1.3029
## factor(sex)2 0.2899 0.5384 -0.60
## Number of obs: 51212, groups: clave, 270
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -13.99965 0.71364 -19.617 < 2e-16 ***
## factor(age.n)2 1.16042 0.73196 1.585 0.112886
## factor(age.n)3 0.93084 0.73196 1.272 0.203478
## factor(age.n)4 2.39438 0.71379 3.354 0.000795 ***
## factor(age.n)5 4.90345 0.70834 6.922 4.44e-12 ***
## factor(age.n)6 5.54318 0.70812 7.828 4.96e-15 ***
## factor(age.n)7 5.70946 0.70810 8.063 7.44e-16 ***
## factor(age.n)8 5.65506 0.70813 7.986 1.40e-15 ***
## factor(age.n)9 5.48699 0.70819 7.748 9.34e-15 ***
## factor(age.n)10 5.22069 0.70841 7.370 1.71e-13 ***
## factor(age.n)11 4.96239 0.70871 7.002 2.52e-12 ***
## factor(age.n)12 4.61725 0.70942 6.508 7.59e-11 ***
## factor(age.n)13 4.39452 0.71036 6.186 6.16e-10 ***
## factor(age.n)14 4.05761 0.71230 5.697 1.22e-08 ***
## factor(age.n)15 3.98089 0.71406 5.575 2.48e-08 ***
## factor(age.n)16 3.49450 0.72183 4.841 1.29e-06 ***
## factor(age.n)17 3.47447 0.72896 4.766 1.88e-06 ***
## factor(age.n)18 2.76843 0.78230 3.539 0.000402 ***
## factor(age.n)19 3.49091 0.75621 4.616 3.91e-06 ***
## factor(sex)2 -2.56511 0.07431 -34.517 < 2e-16 ***
## factor(time)2000 -0.29515 0.04265 -6.920 4.51e-12 ***
## factor(time)2005 -0.25004 0.04144 -6.034 1.60e-09 ***
## factor(time)2010 1.71557 0.03053 56.186 < 2e-16 ***
## factor(time)2015 0.35430 0.03541 10.007 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 24 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
## convergence code: 0
## Model failed to converge with max|grad| = 0.0651098 (tol = 0.001, component 1)
Likelihoo ratio between models and cause of death
# All-causes of death
anova(border, border2)
## Data: border.all
## Models:
## border: cbind(Freq, pop2) ~ factor(age.n) + factor(sex) + factor(time) +
## border: (1 | clave)
## border2: cbind(Freq, pop2) ~ factor(age.n) + factor(sex) + factor(time) +
## border2: (1 + factor(sex) | clave)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## border 25 178282 178504 -89116 178232
## border2 27 177241 177481 -88594 177187 1044.8 2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Homicide as cause of death
anova(borderhom, borderhom2)
## Data: border.hom
## Models:
## borderhom: cbind(Freq, pop2) ~ factor(age.n) + factor(sex) + factor(time) +
## borderhom: (1 | clave)
## borderhom2: cbind(Freq, pop2) ~ factor(age.n) + factor(sex) + factor(time) +
## borderhom2: (1 + factor(sex) | clave)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## borderhom 25 31661 31882 -15806 31611
## borderhom2 27 31514 31753 -15730 31460 151.32 2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Homicides committed with firearms as cause of death
anova(borderfa, borderfa2)
## Data: border.fa
## Models:
## borderfa: cbind(Freq, pop2) ~ factor(age.n) + factor(sex) + factor(time) +
## borderfa: (1 | clave)
## borderfa2: cbind(Freq, pop2) ~ factor(age.n) + factor(sex) + factor(time) +
## borderfa2: (1 + factor(sex) | clave)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## borderfa 25 24082 24303 -12016 24032
## borderfa2 27 24001 24240 -11973 23947 85.035 2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Using the random slope for sex in the models adds some information.
Intra-class correlations suggests that by including the random slope for sex in the second models, we increase the correlation whitin municipalities. Some municipalities are more ‘violent’ than others. (¿?)
library(sjstats)
icc(border); icc(border2)
##
## Generalized linear mixed model
##
## Family : binomial (logit)
## Formula: cbind(Freq, pop2) ~ factor(age.n) + factor(sex) + factor(time) + (1 | clave)
##
## ICC (clave): 0.0387
## Caution! ICC for random-slope-intercept models usually not meaningful. See 'Note' in `?icc`.
##
## Generalized linear mixed model
##
## Family : binomial (logit)
## Formula: cbind(Freq, pop2) ~ factor(age.n) + factor(sex) + factor(time) + (1 + factor(sex) | clave)
##
## ICC (clave): 0.0409
icc(borderhom); icc(borderhom2)
##
## Generalized linear mixed model
##
## Family : binomial (logit)
## Formula: cbind(Freq, pop2) ~ factor(age.n) + factor(sex) + factor(time) + (1 | clave)
##
## ICC (clave): 0.2829
## Caution! ICC for random-slope-intercept models usually not meaningful. See 'Note' in `?icc`.
##
## Generalized linear mixed model
##
## Family : binomial (logit)
## Formula: cbind(Freq, pop2) ~ factor(age.n) + factor(sex) + factor(time) + (1 + factor(sex) | clave)
##
## ICC (clave): 0.2973
icc(borderfa); icc(borderfa2)
##
## Generalized linear mixed model
##
## Family : binomial (logit)
## Formula: cbind(Freq, pop2) ~ factor(age.n) + factor(sex) + factor(time) + (1 | clave)
##
## ICC (clave): 0.3318
## Caution! ICC for random-slope-intercept models usually not meaningful. See 'Note' in `?icc`.
##
## Generalized linear mixed model
##
## Family : binomial (logit)
## Formula: cbind(Freq, pop2) ~ factor(age.n) + factor(sex) + factor(time) + (1 + factor(sex) | clave)
##
## ICC (clave): 0.3404
Finally we plot reggresion lines for each municipality.
# the random effects
rancoefs2<-ranef(border2)
rancoefs22<-ranef(borderhom2)
rancoefs222<-ranef(borderfa2)
Bi-variate plot of each municipality’s random effect for the intercept and the slope for sex. This should show a negative trend for the three models, since in the output the correlation in the random effects was < 0 and stronger for Homicides than for All-causes of death.
par(mfrow=c(1,1))
plot(rancoefs2$clave[,"(Intercept)"],rancoefs2$clave[,"factor(sex)2"], main="Random effects for each municipality", sub = "All causes of death", xlab="Intercept", ylab="Slope - Sex" )
par(mfrow=c(1,1))
plot(rancoefs22$clave[,"(Intercept)"],rancoefs22$clave[,"factor(sex)2"], main="Random effects for each municipality", sub = "Homicides as cause of death", xlab="Intercept", ylab="Slope - Sex" )
par(mfrow=c(1,1))
plot(rancoefs222$clave[,"(Intercept)"],rancoefs222$clave[,"factor(sex)2"], main="Random effects for each municipality", sub = "Homicides with firearm as cause of death", xlab="Intercept", ylab="Slope - Sex" )
Random slopes for each municipality. Intercept of each lines (municipality) is the global intercepto from the second models plus the random intercept component for each municipality. The slope for each line is the average slope - the fixed effect for sex from the second models plus the random slope component for each municipality.
# Border - All-causes of death
plot(NULL, ylim=c(-6, -3), xlim=c(0,1),ylab="ASMR", xlab="sex")
title (main="Regression Lines for each municipality")
cols=sample(rainbow(n=25), size = dim(rancoefs2$clave)[1], replace = T)
for (i in 1: dim(rancoefs2$clave)[1]){
abline(a=fixef(border2)[1]+rancoefs2$clave[[1]][i], b=fixef(border2)[20]+rancoefs2$clave[[2]][i], col=cols[i],lwd=.5 )
}
#this line shows the "average" effect of sex.
abline(a=fixef(border2)[1], b=fixef(border2)[20], col=1, lwd=3)
legend("topright", col=1, lwd=5, legend="Average Effect of sex")
# Border - homicides
plot(NULL, ylim=c(-17, -6), xlim=c(0,1),ylab="ASMR", xlab="sex") #You may have to change the ylim() for your outcome.
title (main="Regression Lines for each municipality")
cols=sample(rainbow(n=25), size = dim(rancoefs22$clave)[1], replace = T)
for (i in 1: dim(rancoefs22$clave)[1]){
abline(a=fixef(borderhom2)[1]+rancoefs22$clave[[1]][i], b=fixef(borderhom2)[20]+rancoefs22$clave[[2]][i], col=cols[i],lwd=.5 )
}
#this line shows the "average" effect of sex
abline(a=fixef(borderhom2)[1], b=fixef(borderhom2)[20], col=1, lwd=3)
legend("topright", col=1, lwd=5, legend="Average Effect of sex")
# Border - homicides with firearms
plot(NULL, ylim=c(-20, -9), xlim=c(0,1),ylab="ASMR", xlab="sex") #You may have to change the ylim() for your outcome.
title (main="Regression Lines for each municipality")
cols=sample(rainbow(n=25), size = dim(rancoefs222$clave)[1], replace = T)
for (i in 1: dim(rancoefs222$clave)[1]){
abline(a=fixef(borderfa2)[1]+rancoefs222$clave[[1]][i], b=fixef(borderfa2)[20]+rancoefs222$clave[[2]][i], col=cols[i],lwd=.5 )
}
#this line shows the "average" effect of sex
abline(a=fixef(borderfa2)[1], b=fixef(borderfa2)[20], col=1, lwd=3)
legend("topright", col=1, lwd=5, legend="Average Effect of sex")
