| Type: | Package |
| Title: | Matrix-Variate Non-Gaussian Linear Regression Models |
| Version: | 0.1.1 |
| Maintainer: | Samuel Soon <samksoon2@gmail.com> |
| Description: | An implementation of the expectation conditional maximization (ECM) algorithm for matrix-variate variance gamma (MVVG) and normal-inverse Gaussian (MVNIG) linear models. These models are designed for settings of multivariate analysis with clustered non-uniform observations and correlated responses. The package includes fitting and prediction functions for both models, and an example dataset from a periodontal on Gullah-speaking African Americans, with responses in 'gaad_res', and covariates in 'gaad_cov'. For more details on the matrix-variate distributions used, see Gallaugher & McNicholas (2019) <doi:10.1016/j.spl.2018.08.012>. |
| License: | MIT + file LICENSE |
| Encoding: | UTF-8 |
| LazyData: | true |
| RoxygenNote: | 7.3.3 |
| Imports: | Bessel, clusterGeneration, DistributionUtils, matlib, maxLik, truncnorm, pracma, matrixcalc, purrr |
| URL: | https://github.com/soonsk-vcu/MVNGmod |
| BugReports: | https://github.com/soonsk-vcu/MVNGmod/issues |
| Suggests: | knitr, rmarkdown |
| VignetteBuilder: | knitr |
| NeedsCompilation: | no |
| Packaged: | 2026-06-04 13:53:19 UTC; soonsk |
| Author: | Samuel Soon [aut, cre], Dipankar Bandyopadhyay [aut], Qingyang Liu [aut] |
| Depends: | R (≥ 3.5.0) |
| Repository: | CRAN |
| Date/Publication: | 2026-06-04 14:10:02 UTC |
ECM 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 ECM 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 (ECM) algorithm, using the density
f(Y_i|M_i,\underline{a},r, \Psi, \tilde\gamma, n_i, p) = \dfrac{2 \exp[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 pcoefficient matrixa:p \times 1skewness vectorrho:Compound symmetry parameter for row correlation matrix
Psi:p \times pcolumn covariance matrixtgamma:Univariate mixing parameter
Value
MVNIGmod returns a list with the following elements:
Iteration:Number of iterations taken to convergence.
Infif 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}|_\inftyat each iteration.AIC:Model AIC
BIC:Model BIC
Author(s)
Samuel Soon
Dipankar Bandyopadhyay
Qingyang Liu
Examples
MVNIGmod(Y,X,theta_mvnig)
ECM Estimation for Matrix-Variate Variance Gamma (MVVG) Models
Description
This function fits MVVG 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
MVVGmod(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 ECM 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, \gamma using the alternating expectation conditional maximization (ECM) algorithm, using the density
f(Y_i| X_i\Theta,\underline{a},r, \Psi, \gamma, n_i, p) = \dfrac{2\gamma^\gamma \exp[tr(\Sigma_i^{-1}(Y_i- X_i\Theta)\Psi^{-1}A_i^T)]}{(2\pi)^{n_ip/2} |\Sigma_i|^{p/2} |\Psi|^{n_i/2} \Gamma(\gamma)} \bigg( \dfrac{\delta(Y_i; X_i\Theta, \Sigma_i, \Psi)}{\rho (A_i, \Sigma_i,\Psi) + 2\gamma} \bigg)^{(\gamma - n_ip/2)/2} \\ \times K_{(\gamma - n_ip/2)} \big( \sqrt{[\rho(A_i, \Sigma_i, \Psi) + 2\gamma][\delta(Y_i; X_i\Theta,\Sigma_i,\Psi)]} \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 pcoefficient matrixa:p \times 1skewness vectorrho:Compound symmetry parameter for row correlation matrix
Psi:p \times pcolumn covariance matrixgamma:Univariate mixing parameter
Value
MVVGmod returns a list with the following elements:
Iteration:Number of iterations taken to convergence.
Infif 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}|_\inftyat each iteration.AIC:Model AIC
BIC:Model BIC
Author(s)
Samuel Soon
Dipankar Bandyopadhyay
Qingyang Liu
Examples
MVVGmod(Y,X,theta_mvvg)
Toy Covariate Matrices
Description
Part of toy dataset for examples.
Usage
X
Format
List of covariate matrices for individual subjects
Examples
X
Toy Response Matrices
Description
Part of toy dataset for examples.
Usage
Y
Format
List of response matrices for individual subjects
Examples
Y
Case Deletion for MVNG Models
Description
This function calculates influence statistics for each element of the data list using a one-step approximation to the generalized Cook's distance.
Usage
case_del(Y, X, mod)
Arguments
Y |
List of |
X |
List of |
mod |
Model object returned by a call to |
Details
Computes influence statistics for the MVVG and MVNIG models using a one-step approximation to the generalized Cook's distance, based on the score and hessian matrices of the complete-data log-likelihood.
GD_i^1(\theta) = \dot Q_{[i]}(\hat\theta|\hat\theta)^T [-\ddot Q(\hat\theta|\hat\theta)]^{-1}\dot Q_{[i]}(\hat\theta|\hat\theta),
where
\dot Q_{[i]}(\hat\theta|\hat\theta) = \frac{\partial Q_{[i]}(\theta|\hat\theta)}{\partial\theta}|_{\theta = \hat\theta}
is the complete-data score function with the ith subject deleted and evaluated in \hat\theta, and
\ddot Q(\hat\theta|\hat\theta) = \frac{\partial^2 Q(\theta|\hat\theta)}{\partial\theta \partial\theta^T}|_{\theta = \hat\theta}
is the complete-data Hessian matrix, as given in Lachos et al. (2015) doi:10.1016/j.jmva.2015.06.014
Value
case_del returns a vector containing the approximate generalized Cook's distance for each corresponding subject in the provided dataset.
Author(s)
Samuel Soon
Dipankar Bandyopadhyay
Qingyang Liu
Examples
mod <- MVNIGmod(Y,X,theta_mvnig)
case_del(Y,X,mod)
GAAD Data
Description
These data sets describe periodontal measurements performed on members of the Gullah-Speaking African American community.
Usage
gaad_cov
Format
Each is a list of matrices, with rows denoting tooth sites and columns denoting CAL/PPD response.
Details
gaad_res and gaad_cov contain the response and covariate matrices of the GAAD data.
Examples
gaad_cov
gaad_res
GAAD Data
Description
These data sets describe periodontal measurements performed on members of the Gullah-Speaking African American community.
Usage
gaad_res
Format
Each is a list of matrices, with rows denoting tooth sites and columns denoting CAL/PPD response.
Details
gaad_res and gaad_cov contain the response and covariate matrices of the GAAD data.
Examples
gaad_cov
gaad_res
MVVG Parameter Format
Description
This is an example of the format of input parameter list theta.
Usage
gaad_theta_mvvg
Format
List of model parameters.
Examples
gaad_theta_mvvg
MVNG Model Prediction
Description
Predicts response values given a list of covariate matrices and a model output from either MVVGmod or MVNIGmod.
Usage
predict(mod, X)
Arguments
mod |
object outputted by either |
X |
Inputted covariate matrix |
Value
Returns a list of predicted response matrices
Author(s)
Samuel Soon
Dipankar Bandyopadhyay
Qingyang Liu
Examples
set.seed(1234)
mvnig_mod <- MVNIGmod(Y,X,theta_mvnig)
predict(mvnig_mod,X)
Toy Response Initial Parameter (MVNIG)
Description
Part of toy dataset for examples.
Usage
theta_mvnig
Format
List of parameters for input to MVNIGmod function
Examples
theta_mvvg
Toy Response Initial Parameter (MVVG)
Description
Part of toy dataset for examples.
Usage
theta_mvvg
Format
List of parameters for input to MVVGmod function
Examples
theta_mvvg