acath
Fernando Lucambio
Última atualização 2022-11-30
Dados
library(Hmisc); library(ggpubr)
Hmisc::getHdata(acath)
head(acath)## sex age cad.dur choleste sigdz tvdlm
## 1 0 73 132 268 1 1
## 2 0 68 85 120 1 1
## 3 0 54 45 NA 1 0
## 4 1 58 86 245 0 0
## 5 1 56 7 269 0 0
## 6 0 64 0 NA 1 0Este conjunto de dados e da Duke University Cardiovascular Disease Databank e consiste em 3504 pacientes e seis variáveis. Os pacientes foram encaminhados para a Duke University Medical Center com dor no peito.
Algumas análises interessantes incluem prever a probabilidade de doenca coronaria, ou seja, >= 75% de estreitamento do diametro em pelo menos uma arteria coronária importante e prever a probabilidade de doença coronária grave, uma vez que é uma doença significativa. A primeira análise usaria sigdz como variável de resposta e a segunda análise usaria tvdlm no subconjunto de pacientes com sigdz = 1.
Doença coronária grave é definida como três vasos doentes ou doença principal na esquerda e é denotada por tvdlm = 1. Temos como covariável o sexo do paciente, denotado como sex = 0 para o sexo masculino e sex = 1 para o sexo feminino.
Referência: http://biostat.mc.vanderbilt.edu/wiki/Main/DataSets
Estudo descritivo
Sexo = factor(acath$sex, levels = c(0,1), labels = c("Masculino","Feminino"))
Idade = acath$age
dados = data.frame(Sexo,Idade)Desta maneira transformamos a variável sex em um fator com dois níveis nominais e na nova variável Idade registramos a idade (age) do paciente.
library(ggplot2)
ggplot(dados, aes(Sexo)) + geom_bar() + ylab("Contagem")
Observamos que a maioria dos pacientes pertencem ao sexo masculino. Temos ainda em cad.dur a duração dos sintomas de doença arterial coronáriana e em choleste o nível de colesterol no sangue.
Prevendo a probabilidade de doença coronária utilizando como resposta sigdz
Resposta = factor(acath$sigdz, levels=c(0,1), labels=c('Não','Sim'))
Colesterol = acath$choleste
Cad.dur = acath$cad.dur
dados1 = data.frame(Resposta, Colesterol, Cad.dur, Sexo, Idade)par(mfrow=c(1,1), mar=c(3,3,1,1), cex=1.0)
mosaicplot(table(Resposta,Sexo), xlab='Doença coronária', main='')
ggstripchart(dados1, y = "Colesterol", x = "Sexo",
color = "Resposta", palette = c("#00AFBB", "#FC4E07"))
ggboxplot(dados1, y = "Colesterol", x = "Sexo",
color = "Resposta", palette = c("#00AFBB", "#FC4E07"))
ggstripchart(dados1, x = "Idade", y = "Colesterol", font.tickslab = c(8, "bold", "black"),
color = "Resposta", palette = c("#00AFBB", "#FC4E07"))
ggstripchart(dados1, x = "Resposta", y = "Cad.dur", color = "Sexo", palette = c("#00AFBB", "#FC4E07"))
Observando as demais variáveis explicativas.
plot(Colesterol, Resposta, pch=19, ylab='Doença coronária', xlab='Colesterol')
plot(Idade, Resposta, pch=19, ylab='Doença coronária', xlab='Idade')
plot(Cad.dur, Resposta, pch=19, ylab='Doença coronária', xlab='Duração dos sintomas')
Percebemos que a doença coronária se manifesta mais no sexo masculino e está clara a influença da idade, mas não está clara a influença do nível de colesterol e nem da duração dos sintomas.
Modelos
modelo00 = glm(Resposta ~ Idade + Sexo + Cad.dur + Colesterol, family=binomial(link='logit'), data = dados1)
modelo01 = update(modelo00, family=binomial(link='probit'))
modelo02 = update(modelo00, family=binomial(link='cloglog'))
modelo03 = update(modelo00, family=binomial(link='cauchit'))Escolhendo a função de ligação
AIC(modelo00,modelo01,modelo02,modelo03)## df AIC
## modelo00 5 2354.626
## modelo01 5 2355.864
## modelo02 5 2353.651
## modelo03 5 2360.281Selecionamos a função de ligação complementar log-log, menor valor de AIC, como a mais adequada.
modelo04 = step(modelo02)## Start: AIC=2353.65
## Resposta ~ Idade + Sexo + Cad.dur + Colesterol
##
## Df Deviance AIC
## - Cad.dur 1 2345.4 2353.4
## <none> 2343.7 2353.7
## - Colesterol 1 2425.2 2433.2
## - Idade 1 2504.3 2512.3
## - Sexo 1 2750.3 2758.3
##
## Step: AIC=2353.43
## Resposta ~ Idade + Sexo + Colesterol
##
## Df Deviance AIC
## <none> 2345.4 2353.4
## - Colesterol 1 2425.5 2431.5
## - Idade 1 2509.8 2515.8
## - Sexo 1 2751.5 2757.5De maneira automática não descartamos nenhuma variável explicativa.
summary(modelo04)##
## Call:
## glm(formula = Resposta ~ Idade + Sexo + Colesterol, family = binomial(link = "cloglog"),
## data = dados1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8908 -0.8981 0.5030 0.7978 2.0801
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.804673 0.228822 -12.257 <2e-16 ***
## Idade 0.040300 0.003337 12.078 <2e-16 ***
## SexoFeminino -1.323691 0.073485 -18.013 <2e-16 ***
## Colesterol 0.005414 0.000601 9.007 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2895.3 on 2257 degrees of freedom
## Residual deviance: 2345.4 on 2254 degrees of freedom
## (1246 observations deleted due to missingness)
## AIC: 2353.4
##
## Number of Fisher Scoring iterations: 5Procurando alternativas.
par(mfrow=c(1,1),mar=c(5,4,1,1),cex=1.0)
Transformada = I(log(Colesterol)/I(Idade^2))
xyplot(Transformada ~ Resposta | Sexo, data = dados1)
modelo05 = glm(Resposta ~ Sexo + Transformada, family=binomial(link='cloglog'))
summary(modelo05)##
## Call:
## glm(formula = Resposta ~ Sexo + Transformada, family = binomial(link = "cloglog"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2174 -0.9963 0.5538 0.7497 2.4785
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.23237 0.08009 15.39 <2e-16 ***
## SexoFeminino -1.19605 0.07000 -17.09 <2e-16 ***
## Transformada -336.01229 31.63797 -10.62 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2895.3 on 2257 degrees of freedom
## Residual deviance: 2443.8 on 2255 degrees of freedom
## (1246 observations deleted due to missingness)
## AIC: 2449.8
##
## Number of Fisher Scoring iterations: 5Residuos
residuos = resid(modelo05, type='pearson')
par(mfrow=c(2,2),mar=c(5,4,1,1),cex=0.8)
plot(residuos, pch=19, cex=0.6)
grid()
hist(residuos)
grid()
box()
qqnorm(residuos, pch=19, cex=0.6)
grid()
plot(residuos~fitted(modelo05), pch=19, cex=0.8)
grid()
Resposta
par(mfrow=c(1,1),mar=c(5,4,1,1),cex=1.0)
plot(modelo05$y~fitted(modelo05), pch=19,cex=0.8, xlab='Valores preditos',
ylab='Valores observados')
grid()
Observamos que no gráfico de valores estimados vs observados há muitos erros nas estimativas da resposta.
PERGUNTA: para quais valores de estimativa da P(Y = 1) deve-se assumir que o valor estimado da resposta é Y = 1 e para quais Y = 0?
mu.preditos=predict(modelo05,type='response')
Y.preditos1=ifelse(mu.preditos>0.5,1,0)
table(modelo05$y,Y.preditos1)## Y.preditos1
## 0 1
## 0 452 316
## 1 268 1222Y.preditos2=ifelse(mu.preditos>0.7,1,0)
table(modelo05$y,Y.preditos2)## Y.preditos2
## 0 1
## 0 534 234
## 1 453 1037Y.preditos3=ifelse(mu.preditos>0.4,1,0)
table(modelo05$y,Y.preditos3)## Y.preditos3
## 0 1
## 0 227 541
## 1 106 1384Y.preditos4=ifelse(mu.preditos>0.9,1,0)
table(modelo05$y,Y.preditos4)## Y.preditos4
## 0 1
## 0 765 3
## 1 1440 50Para decidirmos adequadamente o valor de corte para considerarmos como reposta Y = 1 ou Y = o utilizamos a chamada Curva ROC.
library(pROC)
plot.roc(modelo05$y,fitted(modelo05))
grid()
Y.preditos5=ifelse(mu.preditos>0.7,1,0)
table(modelo05$y,Y.preditos5)## Y.preditos5
## 0 1
## 0 534 234
## 1 453 1037Sensitividade (verdadeiros positivos) 70%
Specificidade (verdadeiros negativos) = 70%