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 |
X |
List of |
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 |
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 matrixa
:p \times 1
skewness vectorrho
:Compound symmetry parameter for row correlation matrix
Psi
:p \times p
column covariance matrixtgamma
: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)