The package provides approximate maximum likelihood estimation of regression coefficients and an unstructured covariance matrix parameterizing dependence between responses of mixed types. A function for approximate likelihood-ratio testing is also included.
The package can be installed from GitHub, using devtools.
## devtools::install_github("koekvall/mmrr")
library(mmrr)Let \(y_1, \dots, y_n\) be independent \(r-\)dimensional response vectors whose elements are possibly of mixed types, some continuous and some discrete for example. Let \(X_1, \dots, X_n\) be corresponding \(r\times p\) design matrices. The package fits the regression model \[ E(y_i \mid w_i) = g^{-1}(w_i), ~~ w_i \sim \mathrm{N}_r(X_i \beta, \Sigma), \] where \(g\) is an \(\mathbb{R}^r\)-valued link function, \(w_i\) is a latent (unobservable) variable, \(\beta \in \mathbb{R}^p\), and \(\Sigma \in \mathbb{R}^{r\times r}\) is a covariance matrix. Details are in the manuscript ``Mixed-type multivariate response regression with covariance estimation" by Ekvall and Molstad.
The package currently supports responses such that, for every \(i = 1, \dots, n\) and \(j = 1, \dots, r\), \(y_{ij} \mid w_i\) is either (a) normal with mean \(w_{ij}\) and variance \(\psi_j\), (b) quasi-Poisson with mean \(\exp(w_{ij})\) and variance \(\psi_j \exp(w_{ij})\), or (c) Bernoulli with success probability \(1 / \{1 + \exp(-w_{ij})\}\). The vector \(\psi = (\psi_1, \dots, \psi_r)\) is assumed known and provided as an argument to the fitting function.
If \(y_{ij}\) is a Bernoulli response, \(\psi_j = 1\) and \(\Sigma_{jj} = 1\) is enforced for identifiability.
We generate \(n = 100\) observations of a \(3\)-dimensional response vector and then fit the model. The matrix of predictors is of size \(nr\times p\), where the first \(r\) rows correspond to the design matrix \(X_1\) for the response vector \(y_1\), the next \(r\) rows correspond to \(X_2\), and so on. That is, \(X\) is the design matrix for the vector of all responses \([y_1', \dots, y_n']' \in \mathbb{R}^{nr}\). We take \(X_i = I_3\), the \(3\times 3\) identity matrix, corresponding to a separate intercept for each of the three responses. The responses are stored in an \(n\times r\) matrix where the \(i\)th row is the \(i\)th response vector. The first response is normal, the second Bernoulli, and the third Poisson (conditionally on the latent vector).
library(mmrr)
# Generate data
n <- 100; r <- 3
set.seed(4)
X <- kronecker(rep(1, n), diag(3)) # Each response has its own intercept
type <- 1:3 # First response normal, second Bernoulli, third Poisson
psi <- rep(1, 3) # psi_j = 1 for all j
Beta0 <- runif(3, -1, 1) # True coefficient vector
Sigma0 <- cov2cor(crossprod(matrix(rnorm(9), 3, 3))) + diag(c(1.5, 0, 0.5)) # True covariance matrix
R0 <- chol(Sigma0) # Function for generating data uses Cholesky root
Y <- generate_mmrr(X = X, Beta = Beta0, R = R0, type = type, psi = psi)
# Fit model
fit <- mmrr(Y = Y, X = X, type = type, psi = psi)
# Get estimates
fit$Beta## [1] 0.1318782 -0.6950476 -0.2211912fit$Sigma## [,1] [,2] [,3]
## [1,] 2.822656 1.6749202 1.0834566
## [2,] 1.674920 1.0000005 0.7222296
## [3,] 1.083457 0.7222296 1.4425031Suppose we want to test whether the Normal and Bernoulli variable are independent. That is equivalent to testing whether \(\Sigma_{12} = \Sigma_{21} = 0\). We fit the null model by incorporating the elementwise restrictions on \(\Sigma\) using the argument \(M\). We also restrict \(\Sigma_{22} = 1\) for identifiability because the second response is Bernoulli (this restriction is automatic when the user does not supply \(M\)).
M <- matrix(NA, 3, 3) # NA means no restriction
M[1, 2] <- M[2, 1] <- 0 # Set covariance restriction
M[2, 2] <- 1 # Set variance restriction
# Inspect M
M## [,1] [,2] [,3]
## [1,] NA 0 NA
## [2,] 0 1 NA
## [3,] NA NA NA# Fit model
fit_null <- mmrr(Y = Y, X = X, type = type, psi = psi, M = M)
# Inspect null estimate
fit_null$Sigma## [,1] [,2] [,3]
## [1,] 2.753878e+00 1.498801e-15 0.9691641
## [2,] 1.554312e-15 1.000000e+00 0.2859969
## [3,] 9.691641e-01 2.859969e-01 1.2902496# Do approximate likelihood ratio test
lrt_approx(fit_null = fit_null, fit_full = fit)## $stat
## [1] 14.5007
##
## $p
## [1] 0.0001401071
##
## $df
## [1] 1Because the predictors (an intercept only) do not depend on \(i\), the marginal moments of the \(y_i \in \mathbb{R}^3\) do not depend on \(i\). With the estimates of \(\beta\) and \(\Sigma\), we can get an estimate of the mean vector and covariance matrix of the \(y_i\) as follows.
# True marginal moments
mean_y <- predict_mmrr(X = diag(3), Beta = Beta0,
sigma = sqrt(diag(Sigma0)),
type = type, num_nodes = 10)
cov_y <- cov_mmrr(X = diag(3), Beta = Beta0, Sigma = Sigma0, psi = psi,
type = type, num_nodes = 10)
# Estimates of marginal moments
mean_y_hat <- predict_mmrr(X = diag(3), Beta = fit$Beta,
sigma = sqrt(diag(fit$Sigma)),
type = type, num_nodes = 10)
cov_y_hat <- cov_mmrr(X = diag(3), Beta = fit$Beta, Sigma = fit$Sigma, psi = psi,
type = type, num_nodes = 10)
# Compare
mean_y; mean_y_hat## [,1] [,2] [,3]
## [1,] 0.1716006 0.3064572 1.40141## [,1] [,2] [,3]
## [1,] 0.1318782 0.3597954 1.648821cov_y; cov_y_hat## [,1] [,2] [,3]
## [1,] 3.5000000 0.1697016 1.0248625
## [2,] 0.1697016 0.2125412 0.1959248
## [3,] 1.0248625 0.1959248 8.2392785## [,1] [,2] [,3]
## [1,] 3.8226563 0.3218452 1.7864258
## [2,] 0.3218452 0.2303427 0.2404325
## [3,] 1.7864258 0.2404325 10.4333959