Movatterモバイル変換

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 modelp(\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 lengthnq 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, theLMN 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 Lysymlysy@uwaterloo.ca

Authors:

• Bryan Yates

See Also

Useful links:

• https://github.com/mlysy/LMN

• Report bugs athttps://github.com/mlysy/LMN/issues

Convert list of MNIW parameter lists to vectorized format.

Description

Converts a list of return values of multiple calls tolmn_prior() orlmn_post() to a single list of MNIW parameters, which can then serve as vectorized arguments to the functions inmniw.

Usage

list2mniw(x)

Arguments

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

Loglikelihood function for LMN models.

Description

Loglikelihood function for LMN models.

Usage

lmn_loglik(Beta, Sigma, suff)

Arguments

Value

Scalar; the value of the loglikelihood.

Examples

# generate datan <- 50q <- 3Y <- matrix(rnorm(n*q),n,q) # response matrixX <- 1 # intercept covariateV <- 0.5 # scalar variance specificationsuff <- lmn_suff(Y, X = X, V = V) # sufficient statistics# calculate loglikelihoodBeta <- 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)

Arguments

Value

The scalar value of the marginal log-posterior.

Examples

# generate datan <- 50q <- 2p <- 3Y <- matrix(rnorm(n*q),n,q) # response matrixX <- matrix(rnorm(n*p),n,p) # covariate matrixV <- .5 * exp(-(1:n)/n) # Toeplitz variance specificationsuff <- 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 parameterslmn_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 (seeDetails).

Usage

lmn_post(suff, prior)

Arguments

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}

where the Matrix-Normal distribution is defined inlmn_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 inprior specifying the parameters of the posterior MNIW distribution. ElementsOmega = NA andnu = NA specify that parametersBeta = 0 andSigma = diag(q), respectively, are known and not to be estimated.

Examples

# generate datan <- 50q <- 2p <- 3Y <- matrix(rnorm(n*q),n,q) # response matrixX <- matrix(rnorm(n*p),n,p) # covariate matrixV <- .5 * exp(-(1:n)/n) # Toeplitz variance specificationsuff <- 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)

Arguments

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}

where the Matrix-Normal distribution is defined inlmn_suff().

Value

A list with elementsLambda,Omega,Psi,nu with the proper dimensions specified above, except possiblyOmega = NA ornu = NA (seeDetails).

Examples

# problem dimensionsp <- 2q <- 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 inlmn_suff() at the MLEBeta = Bhat andSigma = Sigma.hat.

Usage

lmn_prof(suff, noSigma = FALSE)

Arguments

Value

Scalar; the calculated value of the profile loglikelihood.

Examples

# generate datan <- 50q <- 2Y <- matrix(rnorm(n*q),n,q) # response matrixX <- matrix(1,n,1) # covariate matrixV <- exp(-(1:n)/n) # diagonal variance specificationsuff <- lmn_suff(Y, X = X, V = V, Vtype = "diag") # sufficient statistics# profile loglikelihoodlmn_prof(suff)# check that it's the same as loglikelihood at MLElmn_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)

Arguments

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 lengthnq stacking the columns of of\boldsymbol{Y}, and\boldsymbol{\Sigma} \otimes \boldsymbol{V} is the Kronecker product.

The functionlmn_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}).

Whennpred > 0, define the variablesY_star = rbind(Y, y),X_star = rbind(X, x), andV_star = rbind(cbind(V, w), cbind(t(w), v)). Thenlmn_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 tolmn_suff() in this case areY = Y,X = X_star, andV = V_star.

Value

An S3 object of typelmn_suff, consisting of a list with elements:

Bhat

Thep \times q matrix\hat{\boldsymbol{B}} = (\boldsymbol{X}'\boldsymbol{V}^{-1}\boldsymbol{X})^{-1}\boldsymbol{X}'\boldsymbol{V}^{-1}\boldsymbol{Y}.

T

Thep \times p matrix\boldsymbol{T} = \boldsymbol{X}'\boldsymbol{V}^{-1}\boldsymbol{X}.

S

Theq \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 ofV.

n,p,q

The problem dimensions, namelyn = nrow(Y),p = nrow(Beta) (orp = 0 ifX = 0), andq = ncol(Y).

In addition, whennpred > 0 and with\boldsymbol{x},\boldsymbol{w}, andv defined inDetails:

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 scalarV_p = v - \boldsymbol{w}\boldsymbol{V}^{-1}\boldsymbol{w}.

Examples

# Datan <- 50q <- 3Y <- matrix(rnorm(n*q),n,q)# No intercept, diagonal V inputX <- 0V <- 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 <- 2V <- .5lmn_suff(Y, X = X, V = V)# X = dense matrix, Toeplitz variance matrixp <- 2X <- 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")