lmn_suff: Calculate the sufficient statistics of an LMN model.

Description

Calculate the sufficient statistics of an LMN model.

Usage

lmn_suff(Y, X, V, Vtype, npred = 0)

Value

An S3 object of type lmn_suff, consisting of a list with elements:

Bhat: The $p \times q$ matrix $\hat{\boldsymbol{B}} = (\boldsymbol{X}'\boldsymbol{V}^{-1}\boldsymbol{X})^{-1}\boldsymbol{X}'\boldsymbol{V}^{-1}\boldsymbol{Y}$.
T: The $p \times p$ matrix $\boldsymbol{T} = \boldsymbol{X}'\boldsymbol{V}^{-1}\boldsymbol{X}$.
S: The $q \times q$ matrix $\boldsymbol{S} = (\boldsymbol{Y} - \boldsymbol{X} \hat{\boldsymbol{B}})'\boldsymbol{V}^{-1}(\boldsymbol{Y} - \boldsymbol{X} \hat{\boldsymbol{B}})$.
ldV: The scalar log-determinant of V.
n, p, q: The problem dimensions, namely n = nrow(Y), p = nrow(Beta) (or p = 0 if X = 0), and q = ncol(Y).

In addition, when npred > 0 and with $\boldsymbol{x}$, $\boldsymbol{w}$, and $v$ defined in Details:

Ap: The npred x q matrix $\boldsymbol{A}_p = \boldsymbol{w}'\boldsymbol{V}^{-1}\boldsymbol{Y}$.
Xp: The npred x p matrix $\boldsymbol{X}_p = \boldsymbol{x} - \boldsymbol{w}\boldsymbol{V}^{-1}\boldsymbol{X}$.
Vp: The scalar $V_p = v - \boldsymbol{w}\boldsymbol{V}^{-1}\boldsymbol{w}$.

Arguments

Y

An n x q matrix of responses.

X

An N x p matrix of covariates, where N = n + npred (see Details). May also be passed as:

A scalar, in which case the one-column covariate matrix is X = X * matrix(1, N, 1). -X = 0, in which case the mean of Y is known to be zero, i.e., no regression coefficients are estimated.

V, Vtype

The between-observation variance specification. Currently the following options are supported:

Vtype = "full": V is an N x N symmetric positive-definite matrix.
Vtype = "diag": V is a vector of length N such that V = diag(V).
Vtype = "scalar": V is a scalar such that V = V * diag(N).
Vtype = "acf": V is either a vector of length N or an object of class SuperGauss::Toeplitz, such that V = toeplitz(V).

For V specified as a matrix or scalar, Vtype is deduced automatically and need not be specified.

npred

A nonnegative integer. If positive, calculates sufficient statistics to make predictions for new responses. See Details.

Details

The multi-response normal linear regression model is defined as $$ \boldsymbol{Y} \sim \textrm{Matrix-Normal}(\boldsymbol{X}\boldsymbol{B}, \boldsymbol{V}, \boldsymbol{\Sigma}), $$ where $\boldsymbol{Y}_{n \times q}$ is the response matrix, $\boldsymbol{X}_{n \times p}$ is the covariate matrix, $\boldsymbol{B}_{p \times q}$ is the coefficient matrix, $\boldsymbol{V}_{n \times n}$ and $\boldsymbol{\Sigma}_{q \times q}$ are the between-row and between-column variance matrices, and the Matrix-Normal distribution is defined by the multivariate normal distribution $ \textrm{vec}(\boldsymbol{Y}) \sim \mathcal{N}(\textrm{vec}(\boldsymbol{X}\boldsymbol{B}), \boldsymbol{\Sigma} \otimes \boldsymbol{V}), $ where $\textrm{vec}(\boldsymbol{Y})$ is a vector of length $nq$ stacking the columns of of $\boldsymbol{Y}$, and $\boldsymbol{\Sigma} \otimes \boldsymbol{V}$ is the Kronecker product.

The function lmn_suff() returns everything needed to efficiently calculate the likelihood function $$\mathcal{L}(\boldsymbol{B}, \boldsymbol{\Sigma} \mid \boldsymbol{Y}, \boldsymbol{X}, \boldsymbol{V}) = p(\boldsymbol{Y} \mid \boldsymbol{X}, \boldsymbol{V}, \boldsymbol{B}, \boldsymbol{\Sigma}). $$

When npred > 0, define the variables Y_star = rbind(Y, y), X_star = rbind(X, x), and V_star = rbind(cbind(V, w), cbind(t(w), v)). Then lmn_suff() calculates summary statistics required to estimate the conditional distribution $$ p(\boldsymbol{y} \mid \boldsymbol{Y}, \boldsymbol{X}_\star, \boldsymbol{V}_\star, \boldsymbol{B}, \boldsymbol{\Sigma}). $$ The inputs to lmn_suff() in this case are Y = Y, X = X_star, and V = V_star.

Examples

Run this code

# Data
n <- 50
q <- 3
Y <- matrix(rnorm(n*q),n,q)

# No intercept, diagonal V input
X <- 0
V <- exp(-(1:n)/n)
lmn_suff(Y, X = X, V = V, Vtype = "diag")

# X = (scaled) Intercept, scalar V input (no need to specify Vtype)
X <- 2
V <- .5
lmn_suff(Y, X = X, V = V)

# X = dense matrix, Toeplitz variance matrix
p <- 2
X <- matrix(rnorm(n*p), n, p)
Tz <- SuperGauss::Toeplitz$new(acf = 0.5*exp(-seq(1:n)/n))
lmn_suff(Y, X = X, V = Tz, Vtype = "acf")

Run the code above in your browser using DataLab