8 Módulo 6 | Aprendizagem automática com R

Introdução ao R

Autor

Manuel Ribeiro

Data de Publicação

14 de outubro de 2024

8.1 INTRODUÇÃO

9 Programa da aula

O que é, para que serve e como funciona a:

Aprendizagem automática
Aprendizagem supervisionada
- Problemas de regressão
- Problemas de classificação
Aprendizagem não-supervisionada

Vamos usar livrarias open-source. No fim da aula, teremos escrito um programa em R para aprendizagem automática.

Há vários bons recursos sobre os tópicos da aula. Sugiro este:

Gareth James et al (2023) “An Introduction to Statistical Learning with Applications in R”. Springer. Disponivel online: https://www.statlearning.com/ . Bem explicado, muitos exemplos R. Algum (pouco) formalismo estatistico.

10 Aprendizagem

A evolução da indústria da computação e da tecnologia criaram condições impares para acesso a grande volumes de dados na generalidade das áreas de ciências e indústria. Como resultado surgiu a necessidade de desenvolver ferramentas capazes de nos ajudar a processar esses dados para os compreender, detectar padrões, fazer previsões.

Aprendizagem automática

É uma área da computação que permite ao computador aprender sem ser explicitamente programado.

Estas ferramentas são:

aprendizagem supervisionada, que aprende com os dados para fazer previsões
aprendizagem não-supervisionada, que encontra relações e estruturas nos dados
aprendizagem de reforço, que aprende por tentativa e erro para prever melhor

Nesta aula vamos abordar aprendizagem supervisionada e não-supervisionada.

10.1 Supervisionada

Colheita dados

Extração de um dataset:
- dados de input (\(x\)),
- dados de output (\(y\)).

Treinar o algoritmo

Dividir o dataset:
- dados de teste para serem usados na fase de teste,
- dados de treino para serem usados para treino e validação.
Escolhe-se uma família de funções (e.g. árvores de decisão),
Optimizam-se os parâmetros das funções para errarem o menos possível nos dados de treino.
Avalia-se o desempenho com métricas (i.e. funções de custo) nos dados de validação).

Fazer previsões: O erro de teste esperado (na classificação ou previsão de novas observações) é avaliado nos dados de teste com a função aprendida (a que errou menos). Se o resultado é satisfatório, a função é usada para classificar ou fazer previsões.

Pipeline clássico em aprendizagem automática supervisionada

O processo de aprendizagem supervisionada evolui com o treino. Quanto mais treinam, melhor desempenham a tarefa de fazer previsões. O treino faz-se com técnicas de validação cruzada.

10.1.1 Regressão

Quando \(y\) é quantitativa estamos perante um problema de regressão.

Suponha que queremos usar um algoritmo de aprendizagem automática para prever a pressão arterial sistólica esperada de pacientes, a partir da sua idade. Para um algoritmo ser capaz de aprender e prever com exatidão, precisa aprender com exemplos.

A base de dados do hospital fornece um conjunto de exemplos de doentes para os quais se registaram os dados:

\(y\), pressão arterial sistólica.
\(x_1\), idade do paciente.
\(x_2\), peso do paciente.

Os dados destes doentes servirão para aprender uma função para prever a pressão em novos doentes.

id	pressao	idade	peso
1	132	52	78
2	143	59	83
3	153	67	88
4	162	73	96
5	154	64	89
6	168	74	100

Neste dataset, a idade e o peso são preditores \(\mathbf{x}\) (ou features, input features) e a pressão arterial é a resposta \(y\) (ou outcome, output). Poderiamos experimentar treinar uma regressão linear (\(f\)) para mapear a relação entre os preditores e a resposta.

10.1.2 Classificação

Quando \(y\) é qualitativa, categórica (ordinal ou nominal) estamos perante um problema de classificação.

Imagine que o registo oncológico forneceu-nos o dataset seguinte para treinarmos uma função capaz de prever (com muita precisão) se um tumor (\(y\)) é maligno ou benigno a partir do seu tamanho e da idade do paciente (\(\mathbf{x}\)):

id	tamanho	idade	tumor
1	1.5	60.6	Benigno
2	2.2	52.9	Benigno
3	2.8	63.4	Benigno
4	4.9	54.6	Maligno
5	3.5	56.4	Maligno
6	6.1	55.5	Maligno

A resposta (tipo de tumor) é uma variável do tipo categórica (neste caso com duas categorias, “Benigno” ou “Maligno”), estamos perante um problema de classificação.

10.1.3 Bias-variance trade-off

Entre as funções e algoritmos usados em aprendizagem, nenhum é sempre melhor que os outros todos. Num dataset particular, um algoritmo pode funcionar melhor, mas outro algoritmo pode ser melhor noutro dataset semelhante.

Bias-Variance trade-off

É um dos desafios mais importantes em aprendizagem automática. O sucesso de qualquer método de aprendizagem depende da escolha correta de um algoritmo com flexibilidade adequada. Neste contexto é essencial compreender os conceitos de viés e de variância.

Para o processo de aprendizagem especificamos uma função que nos pareça apropriada e treinamos a função com os dados. O treino é acompanhado por uma medida de erro, que reflete o desempenho da função. Vamos falar de duas frequentemente usadas em problemas de regressão e de classificação (respectivamente):

MSE (Mean Squared Error, Erro Quadrático Médio)
ER (Error Rate, Taxa de Erro)

Estas medidas de desempenho são usadas para escolher o algoritmo que erra menos. A escolha deve recair num que capte a informação relevante dos dados e não capte o ruído associado. Encontrar a região ou o ponto certo (sweet spot) passa por encontrar o ponto onde o viés e a variância da função minimizam o valor de MSE|ER - bias-variance tradeoff.

Ajustamento

Quando não capta a informação relevante (variância) dos dados diz-se que a função está sub-ajustada (underfitted). Quando a função é demasiado especifica (ou complexa) e por isso não conseguirá prever com precisão a resposta em dados de teste, diz-se que a função está sobre-ajustada (overfitted).

Nas sub-secções seguintes vamos ver estes conceitos com um exemplo de regressão e outro de classificação.

Regressão

Suponha que queremos prever a intensidade da dor (medida numa escala quantitativa) causada por uma doença especifica, a partir da idades dos doentes. Temos um dataset com 40 doentes de diferentes idades (\(x\)), onde se registou a intensidade da dor (\(y\)) e vamos treinar um função que mapeie a relação entre a idade (preditor) e a intensidade (resposta). Os dados estão projectados no gráfico de dispersão seguinte:

Vamos fingir por momentos que somos Deus, e sabemos que a verdadeira relação entre a idade e a intensidade da dor tem esta forma: \[ f(x)=x^4+ \epsilon \] O gráfico da função é assim:

Vemos que os pontos não se sobrepõem exatamente à linha (função). Isto acontece porque há outros factores que não controlamos e que interferem na relação (ruído).

Vamos agora “voltar à Terra”, imaginar que não sabemos a formula que exprime a relação entre a idade e a intensidade da dor.

Vamos usar os dados para aprender uma função capaz de prever a forma da relação. A primeira coisa que vamos fazer é dividir (aleatoriamente) o dataset em 2 partições, treino e teste (ou validação):

E vamos treinar 3 polinómios com os dados de treino que depois usaremos os dados de teste para medir o erro ou precisão:

A: Polinómio de grau 1
B: Polinomio de grau 3
C: Polinómio de grau 12

Estes polinómios distinguem-se pelo número de parametros a estimar. O que pretendemos é encontrar a função que minimiza o MSE. Vejamos os resultados obtidos:

A primeira função (A) ajusta uma recta aos dados de treino. As segunda (B) e terceira (C) ajustam-se melhor. Indo de A para C observamos que as funções se tornam cada vez mais flexiveis, i.e., “captam” cada vez melhor a variabilidade, aproximando cada vez mais a linha aos dados de treino.

Quiz

Olhando para os gráficos (A,B,C), qual das funções parece melhor para prever?

O erro de treino, \(\varepsilon\), é dado pela distância vertical entre o ponto e a linha. \[ \varepsilon = \text{valor observado}-\text{valor previsto} \] A média do somatório do quadrado dos erros de treino dá-nos o MSE. \[ \text{MSE} = \frac{1}{n}\sum\varepsilon^2 \] Aqui, o MSE é de treino pois é quantificado nos dados de treino. No gráfico seguinte vemos o MSE de treino estimado para as três funções.

O MSE baixa à medida que a função é mais complexa (tem mais parâmetros). O polinómio A nunca conseguirá captar a forma curva da relação entre idade e intensidade da dor (underfitted). Tem um viés elevado. O polinómio C tem o viés menor, pois a aproximação aos pontos é a elevada (entre as 3 funções).

Viés

Viés refere-se ao erro introduzido por aproximarmos um problema complexo com uma função demasiado simples.

Será razoável perguntar nesta altura, porque não escolhemos a função mais flexivel se esta capta melhor a variabilidade? A resposta é simples: não estamos interessados em prever a intensidade da dor nos pacientes para os quais já conhecemos a intensidade da dor. Estamos sim interessados em prever a intensidade da dor em pacientes futuros.

Então como se comportarão então as funções ajustadas com novos dados? Vamos obter uma estimativa com os dados de teste:

O MSE de teste devolve o seguinte resultado:

Com dados novos o polinómio grau três (B) é o que apresenta menor MSE. Consegue captar a informação mais importante, a forma curva da relação e representa um bom compromisso entre tem um viés e a variáncia (bias-variance trade-off). O polinómio de grau um (A) apresenta um MSE superior porque terá um viés elevado. Finalmente, o polinómio de grau 12 (C) apresenta o MSE mais elevado porque terá maior variância.

Variância

A variância refere-se ao quanto teria de mudar a função, se a ajustássemos a um novo dataset de treino.

Regra geral, à medida que a função se torna mais flexivel, o viés diminui e a variância aumenta. Até certo ponto, o viés diminui e o MSE de teste também. Mas a partir de certo ponto, aumentar a flexibilidade já não reduz o viés de forma relevante, mas tende a aumentar muito a variância, e leva a que o MSE de teste comece a subir. O algoritmo entra em sobre-ajustamento. Nos dados de treino o MSE é cada vez mais baixo

Fonte: Hastie, Trevor et al. *The Elements of Statistical Learning*. New York :Springer, 2008

Classificação

No caso de problemas de classificação, lidamos com respostas (\(y\)) qualitativas. Uma medida de erro usual para os classificadores é a taxa de erro: \[ \begin{split} \mathrm{ER} &=\frac{1}{n}\cdot\sum I, \quad I=\begin{cases} 0, & y_i = \hat{f}(x_i) \\ 1, & y_i\ne \hat{f}(x_i) \end{cases} \\ \end{split} \]

Para prever se um tumor é maligno ou benigno (\(y\)) a partir do seu tamanho e da idade do paciente (\(\mathbf{x}\)) vamos usar um dataset (\(n= 120\)) e treinar um algoritmo de classificação simples (e surpreendentemente eficaz), o \(k\)-NN (k-vizinhos mais próximos).

O dataset tem uma distribuição ‘equilibrada’ entre tumores malignos (\(n=74\), 62%) e benignos (\(n=46\), 38%).

A estratégia do algoritmo passa por classificar um novo tumor na classe mais comum entre os \(k\) vizinhos mais próximos. Como em todos os algoritmos de aprendizagem, vamos começar por dividir aleatoriamente o dataset em dados de treino (\(n_{tr}= 60\)) e dados de teste (\(n_{ts}=60\)):

Vamos treinar o algoritmo, com diferentes valores para \(k\) e depois estimamos a exatidão (accuracy) ou o erro (1-accuracy) nos dados de teste. Nos gráficos seguintes, ilustram-se os resultados para \(k=1\), \(k=11\) e \(k=41\).

Hiperparâmetros

Alguns modelos ou parâmetros não podem ser estimados directamente dos dados. Designam-se por hiperparâmetros e são “afinados” (tuned) por métodos de validação cruzada (sub-secção seguinte). O número de vizinhos num modelo \(k\)-NN é um desses parâmetros.

Em cada gráfico, a região azul é a zona onde novas observações (dados de teste) serão classificadas como tumores benignos e a região vermelha a zona onde serão classificadas como tumores malignos. Designam-se por regiões de decisão e são delimitadas pela fronteira de decisão.

O classificador \(k\)-NN com \(k=1\) devolveu o resultado seguinte.

A fronteira de decisão apresenta-se muito recortada, com forma não-linear.

Quiz

Qual a taxa de erro nos dados de treino com classificador \(1\)-nn?

Será que este classificador acertará muitas vezes quando classificar novas observações? Vamos projectar os dados de teste nas regiões definidas e ver quantas vezes os dados são mal classificados:

O classificador cometeu erros de classificação em 6 observações de teste (ER = 0.1). Vamos afinar o classificador com \(k=11\) e \(k=41\) e avaliar os erros com dados de teste.

O gráfico seguinte ilustra o resultado da projecção dos dados de teste (não usados para afinar o valor de \(k\)) nas regiões e fronteiras de decisão obtidas com \(k=11\) e \(k=41\).

Em relação a \(1\)-nn as regiões de decisão com \(11\)-nn e \(41\)-nn tendem a apresentar uma fronteira de decisão menos recortada, porque usam mais observações para classificar cada ponto da região. Nos dados de teste projectados nestas regiões vêem-se algumas observações mal classificadas. A estimativa da taxa de erro de teste está ilustrada no gráfico seguinte. Os classificadores mais e menos flexiveis (\(k= 1\) e \(k=41\), respectivamente) apresentam uma taxa mais elevada em relação ao classificador com \(k\) “intermédio” (\(k=11\)). (ilustrado no gráfico seguinte).

Este último responde melhor ao problema do bias-variance trade-off (compromisso entre o viés e a variância).

Erro de generalização

O valor esperado do erro de teste pode ser decomposto em 3 quantidades fundamentais: a variância e o viés (quadrado) da função, e a variância do erro aleatório (irredutível). A forma em U dos erros de teste que observámos nos dois exemplos são resultado de duas quantidades que competem entre si: o viés e a variância.

A variância da função, Var(\(f\)) refere-se à variação da função quando é ajustada a outro dataset. Quanto mais flexivel for a função, maior a variância esperada de \(f\). O viés da função, Viés(\(f\)), refere-se ao erro introduzido por aproximarmos um problema complexo com um modelo simples. A variância do erro aleatório, Var(\(\epsilon\)), é algo que não controlamos, depende dos dados e não dos algoritmos de aprendizagem que escolhemos.

10.1.4 Validação

Que procedimentos temos para encontrar o ponto de equilíbrio?

O ideal será termos dados teste para estimar o erro de teste.

A partir de um dataset, separamos em 3 partes:

Dados treino,
Dados validação,
Dados teste

Em datasets de menor dimensão, separamos o dataset em 2 partes:

Dados de treino e de validação,
Ajustamos os modelos nos dados de treino,
Medimos o erro nos dados de validação.

O erro nos dados de validação é uma estimativa do erro de teste. Como vimos atrás, escolhe-se uma medida de erro que pode ser uma das seguintes (ER para problemas de classificação e MSE para problemas de regressão):

\[ \begin{split} \mathrm{ER} &=\frac{1}{n}\cdot\sum_{i=1}^n I(y_i\ne \hat{y}_i), \quad I(y_i\ne\hat{y}_i)=\begin{cases} 0, & y_i = \hat{y_i} \\ 1, & y_i\ne \hat{y_i} \end{cases} \\ \mathrm{MSE} &= \frac{1}{n}\cdot\sum_{i=1}^n (y_i- \hat{y}_i)^2 \end{split} \]

E como podemos treinar e avaliar o desempenho de diferentes modelos? Os procedimentos habitu88ais para estimar o erro designam-se por validação simples (que já usámos nos exemplos anteriores) ou validação cruzada:

Conjunto de validação
Validação cruzada (Leave-one-out)
Validação cruzada (k-fold)

A validação simples (ou conjunto de validação) é fácil de implementar. Usámo-la nos toy example.

Esquema ilustrativo da técnica conjunto de validação¹

Vantagens:

é simples e fácil de implementar

Desvantagens:

Estimativa do erro de teste muito variável
Com poucas observações tende a sobre-estimar o erro de teste
‘Desperdício’ de observaçoes (% dados validação)

No procedimento por validação cruzada k-fold dividimos o dataset de treino com \(n\) observações em \(k\) partições de dimensão igual (geralmente com \(k \ll n\)).

Vantagens (em relação ao conjunto de validação):

Estimativa do erro tem menor risco de sobre-estimação

Quando o \(n\) é muito pequeno (e.g. \(n=20\)) podemos usar a validação cruzada com \(k=n\), também designada por Leave-one-out.

Vantagens (em relação à validação cruzada):

O resultado é sempre igual

Desvantagens:

Tempo de processamento com \(n\) grande (por isso só se usa quando \(n\) pequeno).

10.2 Não-supervisionada

Esquema com a visão geral da aprendizagem não-supervisionada⁴

Na aprendizagem supervisionada temos preditores (e.g. x1 e x2) uma resposta (e.g. classe 0 ou classe 1). O objectivo é aprender uma função capaz de prever a resposta com novas observações de x1 e x2.

Na aprendizagem não-supervisionada temos apenas variáveis de input (e.g x1 e x2), não há uma resposta output. O objectivo é encontrar relações ou estruturas nos dados.

Quiz

Olhando para o gráfico seguinte, identifica alguma estrutura nos dados?

Vejamos o seguinte exemplo: imaginemos que temos um dataset como o da figura seguinte, temos dados mas não nos dizem o que fazer com eles. Apesar disso, podemos analisar os dados, procurando relações ou estruturas subjacentes nos dados.

Se usarmos um algoritmo para aprendizagem não-supervisionada conhecido genericamente de clustering, o algoritmo irá identificar dois sub-grupos de dados (ou clusters), e dividir os dados em dois grupos:

Nesta aula vamos ver dois métodos de aprendizagem não-supervisionada:

análise em componentes principais (ACP), usada para visualização de dados ou pre-processamento em aprendizagem supervisionada,
clustering, é uma família de métodos usada para encontrar sub-grupos (clusters) nos dados.

Se quisermos analisar uma tabela de dados com centenas de pacientes e dezenas ou centenas de variáveis ser-nos-á útil usar uma ferramenta capaz de reduzir a quantidade de informação (nr variáveis) sem perder informação relevante. Uma das ferramentas mais usadas para desempenhar essa tarefa é a ACP.

Se quisermos identificar e classificar/agrupar as observações num número pequeno de classes (como no exemplo acima) podemos usar um algoritmo de clustering.

10.2.1 Exemplo ACP

Vejamos um exemplo de aplicação com dados da COVID-19: para cada um dos 278 municipios de Portugal Continental, temos um dataset que contém as taxas de incidência cumulativa a 14 dias por 100 000 habitantes registadas durante 570 dias consecutivos entre Agosto 2020 e Março 2022.

Na figura seguinte ilustra-se com as curvas das séries temporais (cada cor um municipio) a evolução das taxas de incidência ao longo do tempo (a curva preta representa a evolução média nacional).

Séries temporais dos 278 municipios. Adaptado de ‘Understanding spatiotemporal patterns of COVID-19 incidence in Portugal: A functional data analysis from August 2020 to March 2022’, M Ribeiro et al, 2024, PLoS ONE, 19(2): e0297772

A ACP foi realizada com as séries temporais para detetar padrões de evolução ao longo do tempo. Extraímos da ACP três curvas que contêm 72% da variação observada nas 278 curvas(!) correspondentes a três efeitos principais captados no periodo analisado:

efeito global, comum a todos os municipios (PC1)
efeito de prevenção, comum a municipios de incidencia baixa (PC2)
efeito dos surtos, comum aos municipios de menor dimensão onde se registaram surtos (PC3)

Primeiras três componentes principais. Adaptado de ‘Understanding spatiotemporal patterns of COVID-19 incidence in Portugal: A functional data analysis from August 2020 to March 2022’, M Ribeiro et al, 2024, PLoS ONE, 19(2): e0297772

A ACP permitiu caracterizar a evolução temporal das 278 curvas iniciais apenas com 3 curvas que captaram grande parte da variabilidade total (72%) observada nos dados originais.

10.2.2 Exemplo clustering

Nesta mesma análise recorreu-se a um método de clustering (método de classificação hierarquica) para classificar os municipios de acordo com as semelhanças espaço-temporais das curvas.

M Ribeiro et al (2024) ‘Understanding spatiotemporal patterns of COVID-19 incidence in Portugal: A functional data analysis from August 2020 to March 2022’. PLoS ONE, 19(2): e0297772

Os resultados da análise com ACP e clustering permitiram revelar diferenças na dinâmica espaço-temporal entre:

os municipios das regiões norte e litoral vs sul e interior
efeitos ocorridos no outono-inverno 20-21 vs 21-22

Embora possa parecer que estas análises complexas necessitem de muitas linhas de código para extrair as curvas principais ou para classificar todas as curvas, é possível escrever o algoritmo apenas com duas linhas de código. Veremos isso.

11 Aplicações

Vamos recorrer ao dataset HeartDisease (disponível na biblioteca MLDataR) para prever a existência de uma doença cardíaca a partir de dados sobre pressão sanguínea, frequência cardíaca máxima, história de angina de peito e outros preditores. Os exemplos são dados de pacientes, anonimizados, disponibilizados pela British Heart Foundation.

# Ler as bibliotecas de funções para a sessão
library(tidymodels) 
library(MLDataR)
library(rpart.plot) # para plot arvore de decisao
library(kknn)

# Ler o dataset HearDisease
heartd <- MLDataR::heartdisease

O dataset contém 9 preditores para além da variável resposta (HeartDisease) e 918 observações.

dim(heartd)

[1] 918  10

head(heartd)

# A tibble: 6 × 10
    Age Sex   RestingBP Cholesterol FastingBS RestingECG MaxHR Angina
  <dbl> <chr>     <dbl>       <dbl>     <dbl> <chr>      <dbl> <chr> 
1    40 M           140         289         0 Normal       172 N     
2    49 F           160         180         0 Normal       156 N     
3    37 M           130         283         0 ST            98 N     
4    48 F           138         214         0 Normal       108 Y     
5    54 M           150         195         0 Normal       122 N     
6    39 M           120         339         0 Normal       170 N     
# ℹ 2 more variables: HeartPeakReading <dbl>, HeartDisease <dbl>

Cada linha representa um paciente que tem ou não tem doença cardiaca. O dataset tem informação sobre:

Age, idade do paciente em anos;
Sex, sexo do paciente (“F”, “M”);
RestingBP, pressão sanguínea para ritmo cardíaco em repouso;
Cholesterol, nível de colestrol no sangue;
FastingBS, amostra de sangue de glicose após jejum (1=diabetes, 0=normal);
RestingECG, indicador de infartos do miocárdio prévios (“Normal”,“ST”,“LVH”);
MaxHr, batimento cardíaco máximo;
Angina, dor de peito causada por redução do fluxo sanguíneo (“N”,“Y”);
HeartPeakReading, leitura no pico da frequência cardiaca;
HeartDisease, paciente com doença cardíaca (1= sim, 0= não);

Vamos verificar se os tipos de dados das variáveis são adequados para representar os valores.

str(heartd)

tibble [918 × 10] (S3: tbl_df/tbl/data.frame)
 $ Age             : num [1:918] 40 49 37 48 54 39 45 54 37 48 ...
 $ Sex             : chr [1:918] "M" "F" "M" "F" ...
 $ RestingBP       : num [1:918] 140 160 130 138 150 120 130 110 140 120 ...
 $ Cholesterol     : num [1:918] 289 180 283 214 195 339 237 208 207 284 ...
 $ FastingBS       : num [1:918] 0 0 0 0 0 0 0 0 0 0 ...
 $ RestingECG      : chr [1:918] "Normal" "Normal" "ST" "Normal" ...
 $ MaxHR           : num [1:918] 172 156 98 108 122 170 170 142 130 120 ...
 $ Angina          : chr [1:918] "N" "N" "N" "Y" ...
 $ HeartPeakReading: num [1:918] 0 1 0 1.5 0 0 0 0 1.5 0 ...
 $ HeartDisease    : num [1:918] 0 1 0 1 0 0 0 0 1 0 ...

Vamos codificar as variáveis (preditores e resposta) com tipos de dados adequados. Por ex. HearDisease é uma variável do tipo numérico, vamos convertê-la numa variável categorica (no R designa-se por factor). Vamos aplicar o mesmo racional a outros casos semelhantes. Em algumas variáveis adicionamos um label para facilitar a leitura das classes.

heartd$Sex <- factor(heartd$Sex)
heartd$FastingBS  <- factor(heartd$FastingBS, labels = c("Normal","Diabetes"))
heartd$RestingECG  <- factor(heartd$RestingECG)
heartd$Angina <- factor(heartd$Angina)
heartd$HeartDisease <- factor(heartd$HeartDisease, labels = c("N", "Y"))

Nas próximas sub-secções vamos preparar os dados para o processo de aprendizagem, que passará por:

Fase de pre-processamento (recipes)
Especificação dos modelos,
Especificar afinação (tune) e métodos de reamostragem (resamples)

Depois, nas secções seguintes implementaremos o processo de aprendizagem com três familias de classificadores:

Regressão logística,
\(k\)-nearest neighbours (\(k\)-NN),
Random forests

A regressão logística é da familia dos classificadores linear (a fronteira de decisão é linear, isto é, uma linha recta ou, no caso multivariado, um hiperplano), enquanto os restantes são da família dos classificadores não-lineares (nestes, a fronteira de decisão linear não separa as classes, é necessário um modelo não linear para as separar).

Por fim usaremos as funções aprendidas para prever e obter a estimativa do erro de teste.

11.1 Análise exploratória

A análise exploratória é essencial ao processo de modelação estatística. Vamos ver exemplos de como esta etapa nos fornece informação relevante para termos depois melhor desempenho no processo de modelação.

Vamos começar por focar-nos na variável resposta. Podemos criar um gráfico de barras para ilustrar a sua distribuição.

heartd %>%
  count(HeartDisease) %>%
  mutate(freq = round(n / sum(n),2)) %>%
  mutate(cum = cumsum(freq)) %>%
  ggplot(aes(x = "", y = freq, fill = HeartDisease, label = freq)) +
  geom_bar(position = "stack", stat = "identity") +
  scale_fill_viridis_d(end = .5) +
  geom_text(aes(label = freq), position = position_stack(vjust = .95), colour = "white") +
  labs(fill = "Heart disease")

A distribuição de Ye Né ‘equilibrada’. Para modelos lineares como por ex. o modelo de regressão logística, este aspeto é muito relevante, porque se uma das classes fosse muito mais frequente (usual em casos de doenças raras) conduziria a estimativas enviesadas e a um elevado erro de generalização. Neste dataset, no entanto, esse problema não se põe.

Vejamos a distribuição da resposta por sexo:

heartd %>% 
  ggplot(aes(x = HeartDisease, fill = Sex)) +
  geom_bar() +
  scale_fill_viridis_d(end = .5)

Aqui observamos que o peso do sexo feminino F entre pacientes com doença cardíaca é menor que o peso do sexo feminino entre pacientes sem doença cardíaca. Este gráfico sugere que há diferenças na ocorrência da doença cardiaca em função do sexo.

No caso de preditores contínuos podemos criar um histograma. No gráfico seguinte vemos a distribuição do batimento cardíaco máximo. O gráfico ilustra a proporção de pacientes com e sem doença cardiaca ao longo da distribuição. Vemos que até ~ 150 batimentos, a proporção de doentes cardiacos tende a ser superior à dos não cardíacos, e depois disso a tendencia inverte-se.

heartd %>%
  ggplot(aes(x = MaxHR, fill = HeartDisease)) +
  geom_histogram(binwidth = 2) + 
  scale_fill_viridis_d(end = .5)

Verificamos que os preditores analisados parecem estar bem correlacionados com a resposta, dando boas indicações quanto à sua capacidade preditiva.

Uma análise exploratória é essencial antes de qualquer exercicio de modelação, para detetar padrões que aparentam não ser aleatórios, relações entre a resposta e os preditores considerados.

11.2 Pre-processamento

Antes do treino, os dados passam por um pré-processamento para ‘limpar’ e transformar os dados. Estas operações são geralmente muito úteis para preparar os dados e fornecer um dataset mais adequado ao processo de aprendizagem⁵. As necessidades de pré-processamento variam com os modelos usados.

Neste exemplo vamos proceder à implementação de operações frequentemente necessárias e que serão implementadadas em todos os classificadores que vamos usar:

a partição dos dados
transformação de preditores qualitativos em variáveis dummy
normalização de preditores quantitativos

Vamos proceder à partição em dados de treino e dados de teste. Os dados de treino serão usados para estimar parâmetros, afinar hiperparâmetros, comparar modelos. Para assegurar que a partição honra as proporções de pacientes com e sem doença cardíaca, a partição por amostragem aleatória é estratificada. Os dados de teste ficarão “guardados numa gaveta” até ao fim do processo de treino e são usados no final para estimar o erro de generalização.

# sempre que usamos uma função geradora de um output aleatorio, 
# podemos usar o seeding para assegurar sempre o mesmo output (aleatorio) 
set.seed(12399300)

# initial_split() cria as partições e
# arg strata estratifica a amostragem aleatória
# em cada partição
particao <- initial_split(heartd, strata = HeartDisease)

# atribuir nome/objeto aos
# datasets de treino e de teste
heartdtr <- training(particao)
heartdtt  <- testing(particao)

Conforme vimos em Seção 10.1.4 podemos criarmos partições do dataset de treino para estimar o erro de validação. Vamos usar a função vfold_cv() para criar as ditas partições e guardar o resultado num objeto para uso nas fases de treino.

set.seed(4800000)
k_folds <- vfold_cv(heartdtr, v = 10, strata = HeartDisease)

Vamos também normalizar preditores numéricos e criar variáveis dummy para preditores categoricos. Estas operações são feitas com a função recipe() e as funções step_normalize() e step_dummy(). O argumento da função inclui a formula do modelo. O simbolo ~ especifica que é uma formula, onde à sua esquerda está (o nome da) variável resposta e à direita o dos preditores. O simbolo . indica que pretendemos todos os preditores do dataset.

receita <- 
  recipe(HeartDisease ~ ., data = heartdtr) %>% 
  step_normalize(all_numeric(), -all_outcomes()) %>% 
  step_dummy(all_nominal(), -all_outcomes())

Vejamos os resultados do pré-processamentocom as funções prep() e bake(). Estas são usadas respectivamente para 1) estimar os parâmetros dos dados de treino necessários para serem depois aplicados a outros datasets e 2) executar os cálculos aos novos dados.

receita %>% 
  prep(training = heartdtr) %>% 
  bake(new_data = NULL)

# A tibble: 688 × 11
       Age RestingBP Cholesterol    MaxHR HeartPeakReading HeartDisease Sex_M
     <dbl>     <dbl>       <dbl>    <dbl>            <dbl> <fct>        <dbl>
 1 -1.46       0.429      0.835   1.37              -0.812 N                1
 2 -1.78      -0.107      0.780  -1.51              -0.812 N                1
 3 -1.57      -0.644      1.29    1.29              -0.812 N                1
 4 -0.927     -0.107      0.364   1.29              -0.812 N                0
 5  0.0374    -1.18       0.102   0.200             -0.812 N                1
 6 -1.78      -0.107      0.129   0.200             -0.812 N                0
 7 -1.57      -0.644      0.0656  0.317             -0.812 N                1
 8 -1.25      -0.912      0.129   0.00587           -0.812 N                0
 9  0.0374    -0.644      0.690   0.511              0.571 N                0
10 -1.14      -0.644      0.0385  1.09              -0.812 N                0
# ℹ 678 more rows
# ℹ 4 more variables: FastingBS_Diabetes <dbl>, RestingECG_Normal <dbl>,
#   RestingECG_ST <dbl>, Angina_Y <dbl>

11.3 Modelo logístico

11.3.1 Introdução

Com o modelo de regressão logístico o objectivo é prever o valor de uma variável categórica binária (classe 0 ou 1) como uma função linear de preditores. Neste exemplo, vamos usar a regressão logística para prever se um paciente tem (classe 1) ou não (classe 0) uma doença cardíaca.

A variável resposta será a probabilidade de um paciente ter doença cardíaca.

Regressão logística e o logito

A fórmula para calcular a probabilidade \(\Pr(y=1|x)\) com a regressão logística:

\[ \Pr(y=1|x) =\frac{1}{1+\exp^z}=\frac{\exp^z}{1+\exp^z}, \quad z=\beta_0 +\sum\mathrm{\beta} \mathbf{x} \]

onde \(\Pr(y=1|x)\) é a probabilidade de ocorrer a classe de interesse (no nosso exemplo, ter doença cardíaca = 1, não ter doença cardíaca = 0), condicionada pelos preditores \(x\). \(\beta\) são coeficientes associados a \(x\) e \(\beta_0\) a ordenada na origem .

O logito é uma função ligação que transforma a regressão logística num modelo linear: \[ \log\bigg[\frac{\Pr(y=1|x)}{1-\Pr(y=1|x)}\bigg] = z \] O termo \(\frac{\Pr(y=1|x)}{1-\Pr(y=1|x)}\) designa-se usualmente por chances (odds).

Como o modelo devolve a probabilidade de um paciente ter doença cardiaca temos de encontrar um critério (regra de decisão) para convertermos a probabilidade numa das classes.

Regressão logística é um classificador

Embora inclua o termo regressão no nome, a regressão logística é aplicada em problemas de classificação. Para classificar uma observação numa das classes, combinamos o modelo de regressão logistico com um critério ou regra de decisão.

Uma regra de decisão frequentemente aplicada passa por usar o classificador de Bayes (que escolhe a classe que maximiza a probabilidade). Num problema de classificação onde há apenas 2 classes, é usual prever classe 1 se \(\Pr(y=1|x)\ge 0.5\) e na classe 0 em caso contrário. Mas podem escolher-se outras regras de decisão.

11.3.2 Treino

Agora que temos o pré-processamento terminado vamos preparar a aprendizagem. Vamos usar a estrutura tidy para especificar o modelo (model) logístico de regressão, com o package glm (engine) para fazer classificação (mode):

log_spec <-
  logistic_reg() %>% 
  set_engine("glm") %>%
  set_mode("classification")

Vamos usar a função workflow()para integrar de forma encadeada a receita definida na Seção 11.2 e as especificações do modelo.

workflow

Em Python e outras linguagens designa-se uma série de operações encadeadas por pipeline. No tidymodels o termo pipe refere-se ao operador %>%, por isso, para evitar ambiguidades, no tidymodels essa série designa-se por workflow.

A sintaxe seguinte apresenta uma forma de especificar o workflow:

workflow() %>%
  add_recipe(receita) %>%
  add_model(log_spec)

Vamos guardar o workflow num objeto:

log_wflow <- 
  workflow() %>%
  add_recipe(receita) %>%
  add_model(log_spec)

Podemos fazer o ajustamento (fit) aos dados de treino:

log_fit <- 
  log_wflow %>%
  fit(data = heartdtr)

Qual será a classe prevista nos dados de treino? Podemos prever com a função predict(), que devolve uma coluna com a classe prevista.

# predict()
predict(log_fit, new_data = heartdtr)

Usando a função augment() (da biblioteca broom), podemos obter informação adicional sobre os exemplos (e.g. valores previstos)⁶.

# augment()
broom::augment(log_fit, new_data = heartdtr)

# A tibble: 688 × 13
     Age Sex   RestingBP Cholesterol FastingBS RestingECG MaxHR Angina
   <dbl> <fct>     <dbl>       <dbl> <fct>     <fct>      <dbl> <fct> 
 1    40 M           140         289 Normal    Normal       172 N     
 2    37 M           130         283 Normal    ST            98 N     
 3    39 M           120         339 Normal    Normal       170 N     
 4    45 F           130         237 Normal    Normal       170 N     
 5    54 M           110         208 Normal    Normal       142 N     
 6    37 F           130         211 Normal    Normal       142 N     
 7    39 M           120         204 Normal    Normal       145 N     
 8    42 F           115         211 Normal    ST           137 N     
 9    54 F           120         273 Normal    Normal       150 N     
10    43 F           120         201 Normal    Normal       165 N     
# ℹ 678 more rows
# ℹ 5 more variables: HeartPeakReading <dbl>, HeartDisease <fct>,
#   .pred_class <fct>, .pred_N <dbl>, .pred_Y <dbl>

Vamos agora ver como medir o desempenho do classificador e ilustrar resultados:

exatidão (\([\text{TP+TN]}/n\)): proporção de exemplos classificados corretamente;
precisão (\(\text{TP}/(TP+FP)\)): proporção de exemplos previstos na classe 1, que eram da classe 1;
sensibilidade (\(\text{TP}/[\text{TP}+\text{FN}]\)): proporção de exemplos da classe 1, previstos na classe 1;
especificidade (\(\text{TN}/[\text{TN}+\text{FP}]\)): proporção de exemplos da classe 0, previstos na classe 0;
curva ROC (Receiver Operating Characteristic).

Com a tabela expandida podemos ter informação adicional sobre o processo de classificação. Vemos que o classificador usa como regra de decisão atribuir a classe mais provável.

Com a função relocate() podemos mudar a ordem das colunas. Vamos usar a função para ver lado-a-lado as classes observadas e previstas para os dados de treino:

augment(log_fit, new_data = heartdtr) %>%
  relocate(HeartDisease, .pred_class, .pred_Y, .pred_N)

# A tibble: 688 × 13
   HeartDisease .pred_class .pred_Y .pred_N   Age Sex   RestingBP Cholesterol
   <fct>        <fct>         <dbl>   <dbl> <dbl> <fct>     <dbl>       <dbl>
 1 N            N            0.106    0.894    40 M           140         289
 2 N            N            0.296    0.704    37 M           130         283
 3 N            N            0.0941   0.906    39 M           120         339
 4 N            N            0.0541   0.946    45 F           130         237
 5 N            N            0.241    0.759    54 M           110         208
 6 N            N            0.0827   0.917    37 F           130         211
 7 N            N            0.195    0.805    39 M           120         204
 8 N            N            0.0929   0.907    42 F           115         211
 9 N            N            0.199    0.801    54 F           120         273
10 N            N            0.0629   0.937    43 F           120         201
# ℹ 678 more rows
# ℹ 5 more variables: FastingBS <fct>, RestingECG <fct>, MaxHR <dbl>,
#   Angina <fct>, HeartPeakReading <dbl>

Podemos construir a matriz de confusão para calcular diferentes medidas de desempenho (e.g. precisão, exatidão, sensibilidade,..).

augment(log_fit, new_data = heartdtr) %>%
  conf_mat(truth = HeartDisease, estimate = .pred_class)

Para melhorar a visualização da matriz podemos usar a função autoplot()

augment(log_fit, new_data = heartdtr) %>%
  conf_mat(truth = HeartDisease, estimate = .pred_class) %>%
  autoplot(type = "heatmap")

Podemos calcular a exatidão com a função accuracy()

augment(log_fit, new_data = heartdtr) %>%
  accuracy(truth = HeartDisease, estimate = .pred_class)

a sensibilidade com a função sensitivity()

augment(log_fit, new_data = heartdtr) %>%
  sensitivity(truth = HeartDisease, estimate = .pred_class)

a especificidade com a função specificity()

augment(log_fit, new_data = heartdtr) %>%
  specificity(truth = HeartDisease, estimate = .pred_class)

Podemos usar a função yardstick::metric_set()para devolver estas as medidas de uma só vez:

heart_metrics <- metric_set(accuracy, specificity, sensitivity)
augment(log_fit, new_data = heartdtr) %>%
  heart_metrics(truth = HeartDisease, estimate = .pred_class)

# A tibble: 3 × 3
  .metric     .estimator .estimate
  <chr>       <chr>          <dbl>
1 accuracy    binary         0.820
2 specificity binary         0.840
3 sensitivity binary         0.795

Podemos visualizar a curva ROC (Receiver Operating Characteristic):

No argumento da função fit_samples() vamos usar o objecto que criámos na sub-secção Seção 11.2 para criar treinar as partições do dataset de treino e estimar o erro de validação.

log_res <- fit_resamples(log_wflow, k_folds)

O desempenho do modelo nos dados de treino é o seguinte.

log_res %>%
  collect_metrics()

# A tibble: 2 × 6
  .metric  .estimator  mean     n std_err .config             
  <chr>    <chr>      <dbl> <int>   <dbl> <chr>               
1 accuracy binary     0.814    10  0.0184 Preprocessor1_Model1
2 roc_auc  binary     0.872    10  0.0119 Preprocessor1_Model1

Os resultados apresentados nesta seção resultam do ajustamento aos dados de treino. Mas como vimos antes, o que pretendemos é usar este modelo com dados novos. Será que com novos dados este modelo será adequado ? Vamos usar os dados e o modelo especificado, mas agora com procedimentos que nos permitirão ter uma estimativa do erro generalizavel.

11.3.3 Teste

A estimativa do erro de generalização será o seguinte:

# ajustamento nos dados de teste
log_finalfit <- last_fit(log_wflow, particao) 

# guardar resultado final 
log_results <- collect_metrics(log_finalfit)
log_results

# A tibble: 2 × 4
  .metric  .estimator .estimate .config             
  <chr>    <chr>          <dbl> <chr>               
1 accuracy binary         0.804 Preprocessor1_Model1
2 roc_auc  binary         0.876 Preprocessor1_Model1

11.4 \(k\)-NN

11.4.1 Introdução

Na sub-secção Seção 10.1.3.2 vimos como funciona genericamente o algoritmo \(k\)-NN. Nesta sub-secção vamos implementar o algoritmo. Vamos introduzir no processo de treino uma afinação para o valor de \(k\) (número de vizinhos mais próximos) para minimizar o erro (1-accuracy) de classificação.

Sendo um algoritmo de classificação (e regressão) simples, só funciona bem se a distâncias (e.g. euclideanas) refletirem bem a dissemelhança entre os exemplos (quanto mais afastadas estiverem entre si, menor a probabilidade de serem parecidas).

11.4.2 Treino

Especificar o modelo \(k\)-NN, usando a estrutura tidy:

knn_spec <-
  nearest_neighbor(neighbors = tune()) %>% 
  set_engine("kknn") %>% 
  set_mode("classification")

Vamos usar o argumento tune() para “afinar” o valor do hiperparâmetro \(k\) através de validação cruzada. Para tal temos de definir uma grelha (neste caso é apenas um vector).

k_grid <- tibble(neighbors = seq(1,51, by = 2))

Vamos juntar o modelo à receita especificada na Seção 11.2, integrando tudo num workflow:

knn_wflow <- workflow() %>%
  add_recipe(receita) %>%
  add_model(knn_spec)

Agora que temos o processo pronto para o treino adicionamos a afinação do parâmetro \(k\) para fazer previsões:

knn_res <- 
  knn_wflow %>%
  tune_grid(resamples = k_folds, 
            grid = k_grid)

Qual é o valor de \(k\) que mais melhora o desempenho do algortimo?

# top-5
show_best(knn_res, metric="accuracy")

# A tibble: 5 × 7
  neighbors .metric  .estimator  mean     n std_err .config              
      <dbl> <chr>    <chr>      <dbl> <int>   <dbl> <chr>                
1        15 accuracy binary     0.818    10  0.0115 Preprocessor1_Model08
2        13 accuracy binary     0.817    10  0.0122 Preprocessor1_Model07
3        43 accuracy binary     0.815    10  0.0123 Preprocessor1_Model22
4        39 accuracy binary     0.815    10  0.0127 Preprocessor1_Model20
5        19 accuracy binary     0.815    10  0.0117 Preprocessor1_Model10

# top-1
select_best(knn_res, metric="accuracy")

# A tibble: 1 × 2
  neighbors .config              
      <dbl> <chr>                
1        15 Preprocessor1_Model08

Vamos guardar o resultado num objeto.

best_knn <- select_best(knn_res, metric="accuracy")

O último passo é usar a função finalize_workflow() para adicionar ao nosso workflow o modelo com o hiperparâmetro \(k\) melhor afinado:

final_knn_wflow <- 
  knn_wflow %>% 
  finalize_workflow(best_knn)

Este objecto guarda a versão atualizada da função aprendida com os dados de treino/validação.

11.4.3 Teste

Finalmente ajustamos o modelo aos dados de teste:

knn_finalfit <- final_knn_wflow %>% 
           last_fit(split = particao)

A função last_fit() ajusta o modelo final ao conjunto de dados de treino, e avalia o desempenho nos dados de teste. Em baixo, a função collect_predictions() permite visualizar num tibble os resultados da classificação com dados de teste

knn_results <-  
  knn_finalfit %>% 
  collect_metrics()
knn_results

# A tibble: 2 × 4
  .metric  .estimator .estimate .config             
  <chr>    <chr>          <dbl> <chr>               
1 accuracy binary         0.796 Preprocessor1_Model1
2 roc_auc  binary         0.873 Preprocessor1_Model1

11.5 Árvores de decisão

11.5.1 Introdução

As árvores de decisão são métodos de aprendizagem que foram desenvolvidos na década de 1980, e podem ser aplicados em problemas de classificação ou de regressão. São simples de implementar e úteis para interpretar.

Esquema de uma árvore de decisão. Fonte: M.R. Camana et al (2020) ‘Extremely Randomized Trees-Based Scheme for Stealthy Cyber-Attack Detection in Smart Grid Networks’, IEEE Access, vol. 8, pp. 19921-19933

O espaço dos preditores é segmentado/estratificado num conjunto de regiões. As regras usadas para segmentar o espaço de preditores podem ser sintetizadas através de uma árvore. Daí o nome árvores de decisão.

Mais tarde (decada de 1990) apareceram métodos como o bagging ou random forests , que integram árvores de decisão no processo de aprendizagem. Hoje em dia as árvores de decisão são usadas como partes ou módulos das ferramentas de aprendizagem (e.g. ensembles de árvores de decisão).

Vamos ver a intuição antes de vermos uma árvore. No gráfico seguinte temos os dados da Seção 10.1.3.2. São dados de 120 pacientes que têm um tumor benigno (azul) ou tumor maligno (vermelho). Para cada paciente temos os valores de dois preditores: o tamanho do tumor e a idade do paciente.

Queremos criar um modelo (a árvore) para prever se um tumor é maligno ou benigno a partir destes 2 preditores. O gráfico seguinte ilustra a distribuição dos exemplos classificados no espaço dos preditores.

Se quisermos separar os pacientes com tumor maligno dos restantes, o que poderemos fazer? Olhando para o gráfico parece que os pacientes mais novos tendem a ter tumores benignos e os mais velhos tumores malignos.

Vamos tentar separar o espaço dos preditores em 2 regiões, uma com tumores benignos e outra com os tumores malignos.

Esta divisão parece ter separado relativamente bem os tumores benignos dos malignos. Os pacientes mais novos (até \(\approx 44\) anos) tem mais tumores benignos e os mais velhos tem mais tumores malignos (o que faz sentido, pois o risco de ter um tumor maligno aumenta com a idade).

Apesar da separação fazer sentido, ainda há mistura de tumores benignos e malignos, talvez mais evidente na região superior do gráfico. Podemos melhorar esta separação, dividindo a região da seguinte forma:

Ficamos agora com o espaço dos preditores dividido em três regiões:

região do lado superior direito, com a grande parte dos tumores malignos, cujos pacientes são mais velhos, e têm tumores de maior dimensão;
região do lado inferior, com larga maioria de tumores benignos detetados em pacientes mais jovens;
região do lado superior esquerdo, onde estão a grande maioria de tumores benignos, observados em pacientes mais velhos.

A separação nestas três regiões permitiu segmentar o espaço dos preditores em regiões, seguindo o racional das árvores de decisão. A figura seguinte ilustra a árvore de decisão obtida com os mesmos dados:

Vamos compreender e interpretar o gráfico da árvore. A sua construção começa no nó inicial onde estão 100% dos dados. O algoritmo escolhe um preditor e define um valor que melhor separa os exemplos em duas partições - o preditor escolhido é a idade e o valor 44 anos. Do lado esquerdo ficam os exemplos com idade inferior a 44. Do lado direito ficam os restantes. Na partição da direita procede-se a uma nova separação: entre os pacientes com idade \(\ge 44\) anos, os que têm um tumor com menos de 2.5 cm ficam do lado esquerdo e os restantes do lado direito.

Em cada um dos nós estão indicadas a classe prevista e a proporção dos exemplos associada.

A terminologia usada nas árvores de decisão é a seguinte:

o nó inicial onde estão 100% dos dados designa-se por nó inicial (root node);
os nós intermédios designam-se por nós-filhos (child nodes);
os nós finais designam-se por nós-terminais (leaf nodes).

A interpretação dos resultados da árvore é a seguinte:

a primeira separação foi feita com a idade, portanto esse será o preditor mais importante para prever a malignidade do tumor;
sendo um paciente mais jovem, o tamanho do tumor não parece ter um papel muito relevante na malignidade do tumor;
mas entre pacientes mais velhos, o tamanho do tumor parece ter impacto na malignidade: quanto menor o tamanho, maior a probabilidade de ser benigno.

Estas partições especificadas pelo algoritmo correspondem à segmentação que fizemos ‘a olho’ para separar o espaço dos preditores em três regiões. Mas como o oalgoritmo escolheu os preditores? e como sabe que valores usar para separar os dados (os exemplos)? Embora tenhamos visto intuitivamente como responder, falta saber como é que o algortimo o faz de forma automática.

A árvore de decisão é um modelo que:

Divide o espaço de preditores em regiões distintas (não sobreponíveis)
Em cada região a classe prevista é a que tem probabilidade máxima.
A classe prevista para uma observação que caia numa região é a que tem probabilidade máxima nessa região.

Em teoria o espaço de preditores pode ser dividido de muitas maneiras. Nas árvores o espaço é dividido em rectangulos multidimensionais de forma a minimizar o erro de classificação.

As medidas preferenciais para medir esse erro são o índice de Gini ou a entropia e ganho de informação.

Índice de Gini

O indice de Gini é dado por: \[ \text{G} = \sum_{k=1}^K \hat{p}_{rk}(1-\hat{p}_{rk}) \] onde \(\hat{p}_{rk}\) é a proporção de exemplos na região \(r\) que pertencem à classe \(k\). Quanto maior a proporção de exemplos de uma só classe, menor o valor do índice nesse nó.

O índice de Gini é também conhecido como índice de impureza de Gini. Um valor alto do índice indica que não há uma classe predominante no nó, ou seja, é uma região “impura” com mistura de classes entre os exemplos.

A entropia e ganho são alternativas ao índice de Gini mas com interpretação similar (são medidas similares).

Entropia e ganho

A entropia \(E\) do conjunto de exemplos \(s\), é dada por: \[ E(s) = -\sum_{k=1}^K \hat{p}_{rk}\log_2 (\hat{p}_{rk}) \] onde \(\hat{p}_{rk}\) é a proporção de exemplos na região \(r\) que pertencem à classe \(k\).

Quando usamos um nó para fazer uma nova partição, a entropia muda. O ganho de informação é a medida que mede a redução no valor da entropia.

O ganho de informação \(g\) com nova partição usando o preditor \(v\) é dado por: \[ \text{g}(s, v) = E(s) - E(s, v) \] Quanto maior o ganho, mais puro fica o nó.

Infelizmente, em termos computacionais é pesado percorrer todas as partições do espaço possiveis. Por essa razão, o processo de divisão de um nó é feito de cima para baixo (começa no root-node), de forma gananciosa criando em cada nó a partição que melhor resultado dá (de acordo com o indice de Gini ou a entropia e ganho). Este processo de construção designa-se por separação recursiva binária (recursive binary splitting). Repetimos o procedimento até um critério de paragem ser atingido. Podemos impôr limitações à construção da árvore como por ex. parar quando tivermos no máximo 5 exemplos em cada nó-filho (pre-prunning).

O procedimento poderá produzir boas previsões com os dados de treino, mas corre o risco de produzir overfitting. Isto porque uma árvore comprida (i.e. com muitos nós terminais), aprende muito bem os padrões nos dados de treino mas acaba por errar em demasia quando prevê em dados teste. Uma árvore menos complexa poderá produzir melhores resultados (é uma medida de compromisso entre a variância e o viés).

Flexibilidde numa árvore de decisão

Mede-se pelo número de nós terminais (leaf nodes). Quanto maior o número, mais comprida é a árvore e mais complexo/flexivel será o modelo, aumentando o risco de overfitting.

Outra abordagem para procurar atingir um modelo com comprimento ‘ótimo’, i.e., que atinge um ponto equilibrio entre a variància e o viés (entre o overfitting e o underfitting) passa por criar uma árvore longa e depois “podar” alguns ramos da árvore (post-prunning) obtendo uma árvore menor à custa do aumento do viés. No entanto, faze-lo de forma exaustiva é uma solução computacionalmente exigente porque o número de maneiras diferentes de podar árvores pode ser incomportável.

Em alternativa podemos introduzir um hiperparâmetro que adiciona um custo (uma penalização) a árvores mais complexas. Um custo próximo de zero representa uma baixa penalização, i.e., com poucos nós podados (que poderá levar ao overfitting). Um custo elevado irá podar muitos nós da árvore (no limite irá podar todos) e provocar um bias elevado (underfitting).

Este hiperpârametro, designado por custo de complexidade (cost complexity prunning) pode ser escolhido por validação cruzada (ou conjunto de validação).

Custo de complexidade

Seja \(R(T)\) uma medida global de impureza da árvore (e.g. indice de Gini ou Entropia) numa árvore com \(T\) nós terminais. Pretende-se minimizar: \[ R(T)+\alpha|\tilde{T}| \] onde \(\alpha\) é o termo de penalização (estimado por validação cruzada) ou custo e \(|\tilde{T}|\) o número de nós terminais (leaf-nodes) da árvore (que mede a complexidade).

Vamos implementar o algoritmo com o nosso dataset.

11.5.2 Treino

Vamos especificar a árvore de decisão,com um argumento adicional para afinar o hiperparâmetro \(\alpha\) que penaliza modelos mais complexos.

dtree_spec <-
  decision_tree(cost_complexity= tune()) %>% 
  set_engine("rpart") %>%
  set_mode("classification")

Vamos criar o workflow para juntar as tarefas de Seção 11.2 e as especificações para a árvore de decisão.

dtree_wflow <- workflow() %>%
  add_recipe(receita) %>%
  add_model(dtree_spec)

Vamos criar uma grelha de valores para afinar o hiperparâmetro.

dtree_grid <- grid_regular(cost_complexity(), levels = 100)

Agora que temos as especificações todas preparadas, vamos afinar o hiperparâmetro:

dtree_res <- 
  dtree_wflow %>% 
  tune_grid(resamples = k_folds, grid = dtree_grid)

A exatidão obtida para os cinco melhores modelos:

# top-5
show_best(dtree_res, metric="accuracy")

# A tibble: 5 × 7
  cost_complexity .metric  .estimator  mean     n std_err .config               
            <dbl> <chr>    <chr>      <dbl> <int>   <dbl> <chr>                 
1         0.00658 accuracy binary     0.807    10  0.0130 Preprocessor1_Model087
2         0.00152 accuracy binary     0.805    10  0.0119 Preprocessor1_Model080
3         0.00187 accuracy binary     0.804    10  0.0119 Preprocessor1_Model081
4         0.00231 accuracy binary     0.804    10  0.0119 Preprocessor1_Model082
5         0.00285 accuracy binary     0.804    10  0.0119 Preprocessor1_Model083

Vamos guardar num tbl_df (um dataframe com estrutura tibble) e atualizar o modelo especificando o valor do hiperparâmtero que devolveu maior exatidão.

# top-1
select_best(dtree_res, metric="accuracy")

# A tibble: 1 × 2
  cost_complexity .config               
            <dbl> <chr>                 
1         0.00658 Preprocessor1_Model087

best_dtree <- select_best(dtree_res, metric="accuracy")

final_dtree_wflow <- 
  dtree_wflow %>% 
  finalize_workflow(best_dtree)

11.5.3 Teste

Finalmente ajustamos o modelo aos dados de teste:

dtree_finalfit <- final_dtree_wflow %>% 
           last_fit(split = particao)

# guardar resultados de teste
dtree_results <-  
  dtree_finalfit %>% 
  collect_metrics()
dtree_results

# A tibble: 2 × 4
  .metric  .estimator .estimate .config             
  <chr>    <chr>          <dbl> <chr>               
1 accuracy binary         0.809 Preprocessor1_Model1
2 roc_auc  binary         0.832 Preprocessor1_Model1

E visualizamos o modelo obtido com a função rpart.plot():

dtree_finalfit %>%
  extract_fit_engine() %>%
  rpart.plot(roundint = FALSE)

As árvores de decisão tem vantagens em relação a métodos mais clássicos como a regressão logística:

São fáceis de entender e de explicar
Reflectem bem o modo como tomamos decisões
Resultados podem ser visualizados e interpretados (até por um não-especialista)

No entanto tem um ponto fraco não desprezável: é que uma pequena alteração nos dados pode causar uma grande alteração no classificador final.

O desempenho desta familia de classificadores pode no entanto ser significativamente melhorado usando métodos de bagging, boosting ou random forests. Na sub-secções seguintes vamos ver bagging e random forests.

11.6 Bagging

11.6.1 Introdução

Embora seja uma solução interessante para reduzir a variância do modelo, podar as árvores (via pre-prunning ou post-punning) faz aumentar o viés do modelo.

Na década de 1990 aparecem novos métodos de aprendizagem automática que tiram partido das árvores de decisão e usam métodos como o bagging (ou agregação por bootstrapping) e random forests para reduzir a variância do modelo sem necessidade de podar árvores (viés reduzido)!

Bootstrapping

Bootstrapping é um método computacional de reamostragem criado para gerar multiplos subconjuntos de datasets a partir de um dataset original - o dataset de treino. Em cada dataset gerado (bootstrapped), as observações - exemplos - são selecionadas aleatoriamente com reposição. Isto significa que um exemplo pode aparecer mais que uma vez num dataset.

O bagging passa por:

gerar \(n\) datasets independentes a partir do dataset de treino (bootstrapping);
treinar um modelo (árvore) para cada dataset gerado sem a preocupação de o podar;
usar cada modelo treinado para prever o valor de uma nova amostra;
prever a classe da nova amostra atribuindo a classe mais “votada” pelos modelos treinados

Ao conjunto de árvores aprendidas com o conjunto de datasets gerados por bootstrap, designa-se por ensemble de árvores.

Da teoria estatística⁷ sabemos que se tivermos um grande número (\(N\)) datasets de treino independentes com variância \(\sigma^2\), a variância da média será \(\sigma^2/N\). Sendo assim, em vez de podarmos as árvores para reduzir a variância, podemos construir árvores compridas, complexas, que tendem a ter baixo viés e usar métodos como o bagging para reduzir a variância do modelo.

11.6.2 Treino

Vamos implementar este algoritmo, recorrendo à biblioteca baguette que complementa a biblioteca tidymodels para implementação de bagging.

Especificação do modelo de bagging com 50 amostras, ajustar os modelos e ver os resultados:

# biblioteca para bagging com tidymodels
library(baguette)

set.seed(992801000)
bag_spec <- bag_tree() %>%
  set_engine("rpart", times = 50) %>% # 50 amostras bootstrap
  set_mode("classification")
 
bag_wflow <- workflow() %>%
  add_recipe(receita) %>%
  add_model(bag_spec)

bag_fit <- bag_wflow %>% fit(heartdtr)

Dado que o processo de bagging envolve o ajustamento de muitas árvores, deixa de ser possível usar as árvores para interpretar resultados. Podemos no entanto ver outras coisas. Na sub-secção Seção 11.5.1 vimos que podemos estimar a importância de um preditor no modelo com medidas de erro (e.g. o índice de Gini).

No bagging o processo é similar: em cada árvore, em cada partição, calcula-se o índice de Gini para cada preditor, agregam-se os valores de todas as árvores por preditor. Os preditores com índice de Gini médio mais baixo são considerados os mais importantes (maior contribuição para aumentar a pureza das leaf-nodes).

══ Workflow [trained] ══════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: bag_tree()

── Preprocessor ────────────────────────────────────────────────────────────────
2 Recipe Steps

• step_normalize()
• step_dummy()

── Model ───────────────────────────────────────────────────────────────────────
Bagged CART (classification with 50 members)

Variable importance scores include:

# A tibble: 10 × 4
   term               value std.error  used
   <chr>              <dbl>     <dbl> <int>
 1 Cholesterol         94.3     1.97     50
 2 MaxHR               92.9     1.65     50
 3 HeartPeakReading    87.1     1.77     50
 4 Angina_Y            86.6     2.16     50
 5 Age                 81.3     1.37     50
 6 RestingBP           59.3     1.34     50
 7 FastingBS_Diabetes  19.8     1.10     50
 8 Sex_M               11.6     0.741    50
 9 RestingECG_ST       11.0     0.695    50
10 RestingECG_Normal   10.3     0.523    50

Vemos que os mais relevantes preditores para a ocorrência de doença cardiaca são Cholesterol, MaxHR e HearPeakReading.

11.6.3 Teste

Qual será o erro de classificação do algortimo nos dados de teste?

# teste
bag_rs <- heartdtt %>%
  bind_cols(predict(bag_fit, heartdtt)) 

bag_results <- bag_rs %>%
  accuracy(HeartDisease, .pred_class)
bag_results

# A tibble: 1 × 3
  .metric  .estimator .estimate
  <chr>    <chr>          <dbl>
1 accuracy binary         0.796

11.7 Random forests

11.7.1 Introdução

As random forests trazem uma pequena alteração ao algoritmo de bagging que se traduz numa melhoria do algoritmo.

Tal como no bagging modelamos árvores de decisão com datasets gerados por bootstrap a partir do dataset original. Mas ao modelar as árvores, em vez de considerarmos todos os preditores em cada partição, consideramos apenas alguns candidatos (selecionados por amostragem aleatoria simples). Apenas um destes candidatos será usado na partição.

Em cada partição possível, um novo conjunto de preditores candidatos (de dimensão \(m\)) é selecionado por amostragem a partir de todos os candidatos (de dimensão \(p\)). Tipicamente o número de candidatos escolhidos é \(\approx m = \sqrt{p}\).

Escolha dos preditores candidatos

Apesar de usarmos bootstap para gerar datasets independentes, as árvores obtidas com bagging tendem a estar bem correlacionadas, porque todas usam os mesmos preditores. Imagine um caso em que um único preditor tem uma importância muito elevada em relação a todos os outros. Com bagging a grande maioria (ou todas) das árvores iriam usar esse preditor no root-node e consequentemente as árvores seriam parecidas entre si.

Com random forests, procura-se “aliviar” esta correlação entre modelos. Ao não usar todos os preditores em cada partição, as árvores são levadas a usar diferentes preditores a separar em diferentes alturas. Este procedimento reduz a correlação entre os modelos (árvores), tornando-os mais independentes entre si. Como usamos um número grande de árvores neste processo, e depois calculamos a média, este procedimento tende a melhorar o desempenho do modelo em relação ao ensemble de árvores obtido por bagging.

Vamos usar random forests para prever se um paciente tem uma doença cardíaca. Vamos especificar os seguintes hiperparamteros: número de preditores a considerar

11.7.2 Treino

Para além de especificarmos o modelo usando a estrutura habitual, vamos adicionar um um hiperparametro para escolher o número de preditores adequado a ser escolhido aleatoriamente em cada partição do ensemble de árvores - mtry=tune().

# doParallel::registerDoParallel()

# especificacao do modelo
rf_spec <- rand_forest( mtry = tune(), trees = 500) %>%
  set_mode("classification") %>%
  set_engine("ranger")

# workflow 
rf_wflow <- workflow() %>%
  add_recipe(receita) %>%
  add_model(rf_spec)

# grid de afinação
rf_grid <- grid_regular(mtry(range = c(1,4)), levels = 5)

set.seed(33889922)
rf_res <- 
  rf_wflow %>% 
  tune_grid(resamples = k_folds, grid = rf_grid)

# top-5
show_best(rf_res, metric="accuracy")

# A tibble: 4 × 7
   mtry .metric  .estimator  mean     n std_err .config             
  <int> <chr>    <chr>      <dbl> <int>   <dbl> <chr>               
1     2 accuracy binary     0.826    10  0.0135 Preprocessor1_Model2
2     1 accuracy binary     0.818    10  0.0130 Preprocessor1_Model1
3     3 accuracy binary     0.818    10  0.0135 Preprocessor1_Model3
4     4 accuracy binary     0.818    10  0.0120 Preprocessor1_Model4

# top-1
select_best(rf_res, metric="accuracy")

# A tibble: 1 × 2
   mtry .config             
  <int> <chr>               
1     2 Preprocessor1_Model2

best_rf <- select_best(rf_res, metric="accuracy")

# grafico tunning
rf_res %>%
  collect_metrics() %>%
  filter(.metric == "accuracy") %>%
  ggplot(aes(mtry, mean)) +
  geom_line(alpha = 0.5, size = 1.5) +
  geom_point() +
  labs(y = "Accuracy") +
  labs(x = "# preditores")

final_rf_wflow <- 
  rf_wflow %>% 
  finalize_workflow(best_rf)

Ajustamos o modelo aos dados de teste e guardamos os resultados:

rf_finalfit <- final_rf_wflow %>% 
           last_fit(split = particao)

# guardar resultados de teste
rf_results <-  
  rf_finalfit %>% 
  collect_metrics()

# matriz confusao dados teste
rf_finalfit %>%
    collect_predictions() %>%
    conf_mat(HeartDisease, .pred_class)

          Truth
Prediction   N   Y
         N  80  19
         Y  23 108

11.8 Componentes Principais⁸

A análise em componentes principais (ACP) é uma ferramenta de aprendizagem não-supervisionada muito usada em múltiplas áreas (e.g. neurociencias, computação, economia, química,…), pois permite extrair informação relevante a partir de um dataset complexo (i.e. com muitas observações e muitas variáveis).

O método permite resumir a informação relevante do dataset original com um número reduzido de variáveis que coletivamente representam a maior parte da variabilidade encontrada no dataset original.

A ideia-chave é que as variáveis não são todas igualmente ‘interessantes’, não têm todas a mesma quantidade de informação, e que podemos sumarizar a informação relevante com um número (muito) inferior de variáveis (procedendo a uma transformação adequada das variáveis originais)

Componentes principais

As componentes principais \(z_i\) são tranformações obtidas por combinação linear das variáveis originais \(x\). \[ z_i = \phi_{0i} + \phi_{1i} x_1 + \phi_{2i} x_2 + \ldots +\phi_{pi} x_p, \quad i=\{1,\ldots,p\} \] As componentes principais são construídas de tal forma que um sub-conjunto contém aproximadamente a mesma informação contida no dataset original.

Vamos explorar uma análise em componentes principais com o dataset NCI60 disponível na biblioteca ISLR. É um dataset com dados relativos ao genoma de células cancerígenas, que consiste em 6830 medidas de expressão dos genes observadas em 64 linhas celulares de cancro.

library(tidymodels)
library(tidyclust)
library(factoextra)
library(ISLR)

Os dados vêm numa lista com dois objetos, que contém a matriz de dados e os labels das linhas celulares.

data(NCI60, package = "ISLR")
nci60 <- NCI60$data %>%
  as_tibble(.name_repair = ~ paste0("V_", .x)) %>%
  mutate(label = factor(NCI60$labs)) %>%
  relocate(label)

Cada linha está classificada com o tipo de cancro (coluna label). Não vamos usar esta variável (vamos esconder esta variável por agora), pois vamos emular um exercício de aprendizagem não-supervisionada. Mas depois de implementarmos a ACP iremos verificar até que ponto os tipos de cancro são concordantes com os resultados da ACP.

Assim vamos retirar esta variável da análise e usar a função prcomp() para realizar a ACP. Devemos ainda especificar o argumento scale=TRUE para normalizar as variáveis.

nci60_pca <- nci60 %>%
  select(-label) %>%
  prcomp(scale = TRUE)

A quantidade de informação explicada/retida em cada componente principal é dada pelo objecto eigenvalue. Podemos usar a função tidy() with matrix = "eigenvalues" visualizar a proporção da variabilidade (absoluta e acumulada) retida em cada componente (por ordem decrescente):

tidy(nci60_pca, matrix = "eigenvalues")

# A tibble: 64 × 4
      PC std.dev percent cumulative
   <dbl>   <dbl>   <dbl>      <dbl>
 1     1    27.9  0.114       0.114
 2     2    21.5  0.0676      0.181
 3     3    19.8  0.0575      0.239
 4     4    17.0  0.0425      0.281
 5     5    16.0  0.0374      0.318
 6     6    15.7  0.0362      0.355
 7     7    14.5  0.0307      0.385
 8     8    13.5  0.0269      0.412
 9     9    13.1  0.0253      0.438
10    10    12.7  0.0238      0.461
# ℹ 54 more rows

A primeira componente (PC1) retém 11% e as 10 primeiras componentes retêm 46% da variabilidade total dos dados. As 12 primeiras componentes explicam 50% da variabilidade observada nas 64 linhas de celulas tumorais. Podemos visualizar os resultados acumulados no gráfico seguinte:

tidy(nci60_pca, matrix = "eigenvalues") %>%
  ggplot(aes(factor(PC), percent)) +
  geom_bar(stat="identity", aes(fill = percent)) +
  scale_x_discrete(breaks = c(1, 10, 20 ,30,40, 50, 60)) +
  scale_y_continuous(breaks = c(0, 0.05,0.1)) +
  xlab("Componentes principais") +
  ylab("Proporção da variabilidade explicada")

Ou podemos ver a proporção da variabilidade acumulada:

tidy(nci60_pca, matrix = "eigenvalues") %>%
 ggplot(aes(factor(PC), cumulative)) +
  geom_bar(stat="identity", aes(fill = percent)) +
  scale_x_discrete(breaks = c(1, 10, 20 ,30,40, 50, 60)) +
  scale_y_continuous(breaks = c(0, .2,.4, .6,.8,1)) +
  xlab("Componentes principais") +
  ylab("Proporção da variabilidade explicada")

Para fins de visualização, vamos agora juntar o tipo de cancro (label) ao resultado de augment(nci60_pca) para podermos visualizar quão semelhantes são as linhas celulares dos diferentes tipos de cancro.

nci60_pcs <- bind_cols(
  augment(nci60_pca),
  nci60 %>% select(label)
)

Como temos 14 tipos de tumor diferentes, vamos usar uma palete de cores adequada (“Polychrome 36”) que nos ajude a ver as diferenças/semelhanças entre os diferentes tipos de tumor.

# n de categorias
nlevels(nci60$label)

[1] 14

colors <- unname(palette.colors(n = 14, palette = "Polychrome 36"))

Podemos visualizar gráficos de dispersão com pares de componentes principais. É usual visualizar as primeiras componentes para identificar as tendencias. A informação retida nas restantes componentes poderá revelar aspetos particulares que podem ser interessantes analizar. O gráfico seguinte opõe as duas primeiras componentes principais:

nci60_pcs %>%
  ggplot(aes(.fittedPC1, .fittedPC2, color = label)) +
  geom_point(size =3, alpha = 0.5) +
  scale_color_manual(values = colors)

A grande maioria dos melanomas tendem a estar proximos entre si. A mesma observação pode ser feita para os tumores do colon.

O gráfico seguinte opõe a primeira e a terceira componente:

nci60_pcs %>%
  ggplot(aes(.fittedPC1, .fittedPC3, color = label)) +
   geom_point(size =3, alpha = 0.5) +
  scale_color_manual(values = colors)

Vamos experimentar usar um método de clustering, para ver como se agrupam os dados.

11.9 Clustering

Os métodos de clustering têm como objectivo a classificação (não supervisionada) e o reconhecimento de padrões de tal modo que promove:

agrupamento de observações “parecidas” na mesma classe,
assegura que os grupos (clusters) diferentes estão bem “separados”

Há vários algoritmos de clustering, vamos ver a classificação hierárquica.

Antes de começar vamos fazer um pre-processamento de dados, tal como fizemos no passado. Um dos aspetos relevantes para melhorar o desempenho das ferramentas de clustering é normalizar previamente os dados. Podemos usar a função recipe::recipe() para fazer a normalização dos dados e especificar um workflow() para combinar depois este resultado com o modelo de clustering.

nci60_rec <- recipe(~ ., data = nci60) %>%
  step_rm(label) %>%
  step_normalize(all_predictors())

nci60_wf <- workflow() %>%
  add_recipe(nci60_rec)

Vamos experimentar um método de classificação (não-supervisionada) hierarquico. Estes métodos produzem partições em classes sucessivamente mais abrangentes. Existem dois tipos de classificação hirarquica:

classificação hierárquica ascendente,
classificação hierárquica descendente,

Na prática a classificação hierárquica ascendente é mais usada e é a que vamos ver. Estes métodos são bastante usados pois o output gráfico é baseado em dendrogramas (árvores) que permitem uma leitura visual muito intuitiva dos resultados de classificação.

Algoritmo de classificação hierarquica

O algoritmo de classificação hierárquica ascendente começa por:

Calcular a matriz de dissemelhança entre classes,
Aglomerar as classes mais parecidas,
Recalcular a matriz de dissemelhança entre classes
Repetir passos 2. e 3. até restar apenas uma classe

Durante o processo de aglomeração (2.) é necessário definir um critério para aglomerar clusters. Vamos ver três destes critérios:

average
complete
single

Estes processos determinam as regras de agregação entre clusters, e têm geralmente um impacto relevante na classificação, que se manifesta por ex. na configuração do dendrograma.

Vamos usar o package tidyclust e a função tidyclust::hier_clust para parametrizar e executar a classificação:

Podemos visualizar os resultados. A livraria factoextra tem várias funções que permitem visualizar bons gráficos baseados no ggplot2. Vamos ver os resultados com os três critérios de aglomeração. O que pretendemos é uma inspeção visual de qual o que devolve uma separação mais “natural”.

nci60_complete %>%
  extract_fit_engine() %>%
  fviz_dend(main = "Complete")

nci60_average %>% 
  extract_fit_engine() %>%
  fviz_dend(main = "Average")

nci60_single %>%
  extract_fit_engine() %>%
  fviz_dend(main = "Single")

A classificação em clusters é determinada pelo utilizador. Se considerarmos k = 4 obtemos as seguintes separações:

nci60_complete %>%
  extract_fit_engine() %>%
  fviz_dend(k = 4, main = "Classf. hierarquica com k=4 classes")

Em vez de usarmos as 64 linhas de celulas para classificação, poderemos usar o resultado da pca (apenas com as 12 primeiras componentes, que explicam 50% da variabilidade) como dataset de input na classificação. Vamos usar mais uma vez a função recipes() para calcular as componentes principais:

nci60_pca_rec <- recipe(~ ., data = nci60) %>%
  step_rm(label) %>%
  step_normalize(all_predictors()) %>%
  step_pca(all_predictors(), num_comp = 12)

nci60_pca_wf <- workflow() %>%
  add_recipe(nci60_pca_rec)

Podemos agora ajustar o novo workflow:

nci60_pca <- nci60_pca_wf %>%
  add_model(hier_clust(linkage_method = "complete")) %>%
  fit(data = nci60)

Podemos agora visualizar os resultados com este dataset reduzido (64 linhas e 12 colunas) que poderão permitir visualizar algum padrão nos dados agora de menor dimensão.

nci60_pca %>%
  extract_fit_engine() %>%
  fviz_dend(k = 4, main = "Classf. hierarquica com 4 PCs")

11.10 FIM

Fonte: James, Gareth et al. An Introduction to Statistical Learning with Applications in R. New York :Springer, 2013↩︎
Fonte: James, Gareth et al. An Introduction to Statistical Learning with Applications in R. New York :Springer, 2013↩︎
Fonte: James, Gareth et al. An Introduction to Statistical Learning with Applications in R. New York :Springer, 2013↩︎
Adapt. de Andrew Ng 2017↩︎
https://www.tmwr.org/pre-proc-table↩︎
as colunas novas começam com o prefixo . para evitar escrever em cima de variveis originais com o mesmo nome↩︎
Ver e.g. distribuição amostral da média, teorema do limite central↩︎
Fonte: Adaptado de https://github.com/EmilHvitfeldt↩︎

8.1 INTRODUÇÃO

9 Programa da aula

10 Aprendizagem

10.1 Supervisionada

10.1.1 Regressão

10.1.2 Classificação

10.1.3 Bias-variance trade-off

Regressão

Classificação

Erro de generalização

10.1.4 Validação

10.2 Não-supervisionada

10.2.1 Exemplo ACP

10.2.2 Exemplo clustering

11 Aplicações

11.1 Análise exploratória

11.2 Pre-processamento

11.3 Modelo logístico

11.3.1 Introdução

11.3.2 Treino

11.3.3 Teste

11.4 \(k\)-NN

11.4.1 Introdução

11.4.2 Treino

11.4.3 Teste

11.5 Árvores de decisão

11.5.1 Introdução

11.5.2 Treino

11.5.3 Teste

11.6 Bagging

11.6.1 Introdução

11.6.2 Treino

11.6.3 Teste

11.7 Random forests

11.7.1 Introdução

11.7.2 Treino

11.8 Componentes Principais8

11.9 Clustering

11.10 FIM

11.8 Componentes Principais⁸