Modelos Copula Regressão Copula Gaussiana


Referência principal:
Masarotto, G., & Varin, C. (2017).
Gaussian Copula Regression in R. Journal of Statistical Software, 77(8), 1–26.



1. Introdução


Os modelos copula (Joe, 2014) são frequentemente considerados para estender os modelos de regressão univariada assumindo respostas independentes para estruturas mais gerais. Por exemplo, Frees and Valdez (1998); Song (2000); Parsa and Klugman (2011); Kolev and Paiva (2009). O principal mérito da abordagem é que a especificação do modelo de regressão é separada da estrutura de dependência. Este artigo se concentra no método de regressão de copula gaussiana onde a dependência é convenientemente expressa na forma familiar da matriz de correlação de uma distribuição gaussiana multivariada (Song, 2000; Pitt, Chan and Kohn (2006); Masarotto and Varin (2012).

Modelos de regressão copula gaussiana têm sido empregados com sucesso em várias aplicações complexas surgindo, por exemplo, na análise de dados longitudinais, nos trabalhos de Frees and Valdez (2008); Sun, Frees and Rosenberg (2008); Shi and Frees (2011); Song, Li and Zhang (2013), genética em Li, Boehnke, Abecasis and Song (2006); He, Li, Edmondson, Raderand and Li (2012), dados mistos em Song, Li and Yuan (2009); de Leon and Wu (2011); Wu and de Leon (2014); Jiryaie, Withanage, Wu and de Leon (2016), estatísticas espaciais em Kazianka and Pilz (2010); Bai, Kang and Song (2014); Hughes (2015); Nikoloulopoulos (2016), séries temporais em Guolo and Varin (2014).

Vários autores discutiram a inferência por verossimilhança para modelos copula gaussianos, por exemplo, Masarotto and Varin 2012; Song et al. 2013; Nikoloulopoulos 2016. Enquanto os cálculos de verossimilhança para respostas contínuas são diretos, o caso discreto é consideravelmente mais difícil porque a função de verossimilhança envolve integrais gaussianas multidimensionais. Métodos de simulação são frequentemente empregados para aproximar a verossimilhança de modelos copulas gaussianos na presença de respostas discretas de alta dimensão. Pitt et ai. (2006) desenvolveram um algoritmo Monte Carlo Cadeias de Markov para inferência Bayesiana, Masarotto and Varin (2012) adotaram um algoritmo de amostragem de importância sequencial, Nikoloulopoulos (2013) estudou a máxima verossimilhança simulada com base na integração aleatória quasi-Monte Carlo, ver também Nikoloulopoulos (2016). Alternativamente, métodos de verossimilhança composta (Varin, Reid and Firth 2011) têm sido usados para reduzir a dimensionalidade integral e evitar a inversão de matrizes de covariância de alta dimensão, por exemplo, Zhao and Joe 2005; Bai et al. 2014; Hughes 2015.

Limites bem conhecidos da abordagem do copula gaussiano são a impossibilidade de lidar com a dependência assimétrica e a falta de dependência da cauda. Esses limites podem impactar o uso de copulas gaussianos para modelar formas de dependência surgidas, por exemplo, em eventos ambientais extremos ou em dados financeiros. Por outro lado, este artigo se concentra no trabalho de copulas gaussianos usadas para lidar convenientemente com a dependência na análise de regressão, conforme descrito em Masarotto and Varin (2012). Em outros termos, os parâmetros de interesse são os coeficientes de regressão, enquanto a estrutura de dependência identificada pelo copula gaussiano é um componente incômodo.

O pacote gcmr difere dos pacotes acima em termos dos modelos cobertos, o método de ajuste no caso discreto e as funcionalidades para avaliação do modelo ajustado. Os modelos de regressão marginal disponíveis incluem modelos lineares generalizados, regressão binomial negativa para contagens superdispersas e regressão beta para taxas e proporções. As matrizes de correlação de copulas gaussianas implementadas permitem lidar com várias formas de dependência que surgem, por exemplo, na análise de dados longitudinais, séries temporais e geoestatística.

Os modelos são ajustados com o método de máxima verossimilhança no caso contínuo e máxima verossimilhança simulada no caso discreto. Entre as vantagens da inferência de verossimilhança estão a possibilidade de selecionar modelos utilizando critérios de informação como AIC ou BIC e o cálculo da log-verossimilhança perfilada para inferência em um parâmetro de foco. Dadas as preocupações potenciais sobre a copula assumida, vários tipos de erros padrão sanduíche robustos estão disponíveis. Além disso, o pacote gcmr implementa a análise de resíduos para avaliar os desvios das premissas do modelo.


2. Regressão gaussiana copula


Considere um vetor de \(n\) variáveis dependentes \(Y_1,\cdots, Y_n\). A função de distribuição acumulada marginal de uma única variável \(Y_i\) é denotada por \(F(\cdot|x_i)\) e depende de um vetor \(p\)-dimensional de covariáveis \(x_i\). Assumimos que \(F(\cdot|x_i)\) é parametrizado em termos de um parâmetro de localização \(\mu_i\), normalmente correspondente ao valor esperado \(\mbox{E}(Y_i|x_i)\), que depende de \(x_i\) através da relação \[ g_1(\mu_i)=x_i^\top \beta, \]

para uma função de ligação adequada \(g_1(\cdot)\) e um vetor \(p\)-dimensional de coeficientes de regressão \(\beta\). Essa configuração abrange uma variedade de classes de modelos populares como, por exemplo, modelos lineares generalizados, apresentados em McCullagh ande Nelder (1989) ou regressão beta, propostos por Cribari-Neto and Zeileis (2010).

Se a distribuição de \(Y_i\) incluir um parâmetro de dispersão, então o modelo pode ser estendido para permitir a dispersão variável com um segundo modelo de regressão (Cribari-Neto and Zeileis, 2010) \[ g_2(\psi_i)=z_i^\top \gamma, \]

onde \(g_2(\cdot)\) é a função de ligação de dispersão, \(\psi_i\) é o parâmetro de dispersão associado a \(Y_i\), \(z_i\) é o vetor \(q\)-dimensional das covariáveis de dispersão e \(\gamma\) é o vetor correspondente dos coeficientes de regressão. Por uma questão de simplicidade de notação, a partir de então a função de distribuição univariada acumulada marginal de \(Y_i\) será denotada como \(F(\cdot|x_i)\) mesmo no caso de dispersão variável, onde de fato o modelo é \(F(\cdot|x_i, z_i)\).

Na regressão copula gaussiana a dependência entre as variáveis é modelada com uma copula gaussiana de modo que a função de distribuição acumulada de dados conjuntos é dada por \[ P(Y_1\leq y_1,\cdots, Y_n\leq y_n) = \Phi_n(\epsilon_1,\cdots,\epsilon_n; P), \]

onde \(\epsilon_i = \Phi^{-1}\big(F(y_i|x_i)\big)\), com \(\Phi(\cdot)\) denotando a função de distribuição acumulada normal padrão univariada e \(\Phi_n(\cdot;P)\) a função de distribuição acumulada normal multivariada \(n\)-dimensional com matriz de correlação \(P\).

Uma formulação equivalente do modelo de cópula gaussiana que enfatiza a configuração da regressão é descrita em Masarotto and Varin (2012). Considere um modelo de regressão que liga cada variável \(Y_i\) com um vetor de covariáveis \(x_i\) através da relação genérica \[ Y_i = h(x_i,\epsilon_i), \]

onde \(\epsilon_i\) é um erro estaocástico. Entre muitas especificações possíveis da função \(h(\cdot)\) e o erro \(\epsilon_i\), o modelo de regressão copula gaussiano assume que \[ h(x_i,\epsilon_i) = F^{-1}\big( \Phi(\epsilon_i)|x_i\big), \]

e o vetor de termos de erro \(\epsilon = (\epsilon_1,\cdots,\epsilon_n)\) tem uma distribuição normal padrão multivariada com matriz de correlação \(P\). Em outros termos, a copula gaussiana identifica um modelo de regressão construído de forma a:
(i) preservar as distribuições marginais univariadas e
(ii) têm erros normais multivariados.

Uma característica atraente da abordagem copula gaussiana é que várias formas de dependência podem ser expressas através da parametrização adequada da matriz de correlação \(P\). Por exemplo, dados longitudinais podem ser modelados com as matrizes de correlação de trabalho consideradas em equações de estimativa generalizada (Song, 2007), dependência serial em séries temporais com matriz de correlação correspondente a um processo autoregressivo e de média móvel (Guolo and Varin, 2014), dependência espacial com matriz de correlação induzida por um campo aleatório gaussiano (Bai et al., 2014).


2.1. Inferência


O pacote gcmr implementa a inferência por máxima verossimilhança para modelos de regressão copula gaussiana. Seja \(\theta\) o vetor de parâmetros do modelo constituído pelos parâmetros das marginais univariadas e os parâmetros pertencentes à matriz de correlação da copula gaussiana.

A função de verossimilhança para no caso contínuo tem a forma fechada \[ L(\theta)=\phi_n(\epsilon_1,\cdots,\epsilon_n;P)=\prod_{i=1}^n \dfrac{f(y_i|x_i)}{\phi(\epsilon_i)}, \] onde \(\phi(\cdot)\) indica a densidade normal padrão univariada, \(\phi_n(\cdot;P)\) a densidade normal padrão \(n\)-dimensional com matriz de correlação \(P\), \(f(\cdot|x_i)\) a densidade de \(Y_i\) dado \(x_i\) e a dependência das funções de densidade em \(\theta\) é mantido implícito por simplicidade de notação.

O caso discreto é consideravelmente mais complicado porque a verossimilhança é dada pela integral da normal \(n\)-dimensional \[ L(\theta)=\int_{D_1} \cdots \int_{D_n} \phi_n(\epsilon_1,\cdots,\epsilon_n;P)\mbox{d}\epsilon_1\cdots\mbox{d}\epsilon_n, \] onde o domínio integral é o produto cartesiano dos intervalos

\[ D_i=\left[ \Phi^{-1} \big( F(y_i-1|x_i)\big), \; \Phi^{-1} \big( F(y_i|x_i)\big) \right]\cdot \]

Uma quantidade notável de pesquisas tem sido dirigida à aproximação numérica de integrais normais multivariadas. Por exemplo, aproximações quase-Monte Carlo estão disponíveis através do popular pacote R mvtnorm (Genz and Bretz 2009; Genz et al. 2017). Referimo-nos a Nikoloulopoulos (2013, 2016) para estudos numéricos sobre a eficiência da estimativa de máxima verossimilhança simulada com base em mvtnorm em modelos copula gaussianos. Masarotto and Varin (2012) sugerem que aproximações numéricas eficientes de \(L(\theta)\) podem ser obtidas por generalizações adequadas de métodos numéricos para inferência de verossimilhança aproximada em modelos probit multivariados. O mais popular desses métodos é provavelmente o simulador Geweke-Hajivassiliou-Keane (GHK) (Keane 1994). O simulador GHK é um algoritmo de amostragem de importância sequencial amplamente estudado na literatura de econometria computacional, onde é comumente considerado o padrão-ouro para cálculo de verossimilhança em modelos probit multivariados, ver, por exemplo, Train (2003). O pacote gcmr implementa a estimativa de máxima verossimilhança simulada com base em uma variante do algoritmo GHK descrito em Masarotto and Varin (2012).

Existem diferentes opções para avaliar a incerteza do estimador de máxima verossimilhança. Primeiro, a abordagem clássica que avalia a incerteza com o inverso das informações de Fisher observadas. Alternativamente, a variância assintótica do estimador de máxima verossimilhança pode ser estimada com o produto externo das pontuações derivadas da decomposição preditiva da verossimilhança. As opções acima são válidas se o modelo copula gaussiano for especificado corretamente. Dadas as preocupações potenciais sobre a suposição do copula gaussiano é aconselhável comparar erros padrão baseados em modelos com estimadores sanduíche robustos. Divergências significativas entre erros padrão baseados em modelos base e robustos fornecem indicação indireta de especificação incorreta do modelo.


2.2. Resíduos


Masarotto and Varin (2012) sugeriram validar modelos de regressão copula gaussiano para respostas contínuas com resíduos quantílicos preditivos \[ r_i = \Phi^{-1}\big( F(y_i|y_{i-1},\cdots,y_{1};\widehat{\theta})\big), \]

onde \(\widehat{\theta}\) denota o estimador de máxima verossimilhança de \(\theta\). Sob condições de modelo, os resíduos quantílicos \(r_i\) são, aproximadamente, realizações de variáveis normais padrão não correlacionadas e não estão relacionadas às covariáveis \(x_i\).

Os resíduos quantil preditivos para respostas contínuas podem ser expressos na forma familiar de resíduos padronizados \[ r_i=\dfrac{\widehat{\epsilon}_i-\widehat{m}_i}{\widehat{s}_i}, \]

onde \(\widehat{\epsilon}=\Phi^{-1}\big(y_i|x_i;\widehat{\theta} \big)\), \(\widehat{m}_i=\mbox{E}\big(\epsilon_i|\epsilon_{i-1},\cdots,\epsilon_1;\widehat{\theta} \big)\) e \(\widehat{s}_i=\mbox{Var}\big( \epsilon_i|\epsilon_{i-1},\cdots,\epsilon_1;\widehat{\theta}\big)\)

No caso discreto, os resíduos quantílicos \(r_i\) são definidos como qualquer valor arbitrário no intervalo \[ \left[ \Phi^{-1}(a_i), \; \Phi^{-1}(b_i)\right], \] com \(a_i=F(y_i-1|y_{i-1},\cdots,y_1;\widehat{\theta})\) e \(b_i=F(y_i|y_{i-1},\cdots,y_1;\widehat{\theta})\).

A verificação do modelo pode ser baseada nos resíduos quantílicos aleatórios \[ r_i^\mbox{rnd}(u_i)=\Phi^{-1}\big(a_i+u_i(b_i-a_i) \big), \] onde \(u_i\) é gerado a partir de uma variável aleatória uniforme no intervalo unidade (Dunn and Smyth, 1996). Os resíduos de quantis aleatórios \(r_i^\mbox{rnd}(u_i)\) são realizações de variáveis normais padrão não correlacionadas sob condições de modelo e podem ser usados como resíduos ordinários para verificar as premissas do modelo. Como \(r^\mbox{rnd}_i(u_i)\) são randomizados, é oportuno examinar vários conjuntos de resíduos antes de tirar conclusões sobre a qualidade do ajuste do modelo.

Zucchini and MacDonald (2009) sugerem evitar a randomização com os resíduos do quantil de intervalo médio \(r^\mbox{mid}_i = \Phi^{-1}\big((a_i + b_i)/2\big)\), que são, no entanto, nem normalmente distribuídos normal nem não correlacionados.

As quantidades \(r_i\), \(r^\mbox{rnd}_i\) e \(r^\mbox{mid}_i\) são exemplos de resíduos quantílicos condicionais, pois envolvem a distribuição condicional de \(Y_i\) dadas as observações anteriores \(y_{i-1},,\cdots,y_1\). Também é possível considerar versões marginais desses resíduos com base na distribuição marginal univariada de \(Y_i\) obtido, no caso contínuo, com a equação \(r_i = \Phi^{-1}\big( F(y_i|y_{i-1},\cdots,y_{1};\widehat{\theta})\big)\) substituída pela equação \(r^\mbox{marg}_i = \Phi^{-1}\big(F(y_i|x_i;\widehat{\theta})\big)\).

Diferentemente das versões condicionais, os resíduos quantílicos marginais são úteis para verificar as suposições sobre o componente marginal do modelo, mas não são informativos sobre a correção da suposição da copula gaussiana.


3. Implementação no R


A principal função do pacote gcmr é gcmr() que permite ajustar modelos copulas gaussianos por máxima verossimilhança no caso contínuo e por máxima verossimilhança simulada no caso discreto. Os argumentos de gcmr() são os seguintes

\[ \mbox{gcmr(formula, data, subset, offset, marginal, cormat, start,} \\ \mbox{fixed, options = gcmr.options(...), model = TRUE, ...)} \]

A função possui argumentos padrão para especificação de estrutura de modelo (Chambers and Hastie 1993), como uma fórmula formula, a possibilidade de restringir a análise a um subconjunto subset dos dados data, definir um deslocamento offset ou fixar contrastes para fatores contrasts. Os argumentos específicos de gcmr() incluem os dois argumentos chave marginal e cormat, que especificam a parte marginal do modelo e a matriz de correlação copula, respectivamente. Finalmente, existem três argumentos opcionais para fornecer valores iniciais (start), fixar os valores de alguns parâmetros (fixed) e definir as opções de ajuste (options). O restante desta seção descreve os componentes de gcmr() e os métodos relacionados.


3.1. Fórmulas em duas partes


A fórmula básica permitida em gcmr é do tipo \(y\sim x1 + x2\) e especifica o modelo de regressão para a resposta média com a função de ligação \(g_1(\cdot)\) definida no argumento marginal conforme explicado na Seção 3.2 abaixo. Seguindo a implementação da regressão beta no pacote betareg (Cribari-Neto and Zeileis, 2010), o pacote gcmr também permite especificar um segundo modelo de regressão para a dispersão através de uma fórmula em duas partes do tipo \(y\sim x1 + x2 | z1 + z2\) usando funcionalidades herdadas do pacote Formula (Zeileis and Croissant 2010).

No caso de fórmula em duas partes, o modelo tem a mesma expressão de regressão média \(y\sim x1 + x2\), enquanto o parâmetro de dispersão é modelado como uma função do preditor linear \(\mbox{predictor}\sim z1 + z2\). O pacote gcmr assume um modelo log-linear \(g_2(\cdot) = \log(\cdot)\) para o modelo de regressão de dispersão.


3.2. Especificação do modelo marginal


O modelo marginal \(F(\cdot|x_i)\) é especificado através de um objeto da classe marginal.gcmr definido no argumento marginal da função gcmr(). As distribuições marginais disponíveis no gcmr versão 1.0.2 são beta, binomial, gama, Gaussiana, binomial negativa, Poisson e Weibull. Para cada uma dessas distribuições é possível escolher uma função de ligação que relaciona a média da resposta ao preditor linear como nos modelos lineares generalizados tradicionais. Todas as funções de ligação disponíveis na classe link-glm são permitidas.

Figura 1. Modelos marginais disponíveis no gcmr versão 1.0.2 com a função de ligação padrão. A coluna Dispersion (Dispersão) identifica as distribuições com um parâmetro de dispersão.


As marginais gaussianas são incluídas no gcmr para completar, mas não é recomendado usar gcmr para ajustar modelos normais multivariados que surgem trivialmente da combinação de marginais gaussianas com uma copula gaussiana. De fato, o pacote gcmr é projetado para trabalhar com modelos copulas gaussianos com distribuições marginais univariadas genéricas e, portanto, não é numericamente eficiente para inferência em modelos gaussianos lineares multivariados, onde a disponibilidade de resultados analíticos permite uma aceleração significativa dos cálculos.


3.3. Especificação da estrutura de correlação


A matriz de correlação \(P\) da copula gaussiana é especificada através de um objeto da classe cormat.gcmr definido no argumento cormat da função gcmr(). A versão 1.0.2 do gcmr inclui quatro estruturas de correlação de ampla aplicabilidade, consulte na Figura 2 abaixo. A opção de correlação de independência de trabalho é semelhante em espírito à das equações de estimação generalizada. As outras três estruturas de correlação permitem lidar com séries temporais, dados agrupados ou longitudinais e dados espaciais.

Figura 2. Modelos de correlação disponíveis no gcmr versão 1.0.2.


Dados agrupados e longitudinais podem ser analisados com cluster.cormat(id, type) construído sobre funções herdadas do pacote nlme (Pinheiro, Bates, DebRoy, Sarkar and R Core Team, 2017). As entradas de cluster.cormat() são o vetor id dos indivíduos e type, o tipo de correlação com as opções possíveis “independence”, “ar1”, “ma1”, “exchangeable” e “unstructured”.

O identificador da unidade amostral id é um vetor do mesmo comprimento que o número de observações. Supõe-se que os dados sejam classificados de forma que as observações do mesmo assunto ou cluster sejam contíguas, caso contrário, o gcmr para e retorna uma mensagem de erro. A dependência serial em séries temporais pode ser descrita com a função arma.cormat(p, q), que recebe as ordens \(p\) e \(q\) do processo \(ARMA(p, q)\) como entrada. Dados espacialmente correlacionados podem ser modelados assumindo uma função de correlação espacial Mátern definida pela função matern.cormat(D, alpha), onde \(D\) é a matriz das distâncias entre observações e alpha é o parâmetro de forma (Diggle and Ribeiro 2007). A função matern.cormat() é construída sobre a função matern() do pacote geoR (Ribeiro Jr and Diggle, 2016). O valor padrão para o parâmetro alfa é 0.5 e corresponde a um modelo de correlação exponencial.


3.4. Opções de ajuste


As opções de ajuste em gcmr() são definidas por opções de argumento ou por uma chamada direta à função

\[ \mbox{gcmr.options(seed = round(runif(1, 1, 1e+05)), nrep = c(100, 1000),} \\ \mbox{no.se = FALSE, method = c("BFGS", "Nelder-Mead", "CG"), ...} \]

Os argumentos disponíveis são seed, para fixar a semente pseudo-aleatória usada no algoritmo GHK para aproximar a função de verossimilhança com respostas discretas, nrep, para definir o número de replicações Monte Carlo no algoritmo GHK, no.se, para escolher se o cálculo é dos erros padrão ou não e method, para seleção do método de otimização a ser passado para optim().

O algoritmo de otimização padrão é o algoritmo quase-Newton BFGS. É possível fornecer um vetor de replicações Monte Carlo para nrep, de modo que o modelo seja ajustado com uma sequência de diferentes tamanhos de Monte Carlo. Neste caso, os valores iniciais para otimização da verossimilhança são retirados do ajuste anterior. Uma estratégia razoável é ajustar o modelo com um tamanho pequeno de réplicas Monte Carlo para obter valores iniciais sensatos e então reajustar com um tamanho maior de réplicas Monte Carlo. O tamanho padrão de réplicas Monte Carlo é 100 para a primeira otimização e 1.000 para a segunda otimização definitiva. Se as respostas são contínuas, então a função de verossimilhança tem uma expressão de forma fechada e os valores de seed e nrep são ignorados.


3.5. Métodos


O objeto de modelo ajustado retornado da classe gcmr é uma lista que contém, entre outras, as estimativas de máxima verossimilhança, a log-verossimilhança maximizada e as estimativas numéricas da Hessiana e do Jacobiano da verossimilhança calculadas na estimativa de máxima verossimilhança. Um conjunto de métodos padrão está disponível para extrair informações do modelo ajustado, veja a Figura 3. A maioria das funções e métodos tem sintaxe padrão como em outros pacotes R orientados para análise de regressão, veja, por exemplo, betareg (Cribari-Neto and Zeileis, 2010).

Figura 3. Funções e métodos disponíveis para objetos da classe gcmr.


O método plot() produz vários gráficos de diagnóstico do objeto gcmr ajustado que incluem gráficos de dispersão dos resíduos quantílicos em relação aos índices das observações ou em relação ao preditor linear, o gráfico de probabilidade normal com bandas de confiança com base na implementação no pacote car (Fox and Weisberg, 2011), o gráfico de dispersão dos valores previstos contra os valores observados, gráficos de autocorrelação e autocorrelação parcial dos resíduos.

O comportamento padrão do método plot() se adapta ao tipo de matriz de correlação de forma que, por exemplo, gráficos de autocorrelação são exibidos automaticamente para correlação ARMA(p, q) especificada pela função arma.cormat(p, q).

Os resíduos quantis são calculados pelo método \[ \mbox{residuals(object, type = c("conditional", "marginal"),} \\ \mbox{method = c("random", "mid"), ...)} \]

onde o argumento type permite escolher entre resíduos de quantil “conditional” ou “marginal”. O argumento method está ativo apenas no caso discreto para selecionar entre resíduos de quantil “random” ou resíduos de quantil de intervalo médio “mid”.

A log-verossimilhança perfilada pode ser obtida com uma chamada ao método \[ \mbox{profile(fitted, which, low, up, npoints = 10, display = TRUE,} \\ \mbox{alpha = 0.05, progress.bar = TRUE, ...)} \] onde argumento which é o índice do parâmetro a ser perfilado, low e up são os limites inferior e superior usados no cálculo, npoints é o número de pontos usados no cálculo da verossimilhança perfilada, alpha é o nível de significância, display controla se a verossimilhança perfilada deve ser apresentada ou não e progress.bar define uma barra de progresso para visualizar a progressão do cálculo da verossimilhança perfilada. Se os valores dos limites low e up não são fornecidos, então eles são definidos iguais ao parâmetro estimado menos e mais três vezes o erro padrão, respectivamente.


4. Aplicações


Exemplo 1


Como exemplo consideramos os resultados de um estudo inédito de Alzheimer relatados em Hand and Talor (1987). Existem dois grupos de pacientes, um recebendo placebo e outro recebendo uma droga chamada lecitina. Cada paciente foi solicitado a memorizar uma lista de palavras. O resultado foi o número de palavras lembradas corretamente. Este teste foi administrado em cada um dos cinco momentos: 0, 1, 2, 4 e 6 meses. Há 48 pacientes, 22 em drogas e 26 em placebo. As palavras usadas em cada ocasião de teste foram exatamente as mesmas das ocasiões de teste anteriores.

Espera-se que quaisquer efeitos do tratamento possam ser pequenos e lentos para aparecer devido à natureza das demências. A administração do teste de memorização em várias ocasiões também introduz um efeito de prática, isto é, um sujeito pode melhorar com o tempo apenas por causa dos efeitos de transferência de uma ocasião para outra. Isso pode ser devido à exposição repetida às mesmas palavras e/ou exposição repetida à situação de teste. 


alzheim = read.table(file = "http://leg.ufpr.br/~lucambio/CM/alzheim.txt", header = T)
attach(alzheim)
head(alzheim)
##   group time0 time1 time2 time4 time6
## 1     0    20    19    20    20    18
## 2     0    14    15    16     9     6
## 3     0     7     5     8     8     5
## 4     0     6    10     9    10    10
## 5     0     9     7     9     4     6
## 6     0     9    10     9    11    11


As questões de interesse são:
- Existe uma tendência temporal na capacidade de recordação nas cinco ocasiões de teste? Por exemplo, está diminuindo?
- Alguma tendência temporal difere por tratamento?

O conjunto de dados contém o tratamento na coluna 1: 0 para placebo e 1 para lecitina e os resultados do teste nas colunas 2 a 6. Não há outras covariáveis disponíveis.

Fonte de dados: Hand and Taylor (1987). Multivariate analysis of variance and repeated measures. London: Chapman and Hall, Inc.


Definindo os tratamentos:

Tratamento = factor(group, levels=c(0,1), labels=c('Placebo','Lecithin'))
Tratamento = rep(Tratamento, each=5)


Resposta

N.Palabras = matrix(cbind(time0,time1,time2,time4,time6),ncol=48,byrow=T)
N.Palabras = c(N.Palabras)


Controle do tempo de observacao

Meses = factor(rep(c(0,1,2,4,6),48))


Identificador da unidade amostral

ID = rep(1:48, each=5)


Estudo descritivo

library(lattice)
xyplot(N.Palabras~Meses|Tratamento, groups=ID, type='l',
       ylab='No. de palabras corretamente lembradas')


Estrutura de correlacao

corr = cor(cbind(time0,time1,time2,time4,time6))
colnames(corr)=c('Inicio','Mes 1','Mes 2','Mes 4','Mes 6')
rownames(corr)=c('Inicio','Mes 1','Mes 2','Mes 4','Mes 6')
library(graphics)
mosaicplot(corr,main='Correlacoes')


library("ellipse")
plotcorr(corr,main='Correlacoes',col=colors())


Indicam estes graficos que a estrutura mais adequada deve ser constante.
Modelo proposto: Poisson

library(gcmr)
set.seed(123)
mod00.exc = gcmr( N.Palabras ~ Meses + Tratamento,
                marginal = poisson.marg, cormat = cluster.cormat(ID, type = "exc"))
summary(mod00.exc)
## 
## Call:
## gcmr(formula = N.Palabras ~ Meses + Tratamento, marginal = poisson.marg, 
##     cormat = cluster.cormat(ID, type = "exc"))
## 
## 
## Coefficients marginal model:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         2.23315    0.05808  38.451   <2e-16 ***
## Meses1              0.05236    0.04810   1.088    0.276    
## Meses2              0.03929    0.04819   0.815    0.415    
## Meses4              0.01472    0.04854   0.303    0.762    
## Meses6             -0.01445    0.04890  -0.296    0.768    
## TratamentoLecithin -0.04381    0.07265  -0.603    0.546    
## 
## Coefficients Gaussian copula:
##     Estimate Std. Error z value Pr(>|z|)    
## tau  0.49352    0.03168   15.58   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## log likelihood = 673.02,  AIC = 1360


library(ggplot2)
residuos = residuals(mod00.exc)
ggplot(data=as.data.frame(qqnorm( residuos, plot=F)), mapping=aes(x=x, y=y)) + 
  geom_point() + geom_smooth(method="lm", se=FALSE)

Exemplo 2


Como outro exemplo, escolhemos o conjunto de dados worldindices que pode ser lido como indicado: 

exemplo = read.table(file = "http://leg.ufpr.br/~lucambio/MSM/worldindices.txt", header = TRUE)
head(exemplo)
##   ind     X.GSPC      X.N225     X.SSEC     X.GDAXI      X.FCHI      X.FTSE
## 1   1 0.02731297 0.901417531 0.37751835 0.154884969 0.228478606 0.054874507
## 2   2 0.42685972 0.009378426 0.05154599 0.147087136 0.255125709 0.281661523
## 3   3 0.07092167 0.326867186 0.80295656 0.094014688 0.301736635 0.223585323
## 4   4 0.05002729 0.007313291 0.06850590 0.020725559 0.019621662 0.113758341
## 5   5 0.02219325 0.461568811 0.96730505 0.006743126 0.007200171 0.002163702
## 6   6 0.42156859 0.013288659 0.32775430 0.156650766 0.162692380 0.166405342

Este conjunto de dados contém resíduos padronizados transformados de retornos diários dos principais índices de ações mundiais em 2009 e 2010 (396 observações). Os índices considerados são as principais bolsas de valores das seis maiores economias do mundo: o norte-americano S&P 500 (^GSPC) X.GSPC, o japonês Nikkei 225 (^N225) X.N225, o chinês SSE Composite Index (^SSEC) X.SSEC, o alemão DAX (^GDAXI) X.DAXI, o francês CAC 40 (^FCHI) X.FCHI e o britânico FTSE 100 Index (^FTSE) X.FTSE. 


library(VineCopula)
up <- function(x, y) {
  # upper panel: empirical contour plot
  op <- par(usr = c(-3, 3, -3, 3), new = TRUE, pch = 19, cex = 0.6)
  BiCopKDE(x, y,
           levels = c(0.01, 0.05, 0.1, 0.15, 0.2),
           margins = "exp",
           axes = FALSE)
  on.exit(par(op))
}

lp <- function(x, y) {
  # lower panel: scatter plot (copula data) and correlation
  op <- par(usr = c(0, 1, 0, 1), new = TRUE, pch = 19, cex = 0.6)
  points(x, y, pch = 1, col = "black")
  r <- cor(x, y, method = "spearman") # Spearman's rho
  txt <- format(x = r, digits = 3, nsmall = 3)[1]
  text(x = 0.5, y = 0.5, labels = txt, cex = 1 + abs(r) * 2, col = "blue")
  on.exit(par(op))
}

dp <- function(x) {
  # diagonal panel: histograms (copula data)
  op <- par(usr = c(0, 1, 0, 1.5), new = TRUE, pch = 19, cex = 0.6)
  hist(x, freq = FALSE, add = TRUE, col = "brown", border = "black", main = "")
  abline(h = 1, col = "black", lty = 2)
  on.exit(par(op))
}

pairs(exemplo[,2:7], lower.panel = lp, upper.panel = up, diag.panel = dp, gap = 0.5)