plot.ritest: A plot method for ritest objects

In grantmcdermott/ritest: Randomisation Inference Testing

Description

Nice plots of your ritest objects.

Usage

## S3 method for class 'ritest'
plot(
  x,
  type = c("density", "hist"),
  highlight = c("lines", "fill", "both", "none"),
  show_parm = FALSE,
  breaks = "auto",
  family = NULL,
  ...
)

Arguments

x	An ritest object.
type	Character. What type of plot do you want?
highlight	Character. How do you want to highlight the H0 rejection regions in the distribution tails?
show_parm	Logical. Should we highlight the parametric H0 rejection regions too?
breaks	Character. Histogram plot only. What type of breaks do you want? The default method creates more breaks than the standard R behaviour. You can revert to the latter by selecting NULL.
family	Character. The font family. Defaults to 'HersheySans' instead of R's normal Arial plotting font.
...	Other plot arguments. Currently ignored.

Examples

#
## Example 1: Basic functionality
#

# First estimate a simple linaer regression on the base 'npk' dataset. For
# this first example, we won't worry about strata or clusters, or other
# experimental design complications.
est = lm(yield ~ N + P + K, data = npk)

# Conduct RI on the 'N' (i.e. nitrogen) coefficient. We'll do 1,000
# simulations and, just for illustration, limit the number of parallel cores
# to 2 (default is half of the available cores). The 'verbose = TRUE'
# argument simply prints the results upon completion, including the original
# regression model summary.
est_ri = ritest(est, 'N', reps = 1e3, seed = 1234L, verbose = TRUE)

# Result: The RI rejection rate (0.021) is very similar to the parametric
# p-value (0.019).

# We can plot the results and various options are available to customise the appearance.
plot(est_ri)
plot(est_ri, type = 'hist')
# etc

# Aside: By default, ritest() conducts a standard two-sided test against a
# sharp null hypothesis of zero. You can can specify other null hypotheses as
# part of the 'resampvar' string argument. For example, a (left) one-sided
# test...
plot(ritest(est, 'N<=0', reps = 1e3, seed = 1234L, pcores = 2L))
# ... or, null values different from zero.
plot(ritest(est, 'N=2', reps = 1e3, seed = 1234L, pcores = 2L))


#
## Example 2: Real-life example
#

# Now that we've seen the basic functionality, here's a more realistic RI
# example using data from a randomized control trial conducted in Colombia.
# More details on the dataset -- kindly provided by the study authors -- can
# be found in the accompanying helpfile ("?colombia"). The most important
# thing to note is that we need to control for the stratified (aka "blocked")
# and clustered experimental design.

data("colombia")

# We'll use the fixest package to estimate our parametric regression model,
# specifying the strata (here: treatment-control pairs) as fixed-effects and
# clustering the standard errors by location (here: city blocks).
library(fixest)
co_est = feols(dayscorab ~ b_treat + b_dayscorab + miss_b_dayscorab |
                 b_pair + round2 + round3,
               vcov = ~b_block, data = colombia)
co_est

# Run RI on the 'b_treat' variable, specifying the strata and clusters.
co_ri = ritest(co_est, 'b_treat', strata='b_pair', cluster='b_block',
               reps=1e3, seed=123L)
co_ri
plot(co_ri, type = 'hist', highlight = 'fill')

# This time, the RI rejection rate (0.11) is noticeably higher than the
# parametric p-value (0.024) from the regression model.

grantmcdermott/ritest documentation built on Jan. 8, 2022, 12:07 a.m.