aula 10

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

    1. o que é o caret
    1. avaliação dos dados
    1. pré-processamento dos dados
    1. treinamento e avaliação
    1. regressão: rf - random forest
    1. regressão: knn - kth nearest neighbor
    1. avaliação da importância relativa das variáveis
    1. ajuste de hiper-parâmetros de modelos
    1. classificação: preparação dos dados
    1. classificação: SVM (Support Vector Machines)
    1. classificação: árvore de decisão
    1. exercícios
# 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))
Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
logical.return = TRUE, : there is no package called '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 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).

2 avaliação dos dados

os dados são a base do ML e, portanto, deve-se saber se são adequados para o que se pretende fazer: há dados faltantes? os valores estão dentro do esperado?

é importante fazer uma avaliação cuidadosa! lembre-se: garbage in, garbage out!

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 None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
uJAVA_petro 0 1 20.10 1.21 14.95 19.36 20.05 20.76 31.02 ▁▇▁▁▁
F378_petro 0 1 19.85 1.17 12.94 19.15 19.81 20.50 28.69 ▁▃▇▁▁
F395_petro 0 1 19.69 1.22 12.72 18.96 19.66 20.37 29.16 ▁▅▇▁▁
F410_petro 0 1 19.50 1.20 11.89 18.79 19.48 20.19 28.37 ▁▂▇▁▁
F430_petro 0 1 19.25 1.19 11.66 18.54 19.25 19.97 29.13 ▁▃▇▁▁
g_petro 0 1 18.72 1.06 15.10 18.05 18.73 19.45 24.35 ▁▇▇▁▁
F515_petro 0 1 18.42 1.08 10.78 17.74 18.43 19.16 25.06 ▁▁▇▂▁
r_petro 0 1 17.97 1.02 15.01 17.30 18.01 18.76 19.99 ▁▃▇▇▅
F660_petro 0 1 17.84 1.03 10.70 17.16 17.89 18.65 21.11 ▁▁▂▇▂
i_petro 0 1 17.62 1.03 10.15 16.93 17.68 18.43 20.52 ▁▁▂▇▃
F861_petro 0 1 17.46 1.05 10.39 16.75 17.52 18.29 21.00 ▁▁▃▇▁
z_petro 0 1 17.42 1.06 9.83 16.70 17.48 18.26 20.96 ▁▁▂▇▁
z_SDSS 0 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 caret e algumas de suas principais funções, vamos restringir o número de objetos para diminuir o tempo de processamento. Mas note que em ML quanto mais dados para treinamento melhor!

nsample = 10000

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

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

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

# diagrama magnitude redshift
par(mfrow = c(1,2))
plot(zspec,mags$r_petro,xlab='z_SDSS',ylab='r_petro',main='',col=rgb(.2, .2, .2, .4))
smoothScatter(zspec,mags$r_petro,xlab='z_SDSS',ylab='r_petro',main='')

c. vamos visualizar as correlações entre o redshift e as várias magnitudes usando featurePlot:

theme1 <- trellis.par.get()
theme1$plot.symbol$col = rgb(.2, .2, .2, .4)
theme1$plot.symbol$pch = 16
theme1$plot.line$col = rgb(1, 0, 0, .7)
theme1$plot.line$lwd <- 2
trellis.par.set(theme1)
featurePlot(x = mags, y = zspec, plot = "scatter")

Note que cada banda fotométrica tem uma certa correlação com \(z\), contém informação sobre \(z\): os algoritmos de ML buscam capturar simultaneamente a informação contida em todas os atributos para obter a melhor estimativa possível da variável desejada, no caso \(z\).

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 para colocá-los num intervalo para tornar a análise mais eficiente.

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

  1. reescalonar as variáveis entre 0 e 1 (ou entre -1 e 1)
  1. subtrair a média de cada variável e dividir por seu desvio padrão (padronizaçã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

4 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(): ela embaralha os dados e seleciona a fração desejada
numlinhastreinamento =  createDataPartition(zspec, p=0.75, list=FALSE)

# índices das primeiras galáxias selecionadas:
head(numlinhastreinamento)
     Resample1
[1,]         1
[2,]         2
[3,]         4
[4,]         5
[5,]         6
[6,]         7
# conjunto de treinamento:
xtrain=magnorm[numlinhastreinamento,]
ytrain=zspec[numlinhastreinamento]

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

# tamanho dos conjuntos de treinamento e teste:
length(ytrain)
[1] 7502
length(ytest)
[1] 2498

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

Note que cada modelo pode ter hiperparâmetros que têm que ser definidos.

5 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

O algoritmo tem um hyperparâmetro, mtry, que é o número de atributos de cada árvore (cada árvore tem uma seleção aleatória de mtry dos 12 atributos)

# 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:
Trf = Sys.time() - t0
Trf
Time difference of 15.69473 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.05249076  0.6852403  0.03672094
   7    0.05155837  0.6937068  0.03551751
  12    0.05206127  0.6873302  0.03567882

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_rf = 0.7413*(quantile(relz,0.75,names = FALSE) - quantile(relz,0.25,names = FALSE))
sig_G_rf
[1] 0.03470638
#visualização do resultado
# vamos juntar os valores preditos e observados do redshift em um data frame:
my_data = as.data.frame(cbind(observed = ytest,predicted = pred_rf))
head(my_data)                            
   observed  predicted
3     0.086 0.08604740
8     0.159 0.14800790
16    0.057 0.06984820
21    0.078 0.07307017
22    0.187 0.13492240
36    0.077 0.09916043
# 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'

Os outliers são denominados redshifts catastróficos.

6 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
# alguns parâmetros foram introduzidos na train() para encontrar o melhor valor de k: trControl e tuneGrid
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:
Tknn = Sys.time()-t0
Tknn
Time difference of 3.318068 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.07073120  0.4748350  0.04823369
   3  0.05937686  0.5960307  0.04074287
   5  0.05751589  0.6207698  0.03941941
   7  0.05646660  0.6361129  0.03882913
   9  0.05618505  0.6417540  0.03867421
  11  0.05596326  0.6465358  0.03869639
  13  0.05603269  0.6468720  0.03874140
  15  0.05614280  0.6463015  0.03887838
  17  0.05619653  0.6467115  0.03894073
  19  0.05623716  0.6469274  0.03900855
  21  0.05631974  0.6467934  0.03908923
  23  0.05648782  0.6454102  0.03924854
  25  0.05661555  0.6445991  0.03941308
  27  0.05674109  0.6436650  0.03953907
  29  0.05683782  0.6431500  0.03962294
  31  0.05693148  0.6423334  0.03972244
  33  0.05707390  0.6409531  0.03981858
  35  0.05719980  0.6397221  0.03991602
  37  0.05725540  0.6395915  0.03996967
  39  0.05738066  0.6380351  0.04004858
  41  0.05750937  0.6367460  0.04016920

RMSE was used to select the optimal model using the smallest value.
The final value used for the model was k = 11.
# 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_knn = 0.7413*(quantile(relz,0.75,names = FALSE) - quantile(relz,0.25,names = FALSE))
sig_G_knn
[1] 0.04176552
#
# 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'

7 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:

  1. 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
  1. 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
  1. 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 28.71535 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.05476   0.6539 0.03830 0.003548    0.04260 0.002341         
         9 0.05186   0.6906 0.03574 0.002871    0.02800 0.001173         
        12 0.05078   0.7042 0.03482 0.002828    0.02719 0.001133        *

The top 5 variables (out of 12):
   g_petro, F430_petro, F395_petro, F515_petro, z_petro
# vejam quais são, nesse caso, as 5 magnitudes mais importantes

8 regressão: 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 do 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

Para determinar os melhores valores dos hiperparâmetros 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 hiperparâ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)
Tgbm = Sys.time() - t0
Tgbm
Time difference of 5.393274 mins
summary(model_gbm)

                    var   rel.inf
g_petro         g_petro 41.583626
F430_petro   F430_petro 13.186972
F861_petro   F861_petro  6.383008
F515_petro   F515_petro  6.179150
z_petro         z_petro  5.142561
F378_petro   F378_petro  4.713552
i_petro         i_petro  4.587749
uJAVA_petro uJAVA_petro  4.355976
F410_petro   F410_petro  4.088857
F395_petro   F395_petro  3.602628
r_petro         r_petro  3.282497
F660_petro   F660_petro  2.893424
# compare as 5 variáveis mais importantes com as obtidas com random forest

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

relz = (ytest - pred)
sig_G_gbm = 0.7413*(quantile(relz,0.75,names = FALSE) - quantile(relz,0.25,names = FALSE))
sig_G_gbm
[1] 0.03658794
#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'

Vamos comparar os resultados dos modelos rf, knn e gbm:

c(sig_G_rf,sig_G_knn,sig_G_gbm)
[1] 0.03470638 0.04176552 0.03658794

Nesse caso RF teve o melhor desempenho.

Compare os tempos de processamento dos modelos rf, knn e gbm:

c(Trf,Tknn,Tgbm)
Time differences in secs
[1] 941.683740   3.318068 323.596420

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

Vamos agora considerar problemas de classificação, com um conjunto de dados puramente fotométrico.

Em problemas de classificação é recomendável se usar amostras equilibradas, isto é, com o mesmo número de objetos em cada classe.

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
tabela = tabela[1:10000,]

# número de objetos em cada 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])

10 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 14.82921 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 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

O algoritmo dá estatísticas úteis para se avaliar o resultado:

  • acurácia: fração das classificações corretas \[{TP + TN \over TP + TN + FP + FN}\]
  • 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)
  • completeza (ou sensitivity): \[{TP \over TP+FN}\]
  • pureza da classe 1 (Pos Pred Value): \[{TP \over TP+FP}\]
  • pureza da classe 0 (Neg Pred Value): \[{TN \over TN+FN}\]

no exemplo, a acurácia é:

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

Note que se consegue uma excelente acurácia nesse caso.

11 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.4483337 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.125803 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, a relação entre as taxas de positivos verdadeiros e falsos positivos) 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 desta aula:

stopCluster(cl)

Sys.time() - t00
Time difference of 50.23232 mins

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

12 exercícios

  1. 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 5000 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

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