Em muitas áreas de aplicação, quando se trata de modelos de regressão para uma variável de resposta dicotômica, a regressão logística é uma das abordagens de modelagem mais adotadas. A literatura sobre o assunto é imensa e cresce rapidamente; ver, por exemplo, Hosmer et al. (2013).

Um dos problemas que podem ocorrer quando utilizada a regressão logística é o chamado problema de separação ou de verossimilhança monotônica, estudado por Albert & Anderson (1984).

Este problema ocorre no processo de ajuste de um modelo logístico quando a função de verossimilhança converge enquanto pelo menos uma estimativa do coeficiente de regressão diverge para \(-\infty\) ou \(+\infty\), de modo que o estimador de máxima verossimilhança não existe. Tal situação ocorre se os subconjuntos dos dados com variável de resposta igual a 1 e 0 puderem ser perfeitamente separados por uma única covariável ou por uma intrincada combinação linear de covariáveis. Por exemplo, se uma covariável representa um fator de risco com dois níveis, ausência e presença, digamos e para todos os sujeitos com o fator de risco ausente, a resposta é igual a 1, o problema acontece.

Exemplo

Primeiro, carregamos o conjunto de dados de câncer endometrial (Heinze & Schemper, 2002): 

library(hmclearn)
library(MCMCpack)
data(Endometrial)
head(Endometrial)
##   NV PI   EH HG
## 1  0 13 1.64  0
## 2  0 16 2.26  0
## 3  0  8 3.14  0
## 4  0 34 2.68  0
## 5  0 20 1.28  0
## 6  0  5 2.31  0
attach(Endometrial)

Este arquivo de dados contém 79 linhas e 4 variáveis:
NV Indicador de fator de risco de neovascularização (0=Ausente, 1=Presente)
PI Índice de pulsatilidade da artéria uterina
EH Altura do endométrio
HG histologia do paciente (0=Baixo, 1=Alto)

Estudo descritivo

lattice::xyplot(HG ~ PI | NV)

lattice::xyplot(HG ~ EH | NV)

table(HG,NV)
##    NV
## HG   0  1
##   0 49  0
##   1 17 13

Metodologia

Hamiltonian Monte Carlo (HMC) ou Monte Carlo Hamiltoniano é uma ferramenta poderosa para computação Bayesiana. Em comparação com o algoritmo tradicional de Metropolis-Hastings, o HMC oferece maior eficiência computacional, principalmente em situações de modelagem mais dimensionais ou mais complexas.

Cadeia Markov Monte Carlo: O Básico

MCMC é uma ampla classe de ferramentas computacionais para aproximação integral e geração posterior de amostras. Na análise Bayesiana, os algoritmos MCMC são usados principalmente para simular amostras para aproximação da distribuição posterior.

Na análise Bayesiana, a estimação e a inferência dos parâmetros de interesse são feitas com base nos dados observados \(\mathcal{D}\) juntamente com a informação a priori que se tem sobre os parâmetros de interesse \(\theta= (\theta_1,\cdots,\theta_k)^\top \in \mathbb{R}^k\). A distribuição posterior ou a posteriori \(f(\theta |\mathcal{D})\) combina os dados e as informações anteriores de acordo com a fórmula de Bayes e é proporcional ao produto da função de verossimilhança \(f(\mathcal{D}|\theta)\) e a densidade a priori \(\pi(\theta)\) (Carlin & Luís 2008), \[\begin{equation*} f(\theta |\mathcal{D}) = \dfrac{f(\mathcal{D}|\theta)\pi(\theta)}{\displaystyle \int f(\mathcal{D}|\theta)\pi(\theta)\mbox{d}\theta} \propto f(\mathcal{D}|\theta)\pi(\theta)\cdot \end{equation*}\]

A integral no denominador é geralmente difícil de calcular. Mas como o denominador é constante em relação a \(\theta\), pode-se trabalhar com a posteriori não normalizada \(f(\mathcal{D}|\theta)\pi(\theta)\). Na ausência de uma expressão explícita da distribuição posterior, aproximá-la com amostras simuladas seguindo \(f(\theta |\mathcal{D})\) torna-se uma alternativa desejável.

Metropolis-Hastings

O algoritmo Metropolis foi o primeiro método MCMC amplamente utilizado para gerar amostras da Cadeia de Markov seguindo \(f(\theta |\mathcal{D})\). O método originou-se de uma aplicação da física na década de 1950 (Metropolis et al. 1953), e foi ampliado cerca de duas décadas depois por Hastings (1970), dando origem ao nome de algoritmo Metropolis-Hastings (MH). Começamos com uma breve descrição do MH, pois o HMC foi construído sobre um conceito semelhante.

MH gera uma seqüência de valores que formam uma cadeia de Markov, cujos valores podem ser usados para aproximar uma densidade posterior \(f(\theta |\mathcal{D})\) ou a posteriori. Por brevidade, retiramos \(\mathcal{D})\) da expressão e escrevemos a posterior simplesmente como \(f(\theta)\). Os valores na cadeia de Markov \(\theta^{(t)}\) são indexados por \(t = 0,1,\cdots,N\), onde \(\theta^{(0)}\) é um valor inicial especificado pelo usuário ou pelo programa.

MH define uma probabilidade de transição que assegura que a cadeia de Markov é ergódica e satisfaz equilíbrio detalhado e reversibilidade (Chib & Greenberg 1995). Essas condições técnicas são postas em prática para garantir que as amostras da cadeia tenham o suporte de \(\theta\) sem viés.

Em MH, os valores de \(\theta^{(t)}\) na cadeia são definidos em parte por uma densidade de proposta \(\theta^{\mbox{PROP}}\), onde \(\mbox{PROP}\) é uma proposta para o próximo valor na cadeia. Esta densidade proposta está condicionada ao valor anterior \(\theta^{(t-1)}\). Uma variedade de funções de proposta pode ser usada, sendo as propostas de passeio aleatório a escolha mais comum.

Em MH, uma proposta é aceita com probabilidade \[\begin{equation} \alpha = \mbox{min}\left(1, \dfrac{f\big(\theta^\mbox{PROP}\big)\, q\big(\theta^{t-1}|\theta^\mbox{PROP}\big)}{f\big(\theta^{t-1}\big)\, q\big(\theta^\mbox{PROP}|\theta^{t-1}\big)} \right), \end{equation}\] Quando \(q\) é simétrico, ou seja, \(q\big(\theta^{t-1}|\theta^\mbox{PROP}\big)=q\big(\theta^\mbox{PROP}|\theta^{t-1}\big)\), isso simplifica para \[\begin{equation} \alpha = \mbox{min}\left(1, \dfrac{f\big(\theta^\mbox{PROP}\big)}{f\big(\theta^{t-1}\big)} \right), \end{equation}\] que é usado no algoritmo Metropolis original.

O denominador na a posteriori é constante em relação a \(\theta\). Assim, a razão de densidades a posteriori em dois pontos diferentes \(\theta^\mbox{PROP}\) e \(\theta^{(t−1)}\) pode ser comparada mesmo quando o denominador é desconhecido, com os denominadores sendo cancelados. Como uma derivação da distribuição a posteriori completa, numerador e denominador, não é necessária para implementar MH e HMC, como veremos, os analistas de dados têm considerável flexibilidade para selecionar priores de sua preferência.

A taxa de aceitação \(\alpha\) é uma medida importante da eficiência de um algoritmo MH. Um exame cuidadoso dos papéis de \(\alpha\)’s fornece uma compreensão mais intuitiva do algoritmo:

  • 1- Quando \(f\big(\theta^{\mbox{PROP}}\big)\geq f\big(\theta^{(t-1)}\big)\), a proposta \(\theta^{\mbox{PROP}}\) representa um valor mais provável que o valor anterior \(\theta^{(t−1)}\), quantificado pelas funções de densidade. Quando isso ocorre, a proposta é sempre aceita, ou seja, com probabilidade 1.
  • Quando \(f\big(\theta^{\mbox{PROP}}\big)< f\big(\theta^{(t-1)}\big)\), a proposta \(\theta^{\mbox{PROP}}\) tem uma densidade menor em relação ao valor anterior e aceitamos a proposta aleatoriamente com probabilidade \(\alpha\in (0,1)\), o que indica a probabilidade relativa de observar \(\theta^{\mbox{PROP}}\) de \(f\), em comparação com \(\theta^{(t−1)}\). Quanto maior o \(\alpha\), maior a chance de aceitar o \(\theta^{\mbox{PROP}}\). Se a proposta não for aceita, a proposta será descartada e a cadeia permanecerá no lugar \(\theta^{(t)}:=\theta^{(t−1)}\) e iniciaremos com uma nova proposta.

Com tal esquema, o algoritmo freqüenta regiões de maior densidade a posteriori, enquanto ocasionalmente visita as áreas de baixa densidade, por exemplo, caudas em situações unidimensionais. Desde que o algoritmo execute um número suficiente de iterações, a distribuição empírica das amostras da cadeia MCMC deve se aproximar da densidade posterior verdadeira. Os valores simulados podem, portanto, ser usados para estimativa e inferência com base na distribuição a posteriori. Ver Carlin & Louis (2008) Chib & Greenberg (1995) Gelman et al. (2013) para detalhes adicionais sobre MH.

Limitações do algoritmo de Metropolis-Hasting

Os requisitos teóricos para o uso de MH são bastante mínimos, tornando-se uma escolha atraente para inferência Bayesiana. As limitações do MH são principalmente computacionais. Com propostas geradas aleatoriamente, geralmente é necessário um grande número de iterações para entrar em áreas de maior densidade posterior. Mesmo algoritmos de MH eficientes às vezes aceitam menos de 25% das propostas (Roberts et al. 1997). Em situações dimensionais mais baixas, o aumento do poder computacional pode compensar a menor eficiência até certo ponto. Mas em situações de modelagem mais dimensionais e mais complexas, computadores maiores e mais rápidos sozinhos raramente são suficientes para superar o desafio.

O amostrador de Gibbs pode ser uma alternativa viável e mais eficiente ao MH em algumas situações (Geman & Geman 1987). De fato, várias plataformas de software populares, A exigência de Gibbs para densidades a posteriori condicionais expressas explicitamente, no entanto, impediu que ela fosse usada em muitas situações práticas. Além dessa restrição, Gibbs também tem suas próprias limitações de eficiência (Robert 2001). É neste contexto que o HMC surge como uma alternativa preferencial para a análise Bayesiana.

Monte Carlo Hamiltoniano

O Monte Carlo Hamiltoniano melhora a eficiência do MH empregando um esquema de geração de proposta guiada. Mais especificamente, o HMC usa o gradiente do log a posteriori para direcionar a cadeia de Markov para regiões de maior densidade posterior, onde a maioria das amostras é retirada. Como resultado, uma cadeia HMC bem ajustada aceitará propostas a uma taxa muito maior do que o algoritmo MH tradicional (Roberts et al. 1997). Este algoritmo foi proposto e eimplementado por Samuel Thomas & Wanzhu Tu (2021).

A ideia

Os métodos utilizados para gerar propostas influenciam fortemente a eficiência do MCMC. Suponha que \(f(\theta)\) seja uma função de densidade a posteriori unidimensional e \(-\log f(\theta)\) assuma a forma de uma curva inversa em forma de sino. Para gerar em uma região de alta densidade a posteriori, é necessário amostrar na região correspondente aos menores valores de \(-\log f(\theta)\); a região pode ser alcançada com a orientação do gradiente de \(-\log f(\theta)\). Em certo sentido, a abordagem é análoga ao movimento de um objeto hipotético em uma curva sem atrito, onde o objeto atravessa e permanece no fundo do vale enquanto ocasionalmente visita os terrenos mais altos em ambos os lados. Na mecânica clássica, tais movimentos são descritos pelas equações Hamiltonianas, onde as trocas de energia cinética e potencial ditam a localização do objeto em um determinado momento.

Em um sistema Hamiltoniano, as posições horizontal e vertical são dadas por \((\theta,\mbox{p})\). No MCMC, estamos interessados em \(\theta\). O parâmetro \(\mbox{p}\), que é frequentemente chamado de momentum, é uma quantidade auxiliar que usamos para simular \(\theta\) sob as equações Hamiltonianas.

As equações Hamiltonianas

Introduzimos HMC em uma configuração genérica de MCMC, onde \(\theta\) segue a densidade a posteriori \(f(\theta)\) de interesse e o momento \(\mbox{p}\) é gerado a partir de uma distribuição paramétrica. O momentum corresponde à dimensionalidade de \(\theta\) como um vetor de comprimento \(k\).

Escrevemos a função Hamiltoniana como \(H(\theta,\mbox{p})\), que consiste em energia potencial \(U(\theta)\) e energia cinética \[\begin{equation*} K(\mbox{p}): H(\theta,\mbox{p}) = U(\theta)+K(\mbox{p}), \end{equation*}\] onde \(\mbox{p}\) e \(\theta\in\mathbb{R}^k\).

Em aplicações estatísticas de MCMC, estamos principalmente interessados em gerar \(\theta\) a partir de uma dada distribuição \(f(\theta)\). Para isso, fazemos \(U(\theta) := −\log f(\theta)\). Tal designação garantiria que \(\theta\) gerada pela função Hamiltoniana seguisse a distribuição desejada. Para momentum, normalmente assumimos \(\mbox{p}\sim N_k(0,M)\), onde \(M\) é uma matriz de covariância especificada pelo usuário.

Sob esta formulação, temos \[\begin{equation*} H(\theta,\mbox{p})=- \log f(\theta)+\frac{1}{2}\mbox{p}^\top M^{-1}\mbox{p}\cdot \end{equation*}\]

Ao longo do tempo, HMC viaja em trajetórias que são governadas pelas seguintes equações diferenciais de primeira ordem, conhecidas como equações Hamiltonianas \[\begin{equation*} \begin{array}{rcl} \dfrac{\mbox{d}\mbox{p}}{\mbox{d}t} & = & -\dfrac{\partial H(\theta,\mbox{p})}{\partial\theta} \, = \, -\dfrac{\partial U(\theta)}{\partial \theta} \, = \, \nabla_{\theta} \log f(\theta),\\ \dfrac{\mbox{d}\theta}{\mbox{d}t} & = & -\dfrac{\partial H(\theta,\mbox{p})}{\partial\mbox{p}} \, = \, -\dfrac{\partial K(\mbox{p})}{\partial \mbox{p}} \, = \, M^{-1}\mbox{p}, \end{array} \end{equation*}\] onde \(\nabla_{\theta} \log f(\theta)\) é o gradiente do logaritmo da densidade a posteriori. Uma solução para as equações Hamiltonianas é uma função que define o caminho de \((\theta,\mbox{p}\) a partir do qual valores específicos de podem ser amostrados. Dentro de uma iteração MCMC, amostramos um valor de \(\theta\) desse caminho. A aleatoriedade das amostras MCMC vem do momentum \(\mbox{p}\sim N_k(0,M)\) e do valor específico \(\theta\) que escolhemos.

Exemplo do algoritmo Monte Carlo Hamiltoniano

Para resposta binária, queremos estimar \[\begin{equation*} p \, = \, P(Y=1|X) \, = \, \dfrac{1}{1+e^{−X\beta}}. \end{equation*}\]

As observações do vetor de respostas é \(y=(y_1\cdots,y_n)^\top\). Os valores de covariáveis para o sujeito \(i\) são \(x_i^\top=(x_{i0},...,x_{iq})\) para \(q\) covariáveis mais um intercepto. Escrevemos a matriz de projeto completa como \(X=(x_1^\top,\cdots,x_n^\top)^\top \in \mathbb{R}^{n\times (q+1)}\) para \(n\) observações. Os coeficientes de regressão são um vetor de comprimento \(q+1\), \(\beta=(\beta_0,\cdots,\beta_q)^\top\).

Preparação dos dados

Endometrial$PI2 <- with(Endometrial, (PI - mean(PI)) / sd(PI))
Endometrial$EH2 <- with(Endometrial, (EH - mean(EH)) / sd(EH))
Endometrial$NV2 <- Endometrial$NV - 0.5
X <- cbind(1, as.matrix(Endometrial[, which(colnames(Endometrial) %in% c("PI2", "EH2", "NV2"))]))
y <- Endometrial$HG
colnames(X) <- c("(Intercept)", "PI2", "EH2", "NV2")

A preparação anterior dados têm dois motivos: um deles consiste em melhorar a interpretabilidade dos dados e diminuição de possíveis dificuldades para encontrar as estimativas dos parâmetros devido à valores extremos das covariáveis.

O segundo motivo refere-se à necessidade de construir a matriz de projeto \(X\) para a estimação utilizando o algoritmo de Metropolis-Hasting.

Modelo frequentista

ajuste = glm(y ~ X-1, family = binomial(link = "logit"))
(summary(ajuste))
## 
## Call:
## glm(formula = y ~ X - 1, family = binomial(link = "logit"))
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.50137  -0.64108  -0.29432   0.00016   2.72777  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## X(Intercept)    7.8411   857.8755   0.009 0.992707    
## XPI2           -0.4217     0.4432  -0.952 0.341333    
## XEH2           -1.9219     0.5599  -3.433 0.000597 ***
## XNV2           18.1856  1715.7509   0.011 0.991543    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 109.517  on 79  degrees of freedom
## Residual deviance:  55.393  on 75  degrees of freedom
## AIC: 63.393
## 
## Number of Fisher Scoring iterations: 17

Para comparar os resultados, primeiro ajustamos um modelo linear padrão usando a função frequentista glm. Observe as estimativas de erro padrão alto para o \((Intercept)\) (intrecepto) e \(NV2\). Isso ocorre porque todos os 13 indivíduos com neovasculação presente também apresentam \(y_i=1\) (HG=1).

Como resultado dos desafios com o ajuste frequentista, Agresti (2015) na seção 10.3.2 também ajusta este modelo de regressão logística usando um algoritmo Metropolis de passeio aleatório do pacote (Martin, Quinn, Park 2011). Ajustamos o modelo com o mesmo número de iterações MCMC, distribuições anteriores e hiperparâmetros como no texto.

Este modelo também atribui uma priori normal para β onde σ2β=100. utiliza uma parametrização de precisão, onde B0=1/σ2β. O hiperparâmetro médio é 0 conforme especificado em b0.

Observe o alto número de iterações de Metropolis especificadas e o tempo necessário para ajustar, mesmo com código compilado otimizado. Nossa abreviatura MH refere-se ao algoritmo geral de Metropolis-Hastings.

O usuário deve definir fornecer a matriz de design diretamente para uso no hmclearn. Nosso primeiro passo é carregar os dados e armazenar a matriz de projeto X e o vetor de variável dependente y.

library(hmclearn)
N <- 1e4
set.seed(412)
t1.hmc <- Sys.time()
f.hmc <- hmc(N = N, theta.init = rep(1, 4),
            epsilon = 0.1, L = 20,
            logPOSTERIOR = logistic_posterior,
            glogPOSTERIOR = g_logistic_posterior,
            varnames = colnames(X),
            param=list(y = y, X=X, sig2beta=100),
            parallel=FALSE, chains=2)
t2.hmc <- Sys.time()
t2.hmc - t1.hmc
## Time difference of 41.18386 secs
summMH = summary(f.hmc, burnin=3000, probs=c(0.025, 0.5, 0.975))
## Summary of MCMC simulation
summMH
##                  2.5%        50%      97.5%      rhat
## (Intercept) -0.292832  2.9281285  9.2283365 1.0017419
## PI2         -1.444876 -0.4606057  0.3592366 1.0000986
## EH2         -3.423364 -2.1016785 -1.0784341 0.9999496
## NV2          2.192343  8.5423329 21.0647490 1.0016373

Histogramas da distribuição a posteriori mostram que as estimativas dos parâmetros HMC estão relativamente de acordo com a distribuição eperada, ainda observamos que a mediana da distribuição representaria as estimativas adequadas nesta situação dos parâmetros da regressão. 

diagplots(f.hmc, burnin=3000, comparison.theta=summMH[,2], plotfun = 2)
## $histogram

Vejamos, de maneira aproximada, a significância dos coeficientes da regressão. 

est = summMH[,2]
se = est
cf = rep(0,7000); for(i in 1:7000) cf[i] = f.hmc$theta[[1]][[i+3000]][1]; se[1] = sd(cf) 
cf = rep(0,7000); for(i in 1:7000) cf[i] = f.hmc$theta[[1]][[i+3000]][2]; se[2] = sd(cf) 
cf = rep(0,7000); for(i in 1:7000) cf[i] = f.hmc$theta[[1]][[i+3000]][3]; se[3] = sd(cf) 
cf = rep(0,7000); for(i in 1:7000) cf[i] = f.hmc$theta[[1]][[i+3000]][4]; se[4] = sd(cf) 
zval = est/se
coefficients = cbind(Estimate = est, `Std. Error` = se, `Z value` = zval, 
                          `Pr(>|z|)` = 2*pnorm(abs(zval), lower.tail = FALSE))
coefficients
##               Estimate Std. Error    Z value    Pr(>|z|)
## (Intercept)  2.9281285  2.4049814  1.2175265 0.223403974
## PI2         -0.4606057  0.4648396 -0.9908917 0.321738442
## EH2         -2.1016785  0.6136942 -3.4246345 0.000615627
## NV2          8.5423329  4.7970888  1.7807327 0.074956130

Como resultado percebemos que a variável \(PI2\) não é signiicativa para o modelo. 

X1 <- cbind(1, as.matrix(Endometrial[, which(colnames(Endometrial) %in% c("EH2", "NV2"))]))
colnames(X1) <- c("(Intercept)", "EH2", "NV2")
ajuste1 = glm(y ~ X1-1, family = binomial(link = "logit"))
(summary(ajuste1))
## 
## Call:
## glm(formula = y ~ X1 - 1, family = binomial(link = "logit"))
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.76820  -0.59872  -0.30827   0.00016   2.81512  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## X1(Intercept)    7.7802   869.4128   0.009 0.992860    
## X1EH2           -1.8468     0.5447  -3.391 0.000697 ***
## X1NV2           17.9894  1738.8255   0.010 0.991745    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 109.517  on 79  degrees of freedom
## Residual deviance:  56.378  on 76  degrees of freedom
## AIC: 62.378
## 
## Number of Fisher Scoring iterations: 17
t1.hmc <- Sys.time()
f.hmc1 <- hmc(N = N, theta.init = rep(1, 3),
            epsilon = 0.1, L = 20,
            logPOSTERIOR = logistic_posterior,
            glogPOSTERIOR = g_logistic_posterior,
            varnames = colnames(X1),
            param=list(y = y, X=X1, sig2beta=100),
            parallel=FALSE, chains=2)
t2.hmc <- Sys.time()
t2.hmc - t1.hmc
## Time difference of 41.3393 secs
summMH1 = summary(f.hmc1, burnin=3000, probs=c(0.025, 0.5, 0.975))
## Summary of MCMC simulation
summMH1
##                   2.5%       50%     97.5%      rhat
## (Intercept) -0.3858465  2.621532  9.640298 0.9999287
## EH2         -3.1754762 -1.968774 -1.041331 1.0001007
## NV2          1.8849861  7.798243 21.777223 0.9999328
diagplots(f.hmc1, burnin=3000, comparison.theta=summMH1[,2], plotfun = 2)
## $histogram

est = summMH1[,2]
se = est
cf = rep(0,7000); for(i in 1:7000) cf[i] = f.hmc$theta[[1]][[i+3000]][1]; se[1] = sd(cf) 
cf = rep(0,7000); for(i in 1:7000) cf[i] = f.hmc$theta[[1]][[i+3000]][2]; se[2] = sd(cf) 
cf = rep(0,7000); for(i in 1:7000) cf[i] = f.hmc$theta[[1]][[i+3000]][3]; se[3] = sd(cf) 
zval = est/se
coefficients = cbind(Estimate = est, `Std. Error` = se, `Z value` = zval, 
                          `Pr(>|z|)` = 2*pnorm(abs(zval), lower.tail = FALSE))
coefficients
##              Estimate Std. Error   Z value     Pr(>|z|)
## (Intercept)  2.621532  2.4049814  1.090042 2.756945e-01
## EH2         -1.968774  0.4648396 -4.235385 2.281607e-05
## NV2          7.798243  0.6136942 12.707049 5.403564e-37

Modelo estimado

Mostramos nos gráficos a seguir as probabilidades estimadas de \(P(HG=1|X)\) e \(P(HG=0|X)\). Resulta que temos duas vriáveis explicativas 

attach(Endometrial)
x.EH = seq(min(EH2),max(EH2), by=0.01)
# Estimativas P(HG=1|X)
estimativas1 = 1/(1+exp(-cbind(1,x.EH,-0.5)%*%summMH1[,2]))
estimativas2 = 1/(1+exp(-cbind(1,x.EH,0.5)%*%summMH1[,2]))
#
par(mfrow=c(1,2), mar = c(4,5,1,1), pch=19)
plot(x.EH,estimativas1, main="NV=-0.5", type="l", col="red", xlab="EH2", ylab="", ylim=c(0,1))
lines(x.EH,1-estimativas1, col="blue")
points(EH2[NV2==-0.5], HG[NV2==-0.5])
grid()
plot(x.EH,estimativas2, main="NV=0.5", type="l", col="red", xlab="EH2", ylab="", ylim=c(0,1))
lines(x.EH,1-estimativas1, col="blue")
points(EH2[NV2==0.5], HG[NV2==0.5])
grid()
legend("bottom",legend = c("P(HG=1|X)","P(HG=0|X)"), lwd = 2, lty = 1, col = c("red","blue"))

Bibliografia

  1. Albert, A. & Anderson, J. A. (1984). On the existence of maximum likelihood estimates in logistic regression models. Biometrika, 71, 1-10.
  2. Betancourt, M. (2017), ‘A Conceptual Introduction to Hamiltonian Monte Carlo’, arXiv preprint arXiv:1701.02434 pp. 1-60. URL: http://arxiv.org/abs/1701.02434.
  3. Carlin, B.P. & Louis, T.A. (2008), Bayesian Methods for Data Analysis, CRC Press. Google-Books-ID: GTJUt8fcFx8C.
  4. Carpenter, B., Hoffman, M. D., Brubaker, M., Lee, D., Li, P. & Betancourt, M. (2015), ‘The Stan Math Library: Reverse-Mode Automatic Differentiation in C++’, arXiv preprint arXiv:1509.07164.
  5. Chib, S. & Greenberg, E. (1995), ‘Understanding the Metropolis-Hastings Algorithm’, The American Statistician 49(4), 327-335.
  6. Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A. & Rubin, D.B. (2013), Bayesian Data Analysis, CRC Press.
  7. Geman, S. & Geman, D. (1987), Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images, in ‘Readings in Computer Vision’, Elsevier, pp. 564–584.
  8. Hastings, W.K. (1970), ‘Monte Carlo Sampling Methods Using Markov Chains and Their Applications’, Biometrika 57(1), 97–109. URL: https://doi.org/10.1093/biomet/57.1.97.
  9. Heinze, G. & Schemper, M. (2002). A solution to the problem of separation in logistic regression. Statistics in medicine, 21(16), 2409-2419.
  10. Hosmer, D.W., Lemeshow, S. & Sturdivant, R.X. (2013). Applied Logistic Regression. Wiley, Hoboken, NJ, third edition.
  11. Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H. & Teller, E. (1953), ‘Equation of State Calculations by Fast Computing Machines’, The Journal of Chemical Physics 21(6), 1087–1092.
  12. Neal, R.M., Brooks, S., Gelman, A. & Jones, G. (2011), Handbook of Markov Chain Monte Carlo, CRC Press.
  13. Robert, C. (2001), The Bayesian Choice: From Decision-Theoretic Foundations to Computational Implementation, Springer Science & Business Media.
  14. Roberts, G.O., Gelman, A., Gilks, W.R. et al. (1997), ‘Weak convergence and optimal scaling of random walk metropolis algorithms’, The annals of applied probability 7(1), 110- 120.
  15. Ruth, R.D. (1983), ‘A canonical integration technique’, IEEE Trans. Nucl. Sci. 30 (CERNLEP- TH-83-14), 2669–2671.
  16. Thomas, S. & Tu, W. (2021) Learning Hamiltonian Monte Carlo in R, The American Statistician, 75:4, 403-413, DOI: 10.1080/00031305.2020.1865198.