lmn_suff: Calculate the sufficient statistics of an LMN model.
In mlysy/LMN: Inference for Linear Models with Nuisance Parameters

View source: R/lmn_suff.R

lmn_suff

R Documentation

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

Y ~ Matrix-Normal(X B, V, Σ),

where Y_(n x q) is the response matrix, X_(n x p) is the covariate matrix, B_(p x q) is the coefficient matrix, V_(n x n) and Σ_(q x q) are the between-row and between-column variance matrices, and the Matrix-Normal distribution is defined by the multivariate normal distribution vec(Y) ~ N( vec(X B), Σ %x% V ), where vec(Y) is a vector of length nq stacking the columns of of Y, and Σ %x% V is the Kronecker product.

The function lmn_suff() returns everything needed to efficiently calculate the likelihood function

L(B, Σ | Y, X, V) = p(Y | X, V, B, Σ).

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(y | Y, X_star, V_star, B, Σ).

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 x q matrix B_hat = (X'V^{-1}X)^{-1}X'V^{-1}Y.
T: The p x p matrix T = X'V^{-1}X.
S: The q x q matrix S = (Y-X B_hat)'V^{-1}(Y-X B_hat).
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 x, w, and v defined in Details:

Ap: The npred x q matrix A_p = w'V^{-1}Y.
Xp: The npred x p matrix X_p = x - w'V^{-1}X.
Vp: The scalar V_p = v - w'V^{-1}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")

mlysy/LMN documentation built on March 25, 2022, 11:12 a.m.