Package 'LMN' reference manual

Title:	Inference for Linear Models with Nuisance Parameters
Description:	Efficient Frequentist profiling and Bayesian marginalization of parameters for which the conditional likelihood is that of a multivariate linear regression model. Arbitrary inter-observation error correlations are supported, with optimized calculations provided for independent-heteroskedastic and stationary dependence structures.
Authors:	Martin Lysy [aut, cre], Bryan Yates [aut]
Maintainer:	Martin Lysy <[email protected]>
License:	GPL-3
Version:	1.1.2
Built:	2024-11-20 04:20:38 UTC
Source:	https://github.com/mlysy/lmn

Inference for Linear Models with Nuisance Parameters.

Description

Efficient profile likelihood and marginal posteriors when nuisance parameters are those of linear regression models.

Details

Consider a model $p(\boldsymbol{Y} \mid \boldsymbol{B}, \boldsymbol{\Sigma}, \boldsymbol{\theta})$ of the form

$\boldsymbol{Y} \sim \textrm{Matrix-Normal}(\boldsymbol{X}(\boldsymbol{\theta})\boldsymbol{B}, \boldsymbol{V}(\boldsymbol{\theta}), \boldsymbol{\Sigma}),$

where $\boldsymbol{Y}_{n \times q}$ is the response matrix, $\boldsymbol{X}(\theta)_{n \times p}$ is a covariate matrix which depends on $\boldsymbol{\theta}$ , $\boldsymbol{B}_{p \times q}$ is the coefficient matrix, $\boldsymbol{V}(\boldsymbol{\theta})_{n \times n}$ and $\boldsymbol{\Sigma}_{q \times q}$ are the between-row and between-column variance matrices, and (suppressing the dependence on $\boldsymbol{\theta}$ ) 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 model above is referred to as a Linear Model with Nuisance parameters (LMN) $(\boldsymbol{B}, \boldsymbol{\Sigma})$ , with parameters of interest $\boldsymbol{\theta}$ . That is, the LMN package provides tools to efficiently conduct inference on $\boldsymbol{\theta}$ first, and subsequently on $(\boldsymbol{B}, \boldsymbol{\Sigma})$ , by Frequentist profile likelihood or Bayesian marginal inference with a Matrix-Normal Inverse-Wishart (MNIW) conjugate prior on $(\boldsymbol{B}, \boldsymbol{\Sigma})$ .

Author(s)

Maintainer: Martin Lysy [email protected]

Authors:

Bryan Yates

Convert list of MNIW parameter lists to vectorized format.

Description

Converts a list of return values of multiple calls to lmn_prior() or lmn_post() to a single list of MNIW parameters, which can then serve as vectorized arguments to the functions in mniw.

Usage

list2mniw(x)
list2mniw(x)

Arguments

`x`	List of `n` MNIW parameter lists.

Value

A list with the following elements:

Lambda: The mean matrices as an array of size ⁠p x p x n⁠.
Omega: The between-row precision matrices, as an array of size ⁠p x p x ⁠.
Psi: The between-column scale matrices, as an array of size ⁠q x q x n⁠.
nu: The degrees-of-freedom parameters, as a vector of length n.

Loglikelihood function for LMN models.

Description

Loglikelihood function for LMN models.

Usage

lmn_loglik(Beta, Sigma, suff)
lmn_loglik(Beta, Sigma, suff)

Arguments

`Beta`	A `⁠p x q⁠` matrix of regression coefficients (see `lmn_suff()`).
`Sigma`	A `⁠q x q⁠` matrix of error variances (see `lmn_suff()`).
`suff`	An object of class `lmn_suff` (see `lmn_suff()`).

Value

Scalar; the value of the loglikelihood.

Examples

# generate data
n <- 50
q <- 3
Y <- matrix(rnorm(n*q),n,q) # response matrix
X <- 1 # intercept covariate
V <- 0.5 # scalar variance specification
suff <- lmn_suff(Y, X = X, V = V) # sufficient statistics

# calculate loglikelihood
Beta <- matrix(rnorm(q),1,q)
Sigma <- diag(rexp(q))
lmn_loglik(Beta = Beta, Sigma = Sigma, suff = suff)
# generate data
n <- 50
q <- 3
Y <- matrix(rnorm(n*q),n,q) # response matrix
X <- 1 # intercept covariate
V <- 0.5 # scalar variance specification
suff <- lmn_suff(Y, X = X, V = V) # sufficient statistics

# calculate loglikelihood
Beta <- matrix(rnorm(q),1,q)
Sigma <- diag(rexp(q))
lmn_loglik(Beta = Beta, Sigma = Sigma, suff = suff)

Marginal log-posterior for the LMN model.

Description

Marginal log-posterior for the LMN model.

Usage

lmn_marg(suff, prior, post)
lmn_marg(suff, prior, post)

Arguments

`suff`	An object of class `lmn_suff` (see `lmn_suff()`).
`prior`	A list with elements `Lambda`, `Omega`, `Psi`, `nu` corresponding to the parameters of the prior MNIW distribution. See `lmn_prior()`.
`post`	A list with elements `Lambda`, `Omega`, `Psi`, `nu` corresponding to the parameters of the posterior MNIW distribution. See `lmn_post()`.

Value

The scalar value of the marginal log-posterior.

Examples

# generate data
n <- 50
q <- 2
p <- 3
Y <- matrix(rnorm(n*q),n,q) # response matrix
X <- matrix(rnorm(n*p),n,p) # covariate matrix
V <- .5 * exp(-(1:n)/n) # Toeplitz variance specification

suff <- lmn_suff(Y = Y, X = X, V = V, Vtype = "acf") # sufficient statistics
# default noninformative prior pi(Beta, Sigma) ~ |Sigma|^(-(q+1)/2)
prior <- lmn_prior(p = suff$p, q = suff$q)
post <- lmn_post(suff, prior = prior) # posterior MNIW parameters
lmn_marg(suff, prior = prior, post = post)
# generate data
n <- 50
q <- 2
p <- 3
Y <- matrix(rnorm(n*q),n,q) # response matrix
X <- matrix(rnorm(n*p),n,p) # covariate matrix
V <- .5 * exp(-(1:n)/n) # Toeplitz variance specification

suff <- lmn_suff(Y = Y, X = X, V = V, Vtype = "acf") # sufficient statistics
# default noninformative prior pi(Beta, Sigma) ~ |Sigma|^(-(q+1)/2)
prior <- lmn_prior(p = suff$p, q = suff$q)
post <- lmn_post(suff, prior = prior) # posterior MNIW parameters
lmn_marg(suff, prior = prior, post = post)

Parameters of the posterior conditional distribution of an LMN model.

Description

Calculates the parameters of the LMN model's Matrix-Normal Inverse-Wishart (MNIW) conjugate posterior distribution (see Details).

Usage

lmn_post(suff, prior)
lmn_post(suff, prior)

Arguments

`suff`	An object of class `lmn_suff` (see `lmn_suff()`).
`prior`	A list with elements `Lambda`, `Omega`, `Psi`, `nu` as returned by `lmn_prior()`.

Details

The Matrix-Normal Inverse-Wishart (MNIW) distribution $(\boldsymbol{B}, \boldsymbol{\Sigma}) \sim \textrm{MNIW}(\boldsymbol{\Lambda}, \boldsymbol{\Omega}, \boldsymbol{\Psi}, \nu)$ on random matrices $\boldsymbol{X}_{p \times q}$ and symmetric positive-definite $\boldsymbol{\Sigma}_{q \times q}$ is defined as

$\begin{array}{rcl} \boldsymbol{\Sigma} & \sim & \textrm{Inverse-Wishart}(\boldsymbol{\Psi}, \nu) \\ \boldsymbol{B} \mid \boldsymbol{\Sigma} & \sim & \textrm{Matrix-Normal}(\boldsymbol{\Lambda}, \boldsymbol{\Omega}^{-1}, \boldsymbol{\Sigma}), \end{array}$

$\vspace{-1em}$

where the Matrix-Normal distribution is defined in lmn_suff().

The posterior MNIW distribution is required to be a proper distribution, but the prior is not. For example, prior = NULL corresponds to the noninformative prior

$\pi(B, \boldsymbol{\Sigma}) \sim |\boldsymbol{Sigma}|^{-(q+1)/2}.$

Value

A list with elements named as in prior specifying the parameters of the posterior MNIW distribution. Elements Omega = NA and nu = NA specify that parameters Beta = 0 and Sigma = diag(q), respectively, are known and not to be estimated.

Examples

# generate data
n <- 50
q <- 2
p <- 3
Y <- matrix(rnorm(n*q),n,q) # response matrix
X <- matrix(rnorm(n*p),n,p) # covariate matrix
V <- .5 * exp(-(1:n)/n) # Toeplitz variance specification

suff <- lmn_suff(Y = Y, X = X, V = V, Vtype = "acf") # sufficient statistics
# generate data
n <- 50
q <- 2
p <- 3
Y <- matrix(rnorm(n*q),n,q) # response matrix
X <- matrix(rnorm(n*p),n,p) # covariate matrix
V <- .5 * exp(-(1:n)/n) # Toeplitz variance specification

suff <- lmn_suff(Y = Y, X = X, V = V, Vtype = "acf") # sufficient statistics

Conjugate prior specification for LMN models.

Description

The conjugate prior for LMN models is the Matrix-Normal Inverse-Wishart (MNIW) distribution. This convenience function converts a partial MNIW prior specification into a full one.

Usage

lmn_prior(p, q, Lambda, Omega, Psi, nu)
lmn_prior(p, q, Lambda, Omega, Psi, nu)

Arguments

`p`	Integer specifying row dimension of `Beta`. `p = 0` corresponds to no `Beta` in the model, i.e., `X = 0` in `lmn_suff()`.
`q`	Integer specifying the dimension of `Sigma`.
`Lambda`	Mean parameter for `Beta`. Either: A `⁠p x q⁠` matrix. A scalar, in which case `Lambda = matrix(Lambda, p, q)`. Missing, in which case `Lambda = matrix(0, p, q)`.
`Omega`	Row-wise precision parameter for `Beta`. Either: A `⁠p x p⁠` matrix. A scalar, in which case `Omega = diag(rep(Omega,p))`. Missing, in which case `Omega = matrix(0, p, p)`. `NA`, which signifies that `Beta` is known, in which case the prior is purely Inverse-Wishart on `Sigma` (see Details).
`Psi`	Scale parameter for `Sigma`. Either: A `⁠q x q⁠` matrix. A scalar, in which case `Psi = diag(rep(Psi,q))`. Missing, in which case `Psi = matrix(0, q, q)`.
`nu`	Degrees-of-freedom parameter for `Sigma`. Either a scalar, missing (defaults to `nu = 0`), or `NA`, which signifies that `Sigma = diag(q)` is known, in which case the prior is purely Matrix-Normal on `Beta` (see Details).

Details

$\vspace{-1em}$

where the Matrix-Normal distribution is defined in lmn_suff().

Value

A list with elements Lambda, Omega, Psi, nu with the proper dimensions specified above, except possibly Omega = NA or nu = NA (see Details).

Examples

# problem dimensions
p <- 2
q <- 4

# default noninformative prior pi(Beta, Sigma) ~ |Sigma|^(-(q+1)/2)
lmn_prior(p, q)

# pi(Sigma) ~ |Sigma|^(-(q+1)/2)
# Beta | Sigma ~ Matrix-Normal(0, I, Sigma)
lmn_prior(p, q, Lambda = 0, Omega = 1)

# Sigma = diag(q)
# Beta ~ Matrix-Normal(0, I, Sigma = diag(q))
lmn_prior(p, q, Lambda = 0, Omega = 1, nu = NA)
# problem dimensions
p <- 2
q <- 4

# default noninformative prior pi(Beta, Sigma) ~ |Sigma|^(-(q+1)/2)
lmn_prior(p, q)

# pi(Sigma) ~ |Sigma|^(-(q+1)/2)
# Beta | Sigma ~ Matrix-Normal(0, I, Sigma)
lmn_prior(p, q, Lambda = 0, Omega = 1)

# Sigma = diag(q)
# Beta ~ Matrix-Normal(0, I, Sigma = diag(q))
lmn_prior(p, q, Lambda = 0, Omega = 1, nu = NA)

Profile loglikelihood for the LMN model.

Description

Calculate the loglikelihood of the LMN model defined in lmn_suff() at the MLE Beta = Bhat and Sigma = Sigma.hat.

Usage

lmn_prof(suff, noSigma = FALSE)
lmn_prof(suff, noSigma = FALSE)

Arguments

`suff`	An object of class `lmn_suff` (see `lmn_suff()`).
`noSigma`	Logical. If `TRUE` assumes that `Sigma = diag(ncol(Y))` is known and therefore not estimated.

Value

Scalar; the calculated value of the profile loglikelihood.

Examples

# generate data
n <- 50
q <- 2
Y <- matrix(rnorm(n*q),n,q) # response matrix
X <- matrix(1,n,1) # covariate matrix
V <- exp(-(1:n)/n) # diagonal variance specification
suff <- lmn_suff(Y, X = X, V = V, Vtype = "diag") # sufficient statistics

# profile loglikelihood
lmn_prof(suff)

# check that it's the same as loglikelihood at MLE
lmn_loglik(Beta = suff$Bhat, Sigma = suff$S/suff$n, suff = suff)
# generate data
n <- 50
q <- 2
Y <- matrix(rnorm(n*q),n,q) # response matrix
X <- matrix(1,n,1) # covariate matrix
V <- exp(-(1:n)/n) # diagonal variance specification
suff <- lmn_suff(Y, X = X, V = V, Vtype = "diag") # sufficient statistics

# profile loglikelihood
lmn_prof(suff)

# check that it's the same as loglikelihood at MLE
lmn_loglik(Beta = suff$Bhat, Sigma = suff$S/suff$n, suff = 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)
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")

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

Package 'LMN'

Help Index

Inference for Linear Models with Nuisance Parameters.

Description

Details

Author(s)

See Also

Convert list of MNIW parameter lists to vectorized format.

Description

Usage

Arguments

Value

Loglikelihood function for LMN models.

Description

Usage

Arguments

Value

Examples

Marginal log-posterior for the LMN model.

Description

Usage

Arguments

Value

Examples

Parameters of the posterior conditional distribution of an LMN model.

Description

Usage

Arguments

Details

Value

Examples

Conjugate prior specification for LMN models.

Description

Usage

Arguments

Details

Value

Examples

Profile loglikelihood for the LMN model.

Description

Usage

Arguments

Value

Examples

Calculate the sufficient statistics of an LMN model.

Description

Usage

Arguments

Details

Value

Examples