Sabatina,
assunto: Métodos de classificação supervisionados
- Emerson Rigoni
- Henrique Aparecido Laureano [Lattes, GitLab, GitHub, LEG GitLab]
Maio de 2016
Banco de dados
library(MASS)
data("Pima.tr")
data("Pima.te")Uma população de mulheres (532) com no mínimo 21 anos, descendentes da tribo indígina Pima e atualmente vivendo próximo de Phoenix, Arizona, que fizeram o teste para diabetes de acordo o critério da Organização Mundial da Saúde (OMS). Os dados foram coletados pelo Instituto Nacional Norte Americano de Diabetes e Doenças de Digestivas e Renais. No banco de dados de treino (
Pima.tr) temos 200 indivíduos, e no banco de dados de teste (Pima.te) temos os demais 332 indivíduos
npreg- número de gestações
glu- concentração de glicose no plasma em um teste oral de tolerância à glicose
bp- pressão sanguínea diastólica (mm Hg)
skin- espessura da prega cutânea no tríceps (mm)
bmi- índice de massa corporal (\(peso (kg) / altura (m)^{2}\))
ped- função de diabetes pedigree
age- idade, em anos
typeYesouNopara diabetes de acordo com o critério da OMS
summary(Pima.tr) npreg glu bp skin bmi ped age type
Min. : 0.00 Min. : 56.0 Min. : 38.00 Min. : 7.00 Min. :18.20 Min. :0.0850 Min. :21.00 No :132
1st Qu.: 1.00 1st Qu.:100.0 1st Qu.: 64.00 1st Qu.:20.75 1st Qu.:27.57 1st Qu.:0.2535 1st Qu.:23.00 Yes: 68
Median : 2.00 Median :120.5 Median : 70.00 Median :29.00 Median :32.80 Median :0.3725 Median :28.00
Mean : 3.57 Mean :124.0 Mean : 71.26 Mean :29.21 Mean :32.31 Mean :0.4608 Mean :32.11
3rd Qu.: 6.00 3rd Qu.:144.0 3rd Qu.: 78.00 3rd Qu.:36.00 3rd Qu.:36.50 3rd Qu.:0.6160 3rd Qu.:39.25
Max. :14.00 Max. :199.0 Max. :110.00 Max. :99.00 Max. :47.90 Max. :2.2880 Max. :63.00 summary(Pima.te) npreg glu bp skin bmi ped age type
Min. : 0.000 Min. : 65.0 Min. : 24.00 Min. : 7.00 Min. :19.40 Min. :0.0850 Min. :21.00 No :223
1st Qu.: 1.000 1st Qu.: 96.0 1st Qu.: 64.00 1st Qu.:22.00 1st Qu.:28.18 1st Qu.:0.2660 1st Qu.:23.00 Yes:109
Median : 2.000 Median :112.0 Median : 72.00 Median :29.00 Median :32.90 Median :0.4400 Median :27.00
Mean : 3.485 Mean :119.3 Mean : 71.65 Mean :29.16 Mean :33.24 Mean :0.5284 Mean :31.32
3rd Qu.: 5.000 3rd Qu.:136.2 3rd Qu.: 80.00 3rd Qu.:36.00 3rd Qu.:37.20 3rd Qu.:0.6793 3rd Qu.:37.00
Max. :17.000 Max. :197.0 Max. :110.00 Max. :63.00 Max. :67.10 Max. :2.4200 Max. :81.00 Na aplicação dos métodos é utilizada type como resposta e bmi e ped como covariáveis
library(latticeExtra)
print(xyplot(ped ~ bmi
, groups = type
, col = c("#0080ff", "brown")
, pch = 19
, type = c("p", "g")
, main = "Pima.tr"
, key = list(space = "top"
, text = list(c("Diabetes: No", "Diabetes: Yes"))
, points = list(col = c("#0080ff", "brown"), pch = 19, cex = .8)
, columns = 2)
, xlim = c(15, 70)
, ylim = c(0, 2.5)
, Pima.tr)
, position = c(0, 0, .5, 1), more = TRUE)
print(xyplot(ped ~ bmi
, groups = type
, col = c("#0080ff", "brown")
, pch = 19
, type = c("p", "g")
, main = "Pima.te"
, key = list(space = "top"
, text = list(c("Diabetes: No", "Diabetes: Yes"))
, points = list(col = c("#0080ff", "brown"), pch = 19, cex = .8)
, columns = 2)
, xlim = c(15, 70)
, ylim = c(0, 2.5)
, Pima.te)
, position = c(.5, 0, 1, 1))
Regressão Logística
Função logística
\[ h(z) = \frac{1}{1 + e^{-z}} \]
h <- function(z) return(1 / (1 + exp(-z)))Função custo
\[ J(\theta) = \frac{1}{m} \sum_{i=1}^{m} - y_{i} {\rm log}(h_{\theta}(x_{i})) - (1 - y_{i}) {\rm log}(1 - (h_{\theta}(x_{i}))) \]
# covariáveis
X <- as.matrix(cbind(1, Pima.tr[ , c("bmi", "ped")]))
# resposta
y <- ifelse(Pima.tr$type == "No", 0, 1)
# parâmetros
theta <- matrix(0, nrow = 3, ncol = 1)
# custo
custo <- function(theta, X, y){
m <- length(y)
h <- h(X %*% theta)
J <- (1 / m) * (t(- y) %*% log(h) - (1 - t(y)) %*% log(1 - h))
return(J)}Fazendo a estimação pela optim
fit_opt <- optim(c(0, 0, 0), function(theta) custo(theta, X, y))
fit_opt$par[1] -4.37583880 0.09612438 1.17538045Fazendo o ajuste pela glm
fit_glm <- glm(type ~ bmi + ped
, family = "binomial"
, Pima.tr)
fit_glm$coefficients(Intercept) bmi ped
-4.37613317 0.09612468 1.17601545 Comparando os coeficientes estimados
fit_opt$par - fit_glm$coefficients (Intercept) bmi ped
2.943714e-04 -3.062838e-07 -6.350014e-04 Diferenças praticamente nulas (quarta casa decimal em diante)
O que o modelo estimou?
print(xyplot(ped ~ bmi
, groups = type
, col = c("#0080ff", "brown")
, pch = 19
, type = c("p", "g")
, main = "Pima.tr: Observado"
, key = list(space = "top"
, text = list(c("Diabetes: No", "Diabetes: Yes"))
, points = list(col = c("#0080ff", "brown"), pch = 19, cex = .8)
, columns = 2)
, Pima.tr)
, position = c(0, 0, .5, 1), more = TRUE)
print(xyplot(ped ~ bmi
, groups = ifelse(fit_glm$fitted.values < .5, "No", "Yes")
, col = c("#0080ff", "brown")
, pch = 19
, type = c("p", "g")
, main = "Pima.tr: Ajustado"
, key = list(space = "top"
, text = list(c("Diabetes: No", "Diabetes: Yes"))
, points = list(col = c("#0080ff", "brown"), pch = 19, cex = .8)
, columns = 2)
, Pima.tr)
, position = c(.5, 0, 1, 1))
Predição (usando os coeficientes da glm - tínhamos que escolher um deles)
print(xyplot(ped ~ bmi
, groups = type
, col = c("#0080ff", "brown")
, pch = 19
, type = c("p", "g")
, main = "Pima.te: Observado"
, key = list(space = "top"
, text = list(c("Diabetes: No", "Diabetes: Yes"))
, points = list(col = c("#0080ff", "brown"), pch = 19, cex = .8)
, columns = 2)
, Pima.te)
, position = c(0, 0, .5, 1), more = TRUE)
print(xyplot(ped ~ bmi
, groups = ifelse(predict(fit_glm
, data.frame(bmi = Pima.te$bmi, ped = Pima.te$ped)
, type = "response") < .5
, "No", "Yes")
, col = c("#0080ff", "brown")
, pch = 19
, type = c("p", "g")
, main = "Pima.te: Predito"
, key = list(space = "top"
, text = list(c("Diabetes: No", "Diabetes: Yes"))
, points = list(col = c("#0080ff", "brown"), pch = 19, cex = .8)
, columns = 2)
, Pima.te)
, position = c(.5, 0, 1, 1))
Discriminante de Fischer
library(MASS)Linear
Considera que as matrizes de covariância das covariáveis não diferem entre as classes (grupos)
fit_l <- lda(type ~ bmi + ped, Pima.tr)
pred_l <- predict(fit_l, Pima.te)
table(Pima.te$type, pred_l$class, dnn = list("Observado", "Predito")) Predito
Observado No Yes
No 195 28
Yes 65 44223 mulheres foram diagnosticadas com diabetes, e pelo discriminante linear 195 mulheres (87%) foram categorizadas como diabéticas
109 mulheres foram diagnosticadas sem diabetes, e pelo discriminante linear 44 mulheres (40%) foram categorizadas como não diabéticas
print(xyplot(ped ~ bmi
, groups = type
, col = c("#0080ff", "brown")
, pch = 19
, type = c("p", "g")
, main = "Pima.te: Observado"
, key = list(space = "top"
, text = list(c("Diabetes: No", "Diabetes: Yes"))
, points = list(col = c("#0080ff", "brown"), pch = 19, cex = .8)
, columns = 2)
, Pima.te)
, position = c(0, 0, .5, 1), more = TRUE)
print(xyplot(ped ~ bmi
, groups = pred_l$class
, col = c("#0080ff", "brown")
, pch = 19
, type = c("p", "g")
, main = "Pima.te: Predito"
, key = list(space = "top"
, text = list(c("Diabetes: No", "Diabetes: Yes"))
, points = list(col = c("#0080ff", "brown"), pch = 19, cex = .8)
, columns = 2)
, Pima.te)
, position = c(.5, 0, 1, 1))
Quadrático
Considera que as matrizes de covariância das covariáveis diferem entre as classes
fit_q <- qda(type ~ bmi + ped, Pima.tr)
pred_q <- predict(fit_q, Pima.te)
table(Pima.te$type, pred_q$class, dnn = list("Observado", "Predito")) Predito
Observado No Yes
No 204 19
Yes 73 36223 mulheres foram diagnosticadas com diabetes, e pelo discriminante linear 204 mulheres (91%) foram categorizadas como diabéticas
109 mulheres foram diagnosticadas sem diabetes, e pelo discriminante linear 36 mulheres (33%) foram categorizadas como não diabéticas
print(xyplot(ped ~ bmi
, groups = type
, col = c("#0080ff", "brown")
, pch = 19
, type = c("p", "g")
, main = "Pima.te: Observado"
, key = list(space = "top"
, text = list(c("Diabetes: No", "Diabetes: Yes"))
, points = list(col = c("#0080ff", "brown"), pch = 19, cex = .8)
, columns = 2)
, Pima.te)
, position = c(0, 0, .5, 1), more = TRUE)
print(xyplot(ped ~ bmi
, groups = pred_q$class
, col = c("#0080ff", "brown")
, pch = 19
, type = c("p", "g")
, main = "Pima.te: Predito"
, key = list(space = "top"
, text = list(c("Diabetes: No", "Diabetes: Yes"))
, points = list(col = c("#0080ff", "brown"), pch = 19, cex = .8)
, columns = 2)
, Pima.te)
, position = c(.5, 0, 1, 1))
Regularizado
Considera dois parâmetros que flexibilizam a possível diferença entre as matrizes de covariância das covariáveis entre as classes e a dependência entre as mesmas covariáveis
library(klaR)
fit_r <- rda(type ~ bmi + ped, Pima.tr)
pred_r <- predict(fit_r, Pima.te)
table(Pima.te$type, pred_r$class, dnn = list("Observado", "Predito")) Predito
Observado No Yes
No 169 54
Yes 65 44223 mulheres foram diagnosticadas com diabetes, e pelo discriminante linear 169 mulheres (76%) foram categorizadas como diabéticas
109 mulheres foram diagnosticadas sem diabetes, e pelo discriminante linear 44 mulheres (40%) foram categorizadas como não diabéticas
print(xyplot(ped ~ bmi
, groups = type
, col = c("#0080ff", "brown")
, pch = 19
, type = c("p", "g")
, main = "Pima.te: Observado"
, key = list(space = "top"
, text = list(c("Diabetes: No", "Diabetes: Yes"))
, points = list(col = c("#0080ff", "brown"), pch = 19, cex = .8)
, columns = 2)
, Pima.te)
, position = c(0, 0, .5, 1), more = TRUE)
print(xyplot(ped ~ bmi
, groups = pred_r$class
, col = c("#0080ff", "brown")
, pch = 19
, type = c("p", "g")
, main = "Pima.te: Predito"
, key = list(space = "top"
, text = list(c("Diabetes: No", "Diabetes: Yes"))
, points = list(col = c("#0080ff", "brown"), pch = 19, cex = .8)
, columns = 2)
, Pima.te)
, position = c(.5, 0, 1, 1))
Comparando os métodos
print(xyplot(ped ~ bmi
, groups = ifelse(predict(fit_glm
, data.frame(bmi = Pima.te$bmi, ped = Pima.te$ped)
, type = "response") < .5
, "No", "Yes")
, col = c("#0080ff", "brown")
, pch = 19
, type = c("p", "g")
, main = "Regressão Logística"
, key = list(space = "top"
, text = list(c("Diabetes: No", "Diabetes: Yes"))
, points = list(col = c("#0080ff", "brown"), pch = 19, cex = .8)
, columns = 2)
, Pima.te)
, position = c(0, .5, .5, 1), more = TRUE)
print(xyplot(ped ~ bmi
, groups = pred_l$class
, col = c("#0080ff", "brown")
, pch = 19
, type = c("p", "g")
, main = "Discriminante linear"
, key = list(space = "top"
, text = list(c("Diabetes: No", "Diabetes: Yes"))
, points = list(col = c("#0080ff", "brown"), pch = 19, cex = .8)
, columns = 2)
, Pima.te)
, position = c(.5, .5, 1, 1), more = TRUE)
print(xyplot(ped ~ bmi
, groups = pred_q$class
, col = c("#0080ff", "brown")
, pch = 19
, type = c("p", "g")
, main = "Discriminante quadrático"
, key = list(space = "top"
, text = list(c("Diabetes: No", "Diabetes: Yes"))
, points = list(col = c("#0080ff", "brown"), pch = 19, cex = .8)
, columns = 2)
, Pima.te)
, position = c(0, 0, .5, .5), more = TRUE)
print(xyplot(ped ~ bmi
, groups = pred_r$class
, col = c("#0080ff", "brown")
, pch = 19
, type = c("p", "g")
, main = "Discriminante regularizado"
, key = list(space = "top"
, text = list(c("Diabetes: No", "Diabetes: Yes"))
, points = list(col = c("#0080ff", "brown"), pch = 19, cex = .8)
, columns = 2)
, Pima.te)
, position = c(.5, 0, 1, .5), more = TRUE)
Resultados muito próximos são obtidos com a regressão logística e com a análise descriminante linear de Fischer
Nenhum dos métodos chega em resultados próximos a classificação original dos dados, que é bem complicada de ser atingida, já que é muito aleatória. Uma ótima justificativa pra essa aleatoriedade são as covariáveis utilizadas. Talvez com as outras covariáveis a clasificação se torne muito mais clara e simples, consequentemente
As covariáveis utilizadas foram definidas / estabelecidas na proposta da tarefa, sabatina
Curvas ROC
library(pROC)
par(mfrow = c(2, 2))
plot.roc(roc(Pima.te$type, predict(fit_glm
, data.frame(bmi = Pima.te$bmi, ped = Pima.te$ped)
, type = "response"))
, print.auc = TRUE, print.thres = TRUE, las = 1
, xlab = "Especificidade", ylab = "Sensibilidade", main = "Regressão Logística")
plot.roc(roc(Pima.te$type, pred_l$posterior[ , 2]), print.auc = TRUE, print.thres = TRUE, las = 1
, xlab = "Especificidade", ylab = "Sensibilidade", main = "Discriminante linear")
plot.roc(roc(Pima.te$type, pred_q$posterior[ , 2]), print.auc = TRUE, print.thres = TRUE, las = 1
, xlab = "Especificidade", ylab = "Sensibilidade", main = "Discriminante quadrático")
plot.roc(roc(Pima.te$type, pred_r$posterior[ , 2]), print.auc = TRUE, print.thres = TRUE, las = 1
, xlab = "Especificidade", ylab = "Sensibilidade", main = "Discriminante regularizado")
Pelo AUC o modelo com melhor poder preditivo é o discriminante linear, contudo, os valores de área abaixo da curva não diferem muito (maior diferença de 0.051)
Especificidade: Proporção de resultados negativos em indivíduos livres da doença
Sensibilidade: Proporção de resultados positivos em indivíduos com a doença
A doença aqui em questão é a diabetes, e um maior peso pode ser atribuído a sensibilidade, i.e, é mais importante termos uma baixa proporção de indivíduos doentes com resultado negativo no exame (consequentemente não sendo tratados) do que uma baixa proporção de indivíduos livres da doença com resultado positivo no exame (consequentemente sendo tratados)
Olhando para a sensibilidade o melhor modelo é o discriminante linear, e olhando para a sensibilidade e a especificade, juntas, o melhor modelo ainda é o discriminante linear
Os resultados do discriminante linear são muito próximos dos resultados obtidos com a regressão logística