delta_method {switchSelection}

Delta method for mvoprobit and mnprobit functions

Description

This function uses delta method to estimate standard errors of functions of the estimator of the parameters of mnprobit and mvoprobit functions if maximum-likelihood estimator has been used.

Usage

delta_method(
  object,
  fn,
  fn_args = list(),
  eps = max(1e-04, sqrt(.Machine$double.eps) * 10),
  cl = 0.95
)

Arguments

object

an object of class 'mvoprobit' or 'mnprobit'.

fn

function which returns a numeric vector and should depend on the elements of object. This elements should be accessed via coef.mvoprobit and coef.mnprobit functions. Also it is possible to use predict.mvoprobit and predict.mnprobit functions. The first argument of fn should be object. Therefore coef and predict functions in fn should also depend on object.

fn_args

list of additional arguments of fn.

eps

positive numeric value representing the increment used for numeric differentiation of fn.

cl

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

Details

Numeric differentiation is used to estimate derivatives of fn respect to various parameters of mvoprobit and mnprobit functions.

This function may be used only if object$estimator = "ml".

Value

This function returns an object of class delta_method that is a matrix which columns are as follows:

• val - output of the fn function.

• se - numeric vector such that se[i] represents standard error associated with val[i].

• p_value - numeric vector such that p_value[i] represents p-value of the two-sided significance test associated with val[i].

• lwr - realization of the lower (left) bound of the confidence interval.

• upr - realization of the upper (right) bound of the confidence interval.

An object of class delta_method has implementation of summary method summary.delta_method.

Examples

# 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 ordered equation
# mean
li1 <- gamma1[1] * s1 + gamma1[2] * s2
li2 <- gamma2[1] * s1 + gamma2[2] * s3
# variance
li1_het <- gamma1_het[1] * s2 + gamma1_het[2] * s3

# Linear index of 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 ordered 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] <- -1

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

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

# Assign 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)

# Estimation
model <- mvoprobit(list(z1 ~ s1 + s2 | s2 + s3,
                        z2 ~ s1 + s3),
                   list(y ~ s1 + s3 + s4),
                   groups = groups, groups2 = groups2, data = data)

# Use delta method to estimate standard error for each P(z1 = 0, z2 = 2)
prob02_fn <- function(object)
{
   val <- predict(object, group = c(1, 0))
   
   return(val)
}
prob02 <- delta_method(object = model, fn = prob02_fn)
head(prob02)

# Use delta method to estimate standard error for each 
# 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 <- delta_method(object = model, fn = ATE_fn)
summary(ATE)

# Use delta method to estimate standard for the difference
# between beta0 and beta1 coefficients
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_val <- delta_method(object = model, fn = coef_fn)
summary(coef_val)