MLGExemplo04
Fernando L Pérez
09/06/2021
Dados
Este conjunto de dados é 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 doença coronária (>= 75% de estreitamento do diámetro em pelo menos uma artéria coronária importante) e prever a probabilidade de doençca coronária grave, uma vez que é uma doença significativa.
A primeira análise usaria sigdz como variável de resposta e o segundo 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.
Leitura dos dados.
Hmisc::getHdata(acath)
attach(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 0Sexo é denotado como: sex = 0 para o sexo masculino e 1 para o sexo feminino.
Estudo descritivo
Variáveis:
Sexo = factor(sex,labels=c('Masculino','Feminino'))
Idade = agecad.dur = Duração dos sintomas de doença arterial coronariana.
Resposta=factor(sigdz,levels=c(0,1),labels=c('Não','Sim'))
Colesterol = choleste
Cad.dur = cad.dur
rm(acath)
dados0 = data.frame(Resposta,Sexo,Idade,Colesterol,Cad.dur)
dados1 = na.exclude(dados0)
attach(dados1)
rm(dados0)Primeira análise: prevendo a probabilidade de doença coronária.
Resposta: sigdz
Sexo = factor(sex, levels=c(0,1), labels=c('Masculino','Feminino'))
Idade = age
Resposta = factor(sigdz, labels=c('Sim','Não'))
Colesterol = choleste
Cad.dur = cad.dur
dados0 = data.frame(Resposta,Sexo,Idade,Colesterol,Cad.dur)
dados1 = na.exclude(dados0)
attach(dados1)
par(mfrow=c(1,1), mar=c(3,3,1,1), cex=1.0)
mosaicplot(table(dados1$Resposta,dados1$Sexo), xlab='Doença coronária', main='')
ggpubr::ggstripchart(dados1, y = "Colesterol", x = "Sexo",
color = "Resposta", palette = c("#00AFBB", "#FC4E07"))
ggpubr::ggboxplot(dados1, y = "Colesterol", x = "Sexo",
color = "Resposta", palette = c("#00AFBB", "#FC4E07"))
ggpubr::ggstripchart(dados1, x = "Idade", y = "Colesterol", font.tickslab = c(8, "bold", "black"),
color = "Resposta", palette = c("#00AFBB", "#FC4E07"))## Don't know how to automatically pick scale for object of type labelled/integer. Defaulting to continuous.
ggpubr::ggstripchart(dados1, x = "Resposta", y = "Cad.dur", color = "Sexo", palette = c("#00AFBB", "#FC4E07"))## Don't know how to automatically pick scale for object of type labelled/integer. Defaulting to continuous.
plot(density(dados1$Colesterol), pch=19, ylab='Doença coronária', xlab='Colesterol', main="")
grid()
rug(dados1$Colesterol)
plot(dados1$Idade, dados1$Resposta, pch=19, ylab='Doença coronária', xlab='Idade')
plot(dados1$Cad.dur,dados1$Resposta,pch=19, ylab='Doença coronária', xlab='Duração dos sintomas')
plot(density(dados1$Idade[dados1$Resposta=="Não"]), xlab="Idade", ylab="Densidade estimada", main="Idade quando a resposta é negativa")
grid()
rug(dados1$Idade[dados1$Resposta=="Não"])
plot(density(dados1$Colesterol[dados1$Resposta=="Sim"]), xlab="Colesterol", ylab="Densidade estimada", main = "Colesterol quando a resposta é afirmativa")
grid()
rug(dados1$Colesterol[dados1$Resposta=="Sim"])
plot(density(dados1$Colesterol[dados1$Resposta=="Não"]), xlab="Colesterol", ylab="Densidade estimada", main = "Colesterol quando a resposta é negativa")
grid()
rug(dados1$Colesterol[dados1$Resposta=="Não"])
Percebemos que a doença coronária se manifesta mais no sexo masculino e está clara a influença da idade mas não do colesterol e nem da duração dos sintomas.
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 cloglog (menor valor de AIC).
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 descartamos somente a variável explicativa Cad.dur.
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
## AIC: 2353.4
##
## Number of Fisher Scoring iterations: 5par(mfrow=c(1,1),mar=c(5,4,1,1),cex=1.0)
Transformada = I(log(dados1$Colesterol*dados1$Idade))
lattice::xyplot(Transformada ~ Resposta | Sexo, data = dados1)
plot(dados1$Resposta ~ Transformada)
modelo05 = glm(Resposta ~ Idade + Sexo + Colesterol + Transformada, family=binomial(link='cloglog'), data=dados1)Para calcular o R\(^2\) de Nagelkerke devemos utilizar a função lrm{rms}, The rms Package for R : Regression Modeling Strategies.
library(rms)
R2 = lrm(formula(modelo05), x=T, y=T, data=dados1)$stats[10]
R2## R2
## 0.300897O qual e um valor baixissimo!
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()
ResourceSelection::hoslem.test(modelo05$y,fitted(modelo05))##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: modelo05$y, fitted(modelo05)
## X-squared = 13.081, df = 8, p-value = 0.1091O teste Hosmer-Lemeshow de bondade de ajuste aceita o ajuste do modelo (modelo05), porém no gráfico de valores estimados e observados percebemos que há muito erro nas estimativas da resposta.
mu.preditos=predict(modelo05,type='response')PERGUNTA: para quais valores de estimativa da \(P(Y=1)\) assumir que o valor estimado da resposta é Y=1 e para quais Y=0?
Y.preditos1=ifelse(mu.preditos>0.5,1,0)
table(modelo05$y,Y.preditos1)## Y.preditos1
## 0 1
## 0 385 383
## 1 186 1304Y.preditos2=ifelse(mu.preditos>0.7,1,0)
table(modelo05$y,Y.preditos2)## Y.preditos2
## 0 1
## 0 592 176
## 1 496 994Y.preditos3=ifelse(mu.preditos>0.4,1,0)
table(modelo05$y,Y.preditos3)## Y.preditos3
## 0 1
## 0 303 465
## 1 114 1376Y.preditos4=ifelse(mu.preditos>0.9,1,0)
table(modelo05$y,Y.preditos4)## Y.preditos4
## 0 1
## 0 745 23
## 1 1207 283