Modeling choices

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

# 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)
})

# Goodness-of-fit measures to include
gm <- tribble(
  ~raw, ~clean, ~fmt,
  "nobs", "Num.Obs.", 3,
  "icc", "ICC", 3
)

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_lhr <- lmer(latent_hr_mean_lead1 ~ (1 | year) + (1 | gwcode),
               REML = FALSE, data = canary)

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

modelsummary(list(m0_pts, m0_lhr),
             estimate = "{estimate}", statistic = NULL,
             group = group + term ~ model,
             gof_map = gm)

		Model 1	Model 2
	(Intercept)	2.610	0.246
gwcode	SD (Intercept)	0.986	1.266
year		0.058	0.208
Residual	SD (Observations)	0.605	0.486
	Num.Obs.	3779.000	3794.000
	ICC	0.727	0.875

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

performance::icc(m0_pts, by_group = TRUE)

Group	ICC
gwcode	0.7247167
year	0.0024885

performance::icc(m0_lhr, by_group = TRUE)

Group	ICC
gwcode	0.8517684
year	0.0229567

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),
             group = group + term ~ model,
             coef_omit = "sd__", gof_map = gm)

		Null model
	(Intercept)	2.834
		(0.070)
	v2csreprss	−0.265
		(0.017)
gwcode	SD (Intercept)	0.831
year		0.081
Residual	SD (Observations)	0.590
	Num.Obs.	3779.000
	ICC	0.667

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),
             group = group + term ~ model,
             coef_omit = "sd__", gof_map = gm)

	Null model	Within/between model
(Intercept)	2.834	2.978	group
	(0.070)	(0.077)
v2csreprss	−0.265
	(0.017)
v2csreprss_within		−0.242
		(0.018)
v2csreprss_between		−0.435
		(0.048)

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.

BUT for the sake of computational efficiency, we omit the within/between split. It takes 30-120 minutes to run one within/between model, and running dozens can take all day. The predictive results are nearly identical to models without the split, and since we’re more interested in prediction than in individual coefficient effects, we don’t need to wait hours for the same results we’d get with the original non-demeaned values.

Interpretation reminders

NB: We ultimately don’t use this within/between difference since the random effects within/between models take forever to run and provide no improvement over the more basic models, but I’ll keep this here as a reminder to future me in future projects:

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).

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