BHHH (outer product of gradients) information matrix for Mixed Logit
Source:R/RcppExports.R
mxl_bhhh_parallel.RdComputes the BHHH approximation to the observed information matrix for the Mixed Logit model: \(H_{BHHH} = \sum_i w_i \cdot s_i s_i^\top\), where \(s_i\) is the score of likelihood unit i (gradient of its simulated log-likelihood \(\log \bar{P}_i\); a decision maker when Ti is supplied, a choice situation otherwise). This outer product of gradients (OPG) estimator is an alternative to the analytical Hessian for standard errors. It needs no second derivatives, and so costs less. Like the Hessian, it forms each unit's contribution only over the parameters that unit's likelihood involves (beta, mu, L and the free ASCs of the alternatives the unit faced), so the work per unit does not grow with the number of alternatives; the n_params x n_params result and the threads' accumulators do (see acc_bytes).
Usage
mxl_bhhh_parallel(
theta,
X,
W,
alt_idx,
choice_idx,
M,
weights,
eta_draws,
rc_dist,
rc_correlation = TRUE,
rc_mean = FALSE,
use_asc = TRUE,
include_outside_option = FALSE,
gen_seed = -1L,
gen_scramble = 1L,
gen_S = 0L,
Ti = NULL,
draw_batch = 0L,
acc_bytes = 2147483648
)Arguments
- theta
vector collecting model parameters (beta, mu, L, delta (ASCs))
- X
design matrix for covariates with fixed coefficients; sum(M_i) x K_x
- W
design matrix for covariates with random coefficients; sum(M_i) x K_w or J x K_w
- alt_idx
sum(M) x 1 vector with indices of alternatives within each choice set; 1-based indexing
- choice_idx
N x 1 vector with indices of chosen alternatives; 1-based indexing relative to X; 0 is used if include_outside_option=True
- M
N x 1 vector with number of alternatives for each individual
- weights
N x 1 vector with weights for each observation; when Ti is supplied they must be constant within each decision maker
- eta_draws
Array of standard-normal draws, K_w x S x U, where U is the number of decision makers when Ti is supplied and the number of choice situations otherwise
- rc_dist
K_w x 1 integer vector indicating distribution of random coefficients: 0 = normal, 1 = log-normal
- rc_correlation
whether random coefficients should be correlated
- rc_mean
whether to estimate means for random coefficients.
- use_asc
whether to use alternative-specific constants.
- include_outside_option
whether to include outside option normalized to 0 (if so, the outside option is not included in the data)
- gen_seed
Integer master seed for the on-the-fly Halton generator.
< 0(default) uses the materializedeta_drawscube;>= 0generates draws on the fly from this seed.- gen_scramble
Integer scramble mode for on-the-fly generation:
0= identity permutations (plain Halton, compat),1= seeded position-wise digit permutations; other values are an error.- gen_S
Integer number of draws per individual, used only when
gen_seed >= 0.- Ti
Optional integer vector with the number of choice situations of each decision maker (panel likelihood); situations must be sorted by decision maker. NULL (default): every choice situation is its own unit (cross-sectional likelihood).
- draw_batch
Integer;
0(default) forms each decision maker's draws in batches sized to a per-thread memory budget, a positive value caps the number of draws per batch (for tests).- acc_bytes
Bytes for the threads' private upper triangles (2 GiB by default). While the triangles of the T threads in use, T n_params (n_params + 1) / 2 doubles, fit, each thread adds into its own; above, or for a NaN or negative value, the threads add into the result (the continuous parameters' rows through smaller private buffers, the ASC block atomically). One thread adds into the result whatever the value. The two ways differ by rounding only;
0forces the shared result (for tests).
Value
n_params x n_params PSD matrix representing the observed information
matrix estimated by the outer product of gradients (same sign convention
as the negated Hessian returned by mxl_hessian_parallel, so it can
be inverted directly to obtain vcov).
Note
The BHHH/OPG estimator is only asymptotically equivalent to the Hessian-based information matrix at the true MLE. In finite samples it can underestimate standard errors, particularly when the model is mis-specified or away from the optimum.
Examples
# \donttest{
library(data.table)
set.seed(42)
N <- 50; J <- 3
dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N))]
#> id alt x1 w1
#> <int> <int> <num> <num>
#> 1: 1 1 1.3709584 -0.04069848
#> 2: 1 2 -0.5646982 -1.55154482
#> 3: 1 3 0.3631284 1.16716955
#> 4: 2 1 0.6328626 -0.27364570
#> 5: 2 2 0.4042683 -0.46784532
#> ---
#> 146: 49 2 1.1133860 -0.47733551
#> 147: 49 3 -0.4809928 -0.16626149
#> 148: 50 1 -0.4331690 0.86256338
#> 149: 50 2 0.6968626 0.09734049
#> 150: 50 3 -1.0563684 -1.62561674
dt[, choice := 0L]
#> id alt x1 w1 choice
#> <int> <int> <num> <num> <int>
#> 1: 1 1 1.3709584 -0.04069848 0
#> 2: 1 2 -0.5646982 -1.55154482 0
#> 3: 1 3 0.3631284 1.16716955 0
#> 4: 2 1 0.6328626 -0.27364570 0
#> 5: 2 2 0.4042683 -0.46784532 0
#> ---
#> 146: 49 2 1.1133860 -0.47733551 0
#> 147: 49 3 -0.4809928 -0.16626149 0
#> 148: 50 1 -0.4331690 0.86256338 0
#> 149: 50 2 0.6968626 0.09734049 0
#> 150: 50 3 -1.0563684 -1.62561674 0
dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#> id alt x1 w1 choice
#> <int> <int> <num> <num> <int>
#> 1: 1 1 1.3709584 -0.04069848 0
#> 2: 1 2 -0.5646982 -1.55154482 0
#> 3: 1 3 0.3631284 1.16716955 1
#> 4: 2 1 0.6328626 -0.27364570 0
#> 5: 2 2 0.4042683 -0.46784532 0
#> ---
#> 146: 49 2 1.1133860 -0.47733551 0
#> 147: 49 3 -0.4809928 -0.16626149 0
#> 148: 50 1 -0.4331690 0.86256338 1
#> 149: 50 2 0.6968626 0.09734049 0
#> 150: 50 3 -1.0563684 -1.62561674 0
d <- prepare_mxl_data(dt, "id", "alt", "choice", "x1", "w1")
eta <- get_halton_normals(50, d$N, ncol(d$W))
theta <- rep(0, ncol(d$X) + ncol(d$W) + nrow(d$alt_mapping) - 1)
H <- choicer:::mxl_bhhh_parallel(theta, d$X, d$W, d$alt_idx, d$choice_idx,
d$M, d$weights, eta, rc_dist = rep(0L, ncol(d$W)),
rc_correlation = FALSE, rc_mean = FALSE)
dim(H)
#> [1] 4 4
# }