para registrar o tempo total de processamento deste script

t00 = Sys.time()

Este documento foi preparado usando o RStudio, um Integrated Development Environment (IDE) útil para o R.

plano da aula

  1. Introdução ao R
  1. ML: aprendizado não-supervisionado
    Métodos não-paramétricos - exploração de dados
  • 2.1 análise de agrupamento: K-means
  • 2.2 redução de dimensionalidade: PCA, UMAP

3 regressão e classificação

  • 3.1 otimização de modelos por descida do gradiente
  • 3.2 a biblioteca caret
  1. ML: aprendizado supervisionado - regressão
    redshifts fotométricos
  • 4.1 avaliação dos dados
  • 4.2 k-NN
  • 4.3 random forest
  1. ML: aprendizado supervisionado - classificaçãoclassificação estrela/galáxia
  • 5.1 preparação dos dados
  • 5.2 árvores de decisão

I. Introdução ao R

Luis Argerich, professor de Data Science na Universidad de Buenos Aires (UBA):

R is very different than other programming languages, you don’t learn R, you just use it.

This is very important and is something that most programmers will find puzzling. You don’t read a book about R control structures, functions and so on and start doing things to learn the language.

You have a problem, you find how to solve it using R and you learn whatever part of the language you need for it. Repeat this until you can do a lot of things with R. You will find the very curious case of people using R for years and being able to do a lot of things with R and they might not claim they know the language.

Introduções ao R

a introdução do CRAN (The Comprehensive R Archive Network):
https://cran.r-project.org/doc/manuals/R-intro.html

o grupo de Astroestatística da Penn State (https://sites.psu.edu/astrostatistics/) é muito influente - não deixem de olhar

tutorial do Eric Feigelson (Penn State):
https://sites.psu.edu/astrostatistics/files/2021/02/Tutorial_1_Intro_R.html

curso da Carnegie Mellon:
https://www.stat.cmu.edu/~cshalizi/statcomp/14/

CRAN Task Views: informações sobre bibliotecas de vários tipos CRAN task views aim to provide some guidance which packages on CRAN are relevant for tasks related to a certain topic. https://cran.r-project.org/web/views/

ache uma biblioteca ou uma função
it works basically like a filter on google searches, to return only R-relevant content.
https://rseek.org/

novidades no r-bloggers: https://www.r-bloggers.com/

Com R dá para fazer gráficos muito bonitos! Veja exemplos em

https://www.r-graph-gallery.com/index.html

para instalar o R no computador de sua casa: visite http://cran.r-project.org/

CRAN: The Comprehensive R Archive Network

R como linguagem

em R todos os comandos são funções que agem sobre objetos de vários tipos e produzem outros objetos

o usuário pode escrever suas próprias funções

* pode-se rodar R de várias formas:

1 criando um ambiente para rodar R em um terminal

2 usando rstudio (este documento foi feito com rstudio)

3 rodando um script em um terminal

4 como um jupyter notebook

* comando para criar um ambiente para rodar R em um terminal:

R

* comando para sair do ambiente R:

q()

* ajuda sobre uma função (no console):

?nome_da_função

exemplo: ajuda sobre a função plot: ?plot

R como calculadora:

Inicie uma seção em R digitando R no terminal e dando enter; agora você pode começar a entrar comandos:

2+2
## [1] 4
exp(-2)
## [1] 0.1353353
# pi
pi
## [1] 3.141593
# log na base 10
log10(pi)
## [1] 0.4971499
# logaritmo neperiano
log(pi)
## [1] 1.14473
# log na base 2
log(pi, 2) 
## [1] 1.651496
# arredondamento de uma resposta em 3 casas decimais:
round(log(pi),3)
## [1] 1.145

faça ?log para ver as opções da função log

sin(2*pi/3)
## [1] 0.8660254
sqrt(19)
## [1] 4.358899
# expoentes
13^(1/5)
## [1] 1.670278
# atribuindo um valor a uma variável: <- ou =
x <- 10
x = 10

# letras minúsculas e maiúsculas importam: x e' diferente de X
X <- x+x     

# para ver os valores de uma variável:
x
## [1] 10
X
## [1] 20
# nem todos os nomes sao possiveis: T e F são usados para as variáveis lógicas TRUE e FALSE

2*x+X^(1/3)
## [1] 22.71442
y <- x/exp(-x/X)
round(y,3)
## [1] 16.487

R: operadores

sequências e repetições

# sequências
1:5
## [1] 1 2 3 4 5
# colocando uma sequência em um vetor
y <- 1:10

#sequência de 1 a 77 com incrementos de 5 unidades
s <- seq(1,77,5)
s
##  [1]  1  6 11 16 21 26 31 36 41 46 51 56 61 66 71 76
# repetições:
rep(2,10)
##  [1] 2 2 2 2 2 2 2 2 2 2

Calcule a raiz quadrada da sequência de 3 a 9

vetores

# concatenando as componentes com o comando c  (concatenate)
vetor1 <- c(1.,-2.,-7.)
vetor1 <- c(vetor1,3.)
vetor1
## [1]  1 -2 -7  3
# acessando o segundo elemento do vetor
# os índices dos vetores em R começam em 1 (em python eles começam em 0)
vetor1[2]
## [1] -2
# note que, enquanto os elementos de uma função são envolvidos por parênteses, (...), os elementos de um vetor o são por chaves, [...]
# acessando os elementos 2 a 4
vetor1[2:4]
## [1] -2 -7  3
# acessando os elementos > 0
vetor1[vetor1 > 0]
## [1] 1 3
# operadores elemento a elemento: +, -, *, /, ^, etc
vetor2 <- -10+c(12,4,8,11)
vetor2
## [1]  2 -6 -2  1
# R é uma linguagem "vetorizada":
# operações elemento a elemento de cada vetor
vetor1+vetor2
## [1]  3 -8 -9  4
v3 <- vetor1/vetor2^2
v3
## [1]  0.25000000 -0.05555556 -1.75000000  3.00000000
vetor1*vetor2
## [1]  2 12 14  3
# produto escalar: usa-se %*%
vetor1%*%vetor2
##      [,1]
## [1,]   31
# funções podem ser aplicadas a todos os elementos de um vetor de uma única vez:
# criamos um vetor
x<-seq(0,1,0.1)
x
##  [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
# função aplicada ao vetor
y <- x*exp(-x)
y
##  [1] 0.00000000 0.09048374 0.16374615 0.22224547 0.26812802 0.30326533
##  [7] 0.32928698 0.34760971 0.35946317 0.36591269 0.36787944
# máximo e mínimo de um vetor
max(vetor1)
## [1] 3
min(vetor1)
## [1] -7
# dimensão do vetor
length(vetor1)
## [1] 4
# soma dos elementos
sum(vetor1)
## [1] -5
# soma cumulativa
cumsum(vetor1)
## [1]  1 -1 -8 -5
#produto dos elementos
prod(vetor1)
## [1] 42
# produto acumulado
cumprod(vetor1)
## [1]  1 -2 14 42
# ordenação dos elementos de um vetor
# vetor1
vetor1
## [1]  1 -2 -7  3
# ordenação dos elementos
sort(vetor1)
## [1] -7 -2  1  3
# índice da ordenação dos elementos
order(vetor1)
## [1] 3 2 1 4

matrizes

# vetor com 4 elementos
c(1:3, 8) 
## [1] 1 2 3 8
# criando matrizes com matrix()
# Matriz 2x2
matrix(c(1:3, 8), 2)
##      [,1] [,2]
## [1,]    1    3
## [2,]    2    8
# matriz 3x2
A <- matrix(c(1, 1, 1, 2, 4, 20),nrow=3,ncol=2) 
A
##      [,1] [,2]
## [1,]    1    2
## [2,]    1    4
## [3,]    1   20
# acessando os elementos de uma matriz (linha,coluna)
# para acessar o elemento (2,1):
A[2,1]
## [1] 1
# para acessar toda a coluna 1:
A[,1]
## [1] 1 1 1
# para acessar toda a linha 3:
A[3,]
## [1]  1 20
# para acessar os elementos das linhas 2 a 3 e da coluna 2:
A[2:3, 2]
## [1]  4 20
# para acessar todos os elementos maiores que 3
A[A > 3]
## [1]  4 20
# dimensão de uma matriz
dim(A)
## [1] 3 2

Considere a matriz M <- matrix(c(-2,-3,0,1, 3, 1, 2, 4, 20),nrow=3,ncol=3); verifique sua dimensão e o número de valores menores ou igual a zero

operadores matriciais:

Vamos começar criando duas matrizes e depois vamos fazer operações com elas:

A <- matrix(c(1, 1, 2, 4),2) #matriz 2x2
A
##      [,1] [,2]
## [1,]    1    2
## [2,]    1    4
B <- matrix(c(2, 1, 3, 1),2) #matriz 2x2
B
##      [,1] [,2]
## [1,]    2    3
## [2,]    1    1
# multiplicacao de matrizes  %*%
A%*%B 
##      [,1] [,2]
## [1,]    4    5
## [2,]    6    7
# matriz transposta
t(A)
##      [,1] [,2]
## [1,]    1    1
## [2,]    2    4
# matriz inversa
solve(B)
##      [,1] [,2]
## [1,]   -1    3
## [2,]    1   -2
# determinante de uma matriz
det(A)
## [1] 2
# teste de sanidade:
# multiplicação de uma matriz por sua inversa (com %*%)
B%*%solve(B)
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1

funções

definindo funções

sintaxe: nome = function(parâmetros da função){cálculo da função}

vamos considerar alguns exemplos:

duplica = function(u){2*u}
duplica(2)
## [1] 4
# distribuição de Cauchy: ela tem 2 parâmetros, a e b
cauchy = function(x,a,b){
b/(pi*(b^2+(x-a)^2))
}

# exemplo: valor da função para x=-1, a=3 e b=0.1
cauchy(-1,3,0.1)
## [1] 0.001988194
x = seq(-10,13,0.1)
plot(x,cauchy(x,3,1),type='l',ylab='f(x)',col='salmon',lwd=4)
abline(v=3)

plot de uma função

Vamos ver como plotar uma função em R:

# variável x: sequência de -10 a 13 com separação de 0.1
x = seq(-10,13,0.1)
y = cauchy(x,3,1)
# preparando um plot com 1 linha e 2 colunas:
par(mfrow=c(1,2))
# a função plot aceita muitos parâmetros
plot(x,y)
plot(x,y,type='l')

# para fazer uma figura sozinha:
par(mfrow=c(1,1))
# incluindo mais parâmetros:
plot(x,cauchy(x,3,1),type='l',xlab='x',ylab='f(x)',col='salmon',lwd=4)
# incluindo uma linha vertical em x=3:
abline(v=3,lwd=2)
# incluindo uma legenda:
legend('topleft',legend='função de Cauchy',text.col='blue',  bty='n',cex=1.5)

Há vários parâmetros para tornar os gráficos mais informativos:

O comando pch() define o tipo de símbolo que será usado na figura:

símbolos para pch
símbolos para pch

lwd: largura de linha (o default é 1)

type = ‘l’: plota uma linha juntando os pontos
há outras opções: faça ?plot

col: cor dos pontos ou linhas

log = ‘x’,‘y’,‘xy’: escala logarítmica

cex: tamanho dos símbolos

cex.axis, cex.lab, cex.main: magnificação em relação a cex das anotações dos eixos, dos labels x e y e do main (o título do plot)

lty: tipo de linha

lwd: largura da linha (default = 1)

e muito mais… veja:

https://www.statmethods.net/advgraphs/parameters.html

Alguns exemplos:

par(mfrow=c(2,2))
# variável x: sequência de 0 a 10 com separação de 0.1
x = seq(0,10,0.1)
y = exp(x*sin(x)^2)

# tipo de linha
plot(x,y,xlab='x',ylab='y',type='l',col='chocolate1',col.lab='blue',lty=3)
plot(x,y,xlab='x',ylab='y',type='l',col='chocolate1',col.lab='blue',lty=6)

# aumento da largura da linha e dos labels
plot(x,y,xlab='x',ylab='y',type='l',col='chocolate1',lwd=2,cex.lab=1.5,col.lab='blue')

# para por o eixo y em log + mudando o tipo de linha
plot(x,y,xlab='x',ylab='y',type='h',log='y',col='chocolate1',lwd=2,cex.lab=1.5,col.lab='blue')

funções pré-definidas

uma das vantagens do R é que ele tem um grande número de funções pré-definidas

exemplo: simulação de 10 números aleatórios uniformemente distribuidos entre 0 e 7 com a função runif()

reprodutibilidade (importante!):

reprodutibilidade - para garantir que os resultados desta simulação possam ser reproduzidos, deve-se definir uma “semente” para gerar números aleatórios:

# usamos a função set.seed()
set.seed(42)

# a função runif() permite gerar números com distribuição uniforme
# simulando 10 números com distribuição uniforme entre 0 e 7:
x = runif(10,min=0,max=7)  
x
##  [1] 6.4036423 6.5595279 2.0029767 5.8131334 4.4922186 3.6336716 5.1561182
##  [8] 0.9426662 4.5989460 4.9354535
# se você gerar agora uma outra sequência, ela será diferente da anterior: verifique!

# simule duas sequências, estabelecendo antes de cada simulação a mesma seed

loops

loops em R são lentos, e sempre que possível deve-se procurar trabalhar com vetorização, isto é, com expressões que envolvam apenas operações com vetores e matrizes

loops em R: for, while, repeat

# para ilustrar vamos considerar dois vetores, que escreveremos usando *vetorização*:
x = c(1,2,3,4,5)
# dimensão do vetor x:
n = length(x)
n
## [1] 5
# vetor y usando vetorialização
y = 2*x
y
## [1]  2  4  6  8 10
# agora vetor y, usando loop com for()
y = NULL # para definir um objeto y
# loop:
for(i in 1:n){
y[i]=2*x[i]
}
y
## [1]  2  4  6  8 10
# exemplo de um loop com while()
x = 1
while (x < 6) {
x = 2*x
print(x)
}
## [1] 2
## [1] 4
## [1] 8
# exemplo de um loop usando repeat() 
x = 1
repeat {
x = 2*x
print(x)
if (x == 8){break}
}
## [1] 2
## [1] 4
## [1] 8
# sai-se do loop com 'break'

combinação de variáveis em tabelas

já vimos como concatenar variáveis com o comando c(), concatenate. Por exemplo, podemos criar três vetores u, v e z de dimensão 6 como

u = c(3,-4,1,0,1,-1)
v = c(1,-2,5,-1,0,0)
z = c(0,2,2,-2,1,1)

# tamanhos dos vetores:
c(length(u),length(v),length(z))
## [1] 6 6 6

Se quisermos criar uma tabela (matriz) tendo como colunas u, v e z usamos cbind()
(column bind):

tab1 = cbind(u,v,z)
tab1
##       u  v  z
## [1,]  3  1  0
## [2,] -4 -2  2
## [3,]  1  5  2
## [4,]  0 -1 -2
## [5,]  1  0  1
## [6,] -1  0  1
# dimensão da tabela/matriz
dim(tab1)
## [1] 6 3
# nomes das colunas
colnames(tab1)
## [1] "u" "v" "z"
# vamos atribuir novos nomes às colunas: note que os nomes são variáveis nominais
colnames(tab1) = c('var_u','var_v', 'var_z')
tab1
##      var_u var_v var_z
## [1,]     3     1     0
## [2,]    -4    -2     2
## [3,]     1     5     2
## [4,]     0    -1    -2
## [5,]     1     0     1
## [6,]    -1     0     1

Podemos também criar tabelas juntando linhas. Por exemplo, juntando por linha os vetores associados a cada coluna de tab1 (u, v e z):

tab2 = rbind(u,v,z) # row bind
tab2
##   [,1] [,2] [,3] [,4] [,5] [,6]
## u    3   -4    1    0    1   -1
## v    1   -2    5   -1    0    0
## z    0    2    2   -2    1    1
# atribuindo nomes às linhas:
rownames(tab2) = c('var_u','var_v','var_z')
tab2
##       [,1] [,2] [,3] [,4] [,5] [,6]
## var_u    3   -4    1    0    1   -1
## var_v    1   -2    5   -1    0    0
## var_z    0    2    2   -2    1    1

pode-se aplicar funções para variáveis em tabelas com apply():

# vamos ilustrar calculando médias de linhas e colunas de uma tabela
x = 1:5
y = 2^x
# vamos criar uma tabela com os x e y acima
tab = cbind(x,y)  # column bind
tab
##      x  y
## [1,] 1  2
## [2,] 2  4
## [3,] 3  8
## [4,] 4 16
## [5,] 5 32
# média das colunas
apply(tab,2,mean)
##    x    y 
##  3.0 12.4
# média das linhas
apply(tab,1,mean)
## [1]  1.5  3.0  5.5 10.0 18.5
# pode-se usar outras funções
apply(tab,2,duplica)
##       x  y
## [1,]  2  4
## [2,]  4  8
## [3,]  6 16
## [4,]  8 32
## [5,] 10 64

métodos numéricos

Em R existem rotinas e bibliotecas prontas para inúmeros métodos numéricos que usamos em nosso dia a dia como cientistas: vamos ilustrar com a integração numérica de uma função:

f = function(x) 1/sqrt(x)
integrate(f, 0, 1)
## 2 with absolute error < 5.8e-15
# note que a saida da função integrate() é uma lista
# e, para acessar apenas o resultado da integral, pode-se fazer
integrate(f, 0, 1)$value
## [1] 2
# outro exemplo:
fgauss = function(x) exp(-x^2/2)
I = integrate(fgauss, -Inf, Inf)
I
## 2.506628 with absolute error < 0.00023
# I é uma variável tipo lista
# o valor numérico da integral é dado por
I$value
## [1] 2.506628
# note que o valor exato dessa integral é sqrt(2*pi)
I$value/sqrt(2*pi)
## [1] 1

pacotes e bibliotecas

O R tem uma infinidade de pacotes e bibliotecas especializadas. Para instalá-las use o comando install.packages(‘nome da biblioteca’).

Para usar uma biblioteca no R, antes de rodar o comando desejado, use o comando library(nome da biblioteca). Veremos alguns exemplos mais abaixo.

2. Machine Learning:

aprendizado supervisionado x aprendizado não-supervisionado x aprendizado com reforço

aprendizado: determinação dos parâmetros de um modelo a partir dos dados

aprendizado não-supevisionado

útil para exploração dos dados

1 análise de agrupamento: K-means

2 redução de dimensionalidade: PCA, UMAP

2.1 análise de agrupamento: K-means

objetivo: encontrar grupos (clusters) no espaço de dados

pense no conceito de espaço de dados como tendo estruturas que se podem descobrir e explora

grupos: objetos com propriedades similares

técnica não-supervisionada: os grupos não são conhecidos a priori

diferente da classificação, que é um procedimento supervisionado onde os objetos estão associados a grupos pré-definidos

exemplos de métodos:

  • algoritmos combinatórios: K-means, hierarchical clustering
  • modelos probabilísticos: misturas gaussianas,…

para identificar grupos, considera-se as distâncias entre os dados

tipos de distâncias entre 2 objetos (\(\mathbf{x}_i\) e \(\mathbf{x}_j\)):

  • manhattan: \(d_{ij} = \sum_{k=1}^M |\mathbf{x}_i-\mathbf{x}_j|\)
  • euclidiana: \(d_{ij} = \Big[\sum_{k=1}^M |\mathbf{x}_i-\mathbf{x}_j|^2\Big]^{1/2}\)
  • mahalanobis: \(d_{ij} = \Big[(\mathbf{x}_i - \mathbf{x}_j)^T.\mathbf{C}^{-1}.(\mathbf{x}_j-\mathbf{x}_j) \Big]^{1/2}\), onde \(C\) é a matriz de covariância

note que para que estas distâncias tenham sentido é necessário que as variáveis envolvidas sejam de escalas semelhantes: muitas vezes é necessário algum tipo de normalização dos dados

algoritmo K-means:**

  • objetivo: atribuir N objetos a K grupos, K\(< <\)N \(\longrightarrow\) K é dado!
  • custo: o algoritmo procura, iterativamente, minimizar a soma \[ J = \sum_{i=1}^N \sum_{k=1}^K (\mathbf{x}_{i}-\mathbf{\mu}_{k})^2 = \sum_{i=1}^N \sum_{j=1}^M \sum_{k=1}^K (x_{ij}-\mu_{jk})^2\] onde \(\mathbf{\mu}_k\) é o centróide do grupo K

algoritmo:

    1. defina o número de grupos, K
    1. inicialize, escolhendo K objetos ao acaso como centro
    1. atribua cada objeto ao grupo mais próximo e recalcule o centróide do grupo

repita 3 até a convergência (os grupos ‘congelam’)

Note que:

  • K deve ser dado a priori
  • o algoritmo é sensível à inicialização: diferentes ‘rodadas’ podem levar a grupos diferentes
  • os resultados podem depender do tipo de distância escolhida

Vamos considerar inicialmente, para ilustrar, dados simulados em duas dimensões:

# paralelização
library(doParallel)
## Loading required package: foreach
## Loading required package: iterators
## Loading required package: parallel
# vou sar 5 cores de meu laptop para rodar este script
cl <- makePSOCKcluster(5)
registerDoParallel(cl)
# o uso de paralelização pode reduzir significativamente o tempo de processamento 

# reprodutibilidade: MUITO IMPORTANTE
set.seed(123)

# vamos gerar coordenadas (x,y) para dois conjuntos de dados
# grupo 1: mu = (-1,-1)    sd = 0.5 
mu1 = c(-1,-1)
sig = 0.5
grupo1 = cbind(rnorm(50,mu1[1],sig), rnorm(50,mu1[2],sig))

# grupo 2: mu = (1,1)    sd = 0.5 
mu2 = c(1,1)
sig = 0.5
grupo2 = cbind(rnorm(50,mu2[1],sig), rnorm(50,mu2[2],sig))

# vamos juntar os dois conjuntos de dados em um único conjunto:
dados = rbind(grupo1, grupo2)

#visualização
par(mfrow = c(1,1))
plot(dados, pch = 16, xlab = "x", ylab = "y", xlim=c(-5,5),ylim=c(-5,5))

# a função kmeans do R permite rodar o algoritmo com várias inicializações aleatórias diferentes, retornando o resultado de menor custo J
# 2 grupos e 20 reinicializações:
kmres = kmeans(dados,2, nstart = 20)

# centros dos grupos: compare com mu1 e mu2
kmres$centers[1,]
## [1] -0.9827982 -0.9267959
kmres$centers[2,]
## [1] 0.8730498 1.0194034
#visualização
par(mfrow = c(1,1))
plot(dados, pch = 16, xlab = "x", ylab = "y", col = kmres$cluster,xlim=c(-5,5),ylim=c(-5,5))
points(kmres$centers[,1], kmres$centers[,2],pch=4,lwd=4, col="yellow", cex = 1.5)
points(mu1[1],mu1[2],pch=4,lwd=4, col="green", cex = 1.5)
points(mu2[1],mu2[2],pch=4,lwd=4, col="green", cex = 1.5)
legend('topleft', col=c("green", 'yellow'), lwd=2, legend=c("input", "k-means"))

# número de objetos em cada grupo:
kmres$size
## [1] 50 50
# classificação de cada objeto:
kmres$cluster
##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
#######################################################

# vamos ver como fica com sd = 1.5, que mistura os grupos:
sig = 1.5
grupo1 = cbind(rnorm(50,mu1[1],sig), rnorm(50,mu1[2],sig))
grupo2 = cbind(rnorm(50,mu2[1],sig), rnorm(50,mu2[2],sig))
# dados
dados = rbind(grupo1, grupo2)

plot(dados, pch = 16, xlab = "x", ylab = "y", xlim=c(-5,5),ylim=c(-5,5))

# função kmeans
# 2 grupos e 20 reinicializações:
kmres = kmeans(dados,2, nstart = 20)

# centros dos grupos
kmres$centers[1,]
## [1] -1.5149972 -0.5467696
kmres$centers[2,]
## [1] 1.4648223 0.4836374
# plots
plot(dados, pch = 16, xlab = "x", ylab = "y", col = kmres$cluster,xlim=c(-5,5),ylim=c(-5,5))
points(kmres$centers[,1], kmres$centers[,2],pch=4,lwd=4, col="blue", cex = 1.5)
points(mu1[1],mu1[2],pch=4,lwd=4, col="green", cex = 1.5)
points(mu2[1],mu2[2],pch=4,lwd=4, col="green", cex = 1.5)
legend('topleft', col=c("green", 'blue'), lwd=2, legend=c("input", "k-means"))

# número de objetos em cada grupo:
kmres$size
## [1] 47 53
# classificação de cada objeto:
kmres$cluster
##   [1] 1 1 1 2 1 2 1 1 1 1 1 1 2 1 1 1 1 2 1 1 2 1 1 1 2 1 1 1 1 2 1 1 1 1 1 1 1
##  [38] 1 2 2 2 1 1 1 1 2 2 1 1 1 2 2 2 2 2 2 1 2 2 2 1 1 2 2 2 2 2 2 2 2 1 2 2 1
##  [75] 2 2 1 2 2 2 2 2 2 1 2 2 2 2 2 2 1 2 2 2 2 2 1 2 2 2

Aplicação: análise de um conjunto de dados puramente fotométrico.

dados do S-PLUS: fotometria em 12 bandas para 10000 objetos

# vamos verificar o tempo de processamento deste módulo
t0 = Sys.time()

# leitura dod dados:
tabela = 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
# uma função muito útil
summary(tabela)
##     u_petro       J0378_petro     J0395_petro     J0410_petro   
##  Min.   :13.66   Min.   :13.70   Min.   :13.84   Min.   :13.31  
##  1st Qu.:17.51   1st Qu.:17.15   1st Qu.:16.99   1st Qu.:16.57  
##  Median :18.33   Median :18.01   Median :17.75   Median :17.38  
##  Mean   :18.42   Mean   :18.05   Mean   :17.82   Mean   :17.43  
##  3rd Qu.:19.26   3rd Qu.:18.89   3rd Qu.:18.63   3rd Qu.:18.31  
##  Max.   :24.06   Max.   :24.32   Max.   :24.96   Max.   :24.32  
##   J0430_petro       g_petro       J0515_petro       r_petro     
##  Min.   :13.10   Min.   :12.62   Min.   :12.44   Min.   :11.92  
##  1st Qu.:16.43   1st Qu.:16.09   1st Qu.:15.90   1st Qu.:15.54  
##  Median :17.09   Median :16.64   Median :16.43   Median :15.95  
##  Mean   :17.18   Mean   :16.69   Mean   :16.47   Mean   :16.00  
##  3rd Qu.:17.93   3rd Qu.:17.35   3rd Qu.:17.08   3rd Qu.:16.54  
##  Max.   :21.25   Max.   :19.72   Max.   :19.34   Max.   :18.34  
##   J0660_petro       i_petro       J0861_petro       z_petro          classe   
##  Min.   :11.83   Min.   :11.54   Min.   :11.40   Min.   :11.31   Min.   :0.0  
##  1st Qu.:15.49   1st Qu.:15.27   1st Qu.:15.15   1st Qu.:15.11   1st Qu.:0.0  
##  Median :15.89   Median :15.65   Median :15.55   Median :15.51   Median :0.5  
##  Mean   :15.92   Mean   :15.67   Mean   :15.55   Mean   :15.50   Mean   :0.5  
##  3rd Qu.:16.43   3rd Qu.:16.14   3rd Qu.:15.99   3rd Qu.:15.94   3rd Qu.:1.0  
##  Max.   :33.27   Max.   :17.72   Max.   :27.95   Max.   :17.68   Max.   :1.0
classe = tabela[,13]

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

Vamos examinar os dados em um diagrama cor-cor:

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 as variáveis de interesse (as magnitudes) e pré-processá-las:

# definindo as variáveis
mag = as.data.frame(tabela[,1:12])

# magnitudes do primeiro objeto:
mag[1,]
##    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
##    r_petro J0660_petro  i_petro J0861_petro  z_petro
## 1 16.53974    16.47142 16.21496    16.00299 15.91366
# pré-processamento: 
# padronização das variáveis: scale(x, center = TRUE, scale = TRUE) - subtrai a média e divide pelo desvio padrão
# ou
# normalização das magnitudes entre 0 e 1
maxs <- apply(mag, 2, max) 
mins <- apply(mag, 2, min)

# dados normalizados: x_norm = (x-xmin)/(xmax-xmin)
dados <- as.data.frame(scale(mag, center = mins, scale = maxs - mins))

# magnitudes normalizadas do primeiro objeto:
dados[1,]
##     u_petro J0378_petro J0395_petro J0410_petro J0430_petro   g_petro
## 1 0.4725361   0.4038829   0.3765649   0.4351743   0.5920056 0.6565072
##   J0515_petro   r_petro J0660_petro   i_petro J0861_petro  z_petro
## 1    0.660737 0.7196676   0.2165281 0.7568543   0.2780254 0.722876
# função kmeans
# 2 grupos e 20 reinicializações:
kmres = kmeans(dados,2, nstart = 20)

# centros dos grupos: 
kmres$centers[1,]
##     u_petro J0378_petro J0395_petro J0410_petro J0430_petro     g_petro 
##   0.3740121   0.3287923   0.2833960   0.2975459   0.4034623   0.4807753 
## J0515_petro     r_petro J0660_petro     i_petro J0861_petro     z_petro 
##   0.4940441   0.5564046   0.1693070   0.6049207   0.2298577   0.6060853
kmres$centers[2,]
##     u_petro J0378_petro J0395_petro J0410_petro J0430_petro     g_petro 
##   0.5515558   0.4994131   0.4426176   0.4610949   0.6080531   0.6755277 
## J0515_petro     r_petro J0660_petro     i_petro J0861_petro     z_petro 
##   0.6841392   0.7245659   0.2151971   0.7411086   0.2736645   0.7159997
# número de objetos em cada grupo:
kmres$size
## [1] 5282 4718
#visualização no diagrama cor-cor
par(mfrow = c(1,1))
ca = dados$g_petro-dados$r_petro
cb = dados$r_petro-dados$i_petro
smoothScatter(ca,cb,nrpoints=0,add=FALSE,xlab="g-r",ylab="r-i",main='todos')

# centro dos grupos
x1 = kmres$centers[1,6]-kmres$centers[1,8]
y1 = kmres$centers[1,8]-kmres$centers[1,10]
x2 = kmres$centers[2,6]-kmres$centers[2,8]
y2 = kmres$centers[2,8]-kmres$centers[2,10]
points(x1,y1,pch=4,lwd=4, col="yellow", cex = 1.5)
points(x2,y2,pch=4,lwd=4, col="green", cex = 1.5)

# estrelas e galáxias
par(mfrow = c(1,2))
smoothScatter(ca[classe == 0],cb[classe == 0],nrpoints=0,add=FALSE,xlab="g-r",ylab="r-i",main='estrelas')
smoothScatter(ca[classe == 1],cb[classe == 1],nrpoints=0,add=FALSE,xlab="g-r",ylab="r-i",main='galáxias')

# classificação do kmeans dos 10 primeiros objetos:
kmres$cluster[1:10]
##  [1] 2 1 2 1 2 1 1 2 2 1
# classes dos 10 primeiros objetos:
classe[1:10]
##  [1] 1 1 1 1 0 0 0 0 1 0
length(classe[classe == 0 & kmres$cluster == 1])
## [1] 3598
length(classe[classe == 0 & kmres$cluster == 2])
## [1] 1402
length(classe[classe == 1 & kmres$cluster == 1])
## [1] 1684
length(classe[classe == 1 & kmres$cluster == 2])
## [1] 3316
# tempo de processamento deste módulo:
Sys.time()-t0
## Time difference of 0.6501188 secs

2.2 redução de dimensionalidade: PCA, UMAP

útil para compressão de dados e visualização

X: dados em D dimensões

queremos uma nova representação de X, que chamaremos Y, em d \(<<\) D dimensões
\[X = {x_1, x_2, ..., x_D} \longrightarrow Y = {y_1, y_2, ..., y_d}\]

o espaço reduzido é denominado espaço latente

método linear: análise de componentes principais (PCA)

métodos não-lineares: LLE (locally linear embedding), IsoMap, t-SNE (t-distributed stochastic neighbor embedding), UMAP (uniform manifold approximation and projection)

embedding: uma estrutura matemática está contida dentro de outra

pense no conceito de embedding e como ele pode ser útil

PCA

como funciona? PCA cria um novo sistema de coordenadas ortonormal no espaço de dados:

    1. no espaço dos dados identifique a direção que apresenta a maior variância dos dados: esta é a direção da primeira componente principal, \(PC_1\)
    1. no sub-espaço perpendicular a \(PC_1\) identifique a direção que apresenta a maior variância dos dados: esta é a direção da segunda componente principal, \(PC_2\)
    1. no sub-espaço perpendicular a \(PC_1\) e \(PC_2\) identifique…
    1. repita o procedimento até \(PC_D\), onde \(D\) é a dimensão do espaço de dados

em muitos casos, poucas componentes são responsáveis pela maior parte da variância

isso significa que os dados podem ser aproximadamente descritos com um número \(d << D\) de componentes.

matematicamente:

vamos considerar um conjunto de N objetos, cada um descrito por D variáveis:

\(x_{ij}\): variável \(j\) do \(i\)-ésimo objeto

um objeto é um ponto no espaço de dados \(D\)-dimensional

matriz de covariância dos dados: matriz \(D \times D\) definida como \[ C_{jk} = \frac{1}{N} \sum_{i=1}^N (x_{ij}-{\bar x}_j) (x_{ik}-{\bar x}_k) \] onde \[ {\bar x}_j = \frac{1}{N} \sum_{i=1}^N x_{ij} ~~~~~~~~~j=1,...,D \] note que as variâncias dos dados correspondem à diagonal, \(C_{jj}\)

muitas vezes as medidas das variáveis são incompatíveis (como distâncias e magnitudes):

para “resolver” isso padronizam-se as variáveis, convertendo-as em quantidades de média nula e variância unitária: \[ x_{ij}^\prime = \frac{(x_{ij}-{\bar x}_j)}{\sqrt{C_{jj}}} \] e constrói-se a matriz de covariância com as variáveis padronizadas

cada componente principal é descrita como uma combinação linear das variáveis padronizadas: \[ PC_j = \sum_{i=1}^N \sum_{k=1}^D a_{jki} x_{ik}^\prime = \sum_{i=1}^N \mathbf{a}_{ji} \mathbf{x_i^\prime} \]

os \(\mathbf{a}_j\) são os autovetores de \(C\): \[ C\mathbf{a}_j = \lambda_j \mathbf{a}_j\]

o autovetor de \(C\) com o j-ésimo maior autovalor \(\lambda_j\) é a j-ésima componente principal e satisfaz (como as demais componentes):
\[\mathbf{a}_j\mathbf{a}_k^\prime = \delta_{jk}\] (normalização e ortogonalidade; \(\mathbf{a}_k^\prime\) é o vetor transposto)

os autovalores \(\lambda_j\) são proporcionais à variância de cada componente

assim, para descrever os dados pode-se selecionar o número de componentes em termos da fração da variância total que se deseja “explicar”

problemas com o PCA

  • as PC não têm necessariamente significado físico
  • a seleção do número de componentes de interesse é muitas vezes subjetiva
  • a análise depende do tipo de normalização aplicada aos dados
  • é um método linear

Exemplo: dados fotométricos (mag)

# vamos verificar o tempo de processamento deste módulo
t0 = Sys.time()

# para o PCA vamos escalonar as variáveis: média zero e variância unitária
pca = prcomp(mag, retx=TRUE, center=TRUE, scale=TRUE)

# raiz quadrada dos autovalores
sd = pca$sdev

# autovalores
lambda = sd^2

# fraçao da variância explicada
ve = round(cumsum(lambda)/sum(lambda),digits = 3)
ve[1:5]
## [1] 0.812 0.972 0.980 0.987 0.991
# vemos que a primeira componente sozinha explica 81% da variância e as duas primeiras 97%!

# visualização: scree plot - variância associada a cada componente
screeplot(pca, main=" scree plot")

# isso significa que 94% da variância do espaço espectral de 2499 dimensões está contida em um plano definido pelas duas primeiras componentes principais!

# autovetores
autovec = pca$rotation

# componentes principais
pc = pca$x

# visualização: duas primeiras componentes
par(mfrow = c(1,1))
smoothScatter(pc[,1], pc[,2], add=FALSE, xlab="PC 1", ylab="PC 2", main='')

# PCA: visualização do espaço de dados a partir de um referencial especial:
par(mfrow = c(2,2))
smoothScatter(pc[,1], pc[,2], xlab="PC 1", ylab="PC 2",col='blue')
smoothScatter(pc[,1], pc[,3], xlab="PC 1", ylab="PC 3",col='blue')
smoothScatter(pc[,2], pc[,3], xlab="PC 2", ylab="PC 3",col='blue')
smoothScatter(pc[,3], pc[,4], xlab="PC 3", ylab="PC 4",col='blue',ylim=c(-1,1))

# note que estes gráficos não estão em escala: compare os intervalos de valores nos dois eixos

# duração deste módulo
Sys.time()-t0
## Time difference of 0.2199316 secs

Importante: a interpretação das variáveis no espaço latente não é óbvia!

UMAP: Uniform manifold approximation and projection

# biblioteca
library(umap)

# vamos verificar o tempo de processamento deste módulo
t0 = Sys.time()

# UMAP em 2 dimensões
umap_fot = umap(dados, n_components = 2, random_state = 15) 

layout <- umap_fot[["layout"]] 

layout <- data.frame(layout) 

plot(layout)

plot(layout,col=factor(classe),xlab='X1',ylab='X2', pch=16, cex=2, cex.lab=1.5) 
legend("topleft",
       legend = levels(factor(classe)),
       pch = 19,
       col = factor(levels(factor(classe))))

smoothScatter(layout[,1], layout[,2], xlab="X1", ylab="X2",col='blue')

# UMAP em 3 dimensões
umap_fot = umap(dados, n_components = 3, random_state = 15) 

layout <- umap_fot[["layout"]] 

layout <- data.frame(layout) 

plot(layout[,1],layout[,2],col=factor(classe),xlab='X1',ylab='X2', pch=16, cex=2, cex.lab=1.5) 
legend("topleft",
       legend = levels(factor(classe)),
       pch = 19,
       col = factor(levels(factor(classe))))

plot(layout[,1],layout[,3],col=factor(classe),xlab='X1',ylab='X3', pch=16, cex=2, cex.lab=1.5) 
legend("topleft",
       legend = levels(factor(classe)),
       pch = 19,
       col = factor(levels(factor(classe))))

plot(layout[,2],layout[,3],col=factor(classe),xlab='X2',ylab='X3', pch=16, cex=2, cex.lab=1.5) 
legend("topleft",
       legend = levels(factor(classe)),
       pch = 19,
       col = factor(levels(factor(classe))))

# tempo de processamento:
Sys.time()-t0
## Time difference of 1.722094 mins

Note que a morfologia na distribuição dos pontos no espaço latente obtida com diferentes métodos é diferente e pode ser explorada para descobrir relações ocultas entre os dados ou para se fazer classificação de objetos.

vocês conhecem problemas onde essas técnicas poderiam ser úteis?

3. regressão e classificação com a biblioteca caret

caret: Classification and Regression Training

regressão e classificação: exemplos de aprendizado supervisionado

aprendizado supervisionado:

  • temos variáveis de entrada e saída, x e y, e o algoritmo aprende um mapeamento de x em y usando exemplos de x em um conjunto de treinamento para o qual os valores de y são conhecidos

x e y podem ser escalares, vetores, matrizes, tensores, imagens, …

  • classificação- y é uma variável qualitativa ou categórica ou discreta:
    ’estrela’, ‘galáxia Scd’, ’detectado’, ’classe 3’, …
  • regressão: y é um número real:
    redshift, metalicidade, massa …

ML: implementa uma função y=f(x;w), com parâmetros w;
w é determinado usando técnicas de descida do gradiente e o conjunto de treinamento

  • define-se uma função custo:
    p.exemplo, para regressão, \[ l(w) = {\sum_i (y^c_i (x_i|w)-y_i)^2 \over N}\] \(y_i\) são os dados do conjunto de treinamento e \(y^c_i (x_i|w)\) é o valor calculado para o valor de entrada \(x_i\) pelo modelo de ML com parâmetros w
  • treinamento: minimiza-se \(l(w)\) com descida do gradiente
# para monitorar o tempo de processamento:
t0 = Sys.time()

# reprodutibilidade
set.seed(111)

# bibliotecas
suppressMessages(library(caret))
suppressMessages(library(corrplot))
suppressMessages(library(skimr))
suppressMessages(library(ggplot2))
suppressMessages(library(rpart))
suppressMessages(library(rpart.plot))

3.1 otimização de um modelo por gradiente descendente

Este algoritmo é, na essência, o motor dos métodos de aprendizado de máquina, usado para minimizar a função de custo

Inventado por Cauchy (1847) como ferramenta para determinar os 6 parâmetros orbitais que descrevem o movimento dos planetas

Vamos ilustrar este algoritmo encontrando os parâmetros \(\mathbf{w} = \{w_0,w_1\}\) de um modelo linear simples: \(y(x;w)=w_0+w_1 x + \epsilon\), onde \(\epsilon\) se refere aos erros em \(y\)

A função de custo é a soma dos quadrados dos resíduos: \[l(w) \propto \chi^2 = \sum_{i=1}^N \frac{[y_i-(w_0 +w_1 x_i)]^2}{2 \epsilon^2}\] e a otimização é equivalente à minimização do \(\chi^2\)

Na aprendizagem descendo o gradiente, o algoritmo pode ser inicializado com valores aleatórios de \(\mathbf{w}\) e o aprendizado é baseado no gradiente local da função custo

Em cada iteração o algoritmo muda \(\mathbf{w}\) como \[ \mathbf{w} \longleftarrow \mathbf{w} - \lambda \times \nabla l(\mathbf{w})\] \(\lambda\): taxa de aprendizado
\(\lambda\) é um hiperparâmetro, um parâmetro associado ao algoritmo de aprendizagem e não ao modelo propriamente dito

Para um único dado \((x_i,y_i)\), \[\mathbf{w} \longleftarrow \mathbf{w} + 2\lambda [y_i - y(x_i ;\mathbf{w})]\] Note que o termo \([...]\) é o resíduo da medida \(i\) com o valor atual de \(\mathbf{w}\)

Conforme o treinamento prossegue (o número de iterações aumenta), os resíduos diminuem e convergem para um valor estacionário

Vamos ilustrar com um exemplo numérico:

# 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

# reprodutibilidade
set.seed(1)

# número de pontos
npt=100
# valores de x uniformemente distribuidos entre -1 e +1
x <- runif(npt, -1, 1)
# simulação de y
y <- 1+2*x + rnorm(npt,sd=0.2)

# vamos primeiro ajustar um modelo linear usando a função lm():
modelo_linear <- lm( y ~ x )
print(modelo_linear)
## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##      0.9954       2.0312
# visualização:
plot(x,y, col=rgb(0.2,0.4,0.6,0.4), main='regressão linear')
abline(modelo_linear, col='blue')

# função custo: erro quadrático médio/2
# vou usar uma notação matricial
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

# vamos guardar o custo e o valor dos parâmetros em cada iteração:
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) {
  erro <- (X %*% w - y)
  delta <- t(X) %*% erro / length(y)
  w <- w - alpha * delta
  custo_historia[i] <- custo(X, y, w)
  w_historia[[i]] <- w
}

# resultado: valor dos parâmetros
print(w)
##          [,1]
## [1,] 1.001105
## [2,] 1.915090
# visualização do processo de otimizaçãor:
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 evolução da função custo:
plot(custo_historia, type='line', col='blue', lwd=3, main='evolução da função custo', ylab='custo', xlab='Iterações')
## Warning in plot.xy(xy, type, ...): plot type 'line' will be truncated to first
## character

# tempo de processamento deste módulo:
Sys.time()-t0
## Time difference of 1.420139 secs

3.2 o que é o caret

Vamos 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).

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

modelLookup('pcaNNet')
##     model parameter         label forReg forClass probModel
## 1 pcaNNet      size #Hidden Units   TRUE     TRUE      TRUE
## 2 pcaNNet     decay  Weight Decay   TRUE     TRUE      TRUE
# pcaNNet: run PCA on a dataset, then use it in a neural network model

examinem alguns modelos do caret

4. regressão

estimativa de redshifts fotométricos

4.1 avaliação dos dados para regressão

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 pacote e algumas de suas principais funções, vamos restringir o número de objetos:

nsample = 10000

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

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

b. vamos dar uma olhada nos dados:

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

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

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

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

redshift versus cada banda fotométrica:

visualização usando featurePlot:

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

d. 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
  1. subtrair a média de cada variável e dividir por seu desvio padrão

vamos colocar as 12 magnitudes entre 0 e 1:

# valores máximos e mínimos de cada banda
maxs <- apply(mags, 2, max) 
mins <- apply(mags, 2, min)
# reescalonamento:
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

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

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
numlinhastreinamento =  createDataPartition(zspec, p=0.75, list=FALSE)

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

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

4.1 regressão com k-NN

k-NN: algoritmo k-th nearest neighbor

  • y(x) é a média ou a média ponderada dos valores de y dos k pontos mais próximos de x
  • útil para regressão e classificação
  • hiperparâmetro k: pode ser obtido por validação cruzada

implementa a intuição de que valores próximos de x devem ter valores próximos de y

k-NN é um modelo intuitivo e rápido; não é muito preciso mas é bom para se ter uma base line do que se pode ter em um dado problema

# tempo de processamento deste módulo:
t0 = Sys.time()

# modelo:
modelLookup('knn')
##   model parameter      label forReg forClass probModel
## 1   knn         k #Neighbors   TRUE     TRUE      TRUE
# reprodutibilidade
set.seed(400)

# escolha automática do k por validação cruzada
ctrl <- trainControl(method="repeatedcv",repeats = 3)

model_knn <- train(xtrain,ytrain, method = "knn", trControl = ctrl, tuneLength = 20)

#Output of kNN fit
model_knn
## k-Nearest Neighbors 
## 
## 7502 samples
##   12 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 3 times) 
## Summary of sample sizes: 6751, 6752, 6752, 6751, 6753, 6752, ... 
## Resampling results across tuning parameters:
## 
##   k   RMSE        Rsquared   MAE       
##    5  0.05738168  0.6179413  0.03940822
##    7  0.05656514  0.6304058  0.03892341
##    9  0.05614788  0.6372892  0.03874324
##   11  0.05598884  0.6408439  0.03870327
##   13  0.05607594  0.6407082  0.03879974
##   15  0.05597575  0.6432849  0.03875726
##   17  0.05587005  0.6458849  0.03873633
##   19  0.05596962  0.6454175  0.03882718
##   21  0.05611735  0.6439787  0.03896351
##   23  0.05626019  0.6425007  0.03909709
##   25  0.05633778  0.6421681  0.03918249
##   27  0.05645852  0.6411818  0.03928964
##   29  0.05661303  0.6397879  0.03941129
##   31  0.05674481  0.6386296  0.03950853
##   33  0.05684198  0.6377943  0.03961418
##   35  0.05692071  0.6371407  0.03971021
##   37  0.05699808  0.6365043  0.03979948
##   39  0.05708950  0.6357638  0.03988896
##   41  0.05723187  0.6341018  0.04000078
##   43  0.05730753  0.6334136  0.04008218
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was k = 17.
plot(model_knn)

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

# estatística da performance do algoritmo
# sigma equivalente da gaussiana
relz = (ytest - pred_knn)
sig_G = 0.7413*(quantile(relz,0.75,names = FALSE) - quantile(relz,0.25,names = FALSE))
sig_G
## [1] 0.0415019
#visualização do resultado
# vamos juntar os valores preditos e observados do redshift em um fata frame:
my_data = as.data.frame(cbind(observed = ytest,predicted = pred_knn))
head(my_data)                            
##   observed  predicted
## 1    0.082 0.07158824
## 2    0.162 0.11400000
## 3    0.266 0.12682353
## 4    0.197 0.10376471
## 5    0.110 0.11776471
## 6    0.163 0.10111765
# figura com ggplot
ggplot(my_data,aes(predicted, observed)) +
      geom_point(color = "darkred", alpha = 0.5) + 
      geom_smooth(method=lm)+ 
      ggtitle("k-NN: 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'

# tempo de processamento
Sys.time()-t0
## Time difference of 29.089 secs

4.2 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 (as 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

# tempo de processamento deste módulo
t0 = Sys.time()

# 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
model_rf<-train(xtrain,ytrain,method='rf')

# 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.05263951  0.6804789  0.03659620
##    7    0.05158730  0.6904934  0.03538153
##   12    0.05196951  0.6854179  0.03553383
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 7.
plot(model_rf)

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

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

# estatística da performance do algoritmo
# sigma equivalente da gaussiana
relz = (ytest - pred_rf)
sig_G = 0.7413*(quantile(relz,0.75,names = FALSE) - quantile(relz,0.25,names = FALSE))
sig_G
## [1] 0.03541967
# compare esse valor com o obtido por k-NN

#visualização do resultado
# vamos juntar os valores preditos e observados do redshift em um fata frame:
my_data = as.data.frame(cbind(observed = ytest,predicted = pred_rf))
head(my_data)                            
##    observed  predicted
## 2     0.082 0.07098320
## 5     0.162 0.10169323
## 7     0.266 0.12521337
## 9     0.197 0.08415127
## 14    0.110 0.10408917
## 15    0.163 0.13829053
# 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'

# tempo de processamento deste módulo                      
Sys.time() - t0
## Time difference of 22.29201 mins

exercício: verifique se a adoção de outra forma de normalização dos dados afeta os resultados dessa regressão

que tipos de problemas voces conhecem onde regressão pode ser útil?

5. classificação

Vamos agora considerar problemas de classificação, considerando a classificação de um conjunto de dados puramente fotométricos em estrelas ou galáxias.

nesse caso os valores de y são discretos: classes

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

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

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

# alguns detalhes dos dados
dim(tabela)
## [1] 10000    13
head(tabela)
##    u_petro J0378_petro J0395_petro J0410_petro J0430_petro  g_petro J0515_petro
## 1 18.57333    17.99002    18.02439    18.09809    17.92609 17.28233    16.99983
## 2 17.14259    16.78725    16.52434    16.22606    15.83635 15.33351    15.02085
## 3 18.60449    18.16562    17.87111    17.58500    17.17834 16.74320    16.53280
## 4 16.74094    16.38911    16.28852    15.96590    15.66112 15.34018    15.18092
## 5 18.92666    18.72918    18.54957    17.53203    17.32281 16.68435    16.57145
## 6 17.54627    17.25873    17.20477    16.65028    16.57487 16.22518    16.05884
##    r_petro J0660_petro  i_petro J0861_petro  z_petro classe
## 1 16.53974    16.47142 16.21496    16.00299 15.91366      1
## 2 14.59325    14.48909 14.15851    13.96557 13.90085      1
## 3 16.07838    15.98695 15.69374    15.49386 15.44834      1
## 4 14.70637    14.63379 14.39819    14.20232 14.10782      1
## 5 15.87334    15.76929 15.59096    15.50189 15.46906      0
## 6 15.74276    15.70845 15.61446    15.59427 15.56725      0
# número de objetos em cada classe:
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
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])

5.2. classificação com árvores de decisão**

Vamos usar o algoritmo ‘tree’ do pacote rpart

# tempo de processamento:
t0 = Sys.time()

modelLookup('rpart')
##   model parameter                label forReg forClass probModel
## 1 rpart        cp Complexity Parameter   TRUE     TRUE      TRUE
model_tree <- train(xtrain,ytrain,method='rpart',
                      trControl = trainControl(method = "cv"))
# 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.8093333  0.6186667
##   0.05546667  0.7574667  0.5149333
##   0.49813333  0.6968000  0.3936000
## 
## 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               
##

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.1833711 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               
##

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

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

#AUC
AUC=0
for(i in 2:length(a[,1])){
AUC=AUC+(a[i-1,1]+a[i,1])*(a[i,2]-a[i-1,2])/2
}
round(AUC,3)
## [1] 0.831
# tempo de processamento:
Sys.time()-t0
## Time difference of 0.6958725 secs

que problemas de classificação você conhece?

tempo de processamento do script:

stopCluster(cl)

Sys.time() - t00
## Time difference of 24.75258 mins
# comentários sobre o tempo de processamento variando nsample (sec. 4.a) e incluindo ou não paralelização (sec. 1)
# nsample=10000 com paralelização: 30 min
# nsample=10000 sem paralelização: 1.5h
# nsample=1000  com paralelização: 1.6 min
# nsample=1000  sem paralelização: 4.5 min

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

Referências

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

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

Alguns pacotes do R úteis para ML:

Classification And Regression Training: um dos pacotes mais populares. 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.

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

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

classificação e regressão com randomForest

implementa vários algoritmos, como support vector machines

deep learning com TensorFlow

  • e muitos outros!

tempo total de processamento deste script:

Sys.time()-t00
## Time difference of 24.75267 mins