Testes de hipótese em Modelos Multivariados de Covariância Linear Generalizada (McGLM)

Modelos Multivariados de Covariância Linear Generalizada

Lineu Alberto Cavazani de Freitas

Introdução

Os Modelos Linerares Generalizados (GLM), propostos por Nelder and Wedderburn (1972), são uma forma de modelagem univariada para dados de diferentes naturezas, tais como respostas contínuas, binárias e contagens. Tais características tornam essa classe de modelos uma flexível ferramenta de modelagem aplicável a diversos tipos de problema. Contudo, por mais flexível e discutida na literatura, essa classe apresenta duas principais restrições:

  1. A incapacidade de lidar com observações dependentes.

  2. E/ou a incapacidade de lidar múltiplas respostas simultaneamente.

Com o objetivo de contornar estas restrições, foi proposta por Bonat and Jørgensen (2016), uma estrutura geral para análise de dados não gaussianos com múltiplas respostas em que não se faz suposições quanto à independência das observações, os chamados Modelos Multivariados de Covariância Linear Generalizada (McGLM).

Vamos discutir os McGLM como uma extensão dos GLM. Vale ressaltar que é usada uma especificação menos usual de um Modelo Linear Generalizado, porém trata-se de uma notação mais conveniente para chegar à uma especificação melhor construída de um Modelo Multivariado de Covariância Linear Generalizada.

Sendo assim, neste material serão tratados os seguintes assuntos:

  • Especificação de um GLM.

  • Especificação de um cGLM.

  • Especificação de um McGLM.

  • Estimação e inferência dos parâmetros de um McGLM.

GLM

Seja \(\boldsymbol{Y}\) um vetor \(N \times 1\) de valores observados da variável resposta, \(\boldsymbol{X}\) uma matriz de delineamento \(N \times k\) e \(\boldsymbol{\beta}\) um vetor de parâmetros de regressão \(k \times 1\), um GLM pode ser descrito da forma

\[\begin{equation} \begin{aligned} \mathrm{E}(\boldsymbol{Y}) &= \boldsymbol{\mu} = g^{-1}(\boldsymbol{X} \boldsymbol{\beta}), \\ \mathrm{Var}(\boldsymbol{Y}) &= \Sigma = \mathrm{V}\left(\boldsymbol{\mu}; p\right)^{1/2}\left(\tau_0\boldsymbol{I}\right)\mathrm{V}\left(\boldsymbol{\mu}; p\right)^{1/2}, \end{aligned} \end{equation}\]

em que \(g(.)\) é a função de ligação, \(\mathrm{V}\left(\boldsymbol{\mu}; p\right)\) é uma matriz diagonal em que as entradas principais são dadas pela função de variância aplicada ao vetor \(\boldsymbol{\mu}\), \(p\) é o parâmetro de potência, \(\tau_0\) o parâmetro de dispersão e \(\boldsymbol{I}\) é a matriz identidade de ordem \(N\times N\).

Os GLM fazem uso de apenas duas funções, a função de variância e de ligação. Diferentes escolhas de funções de variância implicam em diferentes suposições a respeito da distribuição da variável resposta. Dentre as funções de variância conhecidas, podemos citar:

  1. A função de variância potência, que caracteriza a família Tweedie de distribuições, em que a função de variância é dada por \(\vartheta\left(\boldsymbol{\mu}; p\right) = \mu^p\), na qual destacam-se a distribuições: Normal (\(p\) = 0), Poisson (\(p\) = 1), gama (\(p\) = 2) e Normal inversa (\(p\) = 3) (Jørgensen 1987, 1997).

  2. A função de dispersão Poisson–Tweedie, a qual caracteriza a família Poisson-Tweedie de distribuições, que visa contornar a inflexibilidade da utilização da função de variância potência para respostas discretas. A família Poisson-Tweedie tem função de dispersão dada por \(\vartheta\left(\boldsymbol{\mu}; p\right) = \mu + \mu^p\) e tem como casos particulares os mais famosos modelos para dados de contagem: Hermite (\(p\) = 0), Neyman tipo A (\(p\) = 1), binomial negativa (\(p\) = 2) e Poisson–inversa gaussiana (p = \(3\)) (Jørgensen and Kokonendji 2015).

  3. A função de variância binomial, dada por \(\vartheta(\boldsymbol{\mu}) = \mu(1 - \mu)\), utilizada quando a variável resposta é binária, restrita a um intervalo ou quando tem-se o número de sucessos em um número de tentativas.

Lembre-se que o GLM é uma classe de modelos de regressão univariados em que um dos pressupostos é a indepenência entre as observações. Esta independência é especificada na matriz identidade no centro da equação que define a matriz de variância e covariância. Podemos imaginar que, substituindo esta matriz identidade por uma matriz qualquer que reflita a relação entre os indivíduos da amostra teremos uma extensão do Modelo Linear Generalizado para observações dependentes. É justamente essa a ideia dos Modelos de Covariância Linear Generalizada, o cGLM.

cGLM

Os cGLM são uma alternativa para problemas em que a suposição de independência entre as observações não é atendida. Neste caso, a solução proposta é substituir a matriz identidade \(\boldsymbol{I}\) da equação que descreve a matriz de variância e covariância por uma matriz não diagonal \(\boldsymbol{\Omega({\tau})}\) que descreva adequadamente a estrutura de correlação entre as observações. Trata-se de uma ideia similar à proposta de Liang and Zeger (1986) nos modelos GEE (Generalized Estimation Equation), em que utiliza-se uma matriz de correlação de trabalho para considerar a dependência entre as observações. A matriz \(\boldsymbol{\Omega({\tau})}\) é descrita como uma combinação de matrizes conhecidas tal como nas propostas de Anderson and others (1973) e Pourahmadi (2000), podendo ser escrita da forma

\[\begin{equation} h\left \{ \boldsymbol{\Omega}(\boldsymbol{\tau}) \right \} = \tau_0Z_0 + \ldots + \tau_DZ_D, \end{equation}\]

em que \(h(.)\) é a função de ligação de covariância, \(Z_d\) com \(d\) = 0,\(\ldots\), D são matrizes que representam a estrutura de covariância presente nos dados e \(\boldsymbol{\tau}\) = \((\tau_0, \ldots, \tau_D)\) é um vetor \((D + 1) \times 1\) de parâmetros de dispersão. Tal estrutura pode ser vista como um análogo ao preditor linear para a média e foi nomeado como preditor linear matricial. A especificação da função de ligação de covariância é discutida por Pinheiro and Bates (1996) e é possível selecionar combinações de matrizes para se obter os mais conhecidos modelos da literatura para dados longitudinais, séries temporais, dados espaciais e espaço-temporais. Maiores detalhes são discutidos por Demidenko (2013).

Com isso, substituindo a matriz identidade pela equação do preditor linear matricial, temos uma classe com toda a flexibilidade dos GLM, porém contornando a restrição da independência entre as observações desde que o preditor linear matricial seja adequadamente especificado.

Deste modo, é contornada a primeira restrição dos GLM. A segunda restrição diz respeito às múltiplas respostas e, resolvendo esta restrição, chegamos ao McGLM.

McGLM

O McGLM pode ser entendido como uma extensão multivariada do cGLM e contorna as duas principais restrições presentes nos GLM, pois além de permitir a modelagem de dados com estrutura de covariância, permite modelar múltiplas respostas.

Considre \(\boldsymbol{Y}_{N \times R} = \left \{ \boldsymbol{Y}_1, \dots, \boldsymbol{Y}_R \right \}\) uma matriz de variáveis resposta e \(\boldsymbol{M}_{N \times R} = \left \{ \boldsymbol{\mu}_1, \dots, \boldsymbol{\mu}_R \right \}\) uma matriz de valores esperados. Cada uma das variáveis resposta tem sua própria matriz de variância e covariância, responsável por modelar a covariância dentro de cada resposta, sendo expressa por

\[\begin{equation} \Sigma_r = \mathrm{V}_r\left(\boldsymbol{\mu}_r; p\right)^{1/2}\boldsymbol{\Omega}_r\left(\boldsymbol{\tau}\right)\mathrm{V}_r\left(\boldsymbol{\mu}_r; p\right)^{1/2}. \end{equation}\]

Além disso, é necessária uma matriz de correlação \(\Sigma_b\), de ordem \(R \times R\), que descreve a correlação entre as variáveis resposta. Para a especificação da matriz de variância e covariância conjunta é utilizado o produto Kronecker generalizado, proposto por Martinez-Beneito (2013).

Finalmente, um MCGLM é descrito como

\[ \begin{equation} \label{eq:mcglm} \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)\} \\ \mathrm{Var}(\boldsymbol{Y}) &= \boldsymbol{C} = \boldsymbol{\Sigma}_R \overset{G} \otimes \boldsymbol{\Sigma}_b, \end{aligned} \end{equation} \] em que \(\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)\) é o produto generalizado de Kronecker, a matriz \(\tilde{\boldsymbol{\Sigma}}_r\) denota a matriz triangular inferior da decomposição de Cholesky da matriz \({\boldsymbol{\Sigma}}_r\), o operador \(\mathrm{Bdiag}\) denota a matriz bloco-diagonal e \(\boldsymbol{I}\) uma matriz identidade \(N \times N\).

Toda metodologia do McGLM está implementada no pacote mcglm (Bonat 2018) do software estatístico R (R Core Team 2018).

Estimação e Inferência

Os McGLMs são ajustados baseados no método de funções de estimação descritos em detalhes por Bonat and Jørgensen (2016) e Jørgensen and Knudsen (2004). Nesta seção é apresentada uma visão geral do algoritmo e da distribuição assintótica dos estimadores baseados em funções de estimação.

As suposições de segundo momento dos McGLM permitem a divisão dos parâmetros em dois conjuntos: \(\boldsymbol{\theta} = (\boldsymbol{\beta}^{\top}, \boldsymbol{\lambda}^{\top})^{\top}\). Desta forma, \(\boldsymbol{\beta} = (\boldsymbol{\beta}_1^\top, \ldots, \boldsymbol{\beta}_R^\top)^\top\) é um vetor \(K \times 1\) de parâmetros de regressão e \(\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\) é um vetor \(Q \times 1\) de parâmetros de dispersão. Além disso, \(\mathcal{Y} = (\boldsymbol{Y}_1^\top, \ldots, \boldsymbol{Y}_R^\top)^\top\) denota o vetor empilhado de ordem \(NR \times 1\) da matriz de variáveis resposta \(\boldsymbol{Y}_{N \times R}\) e \(\mathcal{M} = (\boldsymbol{\mu}_1^\top, \ldots, \boldsymbol{\mu}_R^\top)^\top\) denota o vetor empilhado de ordem \(NR \times 1\) da matriz de valores esperados \(\boldsymbol{M}_{N \times R}\).

Para estimação dos parâmetros de regressão é utilizada a função quasi-score (Liang and Zeger 1986), representada por \[\begin{equation} \begin{aligned} \psi_{\boldsymbol{\beta}}(\boldsymbol{\beta}, \boldsymbol{\lambda}) = \boldsymbol{D}^\top \boldsymbol{C}^{-1}(\mathcal{Y} - \mathcal{M}), \end{aligned} \end{equation}\] em que \(\boldsymbol{D} = \nabla_{\boldsymbol{\beta}} \mathcal{M}\) é uma matriz \(NR \times K\), e \(\nabla_{\boldsymbol{\beta}}\) denota o operador gradiente. Utilizando a função quasi-score a matriz \(K \times K\) de sensitividade de \(\psi_{\boldsymbol{\beta}}\) é dada por \[\begin{equation} \begin{aligned} S_{\boldsymbol{\beta}} = E(\nabla_{\boldsymbol{\beta} \psi \boldsymbol{\beta}}) = -\boldsymbol{D}^{\top} \boldsymbol{C}^{-1} \boldsymbol{D}, \end{aligned} \end{equation}\] enquanto que a matriz \(K \times K\) de variabilidade de \(\psi_{\boldsymbol{\beta}}\) é escrita como \[\begin{equation} \begin{aligned} V_{\boldsymbol{\beta}} = VAR(\psi \boldsymbol{\beta}) = \boldsymbol{D}^{\top} \boldsymbol{C}^{-1} \boldsymbol{D}. \end{aligned} \end{equation}\]

Para os parâmetros de dispersão é utilizada a função de estimação de Pearson, definida da forma \[\begin{equation} \begin{aligned} \psi_{\boldsymbol{\lambda}_i}(\boldsymbol{\beta}, \boldsymbol{\lambda}) = \mathrm{tr}(W_{\boldsymbol{\lambda}i} (\boldsymbol{r}^\top\boldsymbol{r} - \boldsymbol{C})), i = 1,.., Q, \end{aligned} \end{equation}\] em que \(W_{\boldsymbol{\lambda}i} = -\frac{\partial \boldsymbol{C}^{-1}}{\partial \boldsymbol{\lambda}_i}\) e \(\boldsymbol{r} = (\mathcal{Y} - \mathcal{M})\). A entrada \((i,j)\) da matriz de sensitividade \(Q \times Q\) de \(\psi_{\boldsymbol{\lambda}}\) é dada por \[\begin{equation} \begin{aligned} S_{\boldsymbol{\lambda_{ij}}} = E \left (\frac{\partial }{\partial \boldsymbol{\lambda_{i}}} \psi \boldsymbol{\lambda_{j}}\right) = -tr(W_{\boldsymbol{\lambda_{i}}} CW_{\boldsymbol{\lambda_{J}}} C). \end{aligned} \end{equation}\] Já a entrada \((i,j)\) da matriz de variabilidade \(Q \times Q\) de \(\psi_{\boldsymbol{\lambda}}\) é definida por \[\begin{equation} \begin{aligned} V_{\boldsymbol{\lambda_{ij}}} = Cov\left ( \psi_{\boldsymbol{\lambda_{i}}}, \psi_{\boldsymbol{\lambda_{j}}} \right) = 2tr(W_{\boldsymbol{\lambda_{i}}} CW_{\boldsymbol{\lambda_{J}}} C) + \sum_{l=1}^{NR} k_{l}^{(4)} (W_{\boldsymbol{\lambda_{i}}})_{ll} (W_{\boldsymbol{\lambda_{j}}})_{ll}, \end{aligned} \end{equation}\] em que \(k_{l}^{(4)}\) denota a quarta cumulante de \(\mathcal{Y}_{l}\). No processo de estimação dos McGLM são usadas as versões empíricas.

Para se levar em conta a covariância entre os vetores \(\boldsymbol{\beta}\) e \(\boldsymbol{\lambda}\), Bonat and Jørgensen (2016) obtiveram as matrizes de sensitividade e variabilidade cruzadas, denotadas por \(S_{\boldsymbol{\lambda \beta}}\), \(S_{\boldsymbol{\beta \lambda}}\) e \(V_{\boldsymbol{\lambda \beta}}\), mais detalhes em Bonat and Jørgensen (2016). As matrizes de sensitividade e variabilidade conjuntas de \(\psi_{\boldsymbol{\beta}}\) e \(\psi_{\boldsymbol{\lambda}}\) são denotados por

\[\begin{equation} \begin{aligned} S_{\boldsymbol{\theta}} = \begin{bmatrix} S_{\boldsymbol{\beta}} & S_{\boldsymbol{\beta\lambda}} \\ S_{\boldsymbol{\lambda\beta}} & S_{\boldsymbol{\lambda}} \end{bmatrix} \text{e } V_{\boldsymbol{\theta}} = \begin{bmatrix} V_{\boldsymbol{\beta}} & V^{\top}_{\boldsymbol{\lambda\beta}} \\ V_{\boldsymbol{\lambda\beta}} & V_{\boldsymbol{\lambda}} \end{bmatrix}. \end{aligned} \end{equation}\]

Seja \(\boldsymbol{\hat{\theta}} = (\boldsymbol{\hat{\beta}^{\top}}, \boldsymbol{\hat{\lambda}^{\top}})^{\top}\) o estimador baseado em funções de estimação de \(\boldsymbol{\theta}\). Então, a distribuição assintótica de \(\boldsymbol{\hat{\theta}}\) é \[\begin{equation} \begin{aligned} \boldsymbol{\hat{\theta}} \sim N(\boldsymbol{\theta}, J_{\boldsymbol{\theta}}^{-1}), \end{aligned} \end{equation}\] em que \(J_{\boldsymbol{\theta}}^{-1}\) é a inversa da matriz de informação de Godambe, dada por \(J_{\boldsymbol{\theta}}^{-1} = S_{\boldsymbol{\theta}}^{-1} V_{\boldsymbol{\theta}} S_{\boldsymbol{\theta}}^{-\top}\), em que \(S_{\boldsymbol{\theta}}^{-\top} = (S_{\boldsymbol{\theta}}^{-1})^{\top}.\)

Para resolver o sistema de equações \(\psi_{\boldsymbol{\beta}} = 0\) e \(\psi_{\boldsymbol{\lambda}} = 0\) faz-se uso do algoritmo Chaser modificado, proposto por Jørgensen and Knudsen (2004), que fica definido como

\[\begin{equation} \begin{aligned} \begin{matrix} \boldsymbol{\beta}^{(i+1)} = \boldsymbol{\beta}^{(i)}- S_{\boldsymbol{\beta}}^{-1} \psi \boldsymbol{\beta} (\boldsymbol{\beta}^{(i)}, \boldsymbol{\lambda}^{(i)}), \\ \boldsymbol{\lambda}^{(i+1)} = \boldsymbol{\lambda}^{(i)}\alpha S_{\boldsymbol{\lambda}}^{-1} \psi \boldsymbol{\lambda} (\boldsymbol{\beta}^{(i+1)}, \boldsymbol{\lambda}^{(i)}). \end{matrix} \end{aligned} \end{equation}\]

Anderson, TW, and others. 1973. “Asymptotically Efficient Estimation of Covariance Matrices with Linear Structure.” The Annals of Statistics 1 (1). Institute of Mathematical Statistics: 135–41.

Bonat, Wagner Hugo. 2018. “Multiple Response Variables Regression Models in R: The mcglm Package.” Journal of Statistical Software 84 (4): 1–30. https://doi.org/10.18637/jss.v084.i04.

Bonat, Wagner Hugo, and Bent Jørgensen. 2016. “Multivariate Covariance Generalized Linear Models.” Journal of the Royal Statistical Society: Series C (Applied Statistics) 65 (5). Wiley Online Library: 649–75.

Demidenko, Eugene. 2013. Mixed Models: Theory and Applications with R. John Wiley & Sons.

Jørgensen, Bent. 1987. “Exponential Dispersion Models.” Journal of the Royal Statistical Society: Series B (Methodological) 49 (2). Wiley Online Library: 127–45.

———. 1997. The Theory of Dispersion Models. CRC Press.

Jørgensen, Bent, and Sven Jesper Knudsen. 2004. “Parameter Orthogonality and Bias Adjustment for Estimating Functions.” Scandinavian Journal of Statistics 31 (1). Wiley Online Library: 93–114.

Jørgensen, Bent, and Célestin C Kokonendji. 2015. “Discrete Dispersion Models and Their Tweedie Asymptotics.” AStA Advances in Statistical Analysis 100 (1). Springer: 43–78.

Liang, Kung-Yee, and Scott L Zeger. 1986. “Longitudinal Data Analysis Using Generalized Linear Models.” Biometrika 73 (1). Oxford University Press: 13–22.

Martinez-Beneito, Miguel A. 2013. “A General Modelling Framework for Multivariate Disease Mapping.” Biometrika 100 (3). Oxford University Press: 539–53.

Nelder, John Ashworth, and Robert William Maclagan Wedderburn. 1972. “Generalized Linear Models.” Journal of the Royal Statistical Society. Series A (General) 135: 370–84.

Pinheiro, José C, and Douglas M Bates. 1996. “Unconstrained Parametrizations for Variance-Covariance Matrices.” Statistics and Computing 6 (3). Springer: 289–96.

Pourahmadi, Mohsen. 2000. “Maximum Likelihood Estimation of Generalised Linear Models for Multivariate Normal Covariance Matrix.” Biometrika 87 (2). Oxford University Press: 425–35.

R Core Team. 2018. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Drawing             Drawing