aula de hoje: regressão e classificação com o pacote caret

objetivos

Parte I

aspectos gerais da regressão

1.1. otimização de um modelo por descida do gradiente

1.2. modelos e generalização

1.3. regressão robusta

1.4. regressão com kernel

Parte II

o que é o caret

avaliação dos dados

pré-processamento dos dados

treinamento e avaliação

regressão: rf - random forest

regressão: knn - kth nearest neighbor

avaliação da importância relativa das variáveis

ajuste de hiper-parâmetros de modelos

classificação: preparação dos dados

classificação: SVM (Support Vector Machines)

classificação: árvore de decisão

exercícios

0 objetivos

parte I: aspectos gerais da regressão

na aula de hoje vou começar discutindo alguns aspectos gerais da regressão; inicio apresentando o algoritmo da descida do gradiente (gradient descent), que é o motor por trás da otimização dos modelos de ML; a seguir discuto alguns aspectos importantes de ML, como a generalização, e termino essa parte inicial com alguns modelos de regressão que podem ser úteis (regressão robusta e com kernel).

parte II: aprendizado de máquina com a biblioteca caret

meu objetivo aqui é apresentar uma receita de bolo para se resolver problemas de regressão e classificação com ML

vou ilustrar a regressão com a estimativa de redshifts fotométricos e depois a classificação entre estrelas e galáxias com dados puramente fotométricos extraídos dos catálogos do S-PLUS (Mendes de Oliveira et al. 2019)

note que para resolver estes problemas no “mundo real” usam-se diversos outros parâmetros estruturais dos objetos, mas estou simplificando para examinar o que se obtém apenas com a fotometria (além de reduzir o tempo de processamento)

1 aspectos gerais da regressão

# para monitorar o tempo de processamento deste script:
t00 = Sys.time()

# bibliotecas
suppressMessages(library(glmnet))
suppressMessages(library(olsrr))
suppressMessages(library(MASS))
suppressMessages(library(quantreg))
suppressMessages(require(robustbase))
suppressMessages(library(ggplot2))
suppressMessages(library(kernlab))
suppressMessages(library(caret))
suppressMessages(library(tidyverse))
suppressMessages(library(readxl))
suppressMessages(library(corrplot))
suppressMessages(library(ranger))
suppressMessages(library(dplyr))
suppressMessages(library(e1071))
suppressMessages(library(skimr))
suppressMessages(library(pROC))
suppressMessages(library("rpart"))
suppressMessages(library("rpart.plot"))
suppressMessages(library(KernSmooth))
suppressMessages(library(doParallel))

# paralelização
cl <- makePSOCKcluster(5)
registerDoParallel(cl)
# o uso de paralelização reduz o tempo de processamento deste script em 1/3

1.1 Otimização de um modelo com o algoritmo descendo o gradiente

Nosso objetivo aqui será ilustrar como funciona o algoritmo do gradient descent, ou descendo o gradiente, para minimizar a função de custo.

Vamos considerar como exemplo a regressão linear ordinária. Temos um conjunto de dados \(\cal D\) que queremos modelar como um modelo linear simples: \[y = w_1 +w_2 x + \epsilon\] onde \(\epsilon\) se refere aos erros em \(y\) e os parâmetros são \(\{w_1 ,w_2 \}\).

A função custo é proporcional ao log da verossimilhança e, supondo erros gaussianos, é dada por: \[l(w) \propto \chi^2 = \sum_{i=1}^N \frac{[y_i-(w_1 +w_2 x_i)]^2}{\epsilon^2}\] e os parâmetros são obtidos por minimização de \(l(w)\) (ou do \(\chi^2\)).

Com o algoritmo de descida do gradiente os parâmetros são obtidos iterativamente como \[w \longleftarrow w - \lambda \times \nabla l(w)\] onde \(0 < \lambda < 1\) é o coeficiente de aprendizado.

O gradiente de uma função em um dado ponto aponta para a direção, na vizinhança do ponto considerado, onde a função cresce mais rapidamente; o sinal negativo implica em ir na direção oposta, que minimiza a função.

Para ilustrar este procedimento vamos simular dados assumindo um modelo linear simples, \(y=1+2x\), com \(-1 < x < 1\) e erros gaussianos de desvio padrão 0.2:

# reprodutibilidade
set.seed(123)

# modelo linear:  y = w1 + w2 x + eps
# modelo linear simples: y=1+2x, com x entre -1 e 1 e com ruído gaussiano com sd = 0.2
# n: número de pontos simulados
n=100
eps=0.2
x <- runif(n, -1, 1)
y <- 1+2*x + rnorm(n,sd=eps)

Vamos, inicialmente, resolver o problema (isto é obter os parâmetros da reta) usando a função lm() do R:

modelo_linear <- lm( y ~ x )
print(modelo_linear)

## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##      0.9892       1.9910

# compare os valores dos parâmetros ajustados com os usados para simular os dados

# visualização:
plot(x,y, col=rgb(0.2,0.4,0.6,0.8), main='regressão linear')
# modelo
abline(1,2,col='red')
# ajuste com lm()
abline(modelo_linear, col='blue')
legend("topleft",legend=c('modelo','ajuste com lm()'),col=c('red','blue'),  text.col = c('red','blue'), bty='n')

Vamos agora ilustrar a minimização de \(l(w)\) com descendo o gradiente. Vamos usar notação matricial neste problema:

temos que o modelo dos dados é \[y_i = w_0+w_1 x_i+\epsilon_i,~~~~~~ i=1,...,n\]

seja \(\mathbf{Y}\) o vetor coluna (\(n \times 1\)) com a variável dependente: \[\begin{align} {\bf Y} &= \begin{bmatrix} y_{1} \\ y_{2} \\ \vdots \\ y_{n} \end{bmatrix} \end{align}\]

seja \({\mathbf X}\) a matriz (\(n \times 2\)) com a variável independende; a coluna com 1 é para considerar o intercepto \(w_0\) na multiplicação matricial (ver mais abaixo): \[\begin{align} {\bf X} &= \begin{bmatrix} 1 & x_{1} \\ 1 & x_{2} \\ \vdots \\ 1 & x_{n} \end{bmatrix} \end{align}\]

os parâmetros do modelo são representados pela matriz \(2 \times 1\)

\[\begin{align} {\bf W} &= \begin{bmatrix} w_{0} \\ w_{1} \end{bmatrix} \end{align}\]

e os erros (matriz \(n \times 1\)):

\[\begin{align} {\bf \epsilon} &= \begin{bmatrix} \epsilon_{1} \\ \epsilon_{2} \\ \vdots \\ \epsilon_{n} \end{bmatrix} \end{align}\]

nesse caso o modelo fica \({\mathbf Y} = {\mathbf X}{\mathbf W}+{\mathbf \epsilon}\), e os parâmetros podem ser estimados como \[{\mathbf W} = \Big({\mathbf X^T} {\mathbf X} \Big)^{-1}{\mathbf X}^T {\mathbf Y}\]

# taxa de aprendizado e número de iterações
lambda <- 0.01
num_iters <- 1000

# y: variável dependente
# ycalc = X %*% w

# função custo: erro quadrático
custo <- function(X, y, w) {
  sum( (X %*% w - y)^2 ) / (eps^2)
}

# vetores auxiliares para guardar a história da função de custo e dos parâmetros:
custo_historia <- rep(0,num_iters)
w_historia <- list(num_iters)

# inicializa os parâmetros (w1,w2)=(0,0)
w <- matrix(c(0,0), nrow=2)

# adiciona a coluna de '1's para os valores de x que multiplicam w1
X <- cbind(1, matrix(x))

# gradiente descendente:
for (i in 1:num_iters) {
  erro <- (X %*% w - y)
  delta <- t(X) %*% erro / length(y)
  w <- w - lambda * delta
  custo_historia[i] <- custo(X, y, w)
  w_historia[[i]] <- w
}

# resultado
print(w)

##           [,1]
## [1,] 0.9888444
## [2,] 1.9114226

# compare os valores dos parâmetros ajustados com os usados para simular os dados

# visualização do resultado:
par(mfrow=c(1,1))
plot(x,y, col=rgb(0.2,0.4,0.6,0.8), main='evolução do ajuste com "descendo o gradiente"')
for (i in c(1,3,6,10,14,seq(20,num_iters,by=10))) {
  abline(coef=w_historia[[i]], col=rgb(0.8,0,0,0.3))
}
abline(coef=w, col='blue')

# visualização da convergência:
plot(custo_historia, type='l', col='blue', lwd=2, main='convergência da função custo', ylab='custo', xlab='Iterações')

1.2 modelos e generalização

Para examinar como a complexidade afeta a generalização, vamos gerar um ‘conjunto de treinamento’ e um ‘conjunto de teste’ com o mesmo polinômio de segundo grau e vamos comparar o ajuste de polinômios de graus variados a esses dois conjuntos:

# reprodutibilidade
set.seed(123)

#modelo: y = 0.6-0.5*x+0.8*x^2
# x entre 0 e 1
x = seq(0,1,0.01)
y = 0.6-0.5*x+0.8*x^2

#conjunto de treinamento: dados com erros de 0.02
ntreino = 10
xtreino = runif(ntreino,0,1)
ytreino = 0.6-0.5*xtreino+0.8*xtreino^2+rnorm(ntreino,0,sd = 0.02)

# conjunto de teste
ntest = 10
xtest = runif(ntest,0,1)
ytest = 0.6-0.5*xtest+0.8*xtest^2+rnorm(ntest,0,sd = 0.02)

plot(x,y,col='blue',type='l',ylim=c(0,1), main='vermelho: treino;  verde: teste',lwd=2)
points(xtreino,ytreino,col='red',pch=20)
points(xtest,ytest,col='green',pch=20)

# vamos examinar modelos polinomiais de várias complexidades
# m: grau do polinômio
# custo: soma do quadrado dos resíduos
custo = rep(0,ntreino-1)
custo_test = rep(0,ntreino-1)
grau =  rep(0,ntreino-1)

for(i in 1:(ntreino-1)){
grau[i]=i
# ajuste de polinômio com lm() e poly()
yf=lm(ytreino ~ poly(xtreino, i, raw=TRUE))
model_treino =  rep(0,ntreino)
for(j in 1:(i+1)){model_treino = model_treino + yf$coefficients[j]*xtreino^(j-1)}
custo[i] = sum((ytreino-model_treino)^2)
model_test = rep(0,ntreino)
for(j in 1:(i+1)){model_test = model_test + yf$coefficients[j]*xtest^(j-1)}
custo_test[i] = sum((ytest-model_test)^2)
}

# custo do conjunto de treinamento: decresce com o grau do polinômio
custo

## [1] 4.669985e-02 2.136856e-03 2.106459e-03 1.564326e-03 6.450436e-04
## [6] 3.614415e-04 3.128014e-04 1.378176e-04 2.080280e-20

# custo do conjunto de teste: decresce e depois cresce
custo_test

## [1]  0.038768122  0.005321278  0.005350283  0.005131396  0.013097903
## [6]  0.027322413  0.187435975 12.365744952 91.947179410

plot(grau,custo,col='red',type='l',main='custo dos conjuntos de treinamento e teste')
points(grau,custo,col='red',pch=19)
lines(grau,custo_test,col='green')
points(grau,custo_test,col='green',pch=19)
legend("topright",legend=c("treino","teste"),pch=19,col=c("red","green"))

# polinômio de grau 9:
yf=lm(ytreino ~ poly(xtreino, 9, raw=TRUE))
xt = seq(0,1,0.001)
yt = rep(0,length(xt))
for(j in 1:10){yt =  yt + yf$coefficients[j]*xt^(j-1)}

plot(x,y,col='blue',type='l',ylim=c(0,1), main='preto: polinômio de grau 9',lwd=2)
points(xtreino,ytreino,col='red',pch=20)
points(xtest,ytest,col='green',pch=20)
lines(xt,yt,lwd=3)

É evidente que, embora o ajuste passe por todos os pontos do conjunto de treinamento, a generalização nesse caso é muito ruim!

1.3 regressão robusta

O R tem várias bibliotecas com métodos robustos, como a MASS, quantreg e robustbase

Os métodos de regressão robusta usam funções de custo menos sensível a outliers

Exemplo: função de Huber \[ L(a) = L_\delta(a) = \left\{ \begin{array}{lr} \frac{1}{2}a^2 & {\rm se}~ |a| \le \delta\\ \delta(|a|- \frac{1}{2}\delta) & {\rm se}~ |a| > \delta \end{array} \right. \] onde \(a = y_i -\mathbf x_i.\mathbf w\) é o resíduo de uma medida

Esta função é quadrática para valores pequenos de \(a\) e linear para valores grandes

Vamos dar uma examinada nos vários métodos usando uma simulação de dados:

# https://rpubs.com/dvallslanaquera/robust_regression
# reprodutibilidade
set.seed(123)

# dados:
n=100
x <- runif(n, -1, 1)
y <- 1+2*x + rnorm(n,sd=0.2)

# introduzindo outliers:
nout = 20
for(i in 1:nout){
x[n+i] = runif(1, -0.8, 0) 
y[n+i] = 2.5 + rnorm(1,mean=2,sd=0.5)
}

# modelo linear
fitLS <- lm(y~x)
summary(fitLS)$r.squared

## [1] 0.1628212

p = predict(fitLS, newx = x)
sd_LS = sd(p-y)
mad_LS = mad(p-y)

# Outliers usando a figura de Cook
# a distância de Cook mede o efeito de se remover um certo ponto da análise e, assim, é útil para identificar outliers
ols_plot_cooksd_bar(fitLS)

# Huber loss
fitH <- rlm(y~x, k2 = 1.345) 
p = predict(fitH, newx = x)
sd_H = sd(p-y)
mad_H = mad(p-y)

# LQS - faz um ajuste aos pontos "bons", sem os outliers
fitLQS <- lqs(y~x, method = "lms") 
p = predict(fitLQS, newx = x)
sd_LQS = sd(p-y)
mad_LQS = mad(p-y)

# LTS: least trimmed squares - minimiza a soma dos quadrados dos resíduos para uma subamostra dos dados
fitLTS <- lqs(y~x, method = "lts") 
p = predict(fitLTS, newx = x)
sd_LTS = sd(p-y)
mad_LTS = mad(p-y)

# LAD: least absolute devition
fitLAD <- rq(y~x, tau = 0.5)
p = predict(fitLAD, newx = x)
sd_LAD = sd(p-y)
mad_LAD = mad(p-y)

# S-estimator
fitS <- lqs(y~x, method = "S")
p = predict(fitS, newx = x)
sd_S = sd(p-y)
mad_S = mad(p-y)

# MM-estimator: combinação dos estimadores M e S
fitMM <- rlm(y~x, method = "MM")
p = predict(fitMM, newx = x)
sd_MM = sd(p-y)
mad_MM = mad(p-y)

# desvio padrão e mínimo desvio absoluto dos vários estimadores
c(sd_LS,sd_H,sd_LQS,sd_LTS,sd_LAD,sd_S,sd_MM)

## [1] 1.524960 1.567875 1.589155 1.582387 1.568842 1.579226 1.579755

c(mad_LS,mad_H,mad_LQS,mad_LTS,mad_LAD,mad_S,mad_MM)

## [1] 0.6686189 0.2624243 0.2457071 0.2476601 0.2615887 0.2406093 0.2419468

# note como o MAD é muito menor para os métodos robustos quando comparado aos da regressão linear ordinária (LS)

# visualização
plot(x,y,
      xlab = "x", ylab = "y", type = "p", 
      pch = 20, cex = .8)
abline(fitLS, col = 1) 
abline(fitH, col = 2) 
abline(fitLAD, col = 3) 
abline(fitLTS, col = 4) 
abline(fitLQS, col = 5) 
abline(fitS, col = 6) 
abline(fitMM, col = 7) 
legend('topright', c("LS", "Huber", "LAD","LTS","LQS",
                  "S-estimador","MM-estimador" ),
           lty = rep(1, 7), bty = "n",
           col = c(1, 2, 3, 4, 5, 6, 7))

Note que, no caso, os métodos “robustos” provêm um ajuste muito melhor que a regressão linear usual (LS), evitando os outliers.

1.4 regressão com kernel

este tipo de regressão é muito usado em interpolação

# Nonparametric Regression and Cross-Validation, Yen-Chi Chen
# reprodutibilidade
set.seed(123)

# simulação dos dados
X = runif(200, min=0, max=4*pi)
Y = sin(X) + rnorm(200, sd=0.3)

plot(X,Y, pch=20)

# vamos usar a função ksmooth para esta regressão
# dois parâmetros principais: o tipo de kernel e a largura de banda (h)
# kernel gaussiano
Kreg = ksmooth(x=X,y=Y,kernel = "normal",bandwidth = 0.5)
plot(X,Y,pch=20)
lines(Kreg, lwd=4, col="purple")

# efeito de variar h
Kreg1 = ksmooth(x=X,y=Y,kernel = "normal",bandwidth =0.1)
Kreg2 = ksmooth(x=X,y=Y,kernel = "normal",bandwidth =0.9)
Kreg3 = ksmooth(x=X,y=Y,kernel = "normal",bandwidth =3.0)
plot(X,Y,pch=20)
lines(Kreg1, lwd=4, col="orange")
lines(Kreg2, lwd=4, col="purple")
lines(Kreg3, lwd=4, col="limegreen")
legend("topright", c("h=0.1","h=0.9","h=3.0"), lwd=6,col=c("orange","purple","limegreen"))

# escolha de h por validação cruzada
#exemplo com 5-CV:
n = length(X)
N_cv = 100
k = 5
cv_lab = sample(n,n,replace=F) %% k
## randomly split all the indices into k numbers
h_seq = seq(from=0.2,to=2.0, by=0.1)
CV_err_h = rep(0,length(h_seq))
for(i_tmp in 1:N_cv){
CV_err_h_tmp = rep(0, length(h_seq))
cv_lab = sample(n,n,replace=F) %% k
for(i in 1:length(h_seq)){
h0 = h_seq[i]
CV_err =0
for(i_cv in 1:k){
w_val = which(cv_lab==(i_cv-1))
X_tr = X[-w_val]
Y_tr = Y[-w_val]
X_val = X[w_val]
Y_val = Y[w_val]
kernel_reg = ksmooth(x = X_tr,y=Y_tr,kernel = "normal",bandwidth=h0,
x.points=X_val)
# WARNING! The ksmooth() function will order the x.points from
# the smallest to the largest!
CV_err = CV_err+mean((Y_val[order(X_val)]-kernel_reg$y)^2,na.rm=T)
# na.rm = T: remove the case of 'NA'
}
CV_err_h_tmp[i] = CV_err/k
}
CV_err_h = CV_err_h+CV_err_h_tmp
}
CV_err_h = CV_err_h/N_cv
plot(h_seq,CV_err_h, type="b", lwd=4, col="blue", xlab="Smoothing Bandwidth",
ylab="5-CV Error")

h_opt = h_seq[which(CV_err_h==min(CV_err_h))]
h_opt

## [1] 0.9

# melhor ajuste
# kernel normal
Kreg = ksmooth(x=X,y=Y,kernel = "normal",bandwidth = h_opt)
plot(X,Y,pch=20)
lines(Kreg, lwd=4, col="orange")

2 o que é o `caret`

caret: Classification And REgression Training

caret - https://cran.r-project.org/web/packages/caret/vignettes/caret.html

http://topepo.github.io/caret/index.html

O caret é uma das bibliotecas mais populares para ML. Tem um conjunto de funções pré-definidas que agilizam a implementação de modelos de ML, como pré-processamento dos dados, seleção de atributos, ajuste de hiperparâmetros, etc, para um número muito grande de modelos.

Existem várias outras bibliotecas do R úteis para ML:

mlr3 - https://mlr3.mlr-org.com/

também provê os blocos essenciais para implementação de projetos de ML

rpart - https://www.rdocumentation.org/packages/rpart/versions/4.1-15

classificação e regressão com códigos de árvore

randomForest - https://cran.r-project.org/web/packages/randomForest/index.html

classificação e regressão com randomForest

e1071 - https://cran.r-project.org/web/packages/e1071/index.html

implementa vários algoritmos, como support vector machines

keras - https://tensorflow.rstudio.com/guide/keras/

deep learning com TensorFlow

torch: https://torch.mlverse.org/

também para deep learning (inclui modelos difusivos)

e muitos outros!

Vamos ilustrar o caret com um problema de regressão: estimativa de redshifts fotométricos para galáxias do S-PLUS
Para saber mais sobre o S-PLUS veja Mendes de Oliveira et al. (2019; arXiv:1907.01567).

3 avaliação dos dados

a. leitura dos dados: vamos considerar uma amostra de galáxias do S-PLUS com fotometria nas 12 bandas

o arquivo splus-mag-z.dat contém, para cada objeto da amostra, as 12 magnitudes e o redshift espectroscópico (do SDSS)

o objetivo será estimar o redshift usando apenas as magnitudes (redshift fotométrico)

vamos selecionar apenas galáxias com magnitudes no intervalo 15 \(<\) r_petro \(<\) 20.

dados = as.data.frame(read.table('splus-mag-z.dat', sep = "", header=TRUE))
# dimensão dos dados:
dim(dados)

## [1] 55803    13

# seleção em magnitudes
sel = rep(0, nrow(dados))
sel[dados$r_petro > 15 & dados$r_petro < 20] = 1
sum(sel)

## [1] 54313

dados = dados[sel == 1,]
dim(dados)

## [1] 54313    13

# topo do arquivo
head(dados)

##   uJAVA_petro F378_petro F395_petro F410_petro F430_petro g_petro F515_petro
## 1       18.71      18.69      18.50      18.12      17.78   17.43      17.09
## 2       18.26      18.12      18.28      17.32      16.95   16.67      16.24
## 3       19.38      19.45      18.76      18.60      18.50   18.07      17.69
## 4       20.30      19.73      20.01      19.18      18.91   18.36      18.00
## 5       19.49      18.98      18.64      18.64      18.02   17.14      16.61
## 6       19.86      19.38      18.93      18.75      18.15   17.44      16.94
##   r_petro F660_petro i_petro F861_petro z_petro z_SDSS
## 1   16.89      16.83   16.56      16.54   16.45  0.111
## 2   15.97      15.90   15.61      15.52   15.41  0.082
## 3   17.55      17.47   17.25      17.31   17.18  0.086
## 4   17.75      17.68   17.51      17.46   17.34  0.111
## 5   16.06      15.93   15.65      15.54   15.38  0.162
## 6   16.58      16.48   16.24      16.11   16.00  0.082

é sempre bom ver se há dados faltantantes ou outros problemas com os dados
a função skim() provê um sumário da estatística descritiva de cada variável:

sumario <- skim(dados)
sumario

Data summary
Name	dados
Number of rows	54313
Number of columns	13
_______________________
Column type frequency:
numeric	13
________________________
Group variables

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
uJAVA_petro	1	20.10	1.21	14.95	19.36	20.05	20.76	31.02	▁▇▁▁▁
F378_petro	1	19.85	1.17	12.94	19.15	19.81	20.50	28.69	▁▃▇▁▁
F395_petro	1	19.69	1.22	12.72	18.96	19.66	20.37	29.16	▁▅▇▁▁
F410_petro	1	19.50	1.20	11.89	18.79	19.48	20.19	28.37	▁▂▇▁▁
F430_petro	1	19.25	1.19	11.66	18.54	19.25	19.97	29.13	▁▃▇▁▁
g_petro	1	18.72	1.06	15.10	18.05	18.73	19.45	24.35	▁▇▇▁▁
F515_petro	1	18.42	1.08	10.78	17.74	18.43	19.16	25.06	▁▁▇▂▁
r_petro	1	17.97	1.02	15.01	17.30	18.01	18.76	19.99	▁▃▇▇▅
F660_petro	1	17.84	1.03	10.70	17.16	17.89	18.65	21.11	▁▁▂▇▂
i_petro	1	17.62	1.03	10.15	16.93	17.68	18.43	20.52	▁▁▂▇▃
F861_petro	1	17.46	1.05	10.39	16.75	17.52	18.29	21.00	▁▁▃▇▁
z_petro	1	17.42	1.06	9.83	16.70	17.48	18.26	20.96	▁▁▂▇▁
z_SDSS	1	0.15	0.09	0.00	0.08	0.13	0.19	1.12	▇▂▁▁▁

Pode-se ver que não há dados faltantes; se houvessem precisaríamos ou removê-los ou atribuir um valor para eles (imputação) pois os algoritmos que vamos usar precisam disso!

Como estamos interessados apenas em explorar o pacote e algumas de suas principais funções, vamos restringir o número de objetos:

nsample = 10000

# magnitudes
mags = dados[1:nsample,1:12]

# redshift
zspec = dados[1:nsample,13]

b. vamos dar uma olhada nos dados:

par(mfrow = c(1,2))
hist(mags$r_petro,xlab='r_petro',main='',col='red')
hist(zspec,xlab='z_SDSS',main='',col='blue')

c. vamos visualizar as correlações entre as várias magnitudes

as variáveis usadas na estimativa (as magnitudes) são denominadas atributos ou “features”

dados.cor = cor(mags)
corrplot(dados.cor, type = "upper", order = "original", tl.col = "black", tl.srt = 45)

redshift versus cada banda fotométrica:

visualização usando featurePlot:

featurePlot(x = mags, y = zspec, plot = "scatter")

4 pré-processamento dos dados

Em geral o desempenho dos algoritmos melhora se os dados tiverem um intervalo de valores “razoável”. Assim, é conveniente pré-processar os dados, isto é, transformá-los para colocá-los num intervalo para tornar a análise mais eficiente.

Há muitos tipos de pré-processamento; os mais comuns são

reescalonar as variáveis entre 0 e 1 (ou entre -1 e 1)

subtrair a média de cada variável e dividir por seu desvio padrão

vamos colocar as 12 magnitudes entre 0 e 1:

# valores máximos e mínimos de cada banda
maxs <- apply(mags, 2, max) 
mins <- apply(mags, 2, min)
# reescalonamento usando a função scale():
magnorm <- as.data.frame(scale(mags, center = mins, scale = maxs - mins))

# exemplo de dados normalizados
summary(magnorm[,1])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.2420  0.2869  0.2889  0.3317  1.0000

#se quiséssemos subtrair a média e dividir pelo desvio padrão:
#magnorm <- as.data.frame(scale(mags, center = TRUE, scale = TRUE))

# não é necessário normalizar z nesta amostra pois ele está num intervalo "razoável":
summary(zspec)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0040  0.0770  0.1150  0.1435  0.1900  0.9420

5 Treinamento e avaliação

Em ML o algoritmo aprende a partir dos dados.

Para avaliar quanto ele aprendeu, comparamos suas predições com resultados conhecidos (aprendizado supervisionado)

Para isso divide-se os dados em conjunto de treinamento, para determinar os parâmetros do modelo, e conjunto de teste, para se determinar a qualidade do resultado.

a. criação de conjuntos de treinamento e teste:

Vamos considerar 75% dos objetos para treinamento e 25% para teste

Note que os objetos no conjunto de treinamento serão usados também para validação

# número das linhas dos objetos selecionados para o conjunto de treinamento
# usamos a função createDataPartition()
numlinhastreinamento =  createDataPartition(zspec, p=0.75, list=FALSE)

# conjunto de treinamento:
xtrain=magnorm[numlinhastreinamento,]
ytrain=zspec[numlinhastreinamento]

# conjunto de teste
xtest=magnorm[-numlinhastreinamento,]
ytest=zspec[-numlinhastreinamento]

b. seleção do algoritmo:

O caret contém um número enorme de algoritmos de ML. Para ver quais são:

paste(names(getModelInfo()), collapse=',  ')

## [1] "ada,  AdaBag,  AdaBoost.M1,  adaboost,  amdai,  ANFIS,  avNNet,  awnb,  awtan,  bag,  bagEarth,  bagEarthGCV,  bagFDA,  bagFDAGCV,  bam,  bartMachine,  bayesglm,  binda,  blackboost,  blasso,  blassoAveraged,  bridge,  brnn,  BstLm,  bstSm,  bstTree,  C5.0,  C5.0Cost,  C5.0Rules,  C5.0Tree,  cforest,  chaid,  CSimca,  ctree,  ctree2,  cubist,  dda,  deepboost,  DENFIS,  dnn,  dwdLinear,  dwdPoly,  dwdRadial,  earth,  elm,  enet,  evtree,  extraTrees,  fda,  FH.GBML,  FIR.DM,  foba,  FRBCS.CHI,  FRBCS.W,  FS.HGD,  gam,  gamboost,  gamLoess,  gamSpline,  gaussprLinear,  gaussprPoly,  gaussprRadial,  gbm_h2o,  gbm,  gcvEarth,  GFS.FR.MOGUL,  GFS.LT.RS,  GFS.THRIFT,  glm.nb,  glm,  glmboost,  glmnet_h2o,  glmnet,  glmStepAIC,  gpls,  hda,  hdda,  hdrda,  HYFIS,  icr,  J48,  JRip,  kernelpls,  kknn,  knn,  krlsPoly,  krlsRadial,  lars,  lars2,  lasso,  lda,  lda2,  leapBackward,  leapForward,  leapSeq,  Linda,  lm,  lmStepAIC,  LMT,  loclda,  logicBag,  LogitBoost,  logreg,  lssvmLinear,  lssvmPoly,  lssvmRadial,  lvq,  M5,  M5Rules,  manb,  mda,  Mlda,  mlp,  mlpKerasDecay,  mlpKerasDecayCost,  mlpKerasDropout,  mlpKerasDropoutCost,  mlpML,  mlpSGD,  mlpWeightDecay,  mlpWeightDecayML,  monmlp,  msaenet,  multinom,  mxnet,  mxnetAdam,  naive_bayes,  nb,  nbDiscrete,  nbSearch,  neuralnet,  nnet,  nnls,  nodeHarvest,  null,  OneR,  ordinalNet,  ordinalRF,  ORFlog,  ORFpls,  ORFridge,  ORFsvm,  ownn,  pam,  parRF,  PART,  partDSA,  pcaNNet,  pcr,  pda,  pda2,  penalized,  PenalizedLDA,  plr,  pls,  plsRglm,  polr,  ppr,  pre,  PRIM,  protoclass,  qda,  QdaCov,  qrf,  qrnn,  randomGLM,  ranger,  rbf,  rbfDDA,  Rborist,  rda,  regLogistic,  relaxo,  rf,  rFerns,  RFlda,  rfRules,  ridge,  rlda,  rlm,  rmda,  rocc,  rotationForest,  rotationForestCp,  rpart,  rpart1SE,  rpart2,  rpartCost,  rpartScore,  rqlasso,  rqnc,  RRF,  RRFglobal,  rrlda,  RSimca,  rvmLinear,  rvmPoly,  rvmRadial,  SBC,  sda,  sdwd,  simpls,  SLAVE,  slda,  smda,  snn,  sparseLDA,  spikeslab,  spls,  stepLDA,  stepQDA,  superpc,  svmBoundrangeString,  svmExpoString,  svmLinear,  svmLinear2,  svmLinear3,  svmLinearWeights,  svmLinearWeights2,  svmPoly,  svmRadial,  svmRadialCost,  svmRadialSigma,  svmRadialWeights,  svmSpectrumString,  tan,  tanSearch,  treebag,  vbmpRadial,  vglmAdjCat,  vglmContRatio,  vglmCumulative,  widekernelpls,  WM,  wsrf,  xgbDART,  xgbLinear,  xgbTree,  xyf"

# você pode saber um pouco mais sobre um modelo com o comando `modelLookup:
modelLookup('lars')

##   model parameter    label forReg forClass probModel
## 1  lars  fraction Fraction   TRUE    FALSE     FALSE

# lars é uma variante da regressão linear
# tem um parâmetro, fraction, pode ser usado para regressão mas não para classificação, e não dá os resultados em probabilidades

6 regressão: rf - random forest

vamos agora usar random forest para estimar redshifts fotométricos a partir das magnitudes: \[ z_{phot} = f_{rf}(\mathbf m) \]

rf é um algoritimo poderoso, usado tanto em classificação como em regressão

rf consiste de um grande número de árvores de decisão que funcionam como um ensamble

árvores de decisão funcionam por partições sucessivas de um conjunto de dados

numa rf, cada árvore de decisão usa apenas um subconjunto dos parâmetros de entrada (isto é, um subconjunto das 12 magnitudes, neste caso)

isso gera um número muito grande de modelos que, juntos, têm um valor preditivo maior que os de cada modelo individual

o resultado é uma combinação dos resultados dos modelos individuais

# reprodutibilidade
set.seed(123)

# modelo:
modelLookup('rf')

##   model parameter                         label forReg forClass probModel
## 1    rf      mtry #Randomly Selected Predictors   TRUE     TRUE      TRUE

# rf pode ser usada em regressão e classificação e seu parâmetro é o tamanho do subconjunto de magnitudes que vai ser usado; vamos deixar ele livre para o próprio algoritmo escolher

# treinando o modelo- vejam a sintaxe do comando
# pode-se mudar o modelo mudando apenas o parâmetro method
# note que o conjunto de teste não participa do treinamento
t0 = Sys.time()

# note como se treina modelos com o caret: com train()
model_rf<-train(xtrain,ytrain,method='rf')

# tempo que levou para rodar o modelo:
Sys.time() - t0

## Time difference of 23.48877 mins

# sumário do modelo
print(model_rf)

## Random Forest 
## 
## 7502 samples
##   12 predictor
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 7502, 7502, 7502, 7502, 7502, 7502, ... 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE        Rsquared   MAE       
##    2    0.05272715  0.6781133  0.03684709
##    7    0.05179765  0.6868900  0.03559428
##   12    0.05226953  0.6808656  0.03573003
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 7.

plot(model_rf)

# o melhor subconjunto tem 7 magnitudes
# a predição será feita com um conjunto de árvores com 7 magnitudes escolhidas aleatoriamente entre as 12

# predição com o conjunto de teste
pred_rf = predict(model_rf, xtest)

# estatística da performance do algoritmo
# sigma equivalente da gaussiana
relz = (ytest - pred_rf)
sig_G = 0.7413*(quantile(relz,0.75,names = FALSE) - quantile(relz,0.25,names = FALSE))
sig_G

## [1] 0.0349661

#visualização do resultado
# vamos juntar os valores preditos e observados do redshift em um fata frame:
my_data = as.data.frame(cbind(observed = ytest,predicted = pred_rf))
head(my_data)

##    observed predicted
## 8     0.159 0.1487884
## 11    0.127 0.1246000
## 16    0.057 0.0687450
## 17    0.103 0.1319060
## 24    0.186 0.1294244
## 29    0.194 0.1515569

# figura com ggplot
ggplot(my_data,aes(predicted, observed)) +
      geom_point(color = "darkred", alpha = 0.5) + 
      geom_smooth(method=lm)+ 
      ggtitle("RF: redshift predito x observado") +
      xlab("z_spec") +
      ylab("z_phot") +
      theme(plot.title = element_text(color="darkgreen",size=18,hjust = 0.5),
                     axis.text.y = element_text(size=12),
            axis.text.x = element_text(size=12,hjust=.5),
                      axis.title.x = element_text(size=14),
                      axis.title.y = element_text(size=14))

## `geom_smooth()` using formula 'y ~ x'

7 regressão: knn - kth nearest neighbor

knn, ou o k-ésimo vizinho mais próximo, faz uma inferência do redshift de uma galáxia no conjunto de teste a partir do redshift das galáxias com magnitudes mais “próximas” a ela no conjunto de treinamento

o parâmetro livre importante é o k, o número de vizinhos a considerar

o valor predito pode ser, por exemplo, o redshift médio desses k vizinhos mais próximos

no exemplo a seguir, k é determinado por validação cruzada, usando-se 5 estimativas para cada valor de k, com k variando de 1 a 41, de 2 em 2

preste atenção na estrutura da função train():

# reprodutibilidade
set.seed(123)

# modelo:
modelLookup('knn')

##   model parameter      label forReg forClass probModel
## 1   knn         k #Neighbors   TRUE     TRUE      TRUE

# na função train(), o trControl
t0 = Sys.time()
model_knn <- train(xtrain,ytrain,method='knn',
  trControl = trainControl(method = "cv", number = 5),
  tuneGrid = expand.grid(k = seq(1, 41, by = 2)))
# tempo que levou para rodar o modelo:
Sys.time()-t0

## Time difference of 4.197748 secs

# compare esse tempo com o algoritmo de RF

# resultado
model_knn

## k-Nearest Neighbors 
## 
## 7502 samples
##   12 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 6001, 6002, 6000, 6003, 6002 
## Resampling results across tuning parameters:
## 
##   k   RMSE        Rsquared   MAE       
##    1  0.07160195  0.4687566  0.04863549
##    3  0.06001804  0.5886704  0.04104289
##    5  0.05801333  0.6147965  0.03995158
##    7  0.05714793  0.6275943  0.03940757
##    9  0.05688358  0.6328385  0.03919225
##   11  0.05673935  0.6362021  0.03922058
##   13  0.05674291  0.6369481  0.03929043
##   15  0.05673483  0.6383155  0.03928524
##   17  0.05688254  0.6368559  0.03940698
##   19  0.05683411  0.6381474  0.03944134
##   21  0.05689367  0.6382290  0.03960099
##   23  0.05701331  0.6375593  0.03965783
##   25  0.05714856  0.6362589  0.03976003
##   27  0.05731816  0.6346990  0.03991523
##   29  0.05739302  0.6340868  0.04000369
##   31  0.05748445  0.6332368  0.04009385
##   33  0.05763726  0.6317552  0.04023181
##   35  0.05779370  0.6298680  0.04035755
##   37  0.05793343  0.6283185  0.04046641
##   39  0.05804766  0.6271098  0.04056960
##   41  0.05820836  0.6252203  0.04069799
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was k = 15.

# veja como a performance do algoritmo varia para diferentes valores de k
plot(model_knn)

# predição com o conjunto de teste
pred_knn = predict(model_knn, xtest)

relz = (ytest - pred_knn)
sig_G = 0.7413*(quantile(relz,0.75,names = FALSE) - quantile(relz,0.25,names = FALSE))
sig_G

## [1] 0.04091976

#
# o desempenho foi pior que o usando rf
#
#visualização do resultado
my_data = as.data.frame(cbind(observed = ytest,
                            predicted = pred_knn))

ggplot(my_data,aes(predicted, observed)) +
      geom_point(color = "darkred", alpha = 0.5) + 
      geom_smooth(method=lm)+ 
      ggtitle("knn: z_spec x z_phot") +
      ylab("z_phot ") +
      xlab("z_spec") +
      theme(plot.title = element_text(color="darkgreen",size=18,hjust = 0.5),
                     axis.text.y = element_text(size=12),
            axis.text.x = element_text(size=12,hjust=.5),
                      axis.title.x = element_text(size=14),
                      axis.title.y = element_text(size=14))

## `geom_smooth()` using formula 'y ~ x'

# comparação dos dois resultados:
my_data = as.data.frame(cbind(pred1 = pred_knn,
                            pred2 = pred_rf))

ggplot(my_data,aes(pred1,pred2)) +
      geom_point(color = "darkred", alpha = 0.5) + 
      geom_smooth(method=lm)+ 
      ggtitle("knn: predições knn & rf ") +
      xlab("z_knn ") +
      ylab("z_rf") +
      theme(plot.title = element_text(color="darkgreen",size=18,hjust = 0.5),
                     axis.text.y = element_text(size=12),
            axis.text.x = element_text(size=12,hjust=.5),
                      axis.title.x = element_text(size=14),
                      axis.title.y = element_text(size=14))

## `geom_smooth()` using formula 'y ~ x'

8 avaliação da importância relativa das variáveis

A determinação da importância relativa das variáveis é denominada feature selection: quais bandas fotométricas são mais relevantes para determinar o redshift?

Vamos considerar o problema de feature selection usando rfe (recursive feature elimination). Este procedimento tem 3 passos:

usa-se parte dos dados como conjunto de treinamento para construir um modelo de ML e estimar a importância relativa de cada variável num conjunto de teste, formado pelo resto dos dados

usa-se um procedimento iterativo para construir modelos com subconjuntos das variáveis, mantendo as mais importantes e ordenando por importância cada variável a cada iteração

seleciona-se os preditores ótimos comparando o desempenho de diferentes conjuntos de variáveis

# reprodutibilidade
set.seed(123)

# além de x e y, rfe tem dois parâmetros importantes: 
# sizes: o núnero de variáveis a ser considerado em cada modelo; no caso vamos considerar modelos com 1, 6 e 12 variáveis
# rfeControl determina o tipo de algoritmo e o método de validação cruzada a ser usado; no caso "random forest"  com 5 validações cruzadas
subsets <- c(6, 9, 12)

ctrl <- rfeControl(functions = rfFuncs,
                   method = "repeatedcv",
                   repeats = 5,
                   verbose = FALSE)
t0 = Sys.time()
importancia <- rfe(x=mags, y=zspec,
                 sizes = subsets,
                 rfeControl = ctrl)
# tempo que levou para rodar o comando
Sys.time() - t0

## Time difference of 39.86736 mins

# pode-se saber mais sobre as opções de cada uma dessas funções fazendo ?nome_da_função

# a saída contém a acurácia, o parâmetro kappa (útil em classificação) e seus desvios padrão para os diferentes modelos. O subconjunto final de variáveis é marcado com uma * na coluna a direita

importancia

## 
## Recursive feature selection
## 
## Outer resampling method: Cross-Validated (10 fold, repeated 5 times) 
## 
## Resampling performance over subset size:
## 
##  Variables    RMSE Rsquared     MAE   RMSESD RsquaredSD    MAESD Selected
##          6 0.05439   0.6588 0.03799 0.003250    0.03837 0.001932         
##          9 0.05183   0.6910 0.03570 0.002829    0.02730 0.001151         
##         12 0.05081   0.7037 0.03485 0.002803    0.02699 0.001126        *
## 
## The top 5 variables (out of 12):
##    g_petro, F430_petro, F395_petro, z_petro, F515_petro

# vejam quais são, nesse caso, as 5 magnitudes mais importantes

9 ajuste de hiperparâmetros de modelos

Vamos ilustrar o ajuste de hiperparâmetros com o método de validação cruzada aplicado ao algoritmo ‘gbm’, Generalized Boosted Regression Modeling.

GBMs combinam duas técnicas: árvores de decisão e métodos de boosting

boosting refere-se a uma família de algoritmos iterativos que convertem “weak learners” (como árvores de decisão) em “strong learners”

o GBM é construído sequencialmente: para cada nova árvore, o modelo considera os erros da última árvore, e a árvore sucessiva é construída sobre os erros cometidos pela árvore anterior, usando um algoritmo tipo descendo o gradiente

Para saber quais são os hiperparâmetros ajustáveis de um modelo :

modelLookup(model='gbm')

##   model         parameter                   label forReg forClass probModel
## 1   gbm           n.trees   # Boosting Iterations   TRUE     TRUE      TRUE
## 2   gbm interaction.depth          Max Tree Depth   TRUE     TRUE      TRUE
## 3   gbm         shrinkage               Shrinkage   TRUE     TRUE      TRUE
## 4   gbm    n.minobsinnode Min. Terminal Node Size   TRUE     TRUE      TRUE

# esses hiperparâmetros são, respectivamente, o número de árvores, o número máximo de nós por árvore, um coeficiente de aprendizagem e o número mínimo de observações nos nós terminais da árvore

Vamos fazer a validação cruzada 3 vezes (divide-se os dados em treinamento e validação 3 vezes) e repetir isso também 3 vezes

fitControl <- trainControl(
  method = "repeatedcv",
  number = 3,
  repeats = 3)

vamos criar uma grade de hiper-parâmetros:

grade = expand.grid(n.trees=c(10,100,1000),shrinkage=c(0.01,0.05,0.1,0.5),n.minobsinnode = c(3,5,10),interaction.depth=c(1,5,10))

treinamos o modelo:

# reprodutibilidade
set.seed(123)

t0 = Sys.time()
model_gbm<-train(xtrain,ytrain,method='gbm',trControl=fitControl,tuneGrid=grade,verbose = FALSE)
Sys.time() - t0

## Time difference of 6.258914 mins

summary(model_gbm)

##                     var   rel.inf
## g_petro         g_petro 42.289841
## F430_petro   F430_petro 12.887905
## z_petro         z_petro  5.563788
## F861_petro   F861_petro  5.496348
## F515_petro   F515_petro  5.409431
## F378_petro   F378_petro  5.162054
## i_petro         i_petro  4.921925
## uJAVA_petro uJAVA_petro  4.584888
## F410_petro   F410_petro  3.838165
## F660_petro   F660_petro  3.410508
## F395_petro   F395_petro  3.234831
## r_petro         r_petro  3.200316

# predição com o conjunto de teste
pred = predict(model_gbm, xtest)

relz = (ytest - pred)
sig_G = 0.7413*(quantile(relz,0.75,names = FALSE) - quantile(relz,0.25,names = FALSE))
sig_G

## [1] 0.03614058

#visualização do resultado
my_data = as.data.frame(cbind(observed = ytest,
                            predicted = pred))

ggplot(my_data,aes(predicted, observed)) +
      geom_point(color = "darkred", alpha = 0.5) + 
      geom_smooth(method=lm)+ 
      ggtitle("GBM: z_spec vs z_phot") +
      ylab("z_phot ") +
      xlab("z_spec") +
      theme(plot.title = element_text(color="darkgreen",size=18,hjust = 0.5),
                     axis.text.y = element_text(size=12),
            axis.text.x = element_text(size=12,hjust=.5),
                      axis.title.x = element_text(size=14),
                      axis.title.y = element_text(size=14))

## `geom_smooth()` using formula 'y ~ x'

Vocês podem comparar os resultados desse modelo com rf e knn.

10 classificação: preparação dos dados

Vamos agora considerar problemas de classificação, com um conjunto de dados puramente fotométricos classificados como estrelas ou galáxias.

Para ilustração, vamos usar 10000 objetos classificados como estrelas (classe 0) ou galáxias (classe 1):

#tabela = as.data.frame(read.table(file="class_estr_gal.dat", header=TRUE))
tabela = as.data.frame(read.table(file="dadosparaestrelagalaxia.dat", header=TRUE))

# alguns detalhes dos dados
dim(tabela)

## [1] 10000    13

head(tabela)

##    u_petro J0378_petro J0395_petro J0410_petro J0430_petro  g_petro J0515_petro
## 1 18.57333    17.99002    18.02439    18.09809    17.92609 17.28233    16.99983
## 2 17.14259    16.78725    16.52434    16.22606    15.83635 15.33351    15.02085
## 3 18.60449    18.16562    17.87111    17.58500    17.17834 16.74320    16.53280
## 4 16.74094    16.38911    16.28852    15.96590    15.66112 15.34018    15.18092
## 5 18.92666    18.72918    18.54957    17.53203    17.32281 16.68435    16.57145
## 6 17.54627    17.25873    17.20477    16.65028    16.57487 16.22518    16.05884
##    r_petro J0660_petro  i_petro J0861_petro  z_petro classe
## 1 16.53974    16.47142 16.21496    16.00299 15.91366      1
## 2 14.59325    14.48909 14.15851    13.96557 13.90085      1
## 3 16.07838    15.98695 15.69374    15.49386 15.44834      1
## 4 14.70637    14.63379 14.39819    14.20232 14.10782      1
## 5 15.87334    15.76929 15.59096    15.50189 15.46906      0
## 6 15.74276    15.70845 15.61446    15.59427 15.56725      0

# número de objetos em cada classe:
# (em classificação é sempre útil se trabalhar com conjuntos de treinamento com o mesmo número de objetos por classe)
length(tabela$classe[tabela$classe == 0])

## [1] 5000

length(tabela$classe[tabela$classe == 1])

## [1] 5000

Vamos examinar os dados em um diagrama cor-magnitude:

par(mfrow = c(1,3))
ca = tabela$g_petro-tabela$r_petro
cb = tabela$r_petro-tabela$i_petro
smoothScatter(ca,cb,nrpoints=0,add=FALSE,xlab="g-r",ylab="r-i",xlim=c(-1,3),ylim=c(-1,3),main='todos')
ca = tabela$g_petro[tabela$classe == 0]-tabela$r_petro[tabela$classe == 0]
cb = tabela$r_petro[tabela$classe == 0]-tabela$i_petro[tabela$classe == 0]
smoothScatter(ca,cb,nrpoints=0,add=FALSE,xlab="g-r",ylab="r-i",xlim=c(-1,3),ylim=c(-1,3),main='classe = 0')
ca = tabela$g_petro[tabela$classe == 1]-tabela$r_petro[tabela$classe == 1]
cb = tabela$r_petro[tabela$classe == 1]-tabela$i_petro[tabela$classe == 1]
smoothScatter(ca,cb,nrpoints=0,add=FALSE,xlab="g-r",ylab="r-i",xlim=c(-1,3),ylim=c(-1,3),main='classe = 1')

Vamos definir os ‘features’ (as magnitudes) e pré-processá-los:

# definindo os features
mag = tabela[,1:12]

# pré-processamento: normalização das magnitudes entre 0 e 1
maxs <- apply(mag, 2, max) 
mins <- apply(mag, 2, min)
normalizacao <- as.data.frame(scale(mag, center = mins, scale = maxs - mins))

Vamos definir os conjuntos de treinamento e teste

o conjunto de treinamento será usado para treinamento e validação; o conjunto de teste não entra no treinamento e serve para aferir o desempenho do modelo

# conjuntos de treinamento (75%) e teste (25%):
ntreino = round(0.75*nrow(tabela))
nteste = nrow(tabela)-ntreino
c(ntreino,nteste)

## [1] 7500 2500

# número das linhas dos objetos selecionados para o conjunto de treinamento
set.seed(123)
indice =  createDataPartition(tabela$classe, p=0.75, list=FALSE)
xtrain = normalizacao[indice,]

# ATENÇÃO: para classificação a variável dependente deve ser tipo 'factor'
ytrain = as.factor(tabela[indice,13])

xtest = normalizacao[-indice,]
ytest = as.factor(tabela[-indice,13])

11 classificação com SVM

SVM, Support Vector Machines, é um algoritmo bastante poderoso e bastante usado em problemas de classificação.

Para ver como o algoritmo funciona, considere a figura abaixo: temos duas classes (vermelho e azul) que queremos separar. SVM ajusta um hiper-plano (uma reta nesse caso) que maximiza a margem entre as classes:

https://www.datacamp.com/community/tutorials/support-vector-machines-r

Assim, objetos de um lado da linha do melhor hiper-plano são classificados em uma classe e os do outro lado da linha na outra classe.

Neste exemplo, as classes são linearmente separáveis, isto é, uma linha (ou hiper-plano), provê uma separação adequada entre as classes.

Considere agora o espaço de dados ilustrado no lado esquerdo da figura abaixo: ele não é linearmente separável- não há como traçar uma reta neste espaço que separe os pontos azuis dos vermelhos.

Mas os pontos azuis e vermelhos estão claramente segregados! Um jeito de torná-los linearmente separáveis é introduzir uma terceira dimensão: \[z = x^2 + y^2.\]

No caso do nosso exemplo isso torna o espaço de dados linearmente separável, como ilustrado no lado direito da figura:

https://www.datacamp.com/community/tutorials/support-vector-machines-r

mapeando este hiper-plano de volta no espaço de dados original obtemos um círculo:

https://www.datacamp.com/community/tutorials/support-vector-machines-r

Mas como SVM faz isso? Pode-se mostrar que a equação que precisa ser minimizada (num espaço de n dimensões) não depende da posição dos pontos mas apenas de seu produto interno (\(\mathbf{x}_i . \mathbf{x}_j\)), que é a distância entre dois pontos no espaço de dados.

Assim, se quisermos transformar os dados para um espaço de maior dimensão, não precisamos calcular a transformação exata dos dados e só precisamos do produto interno dos dados nesse espaço.

kernels permitem isso: por exemplo, um kernel gaussiano mapeia os dados num espaço de dimensão infinita!

Isso é chamado de kernel trick: aumenta-se o espaço de dados para permitir uma separação das classes.

Aqui vamos usar um kernel chamado Radial Basis Function (RBF) \[ K(\mathbf{x}_i,\mathbf{x}_j) = \exp \Bigg[ -\gamma \sum_k^{dim} (x_{ik}-x_{jk})^2 \Bigg]\]

No caret temos o algoritmo svmRadial:

modelLookup('svmRadial')

##       model parameter label forReg forClass probModel
## 1 svmRadial     sigma Sigma   TRUE     TRUE      TRUE
## 2 svmRadial         C  Cost   TRUE     TRUE      TRUE

Este modelo tem dois parâmetros ajustáveis: sigma e C

sigma controla a “suavidade” do limite de decisão do modelo

C penaliza o modelo por classificações erradas: quanto maior C, menor a probabilidade de um erro.

O pacote otimiza esses parâmetros maximizando a acurácia, no caso por validação cruzada:

# note que a sintaxe d a função train() é a mesma em problemas de classificação e regressão
t0 = Sys.time()
model_svm <- train(xtrain,ytrain,method='svmRadial',
                      tuneLength = 10,
                      trControl = trainControl(method = "cv"))
Sys.time()-t0

## Time difference of 25.73528 secs

# modelo ajustado:
model_svm

## Support Vector Machines with Radial Basis Function Kernel 
## 
## 7500 samples
##   12 predictor
##    2 classes: '0', '1' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 6750, 6750, 6750, 6750, 6750, 6750, ... 
## Resampling results across tuning parameters:
## 
##   C       Accuracy   Kappa    
##     0.25  0.9624000  0.9248000
##     0.50  0.9665333  0.9330667
##     1.00  0.9694667  0.9389333
##     2.00  0.9709333  0.9418667
##     4.00  0.9725333  0.9450667
##     8.00  0.9730667  0.9461333
##    16.00  0.9736000  0.9472000
##    32.00  0.9724000  0.9448000
##    64.00  0.9716000  0.9432000
##   128.00  0.9706667  0.9413333
## 
## Tuning parameter 'sigma' was held constant at a value of 0.2589718
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 0.2589718 and C = 16.

predição com o conjunto de teste:

pred = predict(model_svm, xtest)
head(pred)

## [1] 0 0 0 1 0 1
## Levels: 0 1

postResample(pred = pred ,obs= ytest)

## Accuracy    Kappa 
##   0.9724   0.9448

Em problemas de classificação duas estatísticas são muito usadas:

acurácia: fração das classificações corretas

kappa: similar à acurácia mas normalizada na classificação aleatória dos dados; útil quando há “desbalanço” (imbalance) entre as classes (ex.: se a razão entre as classes 0 e 1 é 70:30, tem-se 70% de acurácia prevendo qualquer objeto como classe 0)

Em classificação uma coisa importante é a “matriz de confusão”, onde se compara, no conjunto de teste, as classificações preditas com as “verdadeiras”:

# matriz de confusão
# TPR (True Positive Rate) = TP/(TP+FN)
# FPR (False Positive Rate) = FP/(TN+FP)
print(confusionMatrix(data = as.factor(pred),reference = as.factor(ytest),positive = '1'))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 1238   57
##          1   12 1193
##                                           
##                Accuracy : 0.9724          
##                  95% CI : (0.9652, 0.9785)
##     No Information Rate : 0.5             
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.9448          
##                                           
##  Mcnemar's Test P-Value : 1.177e-07       
##                                           
##             Sensitivity : 0.9544          
##             Specificity : 0.9904          
##          Pos Pred Value : 0.9900          
##          Neg Pred Value : 0.9560          
##              Prevalence : 0.5000          
##          Detection Rate : 0.4772          
##    Detection Prevalence : 0.4820          
##       Balanced Accuracy : 0.9724          
##                                           
##        'Positive' Class : 1               
##

tab = table(pred, ytest)
# essa função calcula a acurácia 
accuracy <- function(x){sum(diag(x)/(sum(rowSums(x)))) * 100}
accuracy(tab)

## [1] 97.24

Nota-se que se consegue uma excelente acurácia nesse caso.

12 classificação com árvores de decisão

Vamos usar o algoritmo ‘tree’ do pacote rpart

# reprodutibilidade
set.seed(123)

t0 = Sys.time()
model_tree <- train(xtrain,ytrain,method='rpart',
                      trControl = trainControl(method = "cv"))
Sys.time()-t0

## Time difference of 0.7247763 secs

# modelo ajustado:
model_tree

## CART 
## 
## 7500 samples
##   12 predictor
##    2 classes: '0', '1' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 6750, 6750, 6750, 6750, 6750, 6750, ... 
## Resampling results across tuning parameters:
## 
##   cp          Accuracy   Kappa    
##   0.03360000  0.8185333  0.6370667
##   0.05546667  0.7648000  0.5296000
##   0.49813333  0.6706667  0.3413333
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was cp = 0.0336.

# predição com o conjunto de teste
pred = predict(model_tree, xtest)
head(pred)

## [1] 0 0 0 1 1 1
## Levels: 0 1

postResample(pred = pred ,obs= ytest)

## Accuracy    Kappa 
##   0.8004   0.6008

# note a diferença na função predict() para devolver as probabilidades, não as classes
pred2 = predict(model_tree, xtest,type='prob')
head(pred2)

##            0          1
## 8  0.9051621 0.09483794
## 11 0.7837743 0.21622575
## 13 0.7837743 0.21622575
## 14 0.1407692 0.85923077
## 15 0.1407692 0.85923077
## 39 0.3311688 0.66883117

predição e acurácia com o conjunto de teste:

pr = rep(0,length(ytest))
pr[pred2[,2] > pred2[,1]] = 1
tab = table(pr, ytest)
print(confusionMatrix(data = as.factor(pr),reference = as.factor(ytest),positive = '1'))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0  989  238
##          1  261 1012
##                                           
##                Accuracy : 0.8004          
##                  95% CI : (0.7842, 0.8159)
##     No Information Rate : 0.5             
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.6008          
##                                           
##  Mcnemar's Test P-Value : 0.3247          
##                                           
##             Sensitivity : 0.8096          
##             Specificity : 0.7912          
##          Pos Pred Value : 0.7950          
##          Neg Pred Value : 0.8060          
##              Prevalence : 0.5000          
##          Detection Rate : 0.4048          
##    Detection Prevalence : 0.5092          
##       Balanced Accuracy : 0.8004          
##                                           
##        'Positive' Class : 1               
##

note que o desempenho deste modelo é muito pior que o da SVM!

Um bom jeito de visualizar a árvore de decisão é diretamente com a função rpart():

# a função rpart() recebe como input uma fórmula:
n <- c(names(mag))
f = as.formula(paste("ytrain ~", paste(n[!n %in% "medv"], collapse = " + ")))
f

## ytrain ~ u_petro + J0378_petro + J0395_petro + J0410_petro + 
##     J0430_petro + g_petro + J0515_petro + r_petro + J0660_petro + 
##     i_petro + J0861_petro + z_petro

# classificação 
t0 = Sys.time()
tree = rpart(f, data = xtrain, method = "class")
Sys.time()-t0

## Time difference of 0.1766489 secs

# matriz de confusão
# TPR (True Positive Rate) = TP/(TP+FN)
# FPR (False Positive Rate) = FP/(TN+FP)
print(confusionMatrix(data = as.factor(pred),reference = as.factor(ytest),positive = '1'))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0  989  238
##          1  261 1012
##                                           
##                Accuracy : 0.8004          
##                  95% CI : (0.7842, 0.8159)
##     No Information Rate : 0.5             
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.6008          
##                                           
##  Mcnemar's Test P-Value : 0.3247          
##                                           
##             Sensitivity : 0.8096          
##             Specificity : 0.7912          
##          Pos Pred Value : 0.7950          
##          Neg Pred Value : 0.8060          
##              Prevalence : 0.5000          
##          Detection Rate : 0.4048          
##    Detection Prevalence : 0.5092          
##       Balanced Accuracy : 0.8004          
##                                           
##        'Positive' Class : 1               
##

# visualização da árvore:
par(mfrow = c(1,1))
rpart.plot(tree)

Nessa figura fica claro como a classificação é feita; note que apenas 6 bandas foram usadas!

predição e acurácia com o conjunto de teste:

pr = predict(tree, xtest)
prl = rep(0,length(ytest))
prl[pr[,2] > pr[,1]] = 1
pr = prl
tab = table(pr, ytest)
print(confusionMatrix(data = as.factor(pr),reference = as.factor(ytest),positive = '1'))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 1078  153
##          1  172 1097
##                                           
##                Accuracy : 0.87            
##                  95% CI : (0.8562, 0.8829)
##     No Information Rate : 0.5             
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.74            
##                                           
##  Mcnemar's Test P-Value : 0.3181          
##                                           
##             Sensitivity : 0.8776          
##             Specificity : 0.8624          
##          Pos Pred Value : 0.8645          
##          Neg Pred Value : 0.8757          
##              Prevalence : 0.5000          
##          Detection Rate : 0.4388          
##    Detection Prevalence : 0.5076          
##       Balanced Accuracy : 0.8700          
##                                           
##        'Positive' Class : 1               
##

# o desempenho dessa árvore é melhor que o da anterior, mas ainda pior que SVM

em classificação binária a curva ROC (Receiver Operating Characteristics) e a estatística AUC (area under the curve) são muito utilizadas:

# ROC: FPR x TPR
# AUC: Area Under The (ROC) Curve
# AUC varia entre 0.5 and 1, com 0.5 aleatório e 1 perfeito
# vamos usar as probabilidades das classes
par(mfrow = c(1,1))
meu_roc <- function(ytest, pred2){
  yt <- ytest[order(pred2, decreasing=TRUE)]
  data.frame(TPR=cumsum(yt)/sum(yt), FPR=cumsum(!yt)/sum(!yt), yt)
}
a = meu_roc(as.integer(ytest)-1, pred2[,2])
plot(a[,2],a[,1],type='l',col='blue',lwd=3,xlab='taxa de FP',ylab='taxa de TP',main='curva ROC')
abline(0,1,lwd=3)

#AUC
AUC=0
for(i in 2:length(a[,1])){
AUC=AUC+(a[i-1,1]+a[i,1])*(a[i,2]-a[i-1,2])/2
}
round(AUC,3)

## [1] 0.831

tempo de processamento do script:

stopCluster(cl)

Sys.time() - t00

## Time difference of 1.174186 hours

Verifique o tempo de processamento dos vários algoritmos para ter uma ideia do custo computacional de cada um.

exercícios

Estime redshifts fotométricos usando cores (em relação à banda r_petro) com o algoritmo random forest da biblioteca caret. Use o arquivo splus-mag-z.dat, que contém magnitudes e redshifts para uma amostra de galáxias.

1.1. inicialize o programa com uma semente aleatória igual ao seu número USP

1.2. usando as magnitudes deste arquivo, calcule as cores em relação à banda r_petro, isto é, para uma dada banda i, calcule \(cor_i = mag_i-mag_r\). Crie um novo arquivo de dados com as cores e excluindo a banda r (que terá todos os valores iguais a zero)

1.3 restrinja a amostra a 10000 objetos para reduzir o tempo de processamento

1.4. pré-processamento dos dados: normalize as cores com padronização: subtraia a média e divida pelo desvio padrão

1.5. crie conjuntos de treinamento e teste com 80% dos objetos para treinamento e 20% para teste

1.6. treine o modelo, faça predições com o conjunto de teste e compare o resultado (usando \(\sigma_G\)) com o obtido com as 12 magnitudes

1.7. faça um gráfico z_spec x z_phot para visualização dos resultados

Resolva o problema da classificação estrela/galáxia acima usando cores e o algoritmo knn. Inicialize o programa com uma semente aleatória igual ao seu número USP. Compare o resultado (acurácia) com os obtidos acima com SVM usando as 12 magnitudes.

Referências

https://www.machinelearningplus.com/machine-learning/caret-package/

Mendes de Oliveira et al. (2019), MNRAS, 489, 241

aula10

lsj

May, 2023

aula de hoje: regressão e classificação com o pacote caret

0 objetivos

1 aspectos gerais da regressão

1.1 Otimização de um modelo com o algoritmo descendo o gradiente

1.2 modelos e generalização

1.3 regressão robusta

1.4 regressão com kernel

2 o que é o `caret`

3 avaliação dos dados

4 pré-processamento dos dados

5 Treinamento e avaliação

6 regressão: rf - random forest

7 regressão: knn - kth nearest neighbor

8 avaliação da importância relativa das variáveis

9 ajuste de hiperparâmetros de modelos

10 classificação: preparação dos dados

11 classificação com SVM

12 classificação com árvores de decisão

exercícios

Referências

aula10

lsj

May, 2023

aula de hoje: regressão e classificação com o pacote caret

0 objetivos

1 aspectos gerais da regressão

1.1 Otimização de um modelo com o algoritmo descendo o gradiente

1.2 modelos e generalização

1.3 regressão robusta

1.4 regressão com kernel

2 o que é o caret

3 avaliação dos dados

4 pré-processamento dos dados

5 Treinamento e avaliação

6 regressão: rf - random forest

7 regressão: knn - kth nearest neighbor

8 avaliação da importância relativa das variáveis

9 ajuste de hiperparâmetros de modelos

10 classificação: preparação dos dados

11 classificação com SVM

12 classificação com árvores de decisão

exercícios

Referências

2 o que é o `caret`