Type: | Package |
Title: | Inference for Linear Models with Nuisance Parameters |
Version: | 1.1.3 |
Date: | 2022-08-11 |
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. |
URL: | https://github.com/mlysy/LMN |
BugReports: | https://github.com/mlysy/LMN/issues |
License: | GPL-3 |
Imports: | Rcpp (≥ 0.12.4.4), SuperGauss, stats |
LinkingTo: | Rcpp, RcppEigen |
Encoding: | UTF-8 |
RoxygenNote: | 7.2.1 |
Suggests: | testthat, numDeriv, mniw, knitr, rmarkdown, bookdown, kableExtra |
VignetteBuilder: | knitr |
NeedsCompilation: | yes |
Packaged: | 2022-08-22 15:54:36 UTC; mlysy |
Author: | Martin Lysy [aut, cre], Bryan Yates [aut] |
Maintainer: | Martin Lysy <mlysy@uwaterloo.ca> |
Repository: | CRAN |
Date/Publication: | 2022-08-22 16:20:02 UTC |
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 mlysy@uwaterloo.ca
Authors:
Bryan Yates
See Also
Useful links:
Report bugs at https://github.com/mlysy/LMN/issues
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)
Arguments
x |
List of |
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)
Arguments
Beta |
A |
Sigma |
A |
suff |
An object of class |
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)
Marginal log-posterior for the LMN model.
Description
Marginal log-posterior for the LMN model.
Usage
lmn_marg(suff, prior, post)
Arguments
suff |
An object of class |
prior |
A list with elements |
post |
A list with elements |
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)
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)
Arguments
suff |
An object of class |
prior |
A list with elements |
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 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
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
p |
Integer specifying row dimension of |
q |
Integer specifying the dimension of |
Lambda |
Mean parameter for
|
Omega |
Row-wise precision parameter for
|
Psi |
Scale parameter for
|
nu |
Degrees-of-freedom parameter for |
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 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)
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)
Arguments
suff |
An object of class |
noSigma |
Logical. If |
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)
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 |
X |
An
|
V , Vtype |
The between-observation variance specification. Currently the following options are supported:
For |
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)
(orp = 0
ifX = 0
), andq = 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")