test_msel function

Tests and confidence intervals for the parameters estimated by the msel function

This function conducts various statistical tests and calculates confidence intervals for the parameters of the model estimated via the msel function.

test_msel(
  object,
  fn,
  fn_args = list(),
  test = "t",
  method = "classic",
  ci = "classic",
  cl = 0.95,
  se_type = "dm",
  trim = 0,
  vcov = object$cov,
  iter = 100,
  generator = rnorm,
  bootstrap = NULL,
  par_ind = 1:object$control_lnL$n_par,
  eps = max(1e-04, sqrt(.Machine$double.eps) * 10),
  n_sim = 1000,
  n_cores = 1
)

Arguments

• object: an object of class 'msel'. It also may be a list of two objects. Then object[[1]] and object[[2]] are supplied to the arguments model1 and model2 of the lrtest_msel function.

• fn: a function which returns a numeric vector and should depend on the elements of object. These elements should be accessed via coef.msel or predict.msel functions. The first argument of fn should be an object. Therefore coef and predict functions in fn should also depend on object.

• fn_args: a list of additional arguments of fn.

• test: a character representing the test to be used. If test = "t" then t-test is used. If test = "wald" then Wald test is applied.

• method: a character representing a method used to conduct a test. If test = "t" or test = "wald" and method = "classic"

  then p-values are calculated by using the quantiles of the standard normal distribution. If test = "t" or test = "wald" and method = "bootstrap"

  then p-values are calculated by using the bootstrap as described in Hansen (2022). If test = "wald" and method = "score" then score bootstrap Wald test of P. Kline and A. Santos (2012) is used.

• ci: a character representing the type of the confidence interval used. Available only if test = "t". If ci = "classic" then quantiles of the standard normal distribution are used to build an asymptotic confidence interval. If ci = "percentile" then percentile bootstrap interval is applied. If ci = "bc" then the function constructs a bias-corrected percentile bootstrap confidence interval of Efron (1982) as described in Hansen (2022).

• cl: a numeric value between 0 and 1 representing a confidence level of the confidence interval.

• se_type: a character representing a method used to estimate the standard errors of the outputs of fn. If se_type = "dm" then delta method is used. If se_type = "bootstrap" then bootstrap is applied.

• trim: a numeric value between 0 and 1 representing the share of bootstrap estimates to be nullified when standard errors are estimated for se_type = "bootstrap".

• vcov: an estimate of the asymptotic covariance matrix of the parameters of the model.

• iter: the number of iterations used by the score bootstrap Wald test.

• generator: function which is used by the score bootstrap to generate random weights. It should have an argument n representing the number of random weights to generate. Other arguments are ignored.

• bootstrap: an object of class 'bootstrap_msel' which is an output of the bootstrap_msel function. This object is used to retrieve the estimates of the bootstrap samples.

• par_ind: a vector of indexes of the model parameters used in the calculation of fn. If only necessary indexes are included then in some cases estimation time may greatly decrease. Set ind = TRUE

  in summary.msel to see the indexes of the model parameters. If eps is a vector then eps[i] determines the increment used to differentiate fn respect to the parameter with par_ind[i]-th index.

• eps: a positive numeric value representing the increment used for the numeric differentiation of fn. It may also be a numeric vector such that eps[i] is an increment used to differentiate the fn respect to the par_ind[i]-th parameter of the model. Set ind = TRUE in summary.msel, to see the indexes of the model parameters. If eps[i] = 0 then derivative of fn respect to par_ind[i]-th parameter is assumed to be zero.

• n_sim: the value passed to the n_sim argument of the msel function.

• n_cores: the value passed to the n_cores argument of the msel function.

Returns

This function returns an object of class 'test_msel' which is a list. It may have the following elements:

• tbl - a list with the elements described below.
• is_bootstrap - a logical value which equals TRUE if bootstrap has been used.
• is_ci - a logical value which equals TRUE if confidence intervals were used.
• test - the same as the input argument test.
• method - the same as the input argument method.
• se_type - the same as the input argument method.
• ci - the same as the input argument ci.
• cl - the same as the input argument cl.
• iter - the same as the input argument iter.
• n_bootstrap - an integer representing the number of the bootstrap iterations used.
• n_val - the length of the vector returned by fn.

A list tbl may have the following elements:

• val - an output of the fn function.
• se - a numeric vector such that se[i] represents a standard error associated with val[i].
• p_value - a numeric vector of p-values.
• lwr - a numeric vector such that lwr[i] is the realization of the lower (left) bound of the confidence interval for the true value of val[i].
• upr - a numeric vector such that upr[i] is the realization of the upper (right) bound of the confidence interval for the true value of val[i].
• stat - a numeric vector of values of the test statistics.

An object of class 'test_msel' has an implementation of the summary method summary.test_msel.

In a special case when object is a list of length 2 the function returns an object of class 'lrtest_msel' since the function lrtest_msel is called internally.

Details

Suppose that $\theta$ is a vector of parameters of the model estimated via the msel function and $g(\theta)$ is a differentiable function representing fn which returns a $m$-dimensional vector of real values:

$$g(\theta) = (g_{1}(\theta),...,g_{m}(\theta))^{T}.$$

Classic and bootstrap t-test

If test = "t" then for each $j\in \{1,...,m\}$ the following hypotheses is tested:

$$H_{0}:g_{j}(\theta) = 0,\qquad H_{1}:g_{j}(\theta)\ne 0.$$

The test statistic is:

$$T = g_{j}(\hat{\theta})/\hat{\sigma}_{j},$$

where $\hat{\sigma}$ is a standard error of $g_{j}(\hat{\theta})$.

If se_type = "dm" then delta method is used to estimate this standard error:

$$\hat{\sigma}_{j} = \sqrt{\nabla g_{j}(\hat{\theta})^{T}\widehat{As.Cov}(\hat{\theta})\nabla g_{j}(\hat{\theta})},$$

where $\nabla g_{j}(\hat{\theta})$ is a gradient as a column vector and the estimate of the asymptotic covariance matrix of the estimates $\widehat{As.Cov}(\hat{\theta})$ is provided via the vcov

argument. Numeric differentiation is used to calculate $\nabla g_{j}(\hat{\theta})$.

If se_type = "bootstrap" then bootstrap is applied to estimate the standard error:

$$\hat{\sigma}_{j} = \sqrt{\frac{1}{B - 1}\sum\limits_{b = 1}^{B}(g_{j}(\hat{\theta}^{(b)}) - \overline{g_{j}(\hat{\theta}^{(b)}}))^2},$$

where $B$ is the number of the bootstrap iterations bootstrap$iter, $\hat{\theta}^{(b)}$ is the estimate associated with the $b$-th of these iterations bootstrap$par[b, ], and $g_{j}(\overline{\hat{\theta}^{(b)}})$ is a sample mean of the bootstrap estimates:

$$\overline{g_{j}(\hat{\theta}^{(b)}}) =\frac{1}{B}\sum\limits_{b = 1}^{B}g_{j}(\hat{\theta}^{(b)}).$$

If method = "classic" it is assumed that if the null hypothesis is true then the asymptotic distribution of the test statistic is standard normal. This distribution is used for the calculation of the p-value:

$$p-value = 2\min(\Phi(T), 1 - \Phi(T)),$$

where $\Phi()$ is a cumulative distribution function of the standard normal distribution.

If method = "bootstrap" then p-value is calculated via the bootstrap as suggested by Hansen (2022):

$$p-value = \frac{1}{B}\sum\limits_{b=1}^{B} I(|T_{b} - T| > |T|),$$

where $T_{b} = g_{j}(\hat{\theta}^{(b)})/\hat{\sigma}_{j}$ is the value of the test statistic estimated on the b-th bootstrap sample and $I(q)$ is an indicator function which equals $1$ when $q$ is a true statement and $0$ - otherwise.

Classic and bootstrap Wald test

Suppose that method = "classic" or method = "bootstrap". If test = "wald" then the null hypothesis is:

$$H_{0}:\begin{cases}g_{1}(\theta) = 0\\g_{2}(\theta) = 0\\\vdots\\g_{m}(\theta) = 0\\\end{cases}.$$

The alternative hypothesis is that there is such $j\in\{1,...,m\}$ that:

$$H_{1}:g_{j}(\theta)\ne 0.$$

The test statistic is:

$$T = g(\hat{\theta})^{T}\widehat{As.Cov}(g(\hat{\theta}))^{-1}g(\hat{\theta}),$$

where $\widehat{As.Cov}(g(\hat{\theta}))$ is the estimate of the asymptotic covariance matrix of $g(\hat{\theta})$.

If se_type = "dm" then delta method is used to estimate this matrix:

$$\widehat{As.Cov}(g(\hat{\theta}))= g'(\hat{\theta})\widehat{As.Cov}(\hat{\theta}) g'(\hat{\theta})^{T},$$

where $g'(\hat{\theta})$ is a Jacobian matrix. A numeric differentiation is used to calculate its elements:

$$g'(\hat{\theta})_{ij} =\frac{\partial g_{i}(\theta)}{\partial \theta_{j}}|_{\theta = \hat{\theta}}.$$

If se_type = "bootstrap" then bootstrap is used to estimate this matrix:

$$\widehat{As.Cov}(g(\hat{\theta}))=\frac{1}{B-1}\sum\limits_{i=1}^{B}q_{b}q_{b}^{T},$$

where:

$$q_{b} = (g(\hat{\theta}^{(b)}) -\overline{g(\hat{\theta}^{(b)})}),$$ $$\overline{g(\hat{\theta}^{(b)})} = \frac{1}{B}\sum\limits_{i=1}^{B}g(\hat{\theta}^{(b)}).$$

If method = "classic" then it is assumed that if the null hypothesis is true then the asymptotic distribution of the test statistic is chi-squared with $m$ degrees of freedom. This distribution is used for the calculation of the p-value:

$$p-value = 1 - F_{m}(T),$$

where $F_{m}$ is a cumulative distribution function of the chi-squared distribution with $m$ degrees of freedom.

If method = "bootstrap" then p-value is calculated via the bootstrap as suggested by Hansen (2022):

$$p-value = \frac{1}{B}\sum\limits_{b=1}^{B} I(T_{b} > T),$$

where:

$$T_{b} = s_{b}^{T}\widehat{As.Cov}(g(\hat{\theta}))^{-1}s_{b},$$ $$s_{b} = g(\hat{\theta}^{(b)}) - g(\hat{\theta}).$$

Score bootstrap Wald test

If method = "score" and test = "Wald" then score bootstrap Wald test of Kline and Santos (2012) is used.

Consider $B$ independent samples of $n$ independent identically distributed random weights with zero mean and unit variance. Let $w_{ib}$ denote the $i$-th weight of the $b$-th sample. Argument generator is used to supply a function which generates these weights $w_{ib}$ and iter argument represents $B$. Also $n$ is the number of observations in the model object$other$n_obs.

Let $J$ denote a matrix of sample scores object$J. Further, denote by $J_{b}$ a matrix such that its $b$-th row is a product of the $w_{ib}$ and the $b$-th row of $J$. Also, denote by $H$ a matrix of mean values of the derivatives of sample scores respect to the estimated parameters object$H.

In addition consider the following notations:

$$A = g'(\theta) H^{-1}, \qquad S_{b} = A J^{(c)}_{b},$$

where $J^{(c)}_{b}$ is a vector of the column sums of $J_{b}$.

The test statistic is as follows:

$$T = g(\hat{\theta})^{T}(A\widehat{Cov}(J)A^{T})^{-1}g(\hat{\theta}) / n,$$

where $\widehat{Cov}(J)$ is a sample covariance matrix of the sample scores of the model cov(object$J).

The test statistic on the $b$-th bootstrap sample is similar:

$$T_{b} = S^{T}(A\widehat{Cov}(J_{b})A^{T})^{-1}S / n.$$

The p-value is estimated as follows:

$$p-value = \frac{1}{B}\sum\limits_{b=1}^{B} I(T_{b} \geq T).$$

Confidence intervals

If test = "t" then the function also returns the realizations of the lower and upper bounds of the $100 \times$cl percent symmetric asymptotic confidence interval of $g_{j}(\theta)$.

If ci = "classic" then classic confidence interval is used which assumes asymptotic normality of $g_{j}(\hat{\theta})$:

$$(g_{j}(\hat{\theta}) + z_{(1 - cl) / 2}\hat{\sigma}_{j},g_{j}(\hat{\theta}) + z_{1 - (1 - cl) / 2}\hat{\sigma}_{j}),$$

where $z_{q}$ is a $q$-th quantile of the standard normal distribution and $cl$ is a confidence level cl. The method used to estimate $\hat{\sigma}_{j}$ depends on the se_type argument as described above.

If ci = "percentile" then percentile bootstrap confidence interval is used. Therefore the sample quantiles of $g_{j}(\hat{\theta}^{(b)})$

are used as the realizations of the lower and upper bounds of the confidence interval.

If ci = "bc" then bias corrected percentile bootstrap confidence interval of Efron (1982) is used as described in Hansen (2022). The default percentile bootstrap confidence interval uses sample quantiles of levels $(1 - cl)/2$ and $1 - (1 - cl)/2$. Bias corrected version uses the sample quantiles of the following levels:

$$(1 - cl)/2 + \Phi(\Phi^{-1}((1 - cl)/2) + s),$$ $$1 - (1 - cl)/2 + \Phi(\Phi^{-1}(1 - (1 - cl)/2) + s),$$

where:

$$s = 2\Phi^{-1}(\frac{1}{B}\sum\limits_{b = 1}^{B}I(g_{j}(\hat{\theta}^{(b)})\leq g_{j}(\hat{\theta}))).$$

Trimming

If se_type = "bootstrap" and trim > 0 then trimming is used as described in Hansen (2022) to estimate $\hat{\sigma}_{j}$ and $\widehat{As.Cov}(g(\hat{\theta}))$. The algorithm is as follows. First, nullify 100trim percent of $g(\hat{\theta}^{(b)})$ with the greatest values of the L2-norm of $q_{b}$ (defined above). Then use this 'trimmed' sample to estimate the standard error and the asymptotic covariance matrix.

Examples

# -------------------------------
# CPS data example
# -------------------------------

# Set seed for reproducibility
set.seed(123)

# Upload the data
data(cps)

# Estimate the employment model
model <- msel(work ~ age + I(age ^ 2) + bachelor + master, data = cps)
summary(model)    

# Use Wald test to test the hypothesis that age has no 
# effect on the conditional probability of employment:
# H0: coef age     = 0
#     coef age ^ 2 = 0
age_fn <- function(object)
{
  lwage_coef <- coef(object, type = "coef")[[1]]
  val        <- c(lwage_coef["age"], lwage_coef["I(age^2)"])
  return(val)
}
age_test <- test_msel(object = model, fn = age_fn, test = "wald")
summary(age_test)

# Use t-test to test for each individual the hypothesis:
# P(work = 1 | x) = 0.8
prob_fn <- function(object)
{
  prob <- predict(object, group = 1, type = "prob")
  val  <- prob - 0.8
  return(val)
}
prob_test <- test_msel(object = model, fn = prob_fn, test = "t")
summary(prob_test)

# -------------------------------
# Simulated data example
# Model with continuous outcome
# and ordinal selection
# -------------------------------

# ---
# Step 1
# Simulation of the data
# ---

# Set seed for reproducibility
set.seed(123)

# Load required package
library("mnorm")

# The number of observations
n <- 10000

# Regressors (covariates)
s1 <- runif(n = n, min = -1, max = 1)
s2 <- runif(n = n, min = -1, max = 1)
s3 <- runif(n = n, min = -1, max = 1)
s4 <- runif(n = n, min = -1, max = 1)

# Random errors
sigma <- matrix(c(1,    0.4,  0.45, 0.7,
                  0.4,  1,    0.54, 0.8,
                  0.45, 0.54, 0.81, 0.81,
                  0.7,  0.8,  0.81, 1), nrow = 4)
errors <- mnorm::rmnorm(n = n, mean = c(0, 0, 0, 0), sigma = sigma)
u1   <- errors[, 1]
u2   <- errors[, 2]
eps0 <- errors[, 3]
eps1 <- errors[, 4]

# Coefficients
gamma1     <- c(-1, 2)
gamma2     <- c(1, 1)
gamma1_het <- c(0.5, -1)
beta0      <- c(1, -1, 1, -1.2)
beta1      <- c(2, -1.5, 0.5, 1.2)
# Linear index of the ordinal equations
# mean part
li1 <- gamma1[1] * s1 + gamma1[2] * s2
li2 <- gamma2[1] * s1 + gamma2[2] * s3
# variance part
li1_het <- gamma1_het[1] * s2 + gamma1_het[2] * s3

# Linear index of the continuous equation
# regime 0
li_y0 <- beta0[1] + beta0[2] * s1 + beta0[3] * s3 + beta0[4] * s4
# regime 1
li_y1 <- beta1[1] + beta1[2] * s1 + beta1[3] * s3 + beta1[4] * s4

# Latent variables
z1_star <- li1 + u1 * exp(li1_het)
z2_star <- li2 + u2
y0_star <- li_y0 + eps0
y1_star <- li_y1 + eps1

# Cuts
cuts1 <- c(-1)
cuts2 <- c(0, 1)

# Observable ordinal outcome
# first
z1                     <- rep(0, n)
z1[z1_star > cuts1[1]] <- 1
# second
z2                                               <- rep(0, n)
z2[(z2_star > cuts2[1]) & (z2_star <= cuts2[2])] <- 1
z2[z2_star > cuts2[2]]                           <- 2
z2[z1 == 0]                                      <- NA

# Observable continuous outcome
y                 <- rep(NA, n)
y[which(z2 == 0)] <- y0_star[which(z2 == 0)]
y[which(z2 != 0)] <- y1_star[which(z2 != 0)]
y[which(z1 == 0)] <- NA

# Data
data <- data.frame(s1 = s1, s2 = s2, s3 = s3, s4 = s4,
                   z1 = z1, z2 = z2, y = y)

# ---
# Step 2
# Estimation of the parameters
# ---

# Assign the groups
groups  <- matrix(c(1, 2,
                    1, 1,
                    1, 0,
                    0, -1), 
                  byrow = TRUE, ncol = 2)
groups2 <- matrix(c(1, 1, 0, -1), ncol = 1)

# Estimate the model
model <- msel(list(z1 ~ s1 + s2 | s2 + s3,
                   z2 ~ s1 + s3),
              list(y  ~ s1 + s3 + s4),
              groups = groups, groups2 = groups2, 
              data   = data)
                   
# ---
# Step 3
# Hypotheses testing
# ---

# Use t-test to test for each observation the hypothesis
# H0: P(z1 = 0, z2 = 2 | Xi) = 0 
prob02_fn <- function(object)
{
   val <- predict(object, group = c(1, 0))
   
   return(val)
}
prob02_test <- test_msel(object = model, fn = prob02_fn, test = "t")
summary(prob02_test)

# Use t-test to test the hypothesis
# H0: E(y1|z1=0, z2=2) - E(y0|z1=0, z2=2)
ATE_fn <- function(object)
{
   val1 <- predict(object, group = c(0, 2), group2 = 1)
   val0 <- predict(object, group = c(0, 2), group2 = 0)
   val <- mean(val1 - val0)
   
   return(val)
}
ATE_test <- test_msel(object = model, fn = ATE_fn)
summary(ATE_test)

# Use Wald to test the hypothesis
# H0: beta1 = beta0
coef_fn <- function(object)
{
   coef1 <- coef(object, regime = 1, type = "coef2")
   coef0 <- coef(object, regime = 0, type = "coef2")
   coef_difference <- coef1 - coef0
   
   return(coef_difference)
}
coef_test <- test_msel(object = model, fn = coef_fn, test = "wald")
summary(coef_test)

# Use t-test to test for each 'k' the hypothesis
# H0: beta1k = beta0k
coef_test2 <- test_msel(object = model, fn = coef_fn, test = "t")
summary(coef_test2)

# Use Wald test to test the hypothesis
# H0: beta11 + beta12 - 0.5 = 0
#     beta11 * beta13 - beta03 = 0
test_fn <- function(object)
{
  coef1 <- coef(object, regime = 1, type = "coef2")
  coef0 <- coef(object, regime = 0, type = "coef2")
  val   <- c(coef1[1] + coef1[2] - 0.5, 
  coef1[1] * coef1[3] - coef0[3])
     
  return(val)
}
# classic Wald test
wald1 <- test_msel(object = model,  fn     = test_fn, 
                   test   = "wald", method = "classic")
summary(wald1)
# score bootstrap Wald test
wald2 <- test_msel(object = model,  fn     = test_fn, 
                   test   = "wald", method = "score")
summary(wald2)

# Replicate the latter test with the 2-step estimator
model2 <- msel(list(z1 ~ s1 + s2 | s2 + s3,
                         z2 ~ s1 + s3),
               list(y ~ s1 + s3 + s4),
               groups    = groups, groups2   = groups2, 
               data      = data,   estimator = "2step")
# classic Wald test
wald1_2step <- test_msel(object = model2, fn     = test_fn, 
                         test   = "wald", method = "classic")
summary(wald1_2step)
# score bootstrap Wald test
wald2_2step <- test_msel(object = model2, fn     = test_fn, 
                         test   = "wald", method = "score")
summary(wald2_2step)

References

B. Efron (1982). The Jackknife, the Bootstrap, and Other Resampling Plans. Society for Industrial and Applied Mathematics.

B. Hansen (2022). Econometrics. Princeton University Press.

P. Kline, A. Santos (2012). A Score Based Approach to Wild Bootstrap Inference. Journal of Econometric Methods, vol. 67, no. 1, pages 23-41.