General linear model

• MANOVA can deal only with Gaussian data.
• Homogeneity of variances and covariances.
• Independent observations within response variables.
• What do we need?
  • Generalized linear models (GLM) to deal with non-Gaussian data.
  • Flexibility to deal with different types of dependence structures.
  • Multivariate response variables, possibly of mixed types.
  • Easy to specify, fit, interpret and teach.

Outline

• Goal 1: Make GLM more flexible to deal with:
  • Non-negative highly right-skewed data and continuous data with probability mass at zero (Bonat e Kokonendji 2017).
  • Under, equi, overdispersed, zero-inflated and heavy-tailed count data (Bonat, Jørgensen, et al. 2017, Bonat, Olivero, et al. (2017)).
  • Bounded data (Bonat e Peterle 2017).
• Goal 2: Extend GLM to cGLM to deal with:
  • Mixed models (Bates et al. 2015).
  • Repeated measures (Diggle et al. 2002).
  • Time series (Box e Jenkins 1994).
  • Spatial and space-time data (Cressie e Huang 1999; Cressie e Wikle 2011).
  • Genetic and Twin data (Sorensen e Gianola 2007).
• Goal 3: Extend cGLM to McGLM to deal with:
  • Multivariate response variables (Bonat e Jørgensen 2016).
  • Response variables of mixed types (Bonat 2017a).
  • Identify the general linear model and MANOVA as special cases.
• Goal 4: Extend multivariate hypotheses tests to non-Gaussian data.
• Goal 5: Provide software implementation in R (Bonat 2017b).

Flexible Generalized Linear Models

• Let \(\boldsymbol{Y}\) be an \(N \times 1\) response vector.
• Let \(\boldsymbol{X}\) be an \(N \times k\) design matrix.
• Let \(\boldsymbol{\beta}\) be a \(k \times 1\) parameter vector.
• Orthodox GLM \[ \begin{aligned} \label{modelGLM} \mathrm{E}(\boldsymbol{Y}) &= \boldsymbol{\mu} = g^{-1}(\boldsymbol{X} \boldsymbol{\beta}) \nonumber \\ \mathrm{Var}(\boldsymbol{Y}) &= \boldsymbol{\Sigma} = \mathrm{V}(\boldsymbol{\mu};p)^{\frac{1}{2}} (\tau_0 \boldsymbol{I}) \mathrm{V}(\boldsymbol{\mu};p)^{\frac{1}{2}}. \quad \quad \quad (1) \end{aligned} \]
• \(g\) is the link function.
• \(\mathrm{V}(\boldsymbol{\mu};p) = \mathrm{diag}(\vartheta(\boldsymbol{\mu};p))\) where \(\vartheta(\boldsymbol{\mu};p)\) is the variance function.
• \(p\) and \(\tau_0\) are the power and dispersion parameters, respectively.
• \(\boldsymbol{I}\) is an identity matrix.
• The standard linear model is obtained for identity link and constant variance functions.

Variance/dispersion functions

• Tweedie family (variance function), characterized by \[ \vartheta(\boldsymbol{\mu};p) = \mu^p. \]
• Gaussian (\(p = 0\)), gamma (\(p = 2\)) and inverse Gaussian (\(p = 3\)).
• Deals with symmetric, skewed and zero-inflated continuous outcomes.
• Poisson-Tweedie family (dispersion function), characterized by \[ \vartheta(\boldsymbol{\mu};p) = \mu + \mu^p. \]
• Hermite (\(p=0\)), Neyman Type A (\(p=1\)), negative binomial (\(p=2\)).
• Special form in (1) \[ \boldsymbol{\Sigma} = \mathrm{diag}(\boldsymbol{\mu}) + \mathrm{V}(\boldsymbol{\mu};p)^{\frac{1}{2}} (\tau_0 \boldsymbol{I}) \mathrm{V}(\boldsymbol{\mu};p)^{\frac{1}{2}}. \]
• Deals with over-/under-dispersion and zero-inflated count outcomes.
• Binomial family, characterized by \[ \vartheta(\boldsymbol{\mu}) = \mu(1-\mu). \]
• Deals with binary, binomial and proportion outcomes.
• Extended binomial variance function, \[ \vartheta(\boldsymbol{\mu}; p, q) = \mu^p(1-\mu)^q. \]
• Extra flexibility to deal with binary, binomial and proportion outcomes.
• Three sets of variance functions.
• Deals with the most common outcomes types.
• Power parameter identifies the distribution.
• Estimation of power parameter (model selection).
• Additional flexibility.
• Deals only with independent observations.

Covariance Generalized Linear Models

• Change the identity matrix in (1) to a non-diagonal matrix \(\boldsymbol{\Omega}(\boldsymbol{\tau})\). \[ \mathrm{Var}(\boldsymbol{Y}) = \boldsymbol{\Sigma} = \mathrm{V}(\boldsymbol{\mu};p)^{\frac{1}{2}} (\boldsymbol{\Omega}(\boldsymbol{\tau})) \mathrm{V}(\boldsymbol{\mu};p)^{\frac{1}{2}}. \]
• Similar to working correlation matrix (GEE) (Liang e Zeger 1986).
• Model \(\boldsymbol{\Omega}(\boldsymbol{\tau})\) as a linear combination of known matrices. \[ h(\boldsymbol{\Omega}(\boldsymbol{\tau})) = \tau_0 Z_0 + \cdots + \tau_D Z_D, \] where \(h\) is a covariance link function (Pourahmadi 2000).
• \(Z_d\) with \(d = 0, \ldots, D\) are known matrices.
• \(\boldsymbol{\tau} = (\tau_0,\cdots, \tau_D)\) is a \(D \times 1\) dispersion parameter vector.
• Identity, inverse, exponential-matrix, modified Cholesky decomposition, etc.
• Deals only with one response variable.
• Examples of covariance linear models:
  • Double generalized linear models.
  • Linear mixed models.
  • Moving average models.
  • Exchangeable or compound symmetry and unstructured models (popular in longitudinal data analysis).
  • Conditional autoregressive models (time series, spatial and space-time data).
  • Models in quantitative genetic and phylogenetic.
  • Models for Twin and family data.
  • and many more.

Example 1: Gaussian linear mixed model

• Gaussian linear mixed model \[ \boldsymbol{Y} = \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{U} \boldsymbol{\gamma} + \boldsymbol{\epsilon} \] where \(\boldsymbol{\gamma} \sim N(\boldsymbol{0},\boldsymbol{D})\) and \(\boldsymbol{\epsilon} \sim N(\boldsymbol{0},\tau_0 \boldsymbol{I})\).
• It is easy to show that \[ \begin{aligned} \mathrm{E}(\boldsymbol{Y}) &= \boldsymbol{\mu} = \boldsymbol{X} \boldsymbol{\beta} \nonumber \\ \mathrm{Var}(\boldsymbol{Y}) &= \boldsymbol{\Omega} = \tau_0 \boldsymbol{I} + \boldsymbol{U}\boldsymbol{D}\boldsymbol{U}^{\top}. \nonumber \end{aligned} \]
• Furthermore, \[ \boldsymbol{U}\boldsymbol{D}\boldsymbol{U}^{\top} = \sum_{j \leq q}^k d_{jq} \boldsymbol{V}_{jq} \] where \[ \boldsymbol{V}_{jj} = \boldsymbol{U_{\cdot j}} \boldsymbol{U_{\cdot j}}^{\top}, \quad \boldsymbol{V}_{jq} = \boldsymbol{U_{\cdot j}} \boldsymbol{U_{\cdot q}}^{\top} + \boldsymbol{U_{\cdot q}} \boldsymbol{U_{\cdot j}}^{\top} \] for \(j \neq q\) and \(\boldsymbol{U_{\cdot j}}\) denotes the j\(th\) column of \(\boldsymbol{U}\).

Example 2: exchangeable and unstructured models

• Exchangeable \[ Z_1 = \begin{bmatrix} 1 & 1 & 1\\ 1 & 1 & 1\\ 1 & 1 & 1 \end{bmatrix}. \]
• Unstructured \[ Z_1 = \begin{bmatrix} 0 & 1 & 0\\ 1 & 0 & 0\\ 0 & 0 & 0 \end{bmatrix}, \quad Z_2 = \begin{bmatrix} 0 & 0 & 1\\ 0 & 0 & 0\\ 1 & 0 & 0 \end{bmatrix} \quad \text{and} \quad Z_3 = \begin{bmatrix} 0 & 0 & 0\\ 0 & 0 & 1\\ 0 & 1 & 0 \end{bmatrix}. \]

Example 3: conditional autoregressive models

• Specify the inverse of the covariance matrix \[ \boldsymbol{\Omega}^{-1}(\tau,\rho) = \tau (\boldsymbol{D} - \rho \boldsymbol{W}), \] where \(\boldsymbol{W}\) is a neighbourhood matrix and \(\boldsymbol{D}\) is a diagonal matrix with the number of neighbours.
• Linear covariance models using the inverse covariance link function \[ \boldsymbol{\Omega}^{-1}(\boldsymbol{\tau}) = \tau_0 \boldsymbol{D} + \tau_1 \boldsymbol{W}, \] where \(\tau_0 = \tau\) and \(\tau_1 = - \tau \rho\).
• \(\boldsymbol{W}\) can describe temporal, spatial and spatial temporal neighborhoods.

Example 4: genetic additive models

• Let \(\boldsymbol{A}\) denote an additive genetic relatedness matrix, \[ \boldsymbol{\Omega}(\boldsymbol{\tau}) = \tau_0 \boldsymbol{I} + \tau_1 \boldsymbol{A}, \] where \(\tau_0\) and \(\tau_1\) denote the environment and genetic variance components, respectively. For details see Bonat (2017a).

Multivariate Covariance Generalized Linear Models

• Let \(\boldsymbol{Y}_{N \times R} = \{\boldsymbol{Y}_1, \ldots, \boldsymbol{Y}_R\}\) be an outcome matrix.
• Let \(\boldsymbol{M}_{N \times R} = \{\boldsymbol{\mu}_1, \ldots, \boldsymbol{\mu}_R\}\) be an expected value matrix.
• Let \(\boldsymbol{\Sigma}_r = \mathrm{V}_r(\boldsymbol{\mu};p)^{\frac{1}{2}} \boldsymbol{\Omega}_r(\boldsymbol{\tau}) \mathrm{V}_r(\boldsymbol{\mu};p)^{\frac{1}{2}}\) be the covariance matrix within outcomes.
• Let \(\boldsymbol{\Sigma}_b\) be the correlation matrix between outcomes.
• We define the McGLM by, \[ \begin{aligned} \mathrm{E}(\boldsymbol{Y}) &= \boldsymbol{M} = \{g_1^{-1}(\boldsymbol{X}_1 \boldsymbol{\beta}_1), \ldots, g_R^{-1}(\boldsymbol{X}_R \boldsymbol{\beta}_R)\} \nonumber \\ \mathrm{Var}(\boldsymbol{Y}) &= \boldsymbol{C} = \boldsymbol{\Sigma}_R \overset{G} \otimes \boldsymbol{\Sigma}_b \nonumber \end{aligned} \] where \[ \boldsymbol{\Sigma}_R \overset{G} \otimes \boldsymbol{\Sigma}_b = \mathrm{Bdiag}(\tilde{\boldsymbol{\Sigma}}_1, \ldots, \tilde{\boldsymbol{\Sigma}}_R) (\boldsymbol{\Sigma}_b \otimes \boldsymbol{I}) \mathrm{Bdiag}(\tilde{\boldsymbol{\Sigma}}_1^\top, \ldots, \tilde{\boldsymbol{\Sigma}}_R^\top). \]
• Generalized Kronecker product proposed by Martinez-Beneito (2013).
• \(\tilde{\boldsymbol{\Sigma}}_r\) is the lower triangular matrix of the Cholesky decomposition of \(\boldsymbol{\Sigma}_r\).
• \(\mathrm{Bdiag}\) denotes a block diagonal matrix.
• General linear model is easily obtained by taking
  • Identity link functions.
  • \(\boldsymbol{X}_r = \boldsymbol{X}\) for all \(r\).
  • Constant variance function.
  • \(\boldsymbol{\Sigma}_r = \tau_{0r} \boldsymbol{I}\) for all \(r\).

McGLM - Parametrization

• Let \(\mathcal{Y} = (\boldsymbol{Y}_1^\top, \ldots, \boldsymbol{Y}_R^\top)^\top\) be the stacked vector \((NR \times 1)\) of the outcome matrix \(\boldsymbol{Y}_{N \times R}\) by columns.
• Let \(\mathcal{M} = (\boldsymbol{\mu}_1^\top, \ldots, \boldsymbol{\mu}_R^\top)^\top\) be the stacked vector \((NR \times 1)\) of the expected value matrix \(\boldsymbol{M}_{N \times R}\) by columns.
• The \(\boldsymbol{C}\) matrix is symmetric and \(NR \times NR\).
• Let \(\boldsymbol{\beta} = (\boldsymbol{\beta}_1^\top, \ldots, \boldsymbol{\beta}_R^\top)^\top\) be a vector \(K \times 1\) of regression parameters.
• Let \(\boldsymbol{\lambda} = (\rho_1, \ldots, \rho_{R(R-1)/2}, p_1, \ldots, p_R, \boldsymbol{\tau}_1^\top, \ldots, \boldsymbol{\tau}_R^\top)^\top\) be a \(Q \times 1\) vector of dispersion parameters.
• McGLM are specified by two sets of parameters \(\boldsymbol{\theta} = (\boldsymbol{\beta}, \boldsymbol{\lambda})\).
• Quasi-score function for regression parameters.
• Pearson estimating function for dispersion parameters.

Estimating functions

• The quasi-score function is defined by \[ \psi_{\boldsymbol{\beta}}(\boldsymbol{\beta}, \boldsymbol{\lambda}) = \boldsymbol{D}^\top \boldsymbol{C}^{-1}(\mathcal{Y} - \mathcal{M}), \] where \(\boldsymbol{D} = \nabla_{\boldsymbol{\beta}} \mathcal{M}\) is an \(NR \times K\) matrix.
• The Pearson estimating function is defined by, \[ \psi_{\boldsymbol{\lambda}_i}(\boldsymbol{\beta}, \boldsymbol{\lambda}) = \mathrm{tr}(W_{\boldsymbol{\lambda}i} (\boldsymbol{r}^\top\boldsymbol{r} - \boldsymbol{C})), \] where \(W_{\boldsymbol{\lambda}i} = -\frac{\partial \boldsymbol{C}^{-1}}{\partial \boldsymbol{\lambda}_i}\) and \(\boldsymbol{r} = (\mathcal{Y} - \mathcal{M})\).

Fitting algorithm

• Modified chaser \[ \begin{aligned} \label{chaser} \boldsymbol{\beta}^{(i+1)} &= \boldsymbol{\beta}^{(i)} - \mathrm{S}_{\boldsymbol{\beta}}^{-1} \psi_{\boldsymbol{\beta}} (\boldsymbol{\beta}^{(i)}, \boldsymbol{\lambda}^{(i)}) \nonumber \\ \boldsymbol{\lambda}^{(i+1)} &= \boldsymbol{\lambda}^{(i)} - \alpha \mathrm{S}_{\boldsymbol{\lambda}}^{-1} \psi_{\boldsymbol{\lambda}}(\boldsymbol{\beta}^{(i+1)}, \boldsymbol{\lambda}^{(i)}). \nonumber \end{aligned} \]
• \(\mathrm{S}_{\boldsymbol{\beta}}\) and \(\mathrm{S}_{\boldsymbol{\lambda}}\) are the sensitivity matrices for \(\boldsymbol{\beta}\) and \(\boldsymbol{\lambda}\), respectively.

Godambe information matrix and asymptotic distribution

• Denote \(\hat{\boldsymbol{\theta}} = (\hat{\boldsymbol{\beta}}, \hat{\boldsymbol{\lambda}})\) the estimating function estimator of \(\boldsymbol{\theta}\).
• The asymptotic distribution of \(\hat{\boldsymbol{\theta}}\) is given by \[ \hat{\boldsymbol{\theta}} \sim \mathrm{N}_{K+Q}(\boldsymbol{\theta}, \mathrm{J}_{\boldsymbol{\theta}}^{-1}) \] where \(\mathrm{J}_{\boldsymbol{\theta}}^{-1}\) is the inverse of the Godambe information matrix, \[ \mathrm{J}_{\boldsymbol{\theta}}^{-1} = \mathrm{S}_{\boldsymbol{\theta}}^{-1} \mathrm{V}_{\boldsymbol{\theta}} \mathrm{S}_{\boldsymbol{\theta}}^{-\top}, \] where \(\mathrm{S}_{\boldsymbol{\theta}}^{-\top} = (\mathrm{S}_{\boldsymbol{\theta}}^{-1})^{\top}\).
• \(\mathrm{V}_{\boldsymbol{\theta}}\) is the variability matrix.
• For details, see Bonat e Jørgensen (2016) and Jørgensen e Knudsen (2004).

Multivariate linear hypotheses tests

• Hypotheses \[ H_{0}: \boldsymbol{L}\boldsymbol{\beta} = \boldsymbol{c}\quad \text{vs} \quad H_{1}:\boldsymbol{L}\boldsymbol{\beta} \neq \boldsymbol{c}. \]
• Wald statistics \[ W_{s} = (\boldsymbol{L}\boldsymbol{\beta})^{\top} (\boldsymbol{L} \mathrm{J}_{\boldsymbol{\beta}}^{-1} \boldsymbol{L})^{-1} (\boldsymbol{L}\boldsymbol{\beta}), \] where \(\mathrm{J}_{\boldsymbol{\beta}} = \boldsymbol{D}^{\top} \boldsymbol{C}^{-1}\boldsymbol{D}.\)
• \(\boldsymbol{L}\) is a \(JR \times K\) known matrix of constant.
• \(J\) denotes the number of regression coefficients under testing.
• Asymptotically \(W_s\) has a \(\chi^2\) distribution with \(JR\) degrees of freedom.
• \(W_{s}\) is a straightforward generalization of the Hotelling-Lawley statistics to multivariate non-Gaussian data.

Overview of the mcglm package.

• Package mcglm available at github and CRAN.
• Main reference Bonat (2017a).
• Instalation

library(devtools)

# From github
install_github("wbonat/mcglm")

# From CRAN
install.packages("mcglm")

• Main function mcglm.

library(mcglm)

## ----------------------------------------------------------------------
##   mcglm: Multivariate Covariance Generalized Linear Models
## 
##   For support, collaboration or bug report, visit: 
##     https://github.com/wbonat/mcglm
## 
##   mcglm version 0.4.0 (build on 2016-10-17) is now loaded.
## ----------------------------------------------------------------------

args(mcglm)

## function (linear_pred, matrix_pred, link, variance, covariance, 
##     offset, Ntrial, power_fixed, data, control_initial = "automatic", 
##     contrasts = NULL, control_algorithm = list()) 
## NULL

• Link functions: logit, probit, cauchit, cloglog, loglog, identity, log, sqrt and inverse.
• Variance functions: constant, tweedie, poisson_tweedie, binomialP and binomialPQ.
• Covariance link functions: identity, inverse and exponential-matrix.

Linear covariance structures

• The users can use any list of symmetric matrices.
• Combine pre-specified structures with new ones.

Methods for mcglm objects

Extra features for mcglm objects

• GOF’s implemented: plogLik, pAIC, pKLIC, pBIC, ESS, GOSHO and RJC.

Discussion

• Universal multivariate statistical modelling framework.
• Flexible class of models based on second-moment assumptions.
• General framework for estimation based on estimating function.
• Efficient algorithms for estimation.
• Asymptotic theory (Godambe Information).
• Multivariate hypotheses tests for non-Gaussian data.

Extensions

• Automatic covariate selection (Score Information Criterion-SIC).
• Automatic covariance selection.
• Measures of goodness-of-fit (Pseudo logLik, Pseudo AIC, Pseudo KLIC and ESS.)
• Simulation from McGLMs based on the NORTA algorithm.

Coming soon

• Residual analysis (improvements).
• Diagnostics (Leverage, DFBETA’s and Cook’s distance).
• Penalized estimating functions (high dimensional data and splines).,
• Prediction (time, space and space-time).
• Special module for Twin data analysis. (NordicTwin cancer project)
• Include parameter constraints to fit structural equation models (SEM).
• Improve package documentation.
• Multivariate hypotheses tests based on the pseudo-likelihood and score statistics.