*! version 1.0.0 7jul2020

Title

MULTIVRS -- Module to conduct multiverse analysis

Syntax

multivrs estimation_command depvar indepvar controlterm [controlterm ...]
[if] [in] [, options]

where estimation_command defines the estimation command to be executed. The following
Stata commands are supported by multivrs:

regress, logit, logistic, probit,
poisson, nbreg, areg, rreg, xtreg

indepvar is the independent variable of interest. All of the estimates reported by multivrs are performed
with respect to indepvar.

Description

multivrs allows users to conduct a ‘many worlds’ or multiverse analysis of modeling assumptions.
Users specify a preferred model, as well as alternative estimation commands, control variables, variable coding,
and standard error calculations that reflect different plausible assumptions. The model specification set is all possible combinations of these inputs.
After estimating each specification, results are reported as a distribution of estimates.

multivrs also conducts an analysis of which model ingredients are most influential for the results.

multivrs can report marginal effects from nonlinear models, allow different treatments for missing data,
and can weight the modeling distribution by different criteria of model fit.

multivrs depends upon the following packages:

moremata (Ben Jann, 2005. "MOREMATA: Stata module (Mata) to provide various functions," Statistical Software Components S455001, Boston College Department of Economics, revised 06 Dec 2020.)

valuesof (Ben Jann, 2006. "VALUESOF: Stata module to return the contents of a variable in a macro," Statistical Software Components S456705, Boston College Department of Economics, revised 26 May 2008.)

Options

Grouping and Either|Or

Grouping is specified such that each controlterm is either varname or (varlist). A varlist in parentheses indicates that this group of variables is to be included or excluded together.
This is useful for categorical variables. varlist can be grouped with indepvar if a priori beliefs require varlist to be included as a control in all models.

Either|Or is specified such that a varlist in the form (d1|d2|...) treats d1, d2. etc. as incompatible terms where both cannot be included in the same model.
This option can be used when there is uncertainty over the functional form of varname (x|x^2|ln(x)) or if there are highly collinear measures of the same latent variable.

Both grouping and either|or syntax can be specified for even more flexibility. For example, ((x1|x1^2) x2) specifies that x1 and x1^2 will not be included in the same model but x2 will be included in all models with either.

Functional Form Robustness

Either|Or syntax can be used within estimation_command to pool results across different functional forms.
An estimation_command of the form (m1|m2|...) will generate pooled results from running multivrs
over each of the model types m1, m2, etc. For example, for binary data one could use (reg | logit | probit),
or for count data (reg | poisson | nbreg). Options specific to one functional form can be specified as (m1(m1_options)|m2(m2_options)|...).

Special Regression Options

vce(vcetype) specifies the standard error calculations. vcetype may be robust or cluster clustervar.
See vce_option.

irr displays the beta coefficients as incidence rate ratios for poisson or negative binomial regression models.

exposure(varname) corrects for unequal exposure times to varname for observations in poisson or negative binomial regression models; specifies that ln(varname) be included in the model with coefficient constrained to 1 in the log-link function

offset(varname) specifies that varname be included in the poisson or negative binomial regression models with the coefficient constrained to be 1.

absorb(varname) for areg, categorical variable to be absorbed. This option is required when areg is in the estimation_command.

Other Estimation Options

margins estimates average marginal effects (dydx) for the variable of interest, saves these estimates (and the associated standard errors and p-values)
rather than the regression coefficients, and summarizes and plots them. See margins.

nolistwise turns off the default listwise delete based on the full list of variables, and instead executes listwise delete for each model separately.

other(other_option) permits any syntax allowed as an option for the estimation command.

Model Space

sample(percent) randomly samples the given percentage of possible models. Percent may be an integer 1-99.

size(#) constrains the number of covariates in the models according to one of the following specifications:

size(#) exactly # control terms

size(#min, #max) between #min and #max control terms inclusive

size(#min, .) greater than or equal to #min control terms

size(., #max) less than or equal to #max control terms

For example, for k = 25, specifying size(4,8) would reduce the model space from 33,554,432 to 1,805,155, specifying size(.,7) would reduce the model space to 726,206,
and specifying size(7) would reduce the model space to 480,700.

Weighting of Estimates

weights(weight_type) specifies the weighting method for summarizing the model results (i.e., calculating means and standard errors).
The default is uniform which assigns
equal weight to each model. Alternatively, weight_type may be bic for BIC weighting (an approximation to the posterior probabilities)
or inf for weighting by model influence, or r2 for weighting by the R2. By default each of these 3 weights is saved in the results file.
If the weight_type is specified as no then no weights are calculated.

Reporting

alpha(#) sets # as the significance threshold. Coefficients with p-value <= # are categorized as statistically significant. Default alpha = .05.

pref(estimate, SE) plots the user's preferred estimate as a red line in the distribution of coefficients,
and calculates model robustness statistics with the SE associated with the preferred estimate.

sig_only suppresses the coefficient summary statistics to only report sign and significance statistics.
This option is often used along with the either|or syntax for estimation_command to pool results across different functional forms.
When sig_only is specified, the influence statistics are reported as marginal effect on probability of significant (or positive) coefficient.

nosig overrides the default switch to sig_only when different functional forms
are specified, allowing full standard output.

noinfluence suppresses the control variable influence statistics.

inf_means for each control variable considered, displays the mean coefficient estimate
from all models including that variable, in addition to the standard influence statistics.

normal overlays a normal distribution on the graph of the empirical modeling distribution.

nozero relaxes the default constraint that the graph includes zero.

noplot suppresses the plot of the coefficient of interest.

Saved Results

saveas(file_name) saves the model results in file_name.dta, and the list of estimated models in file_name.do. Default is for no results to be saved.
The .dta file contains coefficient estimates for the variable of interest from each of the models run. The .do file provides an easy way to see which models the program ran.

replace permits saveas to overwrite existing file_name.dta and/or file_name.do.

Advanced

bs(bs_type) performs resampling of the data or estimate to generate a sampling distribution of size B = 50
for the parameter of interest for each of the J models. These B*J bootstrapped estimates are saved in the results file, composing the
total sampling+modeling distribution, which is used to calculate the total SE and the robustness intervals.
bs_type may be par (parametric) or nonpar (nonparametric).

When bs_type = par, the B estimates for each model come from the normal distribution with
mean = point estimate of the parameter of interest for that model and se = SE of the point estimate.

When bs_type = nonpar(bs_options), Stata's bootstrap command is executed on each model,
resampling the actual data points and re-estimating the model B times to generate B parameter estimates. bs_options
are available to accomodate data-specific resampling requirements such as strata(varlist) or cluster(varlist)), and
are passed in directly to bootstrap on each iteration.

intervals prints the modeling distribution 95% confidence interval cut points and extreme bounds,
and 95% and parametric intervals for the full bootstrapped sampling+modeling distribution.
If bs(bs_type) is not specified,
the total distribution is generated by the parametric bootstrap.

Examples: Basic

Setup

. sysuse nlsw88

Use OLS

. multivrs regress union hours age grade collgrad married south smsa c_city ttl_exp tenure

Use logit model

. multivrs logit union hours age grade collgrad married south smsa c_city ttl_exp tenure

Use either|or syntax to assess functional form robustness

. multivrs (logit|probit|regress) union hours age grade collgrad married south smsa c_city ttl_exp tenure

Treat the two education variables (highest grade completed, and dummy for college graduation) as a single term

. multivrs logit union hours age (grade collgrad) married south smsa c_city ttl_exp tenure

Keep age and tenure as covariates in every model, along with the variable of interest (hours)

. multivrs logit union (hours age tenure) grade collgrad married south smsa c_city ttl_exp

Examples: Special Regression Options

Use robust standard errors

. multivrs regress union hours age grade collgrad married south smsa c_city ttl_exp tenure, vce(robust)

Convert poisson or nbreg coefficients to incidence rate ratios

. multivrs poisson tenure union hours age grade south smsa c_city ttl_exp, irr

Correct for unequal exposure times measured by the variable ttl_exp

. multivrs poisson tenure union hours age grade south smsa c_city, exposure(ttl_exp)

Examples: Other Estimation Options

Estimate and plot marginal effects to directly compare effects across different functional forms

. multivrs (logit|probit|regress) union hours age grade smsa c_city ttl_exp tenure, margins

Execute listwise delete separately within each model, which allows the sample size to vary across models

. multivrs regress wage age hours union grade smsa c_city ttl_exp tenure, nolistwise

For rreg, specify a biweight tuning parameter of 9

. multivrs rreg wage age hours union grade smsa c_city ttl_exp tenure, other(tune(9))

Examples: Weighting of Estimates

Weight each estimate using BIC weights

. multivrs regress union hours south smsa c_city ttl_exp tenure, weights(bic)

Weight each estimate by the R2 value of its model

. multivrs regress union hours south smsa c_city ttl_exp tenure, weights(r2)

Weight each estimate by the absolute value of its influence regression coefficient

. multivrs regress union hours south smsa c_city ttl_exp tenure, weights(inf)

Examples: Model Space

Sample 50 percent of the possible models (usually used for very long lists of control variables)

. multivrs regress union hours age grade collgrad married south smsa c_city ttl_exp tenure, sample(50)

Restrict the model space to models with exactly 4 total covariates

. multivrs regress union hours age grade collgrad married south smsa c_city ttl_exp tenure, size(4)

Restrict the model space to models with at least 4 total covariates

. multivrs regress union hours age grade collgrad married south smsa c_city ttl_exp tenure, size(4, .)

Restrict the model space to models with no more than 4 total covariates

. multivrs regress union hours age grade collgrad married south smsa c_city ttl_exp tenure, size(., 4)

Examples: Reporting

Set the alpha (significance threshold) to 0.1 rather than the default 0.05

. multivrs regress wage age hours union grade smsa c_city ttl_exp tenure, alpha(0.1)

Overlay a normal distribution on the plot

. multivrs regress wage age hours union grade smsa c_city ttl_exp tenure, normal

Compare multiverse results to a preferred estimate (0.003) and SE (0.001)

. multivrs regress union hours age grade collgrad married south smsa c_city ttl_exp tenure, pref(0.003 0.001)

Replay results of significance testing broken down by functional form

. multivrs, more

Replay full results including coefficient summary statistics (regardless of prior sig_only specification)

. multivrs, full

Examples: Multiple Options

Use logit model, robust standard errors, and 50 percent sampling

. multivrs logit union hours age grade collgrad married south smsa c_city ttl_exp tenure, vce(robust) sample(50)

Use either logit or regress, estimate marginal effects, and use BIC weights

. multivrs (logit|regress) union hours age grade collgrad ttl_exp tenure, margins weights(bic)

Restrict to models with at least 3 independent variables and execute listwise delete separately for each model

. multivrs regress wage age hours union grade smsa c_city ttl_exp tenure, size(3,.) nolistwise

Saved results

multivrs saves the following in r():

Scalars

r(N) number of observations

r(Nmin) minimum number of observations across all models (equal to r(N) under the default listwise delete)

r(Nmax) maximum number of observations across all models (equal to r(N) under the default listwise delete)

r(nvars) number of independent variables

r(nterms) number of grouped control terms

r(nmodels) number of models calculated

r(meanr2) simple average of the R-squared values of the models

r(meanb) simple average of the estimates

r(samplingSE) mean of the standard errors for the estimates

r(modelingSE) standard deviation of the estimates

r(totalSE) standard error of the sampling+modeling distribution

r(rratio) mean estimate divided by total standard error

r(skewness) skew of the distribution of estimates

r(kurtosis) kurtosis of the distribution of estimates

r(extremeUB) maximum of the modeling distribution

r(extremeLB) minimum of the modeling distribution

r(modeling95UB) upper bound for 95% interval (97.5 percentile) for modeling distribution

r(modeling95LB) lower bound for 95% interval (2.5 percentile) for modeling distribution

r(total95UB) upper bound for 95% interval (97.5 percentile) for sampling+modeling distribution

r(total95LB) lower bound for 95% interval (2.5 percentile) for sampling+modeling distribution

r(totalParUB) upper bound for parametric interval for sampling + modeling distribution (meanb + 2*totalSE)

r(totalParLB) lower bound for parametric interval for sampling + modeling distribution (meanb - 2*totalSE)

r(sdintvar) standard deviation of the variable of interest (the observations of the variable, not its modeling distribution)

r(mc) multicollinearity of variable of interest

r(prefb) preferred estimate

r(prefse) preferred sampling standard error

r(preftotalse) standard error of the sampling+modeling distribution using prefse as the sampling standard error

r(prefpctile) percentile of modeling distribution for preferred estimate

r(prefUB) upper bound for parametric interval centered on preferred estimate (prefb + 2*preftotalSE)

r(prefLB) lower bound for parametric interval centered on preferred estimate (prefb - 2*preftotalSE)

Macros

r(cmd) multivrs

r(model) estimation command(s) specified

r(modeli_name) the ith estimation command specified

r(modeli_opts) options specified for the ith estimation command

r(depvar) name of the outcome variable(s)

r(varint) name of the interest variable(s)

r(title) type(s) of regression model

r(opts_multivrs) multivrs options specified

r(opts_command) options specified to pass in to every estimation command

Matrices

r(infstats_sig) influence statistics for sign and significance

r(infstats_coef) influence statistics for value of estimate

r(sigrates) significance and sign rates

Additional items not listed here may appear among the saved results.
These results are recorded for use by the multivrs program in case the user chooses to replay the results of the command.

Authors

Cristobal Young (Cornell University) and Katherine Holsteen (Stanford University)

cristobal.young@cornell.edu