MVNIGmod {MVSKmod}R Documentation

AECM Estimation for Matrix-Variate Normal-Inverse Gaussian Models

Description

This function fits MVNIG linear models for matrix-variate skew data with non-uniform data rows between subjects. Exchangeable observation row correlation and skewness structures are imposed to accommodate the varying row counts across matrices. Note that multiple restarts may be needed to account for unstable local maxima.

Usage

MVNIGmod(Y, X, theta_g = NULL, stopping = 0.001, max_iter = 50)

Arguments

Y

List of n_i \times p response matrices. Matrices must have same number of columns.

X

List of n_i \times q design matrices. Matrices must have same number of columns.

theta_g

List of parameters to pass as initial values in the AECM algorithm. If NULL, will be randomly generated. See Details for an in-depth explanation.

stopping

Stopping threshold for the L-infinity norm of differences in consecutive parameter space, evaluated at iteration t+1 as |\hat\theta^{t+1} - \hat\theta^{t}|_\infty. Default is 0.001

max_iter

Maximum number of iterations, default is 50.

Details

Fits the matrix-variate skew regression model

Y_i = X_i \Theta + E_i,

where each response Y_i is a n_i \times p matrix that indexes n_i observations and p response variables. X_i corresponds to a n_i \times q design matrix, and \Theta corresponds to a q \times p coefficient matrix. E_i corresponds to a n_i \times p error matrix, following a matrix-variate variance-gamma distribution.

The model estimates MVVG parameters \Theta, \underline{a}, r, \Psi, \tilde\gamma using the alternating expectation conditional maximization (AECM) algorithm, using the density

f(Y_i|M_i,\underline{a},r, \Psi, \tilde\gamma, n_i, p) = \dfrac{2 \exp[matlib::tr(\Sigma_i^{-1}(Y_i-M_i)\Psi^{-1}A_i^T) + \tilde\gamma]}{(2\pi)^{\frac{n_ip}{2} + 1} |\Sigma_i|^{\frac{p}{2}} |\Psi|^{\frac{n_i}{2}}} \bigg( \dfrac{\delta(Y_i; M_i, \Sigma_i, \Psi) + 1}{\rho (A_i, \Sigma_i,\Psi) + \tilde\gamma^2} \bigg)^{-\frac{(1+n_ip)}{4}} \\ \times K_{-\frac{(1+n_ip)}{2}} \big( \sqrt{[\rho(A_i, \Sigma_i, \Psi) + \tilde\gamma^2][\delta(Y_i;M_i,\Sigma_i,\Psi) + 1]} \big),

where A_i = \underline{1}_{n_i} \times \underline{a}^T, \Sigma_i = I_{n_i} + r(\underline{1}_{n_i}\underline{1}_{n_i}^T - I_{n_i}), \delta(X;M, \Sigma, \Psi) = matlib::tr(\Sigma^{-1}(X-M)\Psi^{-1}(X-M)^T), \rho(A, \Sigma, \Psi) = matlib::tr(\Sigma^{-1}A\Psi^{-1}A^T), and K_{\nu}(x) is the modified Bessel function of the second kind.

The structure of theta_g and parameter estimates returned by the function must be in the form of a list with the following named elements:

Theta:

q \times p coefficient matrix

a:

p \times 1 skewness vector

rho:

Compound symmetry parameter for row correlation matrix

Psi:

p \times p column covariance matrix

tgamma:

Univariate mixing parameter

Value

MVNIGmod returns a list with the following elements:

Iteration:

Number of iterations taken to convergence. Inf if convergence not reached.

⁠Starting Value⁠:

List of initial parameter values.

⁠Final Value⁠:

List of final parameter estimates.

⁠Stopping Criteria⁠:

Vector of |\hat\theta^{t+1} - \hat\theta^{t}|_\infty at each iteration.

AIC:

Model AIC

BIC:

Model BIC

Author(s)

Samuel Soon

Dipankar Bandyopadhyay

Qingyang Liu

Examples


MVNIGmod(Y,X,theta_mvnig)

set.seed(1234)
# num response variables
p <- ncol(gaad_res[[1]])
# num covariates
q <- ncol(gaad_cov[[1]])
# generate initial value to input, then run AECM with MVVG distribution
initial_mvnig_theta <- list(Theta = matrix(stats::rnorm(p*q), nrow = q, ncol = p),
                      A = rep(1,p),
                     rho = 0.3,
                     Psi = diag(p),
                     tgamma = 3)
MVNIGmod(gaad_res[1:30], gaad_cov[1:30], initial_mvnig_theta)


[Package MVSKmod version 0.1.0 Index]