Resumo Técnico: Estimativa de Parâmetros para Distribuições de Tipo-Fase Inhomogêneas com Escala Temporal a partir de Observações Discretas

1. Definição do Problema

As distribuições de Fase-Tipo (PH) clássicas modelam o tempo até a absorção de um Processo de Salto de Markov (MJP) homogêneo de estados finitos. As distribuições de Fase-Tipo Inhomogêneas (IPH) estendem este arcabouço ao permitir que as intensidades de transição variem ao longo do tempo, capturando, assim, comportamentos de cauda pesada, taxas de risco não constantes e dinâmicas dependentes do tempo. Uma subclasse específica e tratável de distribuições IPH assume uma matriz de sub-intensidade escalonada no tempo da forma $\Lambda(t) = h_\beta(t)\Lambda$, onde $\Lambda$ é uma matriz de base constante e $h_\beta(t)$ é uma função de escala não negativa parametrizada por $\beta$.

Embora a inferência estatística para esses modelos tenha sido desenvolvida para casos em que os tempos de absorção ou trajetórias completas são observados exatamente, muitas aplicações práticas envolvem observações discretas e irregularmente espaçadas. Nesses cenários (por exemplo, dados de painel multiestatais em monitoramento médico), os tempos exatos de transição e de absorção não são observados, tornando o problema um desafio de dados faltantes. A estimativa direta de máxima verossimilhança é difícil porque a verossimilhança dos dados observados envolve exponenciais de matrizes complexas e requer otimização não linear restrita sobre as entradas de $\Lambda$ e $\beta$. Este artigo aborda a estimativa de parâmetros $(\pi, \Lambda, \beta)$ para distribuições IPH com escala temporal quando o processo subjacente é observado apenas em pontos de tempo discretos.

2. Metodologia

Os autores propõem um arcabouço de inferência estatística que combina aumento de dados de ponte de Markov (Markov-bridge data augmentation) com um algoritmo de Expectativa-Maximização Estocástica (SEM). A estratégia central consiste em transformar o problema inhomogêneo em um problema homogêneo via uma transformação de tempo, simulando trajetórias contínuas latentes e atualizando os parâmetros com base na verossimilhança de dados completos resultante.

2.1 Transformação de Tempo

Para a subclasse $\Lambda(t) = h_\beta(t)\Lambda$, o processo inhomogêneo $X$ pode ser mapeado para um MJP homogêneo $Y$ via a transformação de tempo $s = g_\beta^{-1}(t) = \int_0^t h_\beta(u) du$. Sob esta transformação, o tempo de absorção $\tau$ do processo inhomogêneo satisfaz $\tau \stackrel{d}{=} g_\beta(\rho)$, onde $\rho$ segue uma distribuição PH homogênea padrão com parâmetros $(\pi, \Lambda)$.

2.2 O Algoritmo SEM

O procedimento de estimativa alterna entre a simulação de dados faltantes e a maximização da verossimilhança de dados completos:

1. Inicialização:

   • A distribuição inicial $\pi$ é estimada via frequências empíricas dos estados iniciais.
   • Uma estimativa inicial para $\beta$ ($\hat{\beta}_0$) é fornecida.
   • Uma matriz de sub-intensidade inicial $\hat{\Lambda}_0$ é obtida tratando as observações discretas como se fossem observações contínuas de um processo homogêneo (ignorando inicialmente a escala temporal).
   • Um refinamento inicial de $\hat{\beta}$ é realizado usando ascensão de gradiente na log-verossimilhança de dados completos derivada de caminhos simulados (sem transformação de tempo).

2. Procedimento Iterativo (Passo SE e Passo M):

   • Passo SE (E-estocástica): Usando a estimativa atual de $\hat{\beta}$, os pontos de tempo discretos observados são transformados para a escala de tempo homogênea. Trajetórias contínuas latentes entre os pontos de observação são simuladas usando pontes de Markov condicionadas aos estados observados. Para caminhos ainda não absorvidos, a trajetória é estendida até o estado absorvente. Esses caminhos simulados constituem uma realização de "dados completos".
   • Passo M (Maximização):
     • Atualização de $\Lambda$: Dados os caminhos homogêneos simulados, a matriz de sub-intensidade de base $\Lambda$ é atualizada usando o Estimador de Máxima Verossimilhança (MLE) de forma fechada para MJPs homogêneos (contagem de transições e tempos de ocupação).
     • Atualização de $\beta$: Os caminhos simulados são mapeados de volta para a escala de tempo inhomogênea. O parâmetro $\beta$ é refinado maximizando a log-verossimilhança de dados completos em relação a $\beta$ usando Ascensão de Gradiente. Este passo envolve avaliar o gradiente da log-verossimilhança, que depende dos tempos de salto e tempos de ocupação escalonados.

3. Terminação: O algoritmo termina quando a ascensão de gradiente interna para $\beta$ converge após uma única atualização, indicando que as estimativas se estabilizaram.

2.3 Identificabilidade

O artigo estabelece condições para a identificabilidade dos parâmetros. Especificamente, a parametrização $(\Lambda, \beta)$ é identificável se, e somente se, as funções de escala $h_{\beta_1}(t)$ e $h_{\beta_2}(t)$ diferirem por mais de uma constante multiplicativa sempre que $\beta_1 \neq \beta_2$. Esta condição garante que a matriz de base $\Lambda$ e a função de escala temporal não sejam confundidas.

3. Principais Resultados

O desempenho do método proposto foi avaliado através de estudos de simulação e uma aplicação de dados reais.

3.1 Estudos de Simulação

• Modelo Matriz-Gompertz: Usando $h_\beta(t) = e^{\beta t}$, os autores simularam 1.000 caminhos com 3 estados transientes. O método recuperou com sucesso os parâmetros verdadeiros $(\pi, \Lambda, \beta)$ através de vários horizontes de observação ($T$). A bondade de ajuste, medida pelo teste Kolmogorov-Smirnov (KS) nos tempos de absorção, permaneceu alta para $T$ maiores, mas degradou significativamente quando $T$ era muito curto para capturar eventos de absorção suficientes (ex: $T=36$ resultou em um p-valor baixo).
• Modelo Matriz-Weibull: Usando $h_\beta(t) = \beta t^{\beta-1}$, o método foi testado com 1.000 caminhos com 2 estados transientes. O algoritmo convergiu para estimativas precisas (ex: $\beta$ verdadeiro $= 3$, estimado $\hat{\beta} \approx 2,97$) após 23 iterações. O teste KS resultou em um p-valor de 0,5004, confirmando a capacidade do modelo de reproduzir a distribuição de tempo de absorção verdadeira.

3.2 Aplicação em Dados Reais

O método foi aplicado a um conjunto de dados de 622 receptores de transplante de coração monitorando a progressão da Vasculopatia de Aloenxerto Coronário (CAV). Os dados consistiam em observações irregularmente espaçadas dos estados dos pacientes (livre de CAV, leve, moderada/grave) e morte (estado absorvente).

• Ajuste do Modelo: Um modelo IPH com escala temporal de Gompertz ($h_\beta(t) = e^{\beta t}$) foi ajustado. Os parâmetros estimados indicaram que as intensidades de transição aumentam exponencialmente ao longo do tempo.
• Comparação: O modelo de escala temporal foi comparado a um modelo PH homogêneo. O modelo IPH proporcionou um ajuste significativamente melhor (p-valor KS = 0,5966) em comparação ao modelo homogêneo (p-valor KS = 0,01066), demonstrando a necessidade de considerar intensidades variantes no tempo neste contexto.
• Interpretação: Os resultados destacaram que pacientes no estado moderado/grave apresentavam o maior risco instantâneo de morte, e o modelo capturou a dinâmica de risco acelerado associada à progressão da doença.

4. Significância e Alegações

O artigo afirma fornecer uma ferramenta flexível e computacionalmente tratável para ajustar modelos IPH de escala temporal a dados multiestatais irregulares. As principais contribuições incluem:

• Tratamento de Dados Faltantes: O arcabouço aborda efetivamente o desafio de tempos de transição e absorção não observados, tratando-os como dados faltantes e reconstruindo-os via pontes de Markov.
• Evitação de Otimização Restrita: Ao utilizar atualizações de dados completos, o método evita a difícil otimização não linear restrita da verossimilhança de dados observados. A atualização para $\Lambda$ é de forma fechada, e a atualização para $\beta$ é uma ascensão de gradiente unidimensional.
• Fundamentação Teórica: Os autores estabelecem as condições de identificabilidade para a subclasse de escala temporal, descartando a confusão de parâmetros entre a matriz de base e a função de escala.
• Utilidade Prática: Através da aplicação em CAV, o artigo demonstra que o método de inferência baseado em SEM proposto produz estimativas de parâmetros confiáveis e captura com precisão a evolução temporal nas intensidades de transição, superando os contrapartes homogêneos em cenários do mundo real com observações irregulares.

Os autores enfatizam que os estimadores relatados são estimadores de máxima verossimilhança avaliados em realizações simuladas de dados completos. Sob as condições padrão de SEM, a sequência de estimadores forma uma cadeia de Markov ergódica que se concentra próximo ao maximizador de dados observados, mantendo a interpretação de estimadores de máxima verossimilhança até o nível da variabilidade de Monte Carlo.