• Introdução
• Definição e construção
• Estimação
• Resultados teóricos
• Gráficos
• Complemento
• Exemplo
• Considerações finais

Neste seminário vamos apresentar o modelo aditivo de Aalen, introduzido em 1980. O presente trabalho está dividido em 5 partes. Primeiramente veremos os conceitos básicos e a construção do modelo, logo depois resultados teóricos, seguiremos com os métodos avaliação e diagnóstico do ajuste, uma seção com alguns exemplos e, por fim, as considerações finais.

O modelo de riscos proporcionais proposto por Cox em 1972 possui vantagens na interpretabilidade dos resultados e é relativamente fácil de construir. Entretanto ele apresenta as seguintes desvantagens:

• A suposição de riscos proporcionais é consideravelmente forte e depende das covariáveis incluídas no modelo.

• Mudanças ao longo do tempo na influência das covariáveis não são facilmente descobertas.

• Os parâmetros \(\beta\) são fixos, implicando que os efeitos são fixos e constantes.

Como forma de contornar estes problemas, pode-se utilizar o Modelo Aditivo de Aalen.

• Modelo não-paramétrico de Aalen é um modelo alternativo robusto ao de Cox.
• É um modelo de risco aditivo para análise de regressão de dados censurados.
• Modelo de regressão que generaliza o estimador de Nelson-Aalen.
• Possibilita o monitoramento do efeito da covariável ao longo do acompanhamento.
• É permitido valores estimados negativos para a função de riscos.

O modelo de Aalen é completamente não-paramétrico no sentido de que as funções são ajustadas e não parâmetros.

A suposição básica do modelo não-paramétrico de Aalen (1980):

\[\lambda(t) = \beta_0(t) + \sum_{j=1}^{p}\beta_j(t)x_j(t)=\beta(t)X(t)\]

i.e, a função de risco é uma combinação linear dos parâmetros e as covariáveis, ambas possivelmente variando no tempo.

A estimação direta de \(\beta_j(t)\) não é viável (como de costume). Sendo assim, a estimação da função acumulada de \(\beta_j(t)\) é realizada (novamente, como de costume).

Aalen propôs o seguinte estimador:

\[\hat{\textbf{B}}(t) = \sum_{t_i \leq t} \textbf{Z}(t_i)\textbf{I}(t_i)\] onde \(\textbf{B}(t)=\left(B_1(t),...,B_p(t)\right)'\) com \(\text{B}_j(t)=\int_0^t\beta_j(s)ds\).

Embora \(\textbf{Z}(t_i)\) possa ser qualquer inversa generalizada, iremos utilizar a forma dos mínimos quadrados:

\[ \textbf{Z}(t_i)= \Big(\textbf{X}'(t_i)\textbf{X}(t_i)\Big)^{-1}\textbf{X}'(t_i) \]

Dado \(\beta\) e considerando todas as covariáveis fixadas no tempo zero podemos estimar a função de risco acumulador por:

\[ \hat{\Lambda}(t|x) = x'\hat{\textbf{B}}(t) = \hat{B_0}(t) + \sum_{j=1}^{p}\hat{B_j}(t)x_j(t)\] Agora, existem duas propostas para estimar a função de sobrevivência:

\[ \tilde{S}(t|x) = exp\{-\tilde{\Lambda}(t|x)\} \] e

\[ \hat{S}(t|x) = \prod_{i: t_i\leq t} \Big[1-\Big(\textbf{Z}(t_i)I(t_i)\Big)'x)\Big]\] Onde este último é baseado no estimador de Kaplan-Meier.

O modelo de Aalen pode ser visto como a expansão de primeira ordem da série de Taylor do modelo de Cox em torno de 0.

Fixando \(t\) e considerando o modelo de Cox para uma covariável:

\[ \begin{align*} \lambda(t, x) &= \lambda(t, 0) + \lambda'(t,0)(x)\\ \lambda(t, x) &= \lambda_0(t)\exp\{0\beta\} +\lambda_0(t)\beta exp\{0\beta\}x\\ \lambda(t, x) &= \lambda_0(t) +\lambda_0(t)x\beta \\ \lambda(t, x) &= \lambda_0(t) +x(t)\beta(t) \\ \lambda(t, x) &= \beta_0(t) + x(t)\beta(t) \end{align*}\]

Ainda ganhamos uma interpretação já conhecida: o termo \(\lambda_0(t)\) é “incorporado” em \(\beta_0(t)\). Deixaremos a generalização a cargo do leitor.

Os resíduos de Cox-Snell podem ser calculados igualmente, ou seja:

\[ \hat{e_i} = \hat{\Lambda}(t_i|x_i, \beta(t_i)) \]

Novamente, espera-se que os resíduos tenham uma distribuição exponencial de parâmetro 1.

Considerando a hipótese de interesse com \(t\in [0,\tau]\):

\[ H_{0j} : \beta_j(t)=0 \]

Podemos usar a seguinte estatística:

\[ U = \sum_{t_i\leq\tau} \textbf{K}(t_i)\textbf{Z}(t_i)\textbf{I}(t_i)\]

com \(K(t)=\Big\{diag\Big[\Big(X(t)'X(t)\Big)^{-1}\Big]\Big\}^{-1}\).

Sendo \(V\) um estimador da matriz de covariância de \(U\), a estatística de teste é então \(U_jV_{jj}^{−1/2}\) que sob \(H_0\) tem uma distribuição normal padrão.

• A importância de uma covariável pode mudar durante o período de acompanhamento.

• As funções de regressão, que estimam a contribuição das covariáveis para a função de risco em cada tempo de falha, podem ser descritas no tempo.

\[\hat{\Lambda}(t)=\hat{B}(t) \ \ vs \ \ t\]

• Inclinações positivas: aumentos nos valores das covariáveis estão associados com aumento na função de risco.

• Inclinações negativas: aumentos nos valores das covariáveis estão associados com decréscimos na função de risco.

• Inclinações aproximadamente iguais a zero: covariáveis não influenciam no risco.

Embora não tenha sido discutido até o presente momento, existem técnicas que permitem quantificar a diferença dos resíduos observados para a real distribuição esperada.

Pode-se utilizar os resultados de Kolmogorov-Smirnov em que a estatística do teste é da forma:

\[ D = \sup_{r_C} |F_n(r_C) - F_0(r_C)|\] Ansin (2015) discorre sobre a avaliação dos resíduos do Cox-Snell utilizando este método. Este teste já está implementado na biblioteca timereg.

A interpretabilidade é simples, porém não é direta. Utilizando as estimações das funções relacionadas e os \(\beta\)’s é possível calcular as quantidades de interesse em qualquer ponto no tempo, sem complicações.

Entretanto duas observações devem ser feitas:

• A modelagem proposta permite que os riscos sejam negativos.

• A função de sobrevivência estimada não é necessariamente monótona.

Vamos construir um exemplo simples e didático. Considere o conjunto de dados:

As matrizes \(\textbf{X}(t_i)\) são dadas por:

\[\textbf{X}(0) =\begin{bmatrix} 1& 1 \\ 1& 0 \\ 1& 1 \\ 1& 1 \end{bmatrix}; \textbf{X}(2^+) =\begin{bmatrix} 0& 0 \\ 1& 0 \\ 1& 1 \\ 1& 1 \end{bmatrix}; \textbf{X}(5^+) =\begin{bmatrix} 0& 0 \\ 0& 0 \\ 1& 1 \\ 1& 1 \end{bmatrix}; \textbf{X}(8^+) =\begin{bmatrix} 0& 0 \\ 0& 0 \\ 0& 0 \\ 1& 1 \end{bmatrix}; \textbf{X}(9^+) =\begin{bmatrix} 0& 0 \\ 0& 0 \\ 0& 0 \\ 0& 0 \end{bmatrix} \]

Neste caso os vetores \(\textbf{I}(t_i)\) são:

\[\begin{align*} \textbf{I}(0) &= \begin{pmatrix} 0& 0& 0& 0 \end{pmatrix} \\ \textbf{I}(2) &= \begin{pmatrix} 1& 0& 0& 0\end{pmatrix} \\ \textbf{I}(5) &= \begin{pmatrix} 0& 1& 0& 0 \end{pmatrix} \\ \textbf{I}(8) &= \begin{pmatrix} 0& 0& 1& 0 \end{pmatrix} \\ \textbf{I}(9) &= \begin{pmatrix} 0& 0& 0& 0 \end{pmatrix} \end{align*} \]

sur<-Surv(Tempo, Cens)
fit1<-aareg(sur~Trat)

summary(fit1)

##             slope   coef se(coef)      z     p
## Intercept  0.1720  0.500    0.500  1.000 0.317
## Trat      -0.0805 -0.294    0.503 -0.585 0.558
## 
## Chisq=0.34 on 1 df, p=0.558; test weights=aalen

fit1$coef

##   Intercept       Trat
## 2         0  0.3333333
## 5         1 -1.0000000

Podemos verificar abrindo as contas manualmente. O estimador para \(\beta(2)\) é dado por:

\[ (\textbf{X}'(2)\textbf{X}(2))^{-1}\textbf{X}'(2)\textbf{I}(2)) =\begin{pmatrix}\begin{pmatrix} 1& 1 \\ 1& 0 \\ 1& 1 \\ 1& 1 \end{pmatrix}'\begin{pmatrix} 1& 1 \\ 1& 0 \\ 1& 1 \\ 1& 1 \end{pmatrix}\end{pmatrix}^{-1}\begin{pmatrix} 1& 1 \\ 1& 0 \\ 1& 1 \\ 1& 1 \end{pmatrix}' \begin{pmatrix}1& 0& 0& 0\end{pmatrix} = \begin{pmatrix} 0 \\ \frac{1}{3} \end{pmatrix}\] De fato:

fit1$coef

##   Intercept       Trat
## 2         0  0.3333333
## 5         1 -1.0000000

Sendo assim, o risco instantâneo de um indivíduo \(i\) falhar no tempo \(t=2\) pode ser calculador por:

\[ \lambda(2) = 0 + 0.334+x_i(2)\beta_1(2)\] Considerando \(x_i(2) = 1\), a sobrevivência do indivíduo é:

\[ S_i(2) = exp\{-0.334\} \approx 0.716 \]

Vamos reproduzir o exemplo dos dados aids.txt, encontrado no livro texto.

A classificação do paciente (soropositivo assintomático, ARC e AIDS) pode mudar ao longo do estudo. Ou seja, alguns pacientes que iniciaram o estudo com a classificação soropositivo assintomático evoluíram para AIDS no final do estudo passando por ARC.

sr<-Surv(data$ti,data$tf,data$cens)
m1.s<-aareg(sr~id +sex+ factor(grp),data=data);summary(m1.s)

##                  slope      coef se(coef)      z        p
## Intercept     2.01e-03  0.024200 0.009180  2.630 0.008490
## id           -5.34e-05 -0.000606 0.000265 -2.290 0.022100
## sex          -2.84e-04 -0.000826 0.005000 -0.165 0.869000
## factor(grp)2 -3.85e-04 -0.003100 0.003890 -0.797 0.426000
## factor(grp)3  3.62e-03  0.029000 0.012800  2.260 0.024000
## factor(grp)4  3.09e-03  0.034500 0.010200  3.370 0.000753
## 
## Chisq=18.48 on 5 df, p=0.0024; test weights=aalen

m2<-aareg(sr~id + factor(grp),data=data);summary(m2)

##                  slope      coef se(coef)      z        p
## Intercept     2.12e-03  0.026700 0.009330  2.870 0.004170
## id           -5.84e-05 -0.000714 0.000279 -2.560 0.010600
## factor(grp)2 -3.75e-04 -0.002660 0.003920 -0.678 0.498000
## factor(grp)3  3.75e-03  0.029700 0.012800  2.320 0.020400
## factor(grp)4  3.45e-03  0.038200 0.010900  3.510 0.000451
## 
## Chisq=19.05 on 4 df, p=0.000769; test weights=aalen

McKeague e Sasieni (1994) propõem um modelo aditivo semiparamétrico. Segundo os autores, embora o modelo aditivo de Aalen tenha uma maior flexibilidade, o número de covariáveis que o modelo pode ser limitado. Sendo assim, eles propõem um modelo semiparamétrico onde somente algumas covariáveis podem variar não parametricamente ao longo do tempo e o restante se mantém fixa.

Este é o modelo está implementado no r no pacote timereg e pela função aalen.

Outra extensão interessante é o modelo aditivo multiplicativo de Cox-Aalen, proposto em 2002 por Zhang e Scheike. Este modelo assume a forma:

\[ \lambda_i(t) = Y_i(t)\big(X_i(t)'\alpha(t)\big) exp\{Z_i(t)' \beta\}\]

m1<-aalen(sr~id +sex+ factor(grp),data=data)
(m1$pval.testBeq0)

##  (Intercept)           id          sex factor(grp)2 factor(grp)3 factor(grp)4 
##        0.076        0.085        0.329        0.704        0.124        0.010

m1$pval.testBeqC

##  (Intercept)           id          sex factor(grp)2 factor(grp)3 factor(grp)4 
##        0.701        0.376        0.396        0.184        0.538        0.350

m2.s<-aalen(sr~id + factor(grp) ,data=data)
pt<-capture.output(m2.s)
print(unname(pt[1:19]))

##  [1] "Additive Aalen Model "                                                 
##  [2] ""                                                                      
##  [3] "Test for nonparametric terms "                                         
##  [4] ""                                                                      
##  [5] "Test for non-significant effects "                                     
##  [6] "             Supremum-test of significance p-value H_0: B(t)=0"        
##  [7] "(Intercept)                           2.60               0.069"        
##  [8] "id                                    2.38               0.113"        
##  [9] "factor(grp)2                          1.24               0.738"        
## [10] "factor(grp)3                          2.38               0.119"        
## [11] "factor(grp)4                          3.22               0.010"        
## [12] ""                                                                      
## [13] "Test for time invariant effects "                                      
## [14] "                   Kolmogorov-Smirnov test p-value H_0:constant effect"
## [15] "(Intercept)                        0.18400                       0.657"
## [16] "id                                 0.00732                       0.486"
## [17] "factor(grp)2                       0.12800                       0.184"
## [18] "factor(grp)3                       0.26100                       0.328"
## [19] "factor(grp)4                       0.32700                       0.403"

m2.s<-aalen(sr~id + factor(grp) ,data=data)
pt<-capture.output(m2.s)
print(unname(pt[20:30]))

##  [1] "                     Cramer von Mises test p-value H_0:constant effect"
##  [2] "(Intercept)                        4.62000                       0.599"
##  [3] "id                                 0.00937                       0.351"
##  [4] "factor(grp)2                       2.26000                       0.237"
##  [5] "factor(grp)3                       4.20000                       0.613"
##  [6] "factor(grp)4                      16.10000                       0.357"
##  [7] ""                                                                      
##  [8] "   "                                                                   
##  [9] "   "                                                                   
## [10] "  Call: "                                                              
## [11] "aalen(formula = sr ~ id + factor(grp), data = data)"

O modelo aditivo de Aalen possui vantagens interessantes, porém a interpretabilidade é muito menos direta. Embora, assim como no livro, chegamos à mesma conclusão que o modelo de Cox com variáveis dependentes do tempo, pouco ganhamos com o novo modelo.

Se por um lado a suposição de taxas proporcionais é forte, a suposição dos componentes aditivos também é. Ainda por cima, temos as desvantagens numéricas das taxas aditivas.

AALEN, O.O. A Model for Nonparametric Regression Analysis of Counting Processes. Klonecki, N., Koesh, A. and Rosinski, J., Eds., Lecture Notes in Statistics, 2th Edition, Springer, New York, 1-25. http://dx.doi.org/10.1007/978-1-4615-7397-5_1

ANSIN. E. An evaluation of the Cox-Snell residuals. Uppsala Universitet, Dissertação de mestrado. 2015.

COLOSIMO, E. A. GIOLO, S. R. Análise de Sobrevivências aplicada. São Paulo: Blucher, 2006.

MCKEAGUE, I. W. SASIENI, P. D. A partly parametric additive risk model, Biometrika, Volume 81, Issue 3, Setembro 1994, pg. 501–514, https://doi.org/10.1093/biomet/81.3.501.

SCHEIKE, T. H. ZHANG, M. An Additive-Multiplicative Cox-Aalen Regression Model, Scandinavian Journal of Statistics, Vol. 29, No. 1, pp. 75-88, 2002, https://www.jstor.org/stable/4616700.