Modeling choices

library(tidyverse)
library(targets)
library(broom)
library(janitor)
library(brms)
library(lme4)
library(lmtest)  # For log ratio tests
library(performance)  # For ICCs
library(scales)
library(modelsummary)
library(here)

CHAINS <- 4
ITER <- 2000
WARMUP <- 1000
BAYES_SEED <- 1234

options(mc.cores = parallel::detectCores(),
        brms.backend = "cmdstanr")

# Generated via random.org
set.seed(9936)

# Load data
# Need to use this withr thing because tar_read() and tar_load() need to see the
# _targets folder in the current directory, but this .Rmd file is in a subfolder
withr::with_dir(here::here(), {
  source(tar_read(plot_funs))
  
  canary <- tar_read(panel_lagged)
})

# canary <- canary %>%
#   mutate(gwcode = as.factor(gwcode), year = as.factor(year)) %>%  # Make a factor
#   mutate(PTS_factor_lead1 = ordered(PTS_factor_lead1)) %>%  # Make this ordered
#   # To do fancy Bell and Jones adjustment (https://doi.org/10.1017/psrm.2014.7)
#   # (aka Mundlak devices), explanatory variables should be split into a meaned
#   # version (\bar{x}) and a de-meaned version (x - \bar{x}) so that they can
#   # explain the within-group and between-group variation. We de-mean three variables here:
#   group_by(gwcode) %>%
#   mutate_at(vars(v2csreprss, v2x_polyarchy, gdpcap_log),
#             list(between = ~mean(., na.rm = TRUE),  # Between, for countries
#                  within = ~. - mean(., na.rm = TRUE))) %>%   # Within, for countries
#   ungroup()

Interpretation reminders

I always forget how to interpret these dang within/between coefficients in REWB models, so here’s a reminder:

• “Within” coefficients show what happens as you move up and down from the average level.
  • Template = increasing 1 unit of x above the average in a typical country is associated with a β change in y
• “Between” coefficients show what happens as you move the whole average up. They’re trickier (and less important) to interpret. This post explains more, along with Bell, Fairbrother, and Jones (2018).

Check if multilevel models are warranted

First we can start with a simple null model to see what the effect of country and year is on our outcome (i.e. how much do country and year alone explain the variation in repression?)

m0_pts <- lmer(PTS_lead1 ~ (1 | year) + (1 | gwcode),
               REML = FALSE, data = canary)
m0_clphy <- lmer(v2x_clphy_lead1 ~ (1 | year) + (1 | gwcode),
               REML = FALSE, data = canary)
m0_clpriv <- lmer(v2x_clpriv_lead1 ~ (1 | year) + (1 | gwcode),
               REML = FALSE, data = canary)

We then calculate the intraclass correlation coefficients (ICC) for year and country:

# This shows just the ICC. modelsummary() helpfully includes it.
# icc(m0_pts)

modelsummary(list(m0_pts, m0_clphy, m0_clpriv),
             estimate = "{estimate}", statistic = NULL,
             gof_map = list(list("raw" = "icc", "clean" = "ICC", "fmt" = 4)))

	Model 1	Model 2	Model 3
(Intercept)	2.610	0.639	0.687
ICC	0.7272	0.8884	0.9315

These values represent the variation in the outcome variable that can be explained by the country/year structure, which is a ton! For the PTS model, for instance, the ICC is 0.7272, meaning that the random effects structure explains 73% of the variation in PTS.

We can also look at how each of the levels are related to the outcome

icc(m0_pts, by_group = TRUE)

Group	ICC
gwcode	0.7247164
year	0.0024886

icc(m0_clphy, by_group = TRUE)

Group	ICC
gwcode	0.8767427
year	0.0116359

icc(m0_clpriv, by_group = TRUE)

Group	ICC
gwcode	0.9207071
year	0.0107510

It looks like pretty much all the variability in PTS is explained at a country level and not at a year level, which makes sense, since there’s probably a ton of autocorrelation in human rights practices within a country, while global human rights levels probably don’t move simultaneously across years.

Just because the bulk of the relative variation is at at the country level doesn’t mean that there isn’t significant variability within countries over time. However, the lagged DV seems to soak up pretty much all year-based variability, so we should really use a multilevel model because the country structure matters a lot when explaining human rights abuses. There aren’t really too many gains in efficiency or predictive power when including (1 | year) though.

Gradually build up the model

Understanding all the different country/year dynamics can get complicated, so John Poe suggests slowly building up the model and inspecting it to see what’s going on.

# First we add the main variable of interest (civil society repression here)
m1 <- lmer(PTS_lead1 ~ v2csreprss + (1 | year) + (1 | gwcode),
           REML = FALSE, data = canary)
modelsummary(list("Null model" = m1),
             coef_omit = "sd__")

	Null model
(Intercept)	2.834
	(0.070)
v2csreprss	-0.265
	(0.017)
Num.Obs.	3779
R2 Marg.	0.120
R2 Cond.	0.707
AIC	7401.4
BIC	7432.6
ICC	0.7

On average, increasing civil society repression by a point leads to a decrease of 0.265 points of human rights repression in the following year.

Finally, we can see if adding v2csreprss is an improvement over the null m0_pts model:

lrtest(m0_pts, m1) %>% as_tibble()

#Df	LogLik	Df	Chisq	Pr(>Chisq)
4	-3806.965
5	-3695.701	1	222.5279	0

It is! m1 is definitely an improvement over m0_pts (p < 0.001).

Next, we can split civil society repression into its within and between versions and see how that performs and how it changes the results

m2 <- lmer(PTS_lead1 ~ v2csreprss_within + v2csreprss_between +
             (1 | year) + (1 | gwcode),
           REML = FALSE, data = canary)

modelsummary(list("Null model" = m1, "Within/between model" = m2),
             coef_omit = "sd__")

	Null model	Within/between model
(Intercept)	2.834	2.978
	(0.070)	(0.077)
v2csreprss	-0.265
	(0.017)
v2csreprss_within		-0.242
		(0.018)
v2csreprss_between		-0.435
		(0.048)
AIC	7401.4	7389.8
BIC	7432.6	7427.2
Log.Lik.	-3695.701	-3688.881

On average, a one-unit increase beyond a country’s average level of civil society repression (i.e. the within-country effect) is associated with a 0.242 point decrease in overall repression the following year

Also, the v2csreprss_between coefficient can be interpreted somehow, but I’m not 100% sure how. (I think something like “as a country’s average civil society repression increases, general repression worsens by 0.435” or something).

And we can check if m2 is an improvement over m1 (it is!):

lrtest(m1, m2) %>% as_tibble()

#Df	LogLik	Df	Chisq	Pr(>Chisq)
5	-3695.701
6	-3688.881	1	13.63967	0.0002215

Thus…

So, using a Mundlak within/between split of time-varying covariates is warranted, as is using a country/year random effects structure. Woot.

Copyright © 2020 Suparna Chaudhry and Andrew Heiss. Licensed under the MIT License. Code available on GitHub.