O presente trabalho tem como finalidade apresentar e discutir aplicações de alguns métodos de classificação, no caso: Regressão Logística, Discriminante Linear e Quadrático de Fisher e extensões.

Pacotes utilizados:

require(MASS)
require(xtable)
require(lattice)
require(hmeasure)
require(klaR)
require(ggplot2)
require(gridExtra)
require(pROC)


1. Regressão Logística


O objetivo é encontrar uma regressão para a chance de ser aprovado em uma universidade dados suas notas em exames anteriores. Os dados referem-se a históricos de candidatos as vagas com suas notas anteriores e o resultado deles na universidade.

Visualização do comportamento das notas nos dois exames.

plot(exemplo[,1], exemplo[,2], las = 1,
     col = ifelse(exemplo$status==0,"red","blue"),
     pch = 20, xlab="Exame 1",ylab="Exame 2")
legend(x=85, y = 104 ,c("Aprovado","Reprovado"),
       pch = 20, col=c("blue","red"), bty ="n")

A função logística é dada por:

\[g(z) = \frac{1}{1+e^{-z}}\]

Com isso teremos valores entre 0 e 1.

sigmoid<- function(z) {
  return(1/(1 + exp(-z)))
}

A função custo para um grupo de paramêtros será

\[J(\theta)= \frac{1}{m} \sum\limits_{i=1}^m [-y^{(i)}log(h_\theta(x^{(i)}))-(1-y^{(i)})log(1-h_{\theta}(x^{(i)}))]\]

e o gradiente do custo será um vetor de tamanho igual a \(\theta\) onde o \(j\) é o ésimo elemento (j=0,1, . , n):

\[\frac{\partial J(\theta)}{\partial \theta_J} = \frac{1}{m} \sum\limits_{i=1}^m (h_{\theta}(x^{(i)})-y^{(i)})x^{(i)}_{j}\]

Com isso geramos a função costFunction. Essa função retornará o custo de usar esse conjunto de parâmetros.

#Parâmetros de entrada
X <- as.matrix( cbind(1, exemplo[, 1:2])) # Matriz das variáveis
theta <- matrix(0, nrow = ncol(exemplo[, 1:2]) + 1, ncol = 1)
y <- exemplo[, 3] # Resposta

#Função de custo
costFunction<- function(theta, X, y) {
  m <- length(y)
  J <- 0
  grad <- matrix(0,nrow=length(theta)+1,ncol=1)
  h <- sigmoid(X %*% theta)
  costPos <- t(-y) %*% log(h)
  costNeg <- (1 - t(y)) %*% log(1 - h)
  J = (1/m) * (costPos - costNeg)
  grad = (1/m) * (t(X) %*% (h - y))
  return(c(J))
}

costFunction(theta, X, y)
## [1] 0.6931472

O objetivo é escolher o grupo de parâmetros, que vai gerar o menor custo, e o maior verossimil. Utilizando optim, que é uma função de otimização de uso geral baseado em algoritmos de gradiente.

#Ajuste da regressão
ajuste<-optim(c(0,0,0),function(theta) 
  {costFunction(theta,X,y)},method="Nelder-Mead")

ajuste
## $par
## [1] -25.1647048   0.2062524   0.2015046
## 
## $value
## [1] 0.2034977
## 
## $counts
## function gradient 
##      156       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
#Custo
costFunction(matrix(ajuste$par, 
                    nrow = ncol(exemplo [, 1:2])+ 1, ncol = 1), X, y)
## [1] 0.2034977

Na saída da função, temos os valores finais de parâmetros no par, o J no value e o total de iterações em counts. Em convergence, 0 indica se a conclusão foi bem sucedida.

Verificando ajuste gerado pela função com o ajuste gerado pela função glm.

ajuste$par #Ajuste pela função
## [1] -25.1647048   0.2062524   0.2015046
glm(status ~ exame1 + exame2,
    data = exemplo ,family = "binomial") #Ajuste glm
## 
## Call:  glm(formula = status ~ exame1 + exame2, family = "binomial", 
##     data = exemplo)
## 
## Coefficients:
## (Intercept)       exame1       exame2  
##    -25.1613       0.2062       0.2015  
## 
## Degrees of Freedom: 99 Total (i.e. Null);  97 Residual
## Null Deviance:       134.6 
## Residual Deviance: 40.7  AIC: 46.7

Observa-se que os valores ajustados pela função coincidem com os valores gerados pelo glm.

No gráfico, temos de um lado temos os indivíduos que esperamos que serão aprovados, e do outro o que esperamos que não serão.

x1 <- c(min(X[,2])-2,max(X[,2])+2)
y1 <- (-1/ajuste$par[3] ) * (ajuste$par[2] * x1 + ajuste$par[1] )
plot(exemplo[,1], exemplo[,2],
     col = ifelse(exemplo$status==0, "red", "blue"),
     pch=19, xlab = "Exame 1", ylab = "Exame 2",
     ylim=c(30,110))
lines(x1, y1, lty = 2, lwd=2 ,col="black")
legend(x = 83, y = 115,c("Aprovado","Reprovado","Margem"),
       pch = c(19,19,NA), lty = c(0,0,2), lwd=2,
       col = c("blue","red","black"), bty = "n")

Com a função encontrada, podemos prever alguns resultados. Se um indivíduo tem suas notas anteriores a 80, ele terá 99% de chance de ser aprovado nessa universidade. E assim por diante para os demais casos.

# Indivíduo com nota 80 e nota 80
sigmoid(ajuste$par[1] + ajuste$par[2] * 80 + ajuste$par[3] * 80)
## [1] 0.9994223

Algumas outras previsões. Para um indivíduo com nota 90 na primeira prova, e 40 na segunda, tem 81% de chance de ser aprovado.

# Indivíduo com nota 90 e nota 40
sigmoid(ajuste$par[1] + ajuste$par[2] * 90 + ajuste$par[3] * 40)
## [1] 0.8112565

Já para o indivíduo com nota 70 na primeira prova, e 30 na segunda, terá 34% de chance de ser aprovado.

# Indivíduo com nota 70 e nota 50
sigmoid(ajuste$par[1] + ajuste$par[2] * 70 + ajuste$par[3] * 50)
## [1] 0.3425826

E aquele com notas anteriores igual a 50, terá 8% de chance de ser aprovado.

# Indivíduo com nota 50 e nota 50
sigmoid(ajuste$par[1] + ajuste$par[2] * 50 + ajuste$par[3] * 50)
## [1] 0.008352106

Retirado Recologia

2. Método Discriminante de Fisher


A análise descriminante basea-se na idéia bayesiana, de que temos como informação \(f(\textbf{X}=\textbf{x}|Y=k)\), com \(\textbf{X}\) um vetor de variáveis aleatórias e \(Y\) uma variável nominal com k níveis. A ideia é obter um classificador a partir de dados de treino, e predizer a classificação de um objeto ou indivíduo a partir de um conjunto de variáveis \(\textbf{X}\). Ou seja, queremos \(f(Y=k|\textbf{X}=x)\), então utiliza-se a regra de Bayes de forma: \(f(Y=k|\textbf{X}=x) \propto f(\textbf{X} \textbf{x}|Y) \pi(Y)\). A análise discriminante de Fisher, tradicionalmente pressupõe que \(\textbf{X}|Y=k \sim NM(\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)\).

\(\pi(Y)\) é uma probabilidade à priori de classificação. Quando não se tem conhecimento desta probabilidade, utiliza-se a proporção de observações contidas em cada classe, que seria equivalente uma priori do tipo flat.

Depois de obtido um classificador com os dados de treino, nós queremos saber, a partir de novos objetos mensurados \(\textbf{x_1, x_2, ..., x_n}\), quais as suas classificações. Deste modo, a ideia é identificar a classe k da variável \(Y\) que maximiza \(f(\textbf{X} \textbf{x}|Y) \pi(Y)\), parecida com a idéia de TRV.

Quando a \(f(\textbf{X} \textbf{x}|Y) \) é assumido distribuição normal, chamamos de discriminante de Fisher. Queremos o \(k\) que maximize,então:

\[\delta_k(k) \propto \text{ argmax } \pi_k(k) f_k(x)\].

Passando o logarítimo tem-se:

\[=\text{ argmax } \left \{ \log{\pi_k} - \frac{1}{2}\log{\mid \boldsymbol{\Sigma}_k \mid} - \frac{1}{2} (\textbf{x} - \boldsymbol{\mu}_k){}' \boldsymbol{\Sigma}_k^{-1} (\textbf{x} - \boldsymbol{\mu}_k)\right \}\]

Um exemplo seria um estudo para saber se o indivíduo tem ou não determinado tipo de câncer, realizado com n indivíduos, na qual estes foram submetidos ao exame padrão ouro. Neste mesmo estudo, foram registradas diversas variáveis a respeito do paciente, como idade, sexo, taxa colesterol, entre outras. A partir desta classificação à priori, com os dados de novos indivíduos na qual não foram submetidos ao exame padrão ouro, utiliza-se as variáveis de mais fácil mensuração para predizer se tem ou não o câncer.


Alguns outros exemplos como:

  1. Classificar vasos de porcelana chinenses conforme o período histórico ou dinastia, a partir de variáveis de medidas como, altura, largura, tipos de desenhos e etc.

  2. Classificar um indivíduo tendo ou não x patologia a partir de medidas de um desenho feito em sessão.

  3. Classificar o indivíduo em bom ou mau pagador a partir de informações socioeconômicas do mesmo.


2.1. Função Linear



Ocorre quando \(\boldsymbol{\Sigma}_i = \boldsymbol{\Sigma}_j \text{, } \forall i \neq j\), com \(i \text{, }j \in k\)

Utilizando a função lda() do pacote MASS para realizar a análise discriminante linear de Fisher. A planilha de dados foi dividida em duas partes, a primeira contendo 50 observações para treino e outras 50 para teste.

set.seed(56567567)

k<- sample(100,50)


treino<- exemplo[k,]

teste<- exemplo[-k,]



fit1<- lda(status~.,data=exemplo,subset = k)

PredVC<- predict(fit1,teste)

A função lda() permite incluir fórmula da variável resposta \(\sim\) as variáveis explanatórias, funciona da mesma maneira que no lm() ou glm(). Existe a opção de mudar as proporções à priori, se não declaradas, são usadas as proporções de observações em cada classe.

Na predição, existem três alternativas: declar por fora dividindo a base de dados de treino e teste e ajustando pela primeira;usando o argumento subset que serve para declarar as linhas da base de dados que serão utilizadas para o ajuste; e uma terceira opção é utilizar o argumento CV=T que realiza validação cruzada pelo método Leave-one-out.



Tabela de contingência dos preditos versus observados:

tc <- table(teste$status, PredVC$class)
dimnames(tc) <- list(Observados = c("Reprovado", "Aprovado"), "Predictos (cv)" = c("Reprovado", "Aprovado"))
print(round(tc))
##            Predictos (cv)
## Observados  Reprovado Aprovado
##   Reprovado        13        4
##   Aprovado          4       29


A tabela de contingência é construída por meio do cruzamento das classes dos valores observados e preditos. Como no exemplo existem duas classes (Aprovado ou Reprovado), então, existem 4 cenários possíveis:

  • Predito como Reprovado sendo que defato o indivíduo Reprovou,

  • Classificar como Aprovado dado que o indivíduo foi Aprovado (as duas anteriores consistem os acertos, sempre na diagonal),

  • Classificar como Reprovado sendo que na verdade o indivíduo foi Aprovado

  • Classificar como Aprovado sendo que na realidade o indivíduo veio a Reprovar.

A partir desta configuração, são dispostas as contagens em uma tabela conforme cada um dos cenários anteriores.


Tabela de proporções de acertos e erros

round(prop.table(tc, 1),3)
##            Predictos (cv)
## Observados  Reprovado Aprovado
##   Reprovado     0.765    0.235
##   Aprovado      0.121    0.879

Disposta na forma de proporção por linha (observados), ou seja, a proporção de classificação correta como Reprovado foi de aproximadamente \(95 \%\), enquanto a porporção de classificação correta de Aprovado foi de \(85 \%\).


Proporção geral de acertos:

sum(diag(prop.table(tc)))
## [1] 0.84

Essa é a proporção geral de acertos, que basicamente é a proporção de elementos na diagonal em relação ao total.


Gráfico dos preditos:

###############################################

densityplot(~PredVC$x, groups=teste$status,bw=1,xlab="Discriminante",ylab="Densidade")

ldahist(PredVC$x[,1],g=teste$status)
title("Discriminante 1")

O primeiro gráfico mostra a densidade estimada dos valores preditos em cada classe, já o segundo, é o gráfico de cima, porém, em abcissas diferentes, e serve para visualizar a intersecção e a forma de divisão feita por cada discriminante.



Estudando a sensibilidade da priori. Foi usada uma sequência de 9 proporções diferentes para declarar como priori, do resultado da predição foram extraídas as proporções de acertos por classe e a proporção geral.


Priori no primeiro exemplo

fit1$prior
##    0    1 
## 0.46 0.54


priorVals<- matrix(rbind(seq(0.1,0.9,by=0.1),1-seq(0.1,0.9,by=0.1)),ncol=2,byrow = T)
results<-list()
tcpriori2<- list()
for(i in 1:9){

fit1priori2<- lda(status~.,data=treino,prior=priorVals[i,])

PredVCpriori2<- predict(fit1priori2,teste)

tcpriori2[[i]] <- table(teste$status, PredVCpriori2$class)
dimnames(tcpriori2[[i]]) <- list(Observados = c("Reprovado", "Aprovado"), "Predictos (cv)" = c("Reprovado", "Aprovado"))

results[[i]]<- data.frame(priori=priorVals[i,],propClass=round(diag(prop.table(tcpriori2[[i]],1)),2),propGeral=sum(diag(prop.table(tcpriori2[[i]]))))

}; results
## [[1]]
##           priori propClass propGeral
## Reprovado    0.1      0.65      0.88
## Aprovado     0.9      1.00      0.88
## 
## [[2]]
##           priori propClass propGeral
## Reprovado    0.2      0.76       0.9
## Aprovado     0.8      0.97       0.9
## 
## [[3]]
##           priori propClass propGeral
## Reprovado    0.3      0.76      0.88
## Aprovado     0.7      0.94      0.88
## 
## [[4]]
##           priori propClass propGeral
## Reprovado    0.4      0.76      0.86
## Aprovado     0.6      0.91      0.86
## 
## [[5]]
##           priori propClass propGeral
## Reprovado    0.5      0.76      0.82
## Aprovado     0.5      0.85      0.82
## 
## [[6]]
##           priori propClass propGeral
## Reprovado    0.6      0.88      0.86
## Aprovado     0.4      0.85      0.86
## 
## [[7]]
##           priori propClass propGeral
## Reprovado    0.7      0.94      0.88
## Aprovado     0.3      0.85      0.88
## 
## [[8]]
##           priori propClass propGeral
## Reprovado    0.8      0.94      0.78
## Aprovado     0.2      0.70      0.78
## 
## [[9]]
##           priori propClass propGeral
## Reprovado    0.9      0.94      0.78
## Aprovado     0.1      0.70      0.78


Pelas tabelas acima, nota-se que conforme muda a informação à priori, as proporções de acertos mudam bastante, devido, também ao conjunto de dados pequeno. Na função qda(), essa mudança das proporções ocorre em uma magnitude muito menor, por estar incluindo mais informações da amostra através do termo quadrático.



2.2. Função quadrática



A função quadrática ocorre quando \(\boldsymbol{\Sigma}_i \neq \boldsymbol{\Sigma}_j\), ou seja, terá o termo quadrático \(\textbf{x}{}' \boldsymbol{\Sigma_k}^{-1} \textbf{x}\), com \(k\) a classe.

fit1Q<- qda(status~.,data=exemplo,subset = k)

PredQVC<-predict(fit1Q,teste)



Tabela de contingência dos preditos versus observados:

tcQ <- table(teste$status, PredQVC$class)
dimnames(tcQ) <- list(Observados = c("Reprovado", "Aprovado"), "Predictos (cv)" = c("Reprovado", "Aprovado"))
print(round(tcQ))
##            Predictos (cv)
## Observados  Reprovado Aprovado
##   Reprovado        14        3
##   Aprovado          5       28


Proporção de acertos em cada grupo

round(prop.table(tcQ, 1),3)
##            Predictos (cv)
## Observados  Reprovado Aprovado
##   Reprovado     0.824    0.176
##   Aprovado      0.152    0.848


Proporção geral de acertos:

sum(diag(prop.table(tcQ)))
## [1] 0.84


Comparando com o resultado obtido pela discriminante linear (lda()), obtivemos um pequeno ganho na proporção de Aprovados ao considerar o termo quadrático, ou seja, matrizes de covariâncias das variáveis diferentes entre os grupos.



2.3. Análise discriminante regularizada


Pode-se dizer que a análise discriminante quadrática é mais geral, no sentido de considerar matrizes de covariâncias das variáveis explanatórias diferentes, tendo assim, o termo quadrático. A análise discriminante linear é um caso particular, onde é pressuposto que estas matrizes de covariâncias sejam iguais para todas as classes. Porém, podemos flexibilizar isto incluindo um parâmetro ponderador:

\[\boldsymbol{\hat{\Sigma}}_k(\lambda) = (1- \lambda)\boldsymbol{\hat{\Sigma}}_k + \lambda \boldsymbol{\hat{\Sigma}} \]

Além disso, se assumisse independência entre as variáveis, ou seja, matrizes de covariâncias apenas com valores nas diagonais, teríamos um caso particular da análise discriminante de Fisher, conhecido como Naive Bayes. Usando também esta idéia, a análise discriminante regularizada também flexibiliza neste sentido, adicionando um parâmetro \(\gamma\), ou seja:

\[\boldsymbol{\hat{\Sigma}}_k(\lambda, \gamma) = (1- \gamma)\boldsymbol{\hat{\Sigma}}_k (\lambda) + \gamma \frac{1}{d} tr \left [\boldsymbol{\hat{\Sigma}}_k (\lambda)\right ]\textbf{I} \].


Os parâmetros \(\lambda\) e \(\gamma\) variam de 0 a 1, o primeiro flexibiliza quanto matrizes de covariâncias diferentes entre as classes, o segundo flexibiliza quanto a dependência entre as variáveis explicativas. No caso, se \(\lambda=0\) e \(\gamma=0\), temos o caso da análise discriminante quadrática, se \(\lambda=0\) e \(\gamma=1\) temos a análise discriminante linear. Na função rda() do pacote klaR, temos os argumentos das funções anteriores, mais a opção de fixar o valor dos parâmetros \(\lambda\) e \(\gamma\), assim como podemos estimá-los, afim de que os dados decidam a melhor ponderação.


Foram considerados três casos, primeiro a análogas da discriminante lineare posteriormente a quadrática, e o terceiro é estimando os parâmetros de ponderação. 

fitRLDA<- rda(status~.,data = treino,gamma = 0,lambda = 1)
PredRLDA<- predict(fitRLDA,teste[,1:2])



fitRQDA<- rda(status~.,data = treino,gamma = 0,lambda = 0)
PredRQDA<- predict(fitRQDA,teste[,1:2])


fitREst<- rda(status~.,data = treino)
PredREst<- predict(fitREst,teste[,1:2])



Tabela de contingência dos preditos versus observados:

tcRLDA <- table(teste$status, PredRLDA$class)
dimnames(tcRLDA) <- list(Observados = c("Reprovado", "Aprovado"), "Predictos (cv)" = c("Reprovado", "Aprovado"))
print(round(tcRLDA))
##            Predictos (cv)
## Observados  Reprovado Aprovado
##   Reprovado        13        4
##   Aprovado          4       29
tcRQDA <- table(teste$status, PredRQDA$class)
dimnames(tcRQDA) <- list(Observados = c("Reprovado", "Aprovado"), "Predictos (cv)" = c("Reprovado", "Aprovado"))
print(round(tcRQDA))
##            Predictos (cv)
## Observados  Reprovado Aprovado
##   Reprovado        14        3
##   Aprovado          5       28
tcREst <- table(teste$status, PredREst$class)
dimnames(tcREst) <- list(Observados = c("Reprovado", "Aprovado"), "Predictos (cv)" = c("Reprovado", "Aprovado"))
print(round(tcREst))
##            Predictos (cv)
## Observados  Reprovado Aprovado
##   Reprovado        14        3
##   Aprovado          4       29


Proporção de acertos em cada grupo

round(prop.table(tcRLDA, 1),3)
##            Predictos (cv)
## Observados  Reprovado Aprovado
##   Reprovado     0.765    0.235
##   Aprovado      0.121    0.879
round(prop.table(tcRQDA, 1),3)
##            Predictos (cv)
## Observados  Reprovado Aprovado
##   Reprovado     0.824    0.176
##   Aprovado      0.152    0.848
round(prop.table(tcREst, 1),3)
##            Predictos (cv)
## Observados  Reprovado Aprovado
##   Reprovado     0.824    0.176
##   Aprovado      0.121    0.879


Proporção geral de acertos:

sum(diag(prop.table(tcRLDA)))
## [1] 0.84
sum(diag(prop.table(tcRQDA)))
## [1] 0.84
sum(diag(prop.table(tcREst)))
## [1] 0.86

Nota-se que os resultados dos modelos análogos ao discriminante linear (lda()) e quadrático (qda()) foram os mesmos, e estimando os parâmetros obtem-se o mesmo resultado da quadrática (neste exemplo).

2.4. Gráficos de comparação de curvas

fit<- glm(status ~ exame1 + exame2,
data = treino ,family = "binomial")
 predGLM<- predict(fit,teste,type="response")
 




contour_data <- expand.grid(exame1 = seq(0, 100, length = 200),
exame2 = seq(0, 100, length = 200))

logistic_predict<- data.frame(contour_data,
status = as.numeric(ifelse(predict(fit, contour_data,type="response")>0.5,1,0)))
lda_predict <- data.frame(contour_data,
status = as.numeric(predict(fit1, contour_data)$class))
qda_predict <- data.frame(contour_data,
status = as.numeric(predict(fit1Q, contour_data)$class))
rda_predict <- data.frame(contour_data,
status = as.numeric(predict(fitREst, contour_data)$class))



p <- ggplot(teste, aes(x = exame1, y = exame2, color = status)) + geom_point()
v1<- p + stat_contour(aes(x = exame1, y = exame2, z = status), data = logistic_predict)+ggtitle("Reg. Logistica")
v2<- p + stat_contour(aes(x = exame1, y = exame2, z = status), data = lda_predict)+ggtitle("LDA")
v3<- p + stat_contour(aes(x = exame1, y = exame2, z = status), data = qda_predict)+ggtitle("QDA")
v4<- p + stat_contour(aes(x = exame1, y = exame2, z = status), data = rda_predict)+ggtitle("RDA")

grid.arrange(v1,v2,v3,v4,ncol=2)


2.5. Curva ROC para a análise discriminante linear, quadrática, regularizada e, regressão logística

HRL<- roc(teste$status,predGLM)
HLDL<- roc(teste$status,PredVC$posterior[,2])
HLDQ<- roc(teste$status,PredQVC$posterior[,2])
HRD<- roc(teste$status,PredREst$posterior[,2])

par(mfrow=c(2,2))
plot.roc(HRL,print.auc = T,auc.polygon.col = "blue",max.auc.polygon = T,print.thres = T,main="Reg. Logística")
## 
## Call:
## roc.default(response = teste$status, predictor = predGLM)
## 
## Data: predGLM in 17 controls (teste$status 0) < 33 cases (teste$status 1).
## Area under the curve: 0.9572
plot.roc(HLDL,print.auc = T,auc.polygon.col = "blue",max.auc.polygon = T,print.thres = T,main="LDA")
## 
## Call:
## roc.default(response = teste$status, predictor = PredVC$posterior[,     2])
## 
## Data: PredVC$posterior[, 2] in 17 controls (teste$status 0) < 33 cases (teste$status 1).
## Area under the curve: 0.9519
plot.roc(HLDQ,print.auc = T,auc.polygon.col = "blue",max.auc.polygon = T,print.thres = T,main="QDA")
## 
## Call:
## roc.default(response = teste$status, predictor = PredQVC$posterior[,     2])
## 
## Data: PredQVC$posterior[, 2] in 17 controls (teste$status 0) < 33 cases (teste$status 1).
## Area under the curve: 0.959
plot.roc(HRD,print.auc = T,auc.polygon.col = "blue",max.auc.polygon = T,print.thres = T,main="RDA")

## 
## Call:
## roc.default(response = teste$status, predictor = PredREst$posterior[,     2])
## 
## Data: PredREst$posterior[, 2] in 17 controls (teste$status 0) < 33 cases (teste$status 1).
## Area under the curve: 0.959



3. Segundo conjunto de dados

3.1. Análise discriminante linear

Aplicando a LDA para o conjunto de dados wine do R. O Conjunto de dados wine do pacote rattle, possui 13 variáveis de características de 178 observações de 3 tipos de cultivos de vinho diferentes. A ideia é achar um classificador em um dos três tipos segundo as características das variáveis apresentadas.

data(wine, package='rattle')

t1<- xtable(head(wine))

print(t1,type = "html")
Type Alcohol Malic Ash Alcalinity Magnesium Phenols Flavanoids Nonflavanoids Proanthocyanins Color Hue Dilution Proline
1 1 14.23 1.71 2.43 15.60 127 2.80 3.06 0.28 2.29 5.64 1.04 3.92 1065
2 1 13.20 1.78 2.14 11.20 100 2.65 2.76 0.26 1.28 4.38 1.05 3.40 1050
3 1 13.16 2.36 2.67 18.60 101 2.80 3.24 0.30 2.81 5.68 1.03 3.17 1185
4 1 14.37 1.95 2.50 16.80 113 3.85 3.49 0.24 2.18 7.80 0.86 3.45 1480
5 1 13.24 2.59 2.87 21.00 118 2.80 2.69 0.39 1.82 4.32 1.04 2.93 735
6 1 14.20 1.76 2.45 15.20 112 3.27 3.39 0.34 1.97 6.75 1.05 2.85 1450


Dividiu-se o conjunto de dados em dois, uma parte contendo 60 observações para treino e as outras 118 observações para teste.

kv<- sample(178,60)

testeWine<- wine[-kv,]

fit2<- lda(Type ~.,data=wine,subset = kv)

PredCV2<- predict(fit2,testeWine)

A partir da predição, serão extraídas informações sobre a qualidade da predição da discriminação obtida.



Tabela de contingência dos preditos versus observados:

tc2 <- table(testeWine$Type, PredCV2$class)
dimnames(tc2) <- list(Observados = as.character(unique(wine$Type)), "Predictos (cv)" = as.character(unique(wine$Type)))
print(round(tc2))
##           Predictos (cv)
## Observados  1  2  3
##          1 39  0  0
##          2  1 44  3
##          3  0  0 31


Proporção de acertos em cada grupo

round(prop.table(tc2, 1),3)
##           Predictos (cv)
## Observados     1     2     3
##          1 1.000 0.000 0.000
##          2 0.021 0.917 0.062
##          3 0.000 0.000 1.000


Proporção geral de acertos:

sum(diag(prop.table(tc2)))
## [1] 0.9661017


Gráfico dos preditos:

###############################################
#Primeira discriminante

ldahist(PredCV2$x[,1],g=testeWine$Type)
title("Discriminante 1")

#Segunda discriminante
ldahist(PredCV2$x[,2],g=testeWine$Type)
title("Discriminante 2")

####################

# Ou de forma sememlhante
densityplot(~PredCV2$x[,1], groups=testeWine$Type,bw=0.4,auto.key = T,xlab="Discriminante 1",ylab="Densidade")

densityplot(~PredCV2$x[,2], groups=testeWine$Type,bw=0.4,auto.key = T,xlab="Discriminante 2",ylab="Densidade")


Pelos gráficos acima, pode-se concluir que a primeira discriminante conseguiu separar muito bem a primeira classe de cultivo das demais e de forma razoável a segunda da terceira, já, a segunda discriminante separa os cultivos das classes 1 e 3 da segunda, conforme os gráficos mostrados acima.



3.2. Análise discriminante quadrática


Para a base de dados wine, temos:

fitQda<- qda(Type ~.,data=wine,subset=kv)

PredQda<- predict(fitQda,testeWine)



Tabela de contingência dos preditos versus observados:

tcQda <- table(testeWine$Type, PredQda$class)
dimnames(tcQda) <- list(Observados = as.character(unique(wine$Type)), "Predictos (cv)" = as.character(unique(wine$Type)))
print(round(tcQda))
##           Predictos (cv)
## Observados  1  2  3
##          1 19 20  0
##          2  1 47  0
##          3  0  7 24


Proporção de acertos em cada grupo

prop.table(tcQda, 1)
##           Predictos (cv)
## Observados          1          2          3
##          1 0.48717949 0.51282051 0.00000000
##          2 0.02083333 0.97916667 0.00000000
##          3 0.00000000 0.22580645 0.77419355


Proporção geral de acertos:

sum(diag(prop.table(tcQda)))
## [1] 0.7627119

No primeiro exemplo, tivemos um pequeno ganho na predição considerando as matrizes de covariâncias das variáveis explicativas diferentes entre os grupos, porém, para este segundo exemplo, observou-se uma perda ao considerar o termo quadrático, significando que seu acréscimo acaba correndo muito atrás dos dados.



3.3. Análise discriminante regularizada


Para o segundo conjunto de dados teremos:

fitRLDA2<- rda(Type~.,data = wine[kv,],gamma = 0,lambda = 1)
PredRLDA2<- predict(fitRLDA2,testeWine)



fitRQDA2<- rda(Type~.,data = wine[kv,],gamma = 0,lambda = 0)
PredRQDA2<- predict(fitRQDA2,testeWine)


fitREst2<- rda(Type~.,data = wine[kv,])
PredREst2<- predict(fitREst2,testeWine)



Tabela de contingência dos preditos versus observados:

tcRLDA2 <- table(testeWine$Type, PredRLDA2$class)
dimnames(tcRLDA2) <- list(Observados = as.character(unique(wine$Type)), "Predictos (cv)" = as.character(unique(wine$Type)))
print(round(tcRLDA2))
##           Predictos (cv)
## Observados  1  2  3
##          1 39  0  0
##          2  1 44  3
##          3  0  0 31
tcRQDA2 <- table(testeWine$Type, PredRQDA2$class)
dimnames(tcRQDA2) <- list(Observados = as.character(unique(wine$Type)), "Predictos (cv)" = as.character(unique(wine$Type)))
print(round(tcRQDA2))
##           Predictos (cv)
## Observados  1  2  3
##          1 19 20  0
##          2  1 47  0
##          3  0  7 24
tcREst2 <- table(testeWine$Type, PredREst2$class)
dimnames(tcREst2) <- list(Observados = as.character(unique(wine$Type)), "Predictos (cv)" = as.character(unique(wine$Type)))
print(round(tcREst2))
##           Predictos (cv)
## Observados  1  2  3
##          1 33  0  6
##          2  1 35 12
##          3  0 17 14

Proporção de acertos em cada grupo

round(prop.table(tcRLDA2, 1),3)
##           Predictos (cv)
## Observados     1     2     3
##          1 1.000 0.000 0.000
##          2 0.021 0.917 0.062
##          3 0.000 0.000 1.000
round(prop.table(tcRQDA2, 1),3)
##           Predictos (cv)
## Observados     1     2     3
##          1 0.487 0.513 0.000
##          2 0.021 0.979 0.000
##          3 0.000 0.226 0.774
round(prop.table(tcREst2, 1),3)
##           Predictos (cv)
## Observados     1     2     3
##          1 0.846 0.000 0.154
##          2 0.021 0.729 0.250
##          3 0.000 0.548 0.452

Proporção geral de acertos:

sum(diag(prop.table(tcRLDA2)))
## [1] 0.9661017
sum(diag(prop.table(tcRQDA2)))
## [1] 0.7627119
sum(diag(prop.table(tcREst2)))
## [1] 0.6949153

No segundo conjunto de dados conclui-se que quanto mais complexo o modelo, maior o erro de predição.

3. Sabatina

Utilizar o banco de dados Pima.tr do pacote MASS com treino e a base de dados Pima.te como teste. A com a resposta type e variáveis explicativas bmi e ped.

  1. Ajustar o modelo pela função de optim e comparar com ajuste do glm.

  2. Ajustar pelo método discriminante linear, e quadrático de Fisher

  3. Ajustar pelo método discriminante regularizado

  4. Obter as curvas ROC pelos métodos acima

  5. Interpretar os resultados