Indo além do PCA, v2

Tarssio Barreto

21 de abril de 2019

Utilizando um pouco do PCA da forma mais clássica:

Objetivo

Vamos utilizar o PCA com o intuíto de realizar a redução de dimensionalidade e posteriormente comparar os resultados obtidos para classificação, através do algoritmo de KNN, para os usos de água domésticos.

Realizando o PCA e analisando os resultados

Carregando os dados:

require(caret)
require(tidyverse)

dados <- load("data_041218_1138.RData")

dados <- features_atual %>% 
  group_by(categoria) %>% 
  dplyr::filter(categoria %in% c("Torneira Interna", "Bacia")) %>% 
  dplyr::sample_n(1000)

dados$categoria <- forcats::fct_drop(dados$categoria)


x <- dados[,c(2:9)]

Realizaremos, primeiro, o PCA através do passo-a-passo a seguir:

  1. Escalonando os dados:
x1 <- scale(x)
  1. Cálculo da covariância:
covariancia <- cov(x1)
  1. Autovetores e valores:
auto <- eigen(covariancia)

auto$vectors # Autovetores
##            [,1]       [,2]         [,3]        [,4]        [,5]
## [1,] -0.1831672 -0.5874009  0.152486675 -0.09862627 -0.19644973
## [2,] -0.1575816 -0.6009885  0.041117060 -0.23839736 -0.35667118
## [3,] -0.3730959 -0.3524196  0.049462290  0.30902317  0.77005447
## [4,] -0.2063803 -0.0776952 -0.954019260  0.14175917 -0.10506422
## [5,] -0.4306120  0.2112778  0.025088756 -0.51528797  0.22871441
## [6,] -0.4438831  0.2233538  0.065592627 -0.07536626 -0.22859947
## [7,] -0.4235998  0.1518765  0.239778137  0.68503151 -0.35101376
## [8,] -0.4448919  0.2143483  0.005704431 -0.27817019 -0.04890946
##             [,6]        [,7]        [,8]
## [1,] -0.11642232  0.71451867 -0.16045897
## [2,]  0.02764619 -0.63261793  0.16514500
## [3,]  0.17440122 -0.11984803 -0.02991271
## [4,] -0.07759416  0.06007287 -0.02130979
## [5,] -0.66325605 -0.07555143 -0.07653110
## [6,]  0.43428249 -0.10408692 -0.70100309
## [7,] -0.33566025 -0.05123078  0.17953275
## [8,]  0.45614385  0.22829505  0.64508668
auto$values # Autovalores
## [1] 4.28442840 2.34829484 0.87554281 0.21013798 0.14628223 0.08280511
## [7] 0.03192991 0.02057873
  1. Predizendo os valores no PCA:
pred <- as.matrix(x) %*% as.matrix(auto$vectors)

head(pred)
##           [,1]        [,2]     [,3]       [,4]       [,5]       [,6]
## [1,] -34.47144  -80.368912 19.65657 -12.176037 -24.313779 -14.572499
## [2,] -32.23194  -72.612858 19.35016 -11.306264 -22.145146 -14.772002
## [3,] -42.22146 -118.250253 29.88610 -19.084559 -37.215254 -21.167016
## [4,] -51.84524 -122.810760 26.82344 -18.324622 -34.644219 -20.335443
## [5,] -19.87064   -7.287271  5.14022  -2.986156  -5.949638  -2.773741
## [6,] -22.36313  -41.918539 10.00206  -7.005171 -13.920866  -8.380698
##           [,7]       [,8]
## [1,]  87.72018 -19.495619
## [2,]  81.71832 -18.186036
## [3,] 128.68782 -28.926727
## [4,] 124.30626 -27.885662
## [5,]  12.72801  -2.588439
## [6,]  47.24702  -9.942864

Exploraremos agora, brevemente, alguns opções de projeção do PCA:

Usaremos para isto o pacote factoextra:

library(factoextra)

Criando nosso objeto PCA:

x_pca <- prcomp(x, scale = TRUE)

summary(x_pca)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6
## Standard deviation     2.0699 1.5324 0.9357 0.45841 0.38247 0.28776
## Proportion of Variance 0.5355 0.2935 0.1094 0.02627 0.01829 0.01035
## Cumulative Proportion  0.5355 0.8291 0.9385 0.96480 0.98309 0.99344
##                            PC7     PC8
## Standard deviation     0.17869 0.14345
## Proportion of Variance 0.00399 0.00257
## Cumulative Proportion  0.99743 1.00000
x_pca$rotation
##                PC1        PC2          PC3         PC4         PC5
## duracao -0.1831672  0.5874009 -0.152486675  0.09862627 -0.19644973
## nmoda   -0.1575816  0.6009885 -0.041117060  0.23839736 -0.35667118
## volume  -0.3730959  0.3524196 -0.049462290 -0.30902317  0.77005447
## inercia -0.2063803  0.0776952  0.954019260 -0.14175917 -0.10506422
## moda    -0.4306120 -0.2112778 -0.025088756  0.51528797  0.22871441
## media   -0.4438831 -0.2233538 -0.065592627  0.07536626 -0.22859947
## pico    -0.4235998 -0.1518765 -0.239778137 -0.68503151 -0.35101376
## mediana -0.4448919 -0.2143483 -0.005704431  0.27817019 -0.04890946
##                 PC6         PC7         PC8
## duracao -0.11642232 -0.71451867 -0.16045897
## nmoda    0.02764619  0.63261793  0.16514500
## volume   0.17440122  0.11984803 -0.02991271
## inercia -0.07759416 -0.06007287 -0.02130979
## moda    -0.66325605  0.07555143 -0.07653110
## media    0.43428249  0.10408692 -0.70100309
## pico    -0.33566025  0.05123078  0.17953275
## mediana  0.45614385 -0.22829505  0.64508668

É interessante, em momento de exercicio, comparar o obtido através do algoritmo e o resultado da função prcomp. É possível perceber, também, que com as três componentes principais podemos explicar cerca de 90% da variância dos dados, reduzindo de forma significativa o numéro de dimensões do nosso problema. Vamos, enfim, explorar a projeção das nossas componentes principais.

factoextra::fviz_pca_biplot(x_pca, repel = FALSE)

Com o gráfico acima verificamos a direção de alguma das variáveis e como os dados se comportam após em função dos dois primeiros componentes gerados pela PCA.

Podemos, ver também como acontece a projeção, tendo em vista as classificações dos usos domésticos de água:

factoextra::fviz_pca_ind(x_pca,
             label = "none", # hide individual labels
             habillage = dados$categoria, # color by groups
             addEllipses = TRUE # Concentration ellipses
             )

Os resultados não apontam um agrupamento conclusivo, testemos, de forma expedita, o ganho de acurácia utilizando o PCA no processamento dos dados:

Definindo grupo de teste e treino:

# Definindo Split

set.seed(17)

split=0.80

trainIndex <- caret::createDataPartition(dados$categoria, p=split, list=FALSE)

# Separando 

data_train <- dados[ trainIndex,c(1:8)]

data_test <- dados[-trainIndex,c(1:8)]

Modelo 1 : Sem pca

model <- caret::train(categoria ~ ., data = data_train,
      method = "knn", preProcess = c("center", "scale"), 
      tuneLength = 10, 
      trControl = trainControl(method="cv", number=10))



pred <- predict(model, data_test)

caret::confusionMatrix(pred, data_test$categoria) 
## Confusion Matrix and Statistics
## 
##                   Reference
## Prediction         Bacia Torneira Interna
##   Bacia              175                7
##   Torneira Interna    25              193
##                                           
##                Accuracy : 0.92            
##                  95% CI : (0.8889, 0.9446)
##     No Information Rate : 0.5             
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.84            
##                                           
##  Mcnemar's Test P-Value : 0.002654        
##                                           
##             Sensitivity : 0.8750          
##             Specificity : 0.9650          
##          Pos Pred Value : 0.9615          
##          Neg Pred Value : 0.8853          
##              Prevalence : 0.5000          
##          Detection Rate : 0.4375          
##    Detection Prevalence : 0.4550          
##       Balanced Accuracy : 0.9200          
##                                           
##        'Positive' Class : Bacia           
## 

Modelo 2 : Com pca

model2 <- caret::train(categoria ~ ., data = data_train,
      method = "knn", preProcess = c("center", "scale", "pca"), 
      tuneLength = 10, 
      trControl = trainControl(method="cv", number=10))



pred <- predict(model2, data_test)

caret::confusionMatrix(pred, data_test$categoria) 
## Confusion Matrix and Statistics
## 
##                   Reference
## Prediction         Bacia Torneira Interna
##   Bacia              175                8
##   Torneira Interna    25              192
##                                           
##                Accuracy : 0.9175          
##                  95% CI : (0.8861, 0.9425)
##     No Information Rate : 0.5             
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.835           
##                                           
##  Mcnemar's Test P-Value : 0.005349        
##                                           
##             Sensitivity : 0.8750          
##             Specificity : 0.9600          
##          Pos Pred Value : 0.9563          
##          Neg Pred Value : 0.8848          
##              Prevalence : 0.5000          
##          Detection Rate : 0.4375          
##    Detection Prevalence : 0.4575          
##       Balanced Accuracy : 0.9175          
##                                           
##        'Positive' Class : Bacia           
## 

Existe um pequeno ganho no PCA, o que pode ser visto é que é preciso uma técnica mais robusta de redução de dimensionalidade para que encontremos uma melhor separação dos dados e consequentemente um resultado melhor. Mais a frente será debatido o t-SNE e seus resultados para este banco de dados.

PCA para dados textuais:

Objetivo:

Um dos claros objetivos desta apresentação é ir além do PCA, porém, pede-se um pouco de paciência para que vejamos de forma breve (sem entrar nas questões de NLP) como esta técnica pode ser utilizada em análise textual.

Utilizaremos os seguintes pacotes para determinar a similaridade entre as palavras utilizados por Albert Camus no seu clássico: “O Estrangeiro”.

knitr::include_graphics("7248407GG.jpg")

“Hoje, a mãe morreu. Ou talvez ontem, não sei bem. Recebi um telegrama do asilo: “Sua mãe falecida: Enterro amanhã. Sentidos pêsames”. Isto não quer dizer nada. Talvez tenha sido ontem.”

Carregando o livro

library(epubr)
library(tm)
library(lexiconPT)
library(broom)
library(tidytext)
library(widyr)



x0 <- epubr::epub("O Estrangeiro - Albert Camus.epub")

estrangeiro <- x0$data[[1]]

Alguns ajustes

Neste momento, vamos separar as palavras utilizadas e remover as “stopwords” em português:

Analisando os “unigrams”:

stop_words <- stopwords(kind = "pt") %>% 
  as.tibble()

unigram_probs <- estrangeiro %>%
  tidytext::unnest_tokens(word, text) %>%
  count(word, sort = TRUE) %>%
  mutate(p = n / sum(n)) 

head(unigram_probs)
## # A tibble: 6 x 3
##   word      n      p
##   <chr> <int>  <dbl>
## 1 que    1227 0.0408
## 2 a      1112 0.0370
## 3 o      1053 0.0350
## 4 e       990 0.0329
## 5 de      905 0.0301
## 6 me      664 0.0221

Vamos escolher uma janela móvel de 15 palavras, pode-se aplicar maiores ou menores, é interessante testar algumas afim de entender qual é o valor que propociona os melhores resultados.

colnames(stop_words) <- c("word")

tidy_skipgrams <- estrangeiro %>%
  tidytext::unnest_tokens(ngram, text, token = "ngrams", n = 15) %>%
  mutate(ngramID = row_number()) %>% 
  unite(skipgramID, ngramID) %>%
  unnest_tokens(word, ngram) %>% 
  anti_join(stop_words, by = "word")

head(tidy_skipgrams)
## # A tibble: 6 x 5
##   section           nword nchar skipgramID word    
##   <chr>             <int> <int> <chr>      <chr>   
## 1 content0002.xhtml     2    14 1          <NA>    
## 2 content0003.xhtml  4167 24617 2          capítulo
## 3 content0003.xhtml  4167 24617 2          i       
## 4 content0003.xhtml  4167 24617 2          hoje    
## 5 content0003.xhtml  4167 24617 2          mãe     
## 6 content0003.xhtml  4167 24617 2          morreu

Calculando as probabilidades:

skipgram_probs <- tidy_skipgrams %>%
  widyr::pairwise_count(word, skipgramID, diag = TRUE, sort = TRUE) %>%
  dplyr::mutate(p = n / sum(n))

Normalizando a probabilidade

Vamos utilizar um indicador para visualizar quais palavras ocorreram juntas com mais frequência do que o esperado, tendo em base a frequência com que elas ocorreram sozinhas.

O quanto maior o resultado, mais estas palavras estão associadas e possuem boa probabilidade de ocorrem juntas, em relação a probabilidade de serem encontradas individualmente.

normalized_prob <- skipgram_probs %>%
  dplyr::filter(n > 20) %>%
  dplyr::rename(word1 = item1, word2 = item2) %>%
  dplyr::left_join(unigram_probs %>%
              select(word1 = word, p1 = p),
            by = "word1") %>%
  dplyr::left_join(unigram_probs %>%
              select(word2 = word, p2 = p),
            by = "word2") %>%
  dplyr::mutate(p_together = p / p1 / p2)


head(normalized_prob)
## # A tibble: 6 x 7
##   word1    word2        n        p      p1      p2 p_together
##   <chr>    <chr>    <dbl>    <dbl>   <dbl>   <dbl>      <dbl>
## 1 disse    disse     2844 0.00157  0.00642 0.00642       38.2
## 2 mim      mim       1514 0.000838 0.00343 0.00343       71.4
## 3 é        é         1448 0.000801 0.00349 0.00349       65.7
## 4 pouco    pouco     1307 0.000723 0.00323 0.00323       69.5
## 5 raimundo raimundo  1303 0.000721 0.00306 0.00306       77.0
## 6 então    então     1227 0.000679 0.00279 0.00279       87.0

Vamos observar quais palavras estão associadas a “praia”:

normalized_prob %>% 
  dplyr::filter(word1 == "praia") %>%
  dplyr::arrange(-p_together)
## # A tibble: 23 x 7
##    word1 word2         n         p      p1        p2 p_together
##    <chr> <chr>     <dbl>     <dbl>   <dbl>     <dbl>      <dbl>
##  1 praia praia       443 0.000245  0.00103 0.00103        231. 
##  2 praia saco         39 0.0000216 0.00103 0.0000998      210. 
##  3 praia brilho       22 0.0000122 0.00103 0.000233        50.7
##  4 praia autocarro    33 0.0000183 0.00103 0.000366        48.4
##  5 praia banho        22 0.0000122 0.00103 0.000266        44.4
##  6 praia longe        44 0.0000243 0.00103 0.000565        41.8
##  7 praia voltei       57 0.0000315 0.00103 0.000765        40.0
##  8 praia longo        24 0.0000133 0.00103 0.000333        38.7
##  9 praia mar          57 0.0000315 0.00103 0.000831        36.8
## 10 praia pequena      24 0.0000133 0.00103 0.000466        27.7
## # ... with 13 more rows

Vamos observar quais palavras estão associadas a “raimundo”:

normalized_prob %>% 
  dplyr::filter(word1 == "raimundo") %>%
  dplyr::arrange(-p_together)
## # A tibble: 108 x 7
##    word1    word2        n         p      p1        p2 p_together
##    <chr>    <chr>    <dbl>     <dbl>   <dbl>     <dbl>      <dbl>
##  1 raimundo raimundo  1303 0.000721  0.00306 0.00306         77.0
##  2 raimundo copo        23 0.0000127 0.00306 0.0000665       62.5
##  3 raimundo bêbedo      22 0.0000122 0.00306 0.0000665       59.8
##  4 raimundo vítima      21 0.0000116 0.00306 0.0000665       57.1
##  5 raimundo ouviu       26 0.0000144 0.00306 0.0000998       47.1
##  6 raimundo ouvimos     22 0.0000122 0.00306 0.0000998       39.9
##  7 raimundo revólver    36 0.0000199 0.00306 0.000200        32.6
##  8 raimundo insultá     24 0.0000133 0.00306 0.000133        32.6
##  9 raimundo navalha     35 0.0000194 0.00306 0.000200        31.7
## 10 raimundo bateu       23 0.0000127 0.00306 0.000133        31.3
## # ... with 98 more rows

Vamos transformar estes dados em uma matriz esparsa:

pmi_matrix <- normalized_prob %>%
    dplyr::mutate(pmi = log10(p_together)) %>%
    tidytext::cast_sparse(word1, word2, pmi)

head(pmi_matrix@x)
## [1] 1.5821239 0.6537674 0.7127086 0.4269799 0.8530843 0.8116916

A transformação de dados textuais em matriz, possuem muitos zeros, esta estrutura utilizada economiza tempo e memória.

Aplicando o PCA:

Neste caso, usaremos a abordagem do PCA afim de reduzir a dimensionalidade dos dados e buscar uma forma mais interessante de representar a similaridade entre as palavras :

library(irlba)

pmi_pca <- irlba::prcomp_irlba(pmi_matrix, n = 256)

word_vectors <- pmi_pca$x

rownames(word_vectors) <- rownames(pmi_matrix)

dim(word_vectors)
## [1] 1884  256

Encontrando similaridades:

Vamos utilizar uma função publicada pela Julia Silge para para varrer o vetor de word_vectors atrás das palavras de maiores similaridades:

search_synonyms <- function(word_vectors, selected_vector) {
  
  similarities <- word_vectors %*% selected_vector %>%
    tidy() %>%
    as_tibble() %>%
    rename(token = .rownames,
           similarity = unrowname.x.)
  
  similarities %>%
    arrange(-similarity)    
}

Testando a função:

praia <- search_synonyms(word_vectors, word_vectors["praia",])

Adicionando sentimentos:

lex <- lexiconPT::oplexicon_v3.0

colnames(lex)[1] <- "word"

unnested_words <- praia %>%
  inner_join(lex, by = c("token" = "word"))

Analisando os sentimentos referentes a palavra praia:

unnested_words %>% 
  top_n(40, abs(similarity)) %>% 
  ggplot(aes(x = reorder(token, similarity), y = similarity, fill = as.factor(polarity))) +
  geom_col() +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        axis.ticks.x = element_blank())

Podemos também realizar de outra forma o PCA:

Recuperemos nosso livro e ajustemos ele para análise:

stop_words <- stopwords(kind = "pt") %>% 
  as.tibble()

colnames(stop_words) <- c("word")

words <- map_df(1,
                ~ unnest_tokens(estrangeiro, word, text, 
                                token = "ngrams", n = .x)) %>%
  anti_join(stop_words, by = "word")

unnested_words <- words %>%
  count(word, sort = TRUE) %>% 
  filter(n > 5)

Filtremos apenas com as palavras encontradas no lexcon:

Criando matriz de tag e transformando em matriz esparsa

sparse_tag_matrix <- unnested_words %>%
  tidytext::cast_sparse(word,word,n)

#PCA -> PCA para matriz esparsa

word_scaled <- scale(sparse_tag_matrix)

set.seed(42)

tags_pca <- irlba::prcomp_irlba(word_scaled, n = 64)

Visualizando

tidied_pca <- bind_cols(Tag = colnames(sparse_tag_matrix),
                        tidy(tags_pca$rotation))

text_pca <- reshape2::melt(tidied_pca)

Melhorando a visualizacao

text_pca %>%
  filter(variable == "PC1") %>%
  top_n(40, abs(value)) %>%  
  mutate(Tag = reorder(Tag, value)) %>%
  ggplot(aes(Tag, value, fill = Tag)) +
  geom_col(show.legend = FALSE, alpha = 0.8) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        axis.ticks.x = element_blank()) +
  labs(x = "Words",
       y = "Relative importance in principle component",
       title = "PC1")

Olhando outros PCS

text_pca %>%
  filter(variable == "PC2") %>%
  top_n(40, abs(value)) %>%  
  mutate(Tag = reorder(Tag, value)) %>%
  ggplot(aes(Tag, value, fill = Tag)) +
  geom_col(show.legend = FALSE, alpha = 0.8) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        axis.ticks.x = element_blank()) +
  labs(x = "Words",
       y = "Relative importance in principle component",
       title = "PC2")

Olhando outros PCS

text_pca %>%
  filter(variable == "PC3") %>%
  top_n(40, abs(value)) %>%
  mutate(Tag = reorder(Tag, value)) %>%
  ggplot(aes(Tag, value, fill = Tag)) +
  geom_col(show.legend = FALSE, alpha = 0.8) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        axis.ticks.x = element_blank()) +
  labs(x = "Words",
       y = "Relative importance in principle component",
       title = "PC3")

Está difícil né? O livro não é muito linear e as situações muitas vezes se repetem e o autor tem sentimentos, por vezes, contraditórios.

# "Pontuação dos PCS"  -> sum(importance * lex)

colnames(text_pca)[1] <- "word"

value_pca <- text_pca %>%
  inner_join(lex) %>%
  group_by(variable) %>%
  summarise(val = sum(value*polarity)) 


value_pca
## # A tibble: 64 x 2
##    variable    val
##    <fct>     <dbl>
##  1 PC1       0.144
##  2 PC2      -0.351
##  3 PC3      -0.539
##  4 PC4       0.683
##  5 PC5      -0.312
##  6 PC6      -0.428
##  7 PC7       0.596
##  8 PC8       0.100
##  9 PC9       0.120
## 10 PC10      0.162
## # ... with 54 more rows

Vamos ver algumas palavras como se comportam nos PCS, para isto teremos que refazer o processo:

stop_words <- stopwords(kind = "pt") %>% 
  as.tibble()

colnames(stop_words) <- c("word")

words <- map_df(1,
                ~ unnest_tokens(estrangeiro, word, text, 
                                token = "ngrams", n = .x)) %>%
  anti_join(stop_words, by = "word")

unnested_words <- words %>%
  count(word, sort = TRUE) %>% 
  filter(n > 1)

### Realizando o pca

sparse_tag_matrix <- unnested_words %>%
  tidytext::cast_sparse(word,word,n)

#PCA -> PCA para matriz esparsa

word_scaled <- scale(sparse_tag_matrix)

set.seed(42)

tags_pca <- irlba::prcomp_irlba(word_scaled, n = 64)

tidied_pca <- bind_cols(Tag = colnames(sparse_tag_matrix),
                        tidy(tags_pca$rotation))

text_pca <- reshape2::melt(tidied_pca)

Vamos ver como a palavra “praia” e “enterro” se comportam nos PCs:

text_pca %>% 
  filter(Tag %in% c("enterro", "praia", "mãe", "maria", "árabe", "juiz")) %>% 
  filter(variable %in% c("PC1","PC2", "PC3", "PC4", "PC5")) %>% 
  ggplot(aes(x = variable, y = value, fill = Tag)) +
  geom_col() + 
  facet_wrap(~Tag, ncol = 3) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        axis.ticks.x = element_blank())