aula de hoje: regressão, classificação e uma introdução ao pacote `caret`

1. regressão

1.1. otimização de um modelo por gradiente descendente

1.2. modelos e generalização

1.3. regularização: LASSO e ridge regression

1.4. regressão robusta

1.5. regressão com kernel

2. introdução ao ‘caret’

2.1 o que é o caret

2.2 avaliação dos dados

2.3 pré-processamento dos dados

2.4 regressão: rf - random forest

2.5 regressão com gbm - Generalized Boosted Regression Modeling - e ajuste de hiper-parâmetros de modelos

3. classificação

3.1 dados

3.2 k-nn

3.3 regressão logística

3.4 SVM - support vector machine

3.5 árvores de decisão

3.6 random forest

3.7 classificação binária com árvores de decisão com a função rpart

0. inicialização: bibliotecas para a aula de hoje

# para monitorar o tempo de processamento:
t00 = Sys.time()
# reprodutibilidade
set.seed(123)
# bibliotecas
library(glmnet)

## Loading required package: Matrix

## Loaded glmnet 3.0-2

require(doMC)

## Loading required package: doMC

## Loading required package: foreach

## Loading required package: iterators

## Loading required package: parallel

library(carData)
library(olsrr)

## 
## Attaching package: 'olsrr'

## The following object is masked from 'package:datasets':
## 
##     rivers

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:olsrr':
## 
##     cement

library(quantreg)

## Loading required package: SparseM

## 
## Attaching package: 'SparseM'

## The following object is masked from 'package:base':
## 
##     backsolve

require(robustbase)

## Loading required package: robustbase

library(ggplot2)
library(kernlab)

## 
## Attaching package: 'kernlab'

## The following object is masked from 'package:ggplot2':
## 
##     alpha

require(plyr)

## Loading required package: plyr

require(reshape2)

## Loading required package: reshape2

library('caret')

## Loading required package: lattice

library(pROC)

## Type 'citation("pROC")' for a citation.

## 
## Attaching package: 'pROC'

## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var

library(ROCR)

## Loading required package: gplots

## 
## Attaching package: 'gplots'

## The following object is masked from 'package:stats':
## 
##     lowess

library("rpart")
library("rpart.plot")
library("randomForest")

## randomForest 4.6-14

## Type rfNews() to see new features/changes/bug fixes.

## 
## Attaching package: 'randomForest'

## The following object is masked from 'package:ggplot2':
## 
##     margin

library(KernSmooth)

## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009

library(klaR)
library(deldir)

## deldir 0.1-23

library(verification)

## Loading required package: fields

## Loading required package: spam

## Loading required package: dotCall64

## Loading required package: grid

## Spam version 2.5-1 (2019-12-12) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction 
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.

## 
## Attaching package: 'spam'

## The following objects are masked from 'package:SparseM':
## 
##     backsolve, forwardsolve

## The following object is masked from 'package:Matrix':
## 
##     det

## The following objects are masked from 'package:base':
## 
##     backsolve, forwardsolve

## Loading required package: maps

## 
## Attaching package: 'maps'

## The following object is masked from 'package:plyr':
## 
##     ozone

## See https://github.com/NCAR/Fields for
##  an extensive vignette, other supplements and source code

## Loading required package: boot

## 
## Attaching package: 'boot'

## The following object is masked from 'package:lattice':
## 
##     melanoma

## The following object is masked from 'package:robustbase':
## 
##     salinity

## Loading required package: CircStats

## 
## Attaching package: 'CircStats'

## The following object is masked from 'package:kernlab':
## 
##     rvm

## Loading required package: dtw

## Loading required package: proxy

## 
## Attaching package: 'proxy'

## The following object is masked from 'package:spam':
## 
##     as.matrix

## The following object is masked from 'package:SparseM':
## 
##     as.matrix

## The following object is masked from 'package:Matrix':
## 
##     as.matrix

## The following objects are masked from 'package:stats':
## 
##     as.dist, dist

## The following object is masked from 'package:base':
## 
##     as.matrix

## Loaded dtw v1.21-3. See ?dtw for help, citation("dtw") for use in publication.

## Registered S3 method overwritten by 'verification':
##   method    from
##   lines.roc pROC

## 
## Attaching package: 'verification'

## The following object is masked from 'package:pROC':
## 
##     lines.roc

1. Regressão

1.1 otimização de um modelo por gradiente descendente

Vamos ilustrar o algoritmo de otimização por gradiente descendente encontrando os parâmetros de um modelo linear simples: \(y=w_0+w_1 x + \epsilon\):

# https://www.r-bloggers.com/linear-regression-by-gradient-descent/
# simulação de um modelo linear: y=1+2x, com x entre -1 e 1 e  ruído gaussiano com sd = 0.2
n=100
x <- runif(n, -1, 1)
y <- 1+2*x + rnorm(n,sd=0.2)

# vamos primeiro ajustar um modelo linear:
modelo_linear <- lm( y ~ x )
print(modelo_linear)

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

# visualização:
plot(x,y, col=rgb(0.2,0.4,0.6,0.4), main='regressão linear por gradiente descendente')
abline(modelo_linear, col='blue')

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

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

# guarda a história:
custo_historia <- double(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 w1
X <- cbind(1, matrix(x))

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

# resultado
print(w)

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

# visualização do resultado:
par(mfrow=c(1,1))
plot(x,y, col=rgb(0.2,0.4,0.6,0.4), main='regressão linear por gradiente descendente')
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='line', col='blue', lwd=2, main='convergência da função custo', ylab='custo', xlab='Iterações')

## Warning in plot.xy(xy, type, ...): plot type 'line' will be truncated to first
## character

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 = 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,pch=20,col='blue',type='l',ylim=c(0,1), main='vermelho: treino;  verde: teste')
points(xtreino,ytreino,col='red')
points(xtest,ytest,col='green')

# 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()
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

## [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_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,pch=20,col='blue',type='l',ylim=c(0,1), main='preto: polinômio de grau 9')
points(xtreino,ytreino,col='red')
points(xtest,ytest,col='green')
lines(xt,yt,lwd=3)

1.3 regularização: LASSO e ridge regression

Nesses métodos um termo de regularização torna os parâmetros mais estáveis

Embora os modelos ajustem os dados pior que mínimos quadrados comum (OLS, ordinary least squares), eles têm melhor capacidade de generalização

Vamos usar o pacote glmnet: ajusta um modelo linear com regularização

esta função tem 2 parâmetros principais: lambda é um parâmetro de regularização (controla a magnitude dos parâmetros da regressão)

alpha controla o método- alpha = 1 corresponde a LASSO, alpha = 0 a ridge regression e alpha entre 0 e 1 interpola entre esses dois métodos

# https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html
# reprodutibilidade
set.seed(123)

# dados: vamos gerar dados com um modelo linear e introduzir alguns outliers
# modelo linear: y=1+2x, com x entre -1 e 1 e  ruído gaussiano com sd = 0.1
n=15
x <- runif(n, -1, 1)
y <- 1+2*x + rnorm(n,sd=0.1)

# outliers: idem com sd = 1
x[n+1]=-0.6
y[n+1]=-2
x[n+2]=-0.4
y[n+2]=-1.5

# vamos começar plotando os dados e o modelo:
ymodelo = 1+2*x

plot(x,y,pch=20)
lines(x,ymodelo,type='l',col='red',main='vermelho: modelo, verde: lm com os outliers')

# modelo linear normal (OLS):
lm.mod = lm(y~x)
summary(lm.mod)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.3846 -0.1214  0.1485  0.3178  0.5787 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.7952     0.1421   5.594 5.12e-05 ***
## x             2.3510     0.2496   9.419 1.09e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5824 on 15 degrees of freedom
## Multiple R-squared:  0.8554, Adjusted R-squared:  0.8457 
## F-statistic: 88.71 on 1 and 15 DF,  p-value: 1.092e-07

abline(lm.mod,col='green')

# ridge regression: 
# vamos considerar inicialmente vários valores para o parâmetro de regularização lambda
lambda <- 10^seq(10, -2, length = 100)

#ridge regression:  alpha = 0
# como x só tem 1 dimensão e o glmnet precisa que x tenha ao menos duas, usamos o truque abaixo:
xx=cbind(x,0)

ridge.mod <- glmnet(xx, y, alpha = 0)

# coeficientes do ajuste
predict(ridge.mod, s = 0, type = 'coefficients')[1:2,]

## (Intercept)           x 
##   0.8079507   2.1519292

# plot
plot(x,y,pch=20)
lines(x,ymodelo,type='l',col='red')
abline(lm.mod,col='green')

a=predict(ridge.mod, s = 0, type = 'coefficients')[1]
b=predict(ridge.mod, s = 0, type = 'coefficients')[2]
lines(x,a+b*x,col='blue')

# valor ótimo de lambda para ridge por validação cruzada
lambdas <- 10^seq(3, -2, by = -.1)
cv_fit <- cv.glmnet(xx, y, alpha = 0, lambda = lambdas, nfolds = 10)

## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per
## fold

# valor ótimo:
opt_lambda <- cv_fit$lambda.min
opt_lambda

## [1] 0.01584893

# modelos ajustados na CV
fit <- cv_fit$glmnet.fit
# coeficientes do ajuste
predict(fit, s = 0, type = 'coefficients')[1:2,]

## (Intercept)           x 
##   0.7962139   2.3347219

# agora LASSO  alpha = 1
lasso.mod <- glmnet(xx, y, alpha = 1)

# coeficientes do ajuste
predict(lasso.mod, s = 0, type = 'coefficients')[1:2,]

## (Intercept)           x 
##   0.7961648   2.3354854

# plot
plot(x,y,pch=20)
lines(x,ymodelo,type='l',col='red')
abline(lm.mod,col='green')
lines(x,a+b*x,col='blue')
a2=predict(lasso.mod, s = 0, type = 'coefficients')[1]
b2=predict(lasso.mod, s = 0, type = 'coefficients')[2]
lines(x,a2+b2*x,col='black')

legend("topright",legend=c("modelo","lm",'ridge','lasso'),pch=19,col=c("red","green",'blue','black'))

1.4 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 os dados (x,y) do exemplo acima

# 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)
# 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

# 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))

1.5 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.9)
plot(X,Y,pch=20)
lines(Kreg, lwd=4, col="purple")

# kernel box
Kreg = ksmooth(x=X,y=Y,kernel = "box",bandwidth = 0.9)
plot(X,Y,pch=20)
lines(Kreg, lwd=4, col="orange")

# 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
# leave-one-out-cross-validation (LOOCV)
n = length(X)
# n: sample size
h_seq = seq(from=0.2,to=2.0, by=0.1)
# smoothing bandwidths we are using
CV_err_h = rep(NA,length(h_seq))
for(j in 1:length(h_seq)){
h_using = h_seq[j]
CV_err = rep(NA, n)
for(i in 1:n){
X_val = X[i]
Y_val = Y[i]
# validation set
X_tr = X[-i]
Y_tr = Y[-i]
# training set
Y_val_predict = ksmooth(x=X_tr,y=Y_tr,kernel = "normal",bandwidth=h_using,x.points = X_val)
CV_err[i] = (Y_val - Y_val_predict$y)^2
# we measure the error in terms of difference square
}
CV_err_h[j] = mean(CV_err)
}
# erro medio da validação cruzada para os vários valores de h
CV_err_h

##  [1] 0.10183202 0.09692202 0.09471400 0.09323028 0.09230145 0.09174197
##  [7] 0.09144213 0.09141794 0.09171727 0.09240507 0.09351454 0.09506037
## [13] 0.09704216 0.09944046 0.10225141 0.10545900 0.10905926 0.11305207
## [19] 0.11744657

plot(x=h_seq, y=CV_err_h, type="b", lwd=3, col="blue",xlab="bandwidth h", ylab="LOOCV prediction error")

which(CV_err_h == min(CV_err_h))

## [1] 8

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

## [1] 0.9

#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

2 a biblioteca `caret`

2.1 o que é o ‘caret’

caret: Classification And REgression Training

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

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

Este pacote do R contém várias funções que facilitam muito várias tarefas em aprendizado de máquina (ML).

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

2.2 avaliação dos dados

leitura dos dados: vamos considerar uma amostra de galáxias do S-PLUS com fotometria nas 12 bandas; o arquivo contém, para cada objeto da amostra, as 12 magnitudes e o redshift espectroscópico (do SDSS).

dados = as.data.frame(read.table('splus-mag-z.dat', sep = "", header=TRUE))

# vamos examinar os dados:
dim(dados)

## [1] 55803    13

# cada linha 
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

Como estamos interessados apenas em explorar o pacote e algumas de suas principais funções, vamos restringir o número de objetos em 1000, para o treinamento dos modelos não demorar muito:

nsample = 1000
dados = dados[1:1000,]

vamos dar uma olhada nos dados:

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

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

vamos visualizar as correlações entre as variáveis (magnitudes no caso); as variáveis usadas na estimativa são denominadas “features”.

correlações entre as variáveis:

library(corrplot)

## corrplot 0.84 loaded

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

visualização usando featurePlot:

redshift versus cada banda fotométrica:

featurePlot(x = dados[, 1:12], y = dados$z_SDSS, plot = "scatter")

2.3 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.

Há muitos tipos de pré-processamento.Os mais comuns são

reescalonar as variáveis entre 0 e 1

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

vamos colocar as variáveis entre 0 e 1:

maxs <- apply(dados, 2, max) 
mins <- apply(dados, 2, min)
dadosnorm <- as.data.frame(scale(dados, center = mins, scale = maxs - mins))

#se quisesse subtrair a média e dividir pela variância:
#dadosnorm <- as.data.frame(scale(x, center = TRUE, scale = TRUE))

Vamos agora considerar a criação dos conjuntos de treinamento e teste e a seleção do algoritmo.

criação de conjuntos de treinamento e teste

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

# número das linhas dos objetos selecionados para o conjunto de treinamento
numlinhastreinamento =  createDataPartition(dados$z_SDSS, p=0.75, list=FALSE)

# conjunto de treinamento:
treinamento =  dadosnorm[numlinhastreinamento,]

# conjunto de teste
teste = dadosnorm[-numlinhastreinamento,]

# vamos definir as variáveis x e y para uso posterior
x = treinamento[,1:12]     # magnitudes
y = treinamento[,13]       # redshift

xtrain = treinamento[,1:12]    # magnitudes
ytrain = treinamento[,13]      # redshift
xtest = teste[,1:12]
ytest = teste[,13]

# o caret requer isso:
ytrain = as.vector(ytrain)
ytest = as.vector(teste[,13])

seleção do algoritmo:

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

nomes_dos_modelos <- paste(names(getModelInfo()), collapse=',  ')
nomes_dos_modelos

## [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,  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('pcaNNet')

##     model parameter         label forReg forClass probModel
## 1 pcaNNet      size #Hidden Units   TRUE     TRUE      TRUE
## 2 pcaNNet     decay  Weight Decay   TRUE     TRUE      TRUE

Vamos considerar um exemplo: regressão com random forest

modelLookup('rf')

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

2.4 regressão com random forest

t0 = Sys.time()
model_rf<-train(xtrain,ytrain,method='rf')
Sys.time() - t0

## Time difference of 2.252486 mins

# sumário do modelo
print(model_rf)

## Random Forest 
## 
## 751 samples
##  12 predictor
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 751, 751, 751, 751, 751, 751, ... 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE        Rsquared   MAE       
##    2    0.10113857  0.5138191  0.07291636
##    7    0.09946456  0.5278380  0.07097013
##   12    0.10013458  0.5211102  0.07131213
## 
## 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)

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

# voltando à escala original para comparação
pred <- pred*(max(y)-min(y))+min(y)
ytest <- ytest*(max(y)-min(y))+min(y)

# estatísticas da performance do algoritmo
# soma dos erros quadráticos
SSE = sum((ytest -pred)^2)   

# erro quadrático médio
RMSE = sqrt(SSE/length(pred))
RMSE

## [1] 0.07843343

# sigma equivalente da gaussiana
relz = (ytest - pred)
sig_G = 0.7413*(quantile(relz,0.75,names = FALSE) - quantile(relz,0.25,names = FALSE))
sig_G

## [1] 0.06489556

#visualização do resultado
my_data = as.data.frame(cbind(predicted = pred,
                            observed = ytest))
# Plot predictions vs test data
ggplot(my_data,aes(predicted, observed)) +
      geom_point(color = "darkred", alpha = 0.5) + 
      geom_smooth(method=lm)+ ggtitle('GLM ') +
      ggtitle("RF: prediction vs test data") +
      xlab("prediction ") +
      ylab("observed data") +
      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))

2.5 regressão com gbm (Generalized Boosted Regression Modeling) e ajuste de hiperparâmetros de modelos

Vamos ilustrar o método de validação cruzada com o algoritmo ‘gbm’, Generalized Boosted Regression Modeling.

Vamos fazer a validação cruzada 3 vezes e repetir isso também 3 vezes

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

Para saber quais são os hiper-parâ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

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:

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

## Time difference of 2.10176 mins

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

# voltando à escala original para comparação
pred <- pred*(max(y)-min(y))+min(y)
ytest <- ytest*(max(y)-min(y))+min(y)

SSE = sum((ytest -pred)^2)    
RMSE = sqrt(SSE/length(pred))
RMSE

## [1] 0.0699529

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

## [1] 0.05740729

#visualização do resultado
my_data = as.data.frame(cbind(predicted = pred,
                            observed = ytest))
# Plot predictions vs test data
ggplot(my_data,aes(predicted, observed)) +
      geom_point(color = "darkred", alpha = 0.5) + 
      geom_smooth(method=lm)+ ggtitle('GBM ') +
      ggtitle("GBM: prediction vs test data for photo-z") +
      xlab("prediction ") +
      ylab("observed data") +
      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))

3. Classificação

3.1 dados

vamos considerar uma amostra de fontes puntuais do SDSS classificadas em 3 classes:

1- QSO, 2- MS+RG, 3- WD

esta amostra é baseada na discutida por Feigelson & Babu em Modern Statistical Methods for Astronomy

(https://astrostatistics.psu.edu/MSMA/MSMA_R_scripts.html)

# leitura dos dados:
tabela = as.data.frame(read.table(file="sdss_point_sources.dat", header=TRUE))

# alguns detalhes dos dados
dim(tabela)

## [1] 9000    5

head(tabela)

##      u_g   g_r   r_i    i_z classe
## 1 -0.079 0.136 0.233  0.046      1
## 2  0.033 0.255 0.454  0.300      1
## 3  0.110 0.425 0.221 -0.158      1
## 4  0.325 0.448 0.114  0.221      1
## 5  0.220 0.049 0.189  0.040      1
## 6  0.171 0.104 0.169  0.188      1

summary(tabela)

##       u_g               g_r               r_i               i_z          
##  Min.   :-0.7380   Min.   :-0.7550   Min.   :-1.0725   Min.   :-2.33300  
##  1st Qu.: 0.3387   1st Qu.: 0.0840   1st Qu.: 0.0290   1st Qu.:-0.01287  
##  Median : 1.0292   Median : 0.3779   Median : 0.1608   Median : 0.06600  
##  Mean   : 1.0081   Mean   : 0.3625   Mean   : 0.1457   Mean   : 0.05700  
##  3rd Qu.: 1.4738   3rd Qu.: 0.5840   3rd Qu.: 0.2544   3rd Qu.: 0.14493  
##  Max.   : 5.3580   Max.   : 4.2650   Max.   : 1.8140   Max.   : 1.31812  
##      classe 
##  Min.   :1  
##  1st Qu.:2  
##  Median :2  
##  Mean   :2  
##  3rd Qu.:2  
##  Max.   :3

# número de objetos em cada classe:
length(tabela$classe[tabela$classe == 1])

## [1] 2000

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

## [1] 5000

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

## [1] 2000

vamos examinar os dados em um diagrama cor-magnitude:

par(mfrow = c(1,2))
ca = tabela$u_g
cb = tabela$g_r
smoothScatter(ca,cb,nrpoints=0,add=FALSE,xlab="u-g",ylab="g-r",xlim=c(-1,3),ylim=c(-1,3),main='SDSS: fontes puntuais')

plot(ca, cb, xlim=c(-1,3), ylim=c(-1,3), pch=20, col=tabela$classe, cex=0.6, cex.lab=1.6, cex.axis=1.6, main='', xlab='u-g', ylab='g-r')
legend('bottomright', c('QSO','MS + RG','WD'), pch=20, col=c('black','red','green'), cex=0.8)

vamos definir os atributos a serem usados na classificação:

# definindo os atributos
cores = tabela[,1:4]

conjuntos de treino e teste:

ntreino = round(0.75*nrow(tabela))
nteste = nrow(tabela)-ntreino
c(ntreino,nteste)

## [1] 6750 2250

# 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 = cores[indice,]

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

xtest = cores[-indice,]
ytest = as.factor(tabela[-indice,5])

3.2 classificação com k-nn

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)

modelLookup('knn')

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

t0 = Sys.time()
set.seed(123)
grid <- data.frame(k=c(1,5,10,20,50))
model_knn <- train(xtrain,ytrain,method='knn',trControl=trainControl(method='cv',number=10), tuneGrid=grid)
Sys.time()-t0

## Time difference of 2.538007 secs

print(model_knn)

## k-Nearest Neighbors 
## 
## 6750 samples
##    4 predictor
##    3 classes: '1', '2', '3' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 6075, 6075, 6074, 6075, 6075, 6075, ... 
## Resampling results across tuning parameters:
## 
##   k   Accuracy   Kappa    
##    1  0.9788144  0.9641938
##    5  0.9805913  0.9671156
##   10  0.9792584  0.9648217
##   20  0.9759987  0.9592099
##   50  0.9669610  0.9435952
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was k = 5.

pred = predict(model_knn, xtest)
print(confusionMatrix(data = as.factor(pred),reference = as.factor(ytest),positive = '1'))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    1    2    3
##          1  478    2   14
##          2   10 1246    0
##          3   13    1  486
## 
## Overall Statistics
##                                           
##                Accuracy : 0.9822          
##                  95% CI : (0.9759, 0.9873)
##     No Information Rate : 0.5551          
##     P-Value [Acc > NIR] : < 2e-16         
##                                           
##                   Kappa : 0.97            
##                                           
##  Mcnemar's Test P-Value : 0.09492         
## 
## Statistics by Class:
## 
##                      Class: 1 Class: 2 Class: 3
## Sensitivity            0.9541   0.9976   0.9720
## Specificity            0.9909   0.9900   0.9920
## Pos Pred Value         0.9676   0.9920   0.9720
## Neg Pred Value         0.9869   0.9970   0.9920
## Prevalence             0.2227   0.5551   0.2222
## Detection Rate         0.2124   0.5538   0.2160
## Detection Prevalence   0.2196   0.5582   0.2222
## Balanced Accuracy      0.9725   0.9938   0.9820

par(mfrow = c(1,2))
plot(xtest$u_g, xtest$g_r, xlim=c(-1,3), ylim=c(-1,3), pch=20, col=pred, cex=0.6, cex.lab=1.6, cex.axis=1.6, main='KNN', xlab='u-g', ylab='g-r')
legend('bottomright', c('QSO','MS + RG','WD'), pch=20, col=c('black','red','green'), cex=0.8)
plot(jitter(as.numeric(pred), factor=0.5), jitter(as.numeric
   (ytest), factor=0.5), pch=20, cex=0.5, xlab='classe do classificador',
   ylab='classe verdadeira', xaxp=c(1,3,2),yaxp=c(1,3,2),col='blue')

tab = table(pred, ytest)

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

## [1] 98.22222

O knn promove uma tesselagem de Voronoi no espaço dos dados:

# tesselagem de Voronoi
#https://letstalkdata.com/2014/05/creating-voronoi-diagrams-with-ggplot/
c1 = xtest$u_g[xtest$u_g < 0 &  xtest$g_r < 0]
c2 =  xtest$g_r[xtest$u_g < 0 &  xtest$g_r < 0]
df <- data.frame(c2,c1)

#This creates the voronoi line segments
voronoi <- deldir(df$c1, df$c2)

## 
##      PLEASE NOTE:  The components "delsgs" and "summary" of the
##  object returned by deldir() are now DATA FRAMES rather than
##  matrices (as they were prior to release 0.0-18).
##  See help("deldir").
##  
##      PLEASE NOTE: The process that deldir() uses for determining
##  duplicated points has changed from that used in version
##  0.0-9 of this package (and previously). See help("deldir").

ggplot(data=df, aes(x=c1,y=c2)) +
  #Plot the voronoi lines
  geom_segment(
    aes(x = x1, y = y1, xend = x2, yend = y2),
    size = 2,
    data = voronoi$dirsgs,
    linetype = 1,
    color= "#FFB958") + 
  #Plot the points
  geom_point(
    fill=rgb(70,130,180,255,maxColorValue=255),
    pch=21,
    size = 4,
    color="#333333") +
    labs(x = "u-g") +  labs(y = "g-r") + labs(title = "tesselagem de Voronoi")

3.3 regressão logística

função sigmóide ou logística:

par(mfrow = c(1,1))
xx = seq(-3,3,0.01)
flogistica = function(x,a){1/(1+exp(-a*x))}
plot(xx,flogistica(xx,a=1.702),xlab="x",ylab="Sigmoide(ax)",type='l',col='red', cex=0.6, cex.lab=1.6, cex.axis=1.6,lwd=2,main='azul: a=1  vermelho: a=1.702')
lines(xx,flogistica(xx,a=1),col='blue',lwd=2)

 #note que sigmóide com  a=1.702 é quase o mesmo que pnorm

para ilustrar a classificação logística vamos considerar uma classificação binária: 1- QSO, 2- estrela

ytrain2 = as.numeric(ytrain)
ytrain2[ytrain != 1] = 2
ytrain2 = as.factor(ytrain2)
 
ytest2 = as.numeric(ytest)
ytest2[ytest != 1] = 2
ytest2 = as.factor(ytest2)

# regressão logística:
t0 = Sys.time()
set.seed(123)
model_log <- train(xtrain,ytrain2, method="glm", family="binomial")
Sys.time()-t0

## Time difference of 1.372689 secs

print(model_log)

## Generalized Linear Model 
## 
## 6750 samples
##    4 predictor
##    2 classes: '1', '2' 
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 6750, 6750, 6750, 6750, 6750, 6750, ... 
## Resampling results:
## 
##   Accuracy   Kappa    
##   0.8941216  0.6675203

pred = predict(model_log, xtest)
print(confusionMatrix(data = as.factor(pred),reference = as.factor(ytest2),positive = '1'))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    1    2
##          1  321   53
##          2  180 1696
##                                           
##                Accuracy : 0.8964          
##                  95% CI : (0.8831, 0.9087)
##     No Information Rate : 0.7773          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.6711          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.6407          
##             Specificity : 0.9697          
##          Pos Pred Value : 0.8583          
##          Neg Pred Value : 0.9041          
##              Prevalence : 0.2227          
##          Detection Rate : 0.1427          
##    Detection Prevalence : 0.1662          
##       Balanced Accuracy : 0.8052          
##                                           
##        'Positive' Class : 1               
##

par(mfrow = c(1,2))
plot(xtest$u_g, xtest$g_r, xlim=c(-1,3), ylim=c(-1,3), pch=20, col=pred, cex=0.6, cex.lab=1.6, cex.axis=1.6, main='reg. logística', xlab='u-g', ylab='g-r')
legend('bottomright', c('QSO','MS + RG + WD'), pch=20, col=c('red','black'), cex=0.8)
plot(jitter(as.numeric(pred), factor=0.5), jitter(as.numeric
   (ytest2), factor=0.5), pch=20, cex=0.5, xlab='classe do classificador',
   ylab='classe verdadeira', xaxp=c(1,3,2),yaxp=c(1,3,2),col='blue')

tab = table(pred, ytest2)
ac_log=accuracy(tab)

3.4 classificação com svm: support vector machines

modelLookup('svmLinear')

##       model parameter label forReg forClass probModel
## 1 svmLinear         C  Cost   TRUE     TRUE      TRUE

t0 = Sys.time()
set.seed(123)
model_svmLinear <- train(xtrain,ytrain,method='svmLinear',trControl=trainControl(method='cv',number=10))
Sys.time()-t0

## Time difference of 5.408498 secs

print(model_svmLinear)

## Support Vector Machines with Linear Kernel 
## 
## 6750 samples
##    4 predictor
##    3 classes: '1', '2', '3' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 6075, 6075, 6074, 6075, 6075, 6075, ... 
## Resampling results:
## 
##   Accuracy   Kappa    
##   0.9250352  0.8720519
## 
## Tuning parameter 'C' was held constant at a value of 1

pred = predict(model_svmLinear, xtest)
print(confusionMatrix(data = as.factor(pred),reference = as.factor(ytest),positive = '1'))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    1    2    3
##          1  460    2  117
##          2   34 1247    2
##          3    7    0  381
## 
## Overall Statistics
##                                           
##                Accuracy : 0.928           
##                  95% CI : (0.9165, 0.9383)
##     No Information Rate : 0.5551          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.8775          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
## 
## Statistics by Class:
## 
##                      Class: 1 Class: 2 Class: 3
## Sensitivity            0.9182   0.9984   0.7620
## Specificity            0.9320   0.9640   0.9960
## Pos Pred Value         0.7945   0.9719   0.9820
## Neg Pred Value         0.9755   0.9979   0.9361
## Prevalence             0.2227   0.5551   0.2222
## Detection Rate         0.2044   0.5542   0.1693
## Detection Prevalence   0.2573   0.5702   0.1724
## Balanced Accuracy      0.9251   0.9812   0.8790

par(mfrow = c(1,2))
plot(xtest$u_g, xtest$g_r, xlim=c(-1,3), ylim=c(-1,3), pch=20, col=pred, cex=0.6, cex.lab=1.6, cex.axis=1.6, main='SVMLINEAR', xlab='u-g', ylab='g-r')
legend('bottomright', c('QSO','MS + RG','WD'), pch=20, col=c('black','red','green'), cex=0.8)
plot(jitter(as.numeric(pred), factor=0.5), jitter(as.numeric
   (ytest), factor=0.5), pch=20, cex=0.5, xlab='classe do classificador',
   ylab='classe verdadeira', xaxp=c(1,3,2),yaxp=c(1,3,2),col='blue')

tab = table(pred, ytest)
ac_svmLinear=accuracy(tab)

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

modelLookup('rpart')

##   model parameter                label forReg forClass probModel
## 1 rpart        cp Complexity Parameter   TRUE     TRUE      TRUE

t0 = Sys.time()
set.seed(123)
model_rpart <- train(xtrain,ytrain,method='rpart',trControl=trainControl(method='cv',number=10))
Sys.time()-t0

## Time difference of 0.9364991 secs

print(model_rpart)

## CART 
## 
## 6750 samples
##    4 predictor
##    3 classes: '1', '2', '3' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 6075, 6075, 6074, 6075, 6075, 6075, ... 
## Resampling results across tuning parameters:
## 
##   cp         Accuracy   Kappa    
##   0.0390130  0.8998500  0.8283582
##   0.2954318  0.8208871  0.6933862
##   0.4888296  0.6415136  0.2423872
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was cp = 0.039013.

pred = predict(model_rpart, xtest)
print(confusionMatrix(data = as.factor(pred),reference = as.factor(ytest),positive = '1'))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    1    2    3
##          1  393    0   94
##          2   50 1248   14
##          3   58    1  392
## 
## Overall Statistics
##                                           
##                Accuracy : 0.9036          
##                  95% CI : (0.8906, 0.9154)
##     No Information Rate : 0.5551          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.8347          
##                                           
##  Mcnemar's Test P-Value : 4.727e-15       
## 
## Statistics by Class:
## 
##                      Class: 1 Class: 2 Class: 3
## Sensitivity            0.7844   0.9992   0.7840
## Specificity            0.9463   0.9361   0.9663
## Pos Pred Value         0.8070   0.9512   0.8692
## Neg Pred Value         0.9387   0.9989   0.9400
## Prevalence             0.2227   0.5551   0.2222
## Detection Rate         0.1747   0.5547   0.1742
## Detection Prevalence   0.2164   0.5831   0.2004
## Balanced Accuracy      0.8653   0.9676   0.8751

par(mfrow = c(1,2))
plot(xtest$u_g, xtest$g_r, xlim=c(-1,3), ylim=c(-1,3), pch=20, col=pred, cex=0.6, cex.lab=1.6, cex.axis=1.6, main='RPART', xlab='u-g', ylab='g-r')
legend('bottomright', c('QSO','MS + RG','WD'), pch=20, col=c('black','red','green'), cex=0.8)
plot(jitter(as.numeric(pred), factor=0.5), jitter(as.numeric
   (ytest), factor=0.5), pch=20, cex=0.5, xlab='classe do classificador',
   ylab='classe verdadeira', xaxp=c(1,3,2),yaxp=c(1,3,2),col='blue')

tab = table(pred, ytest)
ac_rpart=accuracy(tab)

3.6 classificação com random forest

modelLookup('rf')

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

t0 = Sys.time()
set.seed(123)
model_rf <- train(xtrain,ytrain,method='rf',trControl=trainControl(method='cv',number=10))
Sys.time()-t0

## Time difference of 56.18436 secs

print(model_rf)

## Random Forest 
## 
## 6750 samples
##    4 predictor
##    3 classes: '1', '2', '3' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 6075, 6075, 6074, 6075, 6075, 6075, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##   2     0.9808871  0.9676589
##   3     0.9810355  0.9679088
##   4     0.9798503  0.9659057
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 3.

pred = predict(model_rf, xtest)
print(confusionMatrix(data = as.factor(pred),reference = as.factor(ytest),positive = '1'))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    1    2    3
##          1  484    3   13
##          2    8 1245    0
##          3    9    1  487
## 
## Overall Statistics
##                                           
##                Accuracy : 0.9849          
##                  95% CI : (0.9789, 0.9895)
##     No Information Rate : 0.5551          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.9745          
##                                           
##  Mcnemar's Test P-Value : 0.2615          
## 
## Statistics by Class:
## 
##                      Class: 1 Class: 2 Class: 3
## Sensitivity            0.9661   0.9968   0.9740
## Specificity            0.9909   0.9920   0.9943
## Pos Pred Value         0.9680   0.9936   0.9799
## Neg Pred Value         0.9903   0.9960   0.9926
## Prevalence             0.2227   0.5551   0.2222
## Detection Rate         0.2151   0.5533   0.2164
## Detection Prevalence   0.2222   0.5569   0.2209
## Balanced Accuracy      0.9785   0.9944   0.9841

par(mfrow = c(1,1))
plot(model_rf)

par(mfrow = c(1,2))
plot(xtest$u_g, xtest$g_r, xlim=c(-1,3), ylim=c(-1,3), pch=20, col=pred, cex=0.6, cex.lab=1.6, cex.axis=1.6, main='RF', xlab='u-g', ylab='g-r')
legend('bottomright', c('QSO','MS + RG','WD'), pch=20, col=c('black','red','green'), cex=0.8)
plot(jitter(as.numeric(pred), factor=0.5), jitter(as.numeric
   (ytest), factor=0.5), pch=20, cex=0.5, xlab='classe do classificador',
   ylab='classe verdadeira', xaxp=c(1,3,2),yaxp=c(1,3,2),col='blue')

tab = table(pred, ytest)
ac_rf=accuracy(tab)

3.7 classificação binária com árvores de decisão usando a função rpart

a função rpart também é útil para árvores de decisão pois permite fazer um gráfico bonito da árvore!

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

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

# alguns detalhes dos dados
dim(tabela)

## [1] 1000   13

head(tabela)

##   uJAVA_petro F378_petro F395_petro F410_petro F430_petro g_petro F515_petro
## 1       23.23      24.09      21.06      22.11      21.36   21.76      20.69
## 2       21.92      20.95      21.30      20.00      21.24   20.73      19.86
## 3       20.44      20.78      20.70      20.92      20.05   20.15      19.55
## 4       21.55      22.96      20.22      19.63      19.83   19.87      19.39
## 5       22.32      21.41      26.56      20.28      20.52   20.28      20.02
## 6       21.26      20.43      20.60      19.44      19.16   18.52      18.34
##   r_petro F660_petro i_petro F861_petro z_petro classe
## 1   20.87      20.53   20.14      19.83   19.71      0
## 2   20.01      20.18   20.36      20.53   20.00      0
## 3   19.49      19.36   19.08      19.00   19.07      0
## 4   19.02      18.78   18.55      18.02   18.17      0
## 5   19.60      19.53   19.34      19.47   18.99      1
## 6   17.65      17.56   17.34      17.30   17.19      0

# número de objetos em cada classe:
length(tabela$classe[tabela$classe == 0])

## [1] 493

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

## [1] 507

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

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

# pré-processamento
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:

# conjuntos de treino e teste:
ntreino = round(0.75*nrow(tabela))
nteste = nrow(tabela)-ntreino
c(ntreino,nteste)

## [1] 750 250

# 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])

Vamos definir ‘f’ como a fórmula a ser passada para rpart:

n <- c(names(mag))
f = as.formula(paste("ytrain ~", paste(n[!n %in% "medv"], collapse = " + ")))
f

## ytrain ~ uJAVA_petro + F378_petro + F395_petro + F410_petro + 
##     F430_petro + g_petro + F515_petro + r_petro + F660_petro + 
##     i_petro + F861_petro + z_petro

tree = rpart(f, data = xtrain, method = "class")

visualização da árvore:

rpart.plot(tree)

predição e acurácia:

pred = predict(tree, xtest)
prl = rep(0,length(ytest))
prl[pred[,2] > pred[,1]] = 1
pr = prl
tab = table(pr, ytest)
# tab é a matriz de confusão
tab

##    ytest
## pr   0  1
##   0 89 39
##   1 37 85

length(pr[pr == 1 & ytest == 0])/length(pr)

## [1] 0.148

length(pr[pr == 0 & ytest == 1])/length(pr)

## [1] 0.156

print(confusionMatrix(data = as.factor(pr),reference = as.factor(ytest),positive = '1'))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  0  1
##          0 89 39
##          1 37 85
##                                           
##                Accuracy : 0.696           
##                  95% CI : (0.6349, 0.7524)
##     No Information Rate : 0.504           
##     P-Value [Acc > NIR] : 5.673e-10       
##                                           
##                   Kappa : 0.3919          
##                                           
##  Mcnemar's Test P-Value : 0.9087          
##                                           
##             Sensitivity : 0.6855          
##             Specificity : 0.7063          
##          Pos Pred Value : 0.6967          
##          Neg Pred Value : 0.6953          
##              Prevalence : 0.4960          
##          Detection Rate : 0.3400          
##    Detection Prevalence : 0.4880          
##       Balanced Accuracy : 0.6959          
##                                           
##        'Positive' Class : 1               
##

accuracy(tab)

## [1] 69.6

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

curva ROC.

# 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

par(mfrow = c(1,1))
meu_roc <- function(ytest, pred){
  yt <- ytest[order(pred, decreasing=TRUE)]
  data.frame(TPR=cumsum(yt)/sum(yt), FPR=cumsum(!yt)/sum(!yt), yt)
}
a = meu_roc(as.integer(ytest)-1, pred[,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.677

# tempo de processamento deste script
Sys.time()-t00

## Time difference of 5.698458 mins

aula12

lsj

June 06, 2020

aula de hoje: regressão, classificação e uma introdução ao pacote `caret`

1. regressão

2. introdução ao ‘caret’

3. classificação

0. inicialização: bibliotecas para a aula de hoje

1. Regressão

1.1 otimização de um modelo por gradiente descendente

1.2 modelos e generalização

1.3 regularização: LASSO e ridge regression

1.4 regressão robusta

1.5 regressão com kernel

2 a biblioteca `caret`

2.1 o que é o ‘caret’

2.2 avaliação dos dados

2.3 pré-processamento dos dados

2.4 regressão com random forest

2.5 regressão com gbm (Generalized Boosted Regression Modeling) e ajuste de hiperparâmetros de modelos

3. Classificação

3.1 dados

3.2 classificação com k-nn

3.3 regressão logística

3.4 classificação com svm: support vector machines

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

3.6 classificação com random forest

3.7 classificação binária com árvores de decisão usando a função rpart

aula12

lsj

June 06, 2020

aula de hoje: regressão, classificação e uma introdução ao pacote caret

1. regressão

2. introdução ao ‘caret’

3. classificação

0. inicialização: bibliotecas para a aula de hoje

1. Regressão

1.1 otimização de um modelo por gradiente descendente

1.2 modelos e generalização

1.3 regularização: LASSO e ridge regression

1.4 regressão robusta

1.5 regressão com kernel

2 a biblioteca caret

2.1 o que é o ‘caret’

2.2 avaliação dos dados

2.3 pré-processamento dos dados

2.4 regressão com random forest

2.5 regressão com gbm (Generalized Boosted Regression Modeling) e ajuste de hiperparâmetros de modelos

3. Classificação

3.1 dados

3.2 classificação com k-nn

3.3 regressão logística

3.4 classificação com svm: support vector machines

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

3.6 classificação com random forest

3.7 classificação binária com árvores de decisão usando a função rpart

aula de hoje: regressão, classificação e uma introdução ao pacote `caret`

2 a biblioteca `caret`