lmn_suff {LMN}

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)

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.

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}.

Examples

# 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")