In mwt/rART: Approximate Randomization Tests with a Small Number of Clusters

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(rART)
set.seed(123)
data("mqy2015")

This example is from Section 4.1 of "A User's Guide to Approximate Randomization Tests with a Small Number of Clusters". It is a replication of Meng, Qian and Yared (2015) which studies the effect of changes in the grain supply and mortality in China during the great famine.

Create dummy variables

We start by making two dummy variables. The first is for the years of the famine, which spanned from 1958-1960. The second is a subset of years used in the original paper to restrict the regression.

# famine years
mqy2015$famine <- as.integer((mqy2015$year >= 1958) & (mqy2015$year <= 1960))

# subset of years for analysis 2, 4, and 6
mqy2015$a2dumm <- (mqy2015$year >= 1953) & (mqy2015$year <= 1965)

Regressions

We start with the first specification as it is the simplest.

Specification 1

We run the first specification using the artlm function.

spreg1 <- artlm(
  ldeath_b ~ lgrain_pred:famine + lgrain_pred + lurbpop + ltotpop + year,
  cluster = province,
  data = mqy2015,
  subset = mqy2015$main,
  select = "lgrain_pred:famine"
)

We can get the p-value using summary

summary(spreg1)

and the confidence intervals using confint.

confint(spreg1)

Note that the numbers are exactly the same as in Table 1 of the paper.

Run all six specifications

spregs <- list(
  # Specification 1
  artlm(
    ldeath_b ~ lgrain_pred:famine + lgrain_pred + lurbpop + ltotpop + year,
    cluster = province,
    data = mqy2015,
    subset = mqy2015$main,
    select = "lgrain_pred:famine"
  ),
  # Specification 2
  artlm(
    ldeath_b ~ lgrain_pred:famine + lgrain_pred + lurbpop + ltotpop + year,
    cluster = province,
    data = mqy2015,
    subset = mqy2015$main & mqy2015$a2dumm,
    select = "lgrain_pred:famine"
  ),
  # Specification 3
  artlm(
    ldeath_b ~ lgrain_pred:famine + lgrain_pred + lurbpop + ltotpop + year,
    cluster = province,
    data = mqy2015,
    select = "lgrain_pred:famine"
  ),
  # Specification 4
  artlm(
    ldeath_b ~ lgrain_pred:famine + lgrain_pred + lurbpop + ltotpop + year,
    cluster = province,
    data = mqy2015,
    subset = mqy2015$a2dumm,
    select = "lgrain_pred:famine"
  ),
  # Specification 5
  artlm(
    ldeath_b ~ lgrain:famine + lgrain + lurbpop + ltotpop + year,
    cluster = province,
    data = mqy2015,
    subset = mqy2015$main,
    select = "lgrain:famine"
  ),
  # Specification 6
  artlm(
    ldeath_b ~ lgrain:famine + lgrain + lurbpop + ltotpop + year,
    cluster = province,
    data = mqy2015,
    subset = mqy2015$main & mqy2015$a2dumm,
    select = "lgrain:famine"
  )
)

Check the values

# get p-values (fourth column of summary)
p_values <- vapply(spregs, function(z) {summary(z)$coefficients[,4L]}, 0)

# get confidence intervals
confidence_intervals <- vapply(spregs, confint, numeric(2L))

knitr::kable(cbind("Pr(>|t|)" = p_values, t(confidence_intervals)))