*! version 3.0.0 July 2023
help randcmdci

Title

randcmdci -- Randomization-t p-values and confidence intervals for OLS regression coefficients that are asymptotically robust to departures from the sharp null.

Syntax

randcmdci cmd depvar indepvars [if] [in] [weight], treatvars(varlist) testvars(varlist) [, options]

Options

Estimation and Conventional Standard Error

noconstant suppress constant term.
absorb(varname) fixed effects for the categories defined by varname.
robust robust conventional standard errors, the default.
cluster(varname) clustered conventional standard errors.

Randomization Inference

strata(varlist) variables identifying randomization strata.
groupvar(varlist) variables identifying treatment groups.
seed(#) set random-number seed to #, default is 1.
reps(#) integer number of attempted randomization iterations, default is 999.
calc1(string) execute stata command after each treatment randomization.
calc2(string) execute stata command after each treatment randomization.
...
calc20(string) execute stata command after each treatment randomization.
null1(string) test a particular null
null2(string) test a particular null
...
null10(string) test a particular null
max(varlist) calculate the maximum p-value for each treatment effect in varlist across all possible nulls for other testvars.
maxlevel(#) the intensity of the search used to find such maxima, default is 100, but higher values should be used to confirm results.
maxcoef(#) maximum bound on the null for other treatment effects based upon deviation from estimated values.
maxwald(#) maximum bound on the null for other treatment effects based upon conventional wald statistic.

cmd must be Stata's regress or areg command.

treatvars indicates the base treatment variables that are permuted/rerandomized on each iteration and is not optional.
treatvars cannot vary within groups identified by groupvar.
Observations for which treatvars are missing are dropped.
testvars indicates the treatment based measures within indepvars in the estimating equation whose effects are tested.
testvars may differ from treatvars when, for example, treatvars are interacted with participant characteristics in the estimating equation (see example 2 below).

To update treatment measures that involve interactions with participant characteristics, up to 20 post-treatment permutation calculations (calc) are allowed and are executed in numerical order.
Permutation of treatment and post-permutation calculations are based on the regression sample alone.

In addition to providing p-values and confidence intervals for individual coefficients, the program also allows the user to call for the testing of up to 10 specific nulls (null).
Numeric values, separated by commas, must be listed for each parameter in testvars (see example 6 below). For equations with multiple treatment measures this is a test
of a joint null.

In equations with multiple testvars, randcmdci reports baseline confidence intervals and p-values for individual treatment effects setting the nulls on
other treatment effects equal to estimated values. The program also allows the user to call for p-values that do not depend upon the nulls for other treatment
effects by calculating the maximum p-value for each individual treatment effect across all possible nulls for other (untested) treatment effects. As this is time
consuming, such maxima are only calculated for the treatment effects requested in max(varlist). maxlevel(#) modifies the intensity of
the search for maxima and either maxcoef(#) or maxwald(#) may be used to put bounds on the nulls considered (further detail below).

If the user does not provide a random number seed, the seed is set to 1 (to ensure replicability, provided the data has not been reordered, on subsequent runs). The seed is restored to its pre-program value on termination.

Description

randcmdci computes randomization confidence intervals and p-values for experiments where permutations of treatvars,
possibly stratified and/or by groups of observations, are valid counterfactual experimental outcomes.
Randomization inference of all forms is exact in finite samples when the sharp null is true (i.e. the treatment effect is the same for each and every observation).
When the sharp null is not true and coefficients estimate the average of heterogeneous treatment effects,
randomization inference based upon studentized test statistics, the randomization-t, has the same asymptotic validity as the clustered/robust covariance estimate provided
(1) treatment is iid (by group if grouping); (2) covariates interacted with treatment are included separately as regressors in their own right;
(3) treatment, regressors and errors have sufficiently high moments;
(4) the maximum number of observations within which treatment, regressors or errors are correlated is bounded;
(5) if stratifying, the first and second moments of treatment are asymptotically identical across strata. For details, see:

Young, Alwyn.
Asymptotically Robust Permutation-Based Randomization Confidence Intervals for Parametric OLS Regression.

randcmdci begins by reporting the conventional Stata estimation results so that the user may confirm the estimating equation is correct.
The sample is automatically restricted to observations for which treatvars are non-missing, but additional sample qualifiers and estimation options are allowed.
randcmdci then rerandomizes (i.e. permutes) treatvars, following any strata and grouping specified by the user, within the regression sample and performs optional calc1 ... calc20 (see examples below).
This sampling from the distribution of potential treatment outcomes is used to construct confidence intervals whose probability of covering the true parameter value
is either exact in finite samples (when the sharp null is true) or otherwise asymptotically accurate (when it is not).

When the regression contains multiple treatment measures (testvars), randomization confidence intervals for each coefficient depend upon the null for others,
although this characteristic vanishes asymptotically. In such cases, randcmdci calculates the confidence interval and p-value for each individual treatment effect under the null
that other treatment effects equal their estimated values, as asymptotically those estimated values typically converge to the true parameters. This attempts to mimic the (impractical) exact
test that sets the nulls on untested measures equal to their true values, but because estimated values are not equal to true parameter values will not necessarily be exact in finite samples.
However, to provide a conservative test that is guaranteed to control size at the desired level if treatment effects are sharp, randcmdci allows the user to call for the calculation of the maximum p-value
for each individual treatment effect across all possible nulls for untested treatment measures. For details on these issues, see the paper above.

In the case of equations with k = 2 treatment effects, the maximum p-value for each treatment effect across all possible nulls for the other can be calculated analytically. In the case
of equations with k = 3 treatment effects, calculating the maximum p-value for each individual treatment effect across all possible nulls for the others can be reduced to searching on the bounded space [0,pi]. randcmdci implements a
maxlevel(#) grid search on this space. In the case of equations with k > 3 treatment effects, finding the maximum can be reduced to searching on the bounded space [0,pi]^(k-2). randcmdci
implements maxlevel(#) Nelder-Mead simplex/amoeba searches on this space, taking the maximum of them all. The program uses a baseline maxlevel(#) of 100 to speed execution.
Before reporting final results, the user should call for a more intensive search for a maximum, e.g. setting maxlevel to 10000 for k = 3 and 1000 for k > 3.

In calculating the maximum p-value of a treatment effect of zero for individual treatment effects across possible nulls for untested treatment measures, some extreme values for other
treatment measures may be considered implausible. The user may use either maxcoef(#) or maxwald(#) to limit the space of possible nulls for other treatment measures across which the maximum p-value
of a treatment effect of zero for an individual treatment effect is calculated (see example 7 below).

Examples

Example 1: One treatment variable, with estimation options.

randcmdci areg outcome treatment covariates [if] [in] [weights], absorb(absorb) robust treatvars(treatment) testvars(treatment)

Example 2: Treatment variables include interactions with participant characteristics which are not part of treatment and have to be recalculated after treatment is rerandomized.

randcmdci reg outcome treatment treatage covariates, treatvars(treatment) testvars(treatment treatage) calc1(replace treatage = treatment*age)

Example 3: Interactions between treatment measures can be coded with calc or without.

randcmdci reg outcome treat1 treat2 treat12 covariates, treatvars(treat1 treat2) testvars(treat1 treat2 treat12) calc1(replace treat12 = treat1*treat2)

randcmdci reg outcome treat1 treat2 treat12 covariates, treatvars(treat1 treat2 treat12) testvars(treat1 treat2 treat12)

Comment: Since treat12 is the product of randomized treatment variables, and does not involve interaction with non-randomized participant characteristics, reallocating pairs of (treat1,treat2) to participants
and then calculating treat12 is equivalent to reallocating triplets of (treat1,treat2,treat12) across participants.

Example 4: Treatment is stratified and applied in groups.

randcmdci areg outcome treat1 treat2 covariates, absorb(village) cluster(village) treatvars(treat1 treat2) testvars(treat1 treat2) strata(strata) groupvar(village)

Comment: treat1 and treat2 cannot vary within groupvar. The pair (treat1, treat2) is rerandomized across villages (within strata) with common values applied to everyone in each village.

Example 5: Testing of specific joint nulls. In example 4, we want to test the joint null that treat1 = 1.5 and treat2 = 1,
as well as the joint null that treat1 = 2 and treat2 = 1:

randcmdci areg outcome treat1 treat2 covariates, absorb(village) cluster(village) treatvars(treat1 treat2) testvars(treat1 treat2) strata(strata) groupvar(village) null1(1.5,1) null2(2,1)

Comment: Numeric values, separated by commas, for the average effect of each variable in testvars must be given. Results report the p-value for the joint test
of both treatment effects at the values given.

Example 6: Maximum p-values across all possible nulls for untested treatment meausures. In example 4, we want to know the maximum p-value of a treatment effect treat1 = 0 across all possible nulls for
treat2 and the maximum pvalue of a treatment effect treat2 = 0 across all possible nulls for treat1. This test is guaranteed to have size less than nominal level (e.g. .05) if treatment
effects are sharp.

randcmdci areg outcome treat1 treat2 covariates, absorb(village) cluster(village) treatvars(treat1 treat2) testvars(treat1 treat2) strata(strata) groupvar(village) max(treat1 treat2)

Comment: Since these calculations are time consuming, it is possible to call for maximum p-values for only a subset of testvars, e.g. max(treat1) or max(treat2).

Example 7: The user wants to limit the range of nulls considered for untested treatment measures in calculating maximum p-values for individual treatment measures.

randcmdci areg outcome treat1 treat2 treat12 covariates, absorb(village) cluster(village) treatvars(treat1 treat2 treat12) testvars(treat1 treat2 treat12) strata(strata) groupvar(village) max(treat1) maxcoef(5)

local c = sqrt(invchi2tail(2,1e-10))

randcmdci areg outcome treat1 treat2 covariates, absorb(village) cluster(village) treatvars(treat1 treat2) testvars(treat1 treat2) strata(strata) groupvar(village) max(treat1 treat2) maxwald(`c')

Comment: Either maxcoef(#) or maxwald(#) may be called, but not both.

When the user selects maxcoef, the range of nulls on untested measures considered is limited to those where the
sum of squared deviations of those treatment measures from estimated values is less that #^2, i.e. restricts the null to lie within a sphere of radius # around the regression point estimates. In the
example above, in calculating the maximum p-value for treat1 the sum of squared deviations of the nulls for treat2 and treat12 from estimated values must be less than 25, while in calculating the
maximum p-value for treat2 the sum of squared deviations of the nulls for treat1 and treat12 from estimated values must be less than 25.

When the user selects maxwald, the wald statistic for the deviation of the nulls on untested measures from estimated values is limited to less than #^2. In the example above, the user
selects a bound such that the p-value of the conventional test of the null for untested measures must be greater than or equal to 10^(-10). This is the default bound when neither
maxcoef or maxwald is specfied by the user.

Reported results

In comparing potential outcomes, randomization inference always produces "ties", as the original experimental outcome (at a minimum) must be considered a tie with itself.
P-values are randomly distributed between minimum and maximum values determined by these ties and this randomization (correctly) determines the boundaries of the confidence interval (see the paper referenced above).
If (1) the number of potential experimental outcomes is large so that there are no ties beyond that of the original experimental outcome and (2) the number of iterations plus 1 times the nominal coverage probability equals an integer,
ties play no role in determining the boundaries of the confidence interval or the significance of coefficients at the given nominal level.
The default number of 999 iterations meets the second of these conditions for the .01, .05, and .10 levels (.99, .95 and .90 confidence intervals).

Randomization confidence intervals may be unbounded (especially when the number of ties is large) and non-convex.
Unbounded upper or lower bounds to confidence intervals are indicated using Stata's symbol for missing (.).
If the confidence interval is not-convex this fact is reported to the user, along with the convex cover of the disjoint confidence regions.

Only "successful" iterations, where coefficients and covariance estimates are both identified, are used in calculating p-values and confidence intervals,
as the reporting of the original experimental result was conditional on the allocation of treatment producing these.

randcmdci stores the following matrices in e():

e(HB) Confidence intervals & p-values under the null that other treatment effects equal their estimated value.
e(JTest) P-values for joint tests, if requested.
e(Pmax) Maximum p-values for tests of individual zero treatment effects across possible nulls for other treatment effects, if requested.

Author

Alwyn Young, London School of Economics, a.young@lse.ac.uk.