Análise Multivariada - Atividade 01

1 Introdução

library(corrplot)
library(ggcorrplot)
library(ggplot2)
library(GGally)
library(dplyr)
library(factoextra)
library(viridis)
library(knitr)
library(pacman)
library(rstatix)
library(e1071)
library(RColorBrewer)
library(ggmap)
library(maps)
library(ggrepel)
library(MASS)
library(magrittr)
library(ggpubr)
library(Rtsne)
library(bigmds)
library(tictoc)
pacman::p_load(knitr, captioner, bundesligR, stringr)

A concentração média anual de dióxido de enxofre (\(SO_2\)), em microgramas por metro cúbico, é uma medida da poluição do ar em cidades. A questão de interesse é:

Quais ou como os aspectos do clima e da ecologia humana medidos por outras variáveis influenciam a poluição do ar?

A fim de responder esse questionamento, serão realizadas diversas análises no dataset usairpollution.

2 Análise Descritiva

A Análise Descritiva é a fase inicial do processo de estudo dos dados coletados, em que são utilizados métodos de Estatística Descritiva para organizar, resumir, descrever e comparar aspectos importantes de um conjunto de variáveis. A identificação de anomalias e e dados dispersos também faz parte desse tipo de análise.

2.1 Características básicas

table_nums <- captioner::captioner(prefix = "Tabela")
tab.1_cap <- table_nums(name = "tab_1", 
                         caption = "Características básicas do dataset usairpollution")
tab.2_cap <- table_nums(name = "tab_2", 
                         caption = "Tipos e descrições das variáveis do dataset usairpollution")
tab.3_cap <- table_nums(name = "tab_3", 
                         caption = "Amostra de 6 observações de usairpollution")
tab.4_cap <- table_nums(name = "tab_4", 
                         caption = "Estatísticas resumo para as variáveis quantitativas de usairpollution")
tab.5_cap <- table_nums(name = "tab_5", 
                         caption = "Coeficiente de assimetria para as variáveis quantitativas de usairpollution")

tab.6_cap <- table_nums(name = "tab_6", 
                         caption = "Outliers para a variável SO2")
tab.7_cap <- table_nums(name = "tab_7", 
                         caption = "Outliers para a variável temp")
tab.8_cap <- table_nums(name = "tab_8", 
                         caption = "Outliers para a variável manu")
tab.9_cap <- table_nums(name = "tab_9", 
                         caption = "Outliers para a variável popul")
tab.10_cap <- table_nums(name = "tab_10", 
                         caption = "Outliers para a variável precip")
tab.11_cap <- table_nums(name = "tab_11", 
                         caption = "Outliers para a variável predays")
tab.12_cap <- table_nums(name = "tab_12", 
                         caption = "Vetores de carga")
tab.13_cap <- table_nums(name = "tab_13", 
                         caption = "Resumo PCA")

A Tabela 1 contém algumas características básicas do dataset, quais sejam:

• Existência de NAs;
• Total de observações;
• Total de variáveis;
• Existência de linhas duplicadas.

De acordo com a Tabela 1:

• O dataset possui 41 observações e 8 variáveis;
• Não possui NAs e nem linhas duplicadas.

Tabela 1: Características básicas do dataset usairpollution

usairpollution <- read.csv2("E:/Multivariada/Atividade 01/usairpollution.csv", header= TRUE)
basic1_charact <- data.frame(NAs = sum(is.na(usairpollution)), Observações = nrow(usairpollution), Variáveis = ncol(usairpollution), Duplicidades = anyDuplicated.data.frame(usairpollution))
colnames(usairpollution)[1] <- "City"                            

knitr::kable(basic1_charact)

NAs	Observações	Variáveis	Duplicidades
0	41	8	0

Em seguida, a Tabela 2 contém os tipos e descrições das 8 variáveis de usairpollution.

De acordo com a Tabela 2, o dataset possui:

• 4 variáveis quantitativas discretas (SO2,manu,popul e predays);
• 3 variáveis quantitativas contínuas (temp, wind, precip);
• 1 variável qualitativa nominal (City).

Tabela 2: Tipos e descrições das variáveis do dataset usairpollution

tipos <- NULL
for(i in 1:ncol(usairpollution)){
 tipos<- rbind(tipos,(typeof(usairpollution[,i])))
}

basic2_charact <- data.frame(Variável = colnames(usairpollution),
                            Tipo = tipos, 
                            Descrição = c("Cidade dos E.U.A","Teor de $SO_2$ (dióxido de enxofre) do ar em microgramas por metro cúbico","Temperatura média anual em Fahrenheit","Nº de empresas manufatureiras que empregam 20 ou mais trabalhadores","Tamanho da população em milhares (Censo de 1970)", "Velocidade média anual do vento em milhas por hora", "Precipitação média anual em polegadas", "Nº médio de dias com precipitação por ano"))

knitr::kable(basic2_charact)

Variável	Tipo	Descrição
City	character	Cidade dos E.U.A
SO2	integer	Teor de \(SO_2\) (dióxido de enxofre) do ar em microgramas por metro cúbico
temp	double	Temperatura média anual em Fahrenheit
manu	integer	Nº de empresas manufatureiras que empregam 20 ou mais trabalhadores
popul	integer	Tamanho da população em milhares (Censo de 1970)
wind	double	Velocidade média anual do vento em milhas por hora
precip	double	Precipitação média anual em polegadas
predays	integer	Nº médio de dias com precipitação por ano

Por fim, ao retirarmos uma pequena amostra do dataset, podemos observar que cada linha traz informação sobre uma cidade.

Tabela 3: Amostra de 6 observações de usairpollution

knitr::kable(head(usairpollution))

City	SO2	temp	manu	popul	wind	precip	predays
Albany	46	47.6	44	116	8.8	33.36	135
Albuquerque	11	56.8	46	244	8.9	7.77	58
Atlanta	24	61.5	368	497	9.1	48.34	115
Baltimore	47	55.0	625	905	9.6	41.31	111
Buffalo	11	47.1	391	463	12.4	36.11	166
Charleston	31	55.2	35	71	6.5	40.75	148

2.2 Estatísticas resumo

A Tabela 4 contém algumas estatísticas resumo para as variáveis quantitativas do dataset usairpollution.

De acordo com a Tabela 4, para a variável de interesse teor de dióxido de enxofre em microgramas por metro cúbico (\(SO_2/m^3\)), temos que:

• Os teores mínimo e máximo observados são 8µg/m³ e 110µg/m³;
• De acordo com o primeiro quartil, 25% das cidades do conjunto de dados possuem teor \(SO_2\) igual ou inferior a 13µg/m³;
• De acordo com a mediana, 50% das cidades do conjunto de dados possuem teor de \(SO_2\) igual ou acima de 26µg/m³;
• De acordo com o terceiro quartil, 25% das cidades do conjunto de dados possuem teor \(SO_2\) igual ou acima de 35µg/m³;
• De acordo com o intervalo interquartil, 50% das cidades possui teor de \(SO_2\) entre 13µg/m³ e 35µg/m³.
• O teor médio de \(SO_2/m^3\) é 30.049µg/m³;
• O desvio padrão do teor de \(SO_2/m^3\) é 23.472µg/m³.

Tabela 4: Estatísticas resumo para as variáveis quantitativas de usairpollution

knitr::kable(get_summary_stats(usairpollution[,-1])[,-c(2,8,9,12,13)])

variable	min	max	median	q1	q3	mean	sd
manu	35.00	3344.0	347.00	181.00	462.00	463.098	563.474
popul	71.00	3369.0	515.00	299.00	717.00	608.610	579.113
precip	7.05	59.8	38.74	30.96	43.11	36.769	11.772
predays	36.00	166.0	115.00	103.00	128.00	113.902	26.506
SO2	8.00	110.0	26.00	13.00	35.00	30.049	23.472
temp	43.50	75.5	54.60	50.60	59.30	55.763	7.228
wind	6.00	12.7	9.30	8.70	10.60	9.444	1.429

2.3 Assimetria e Densidades

A Tabela 5 contém os coeficientes de assimetria para as variáveis quantitativas do dataset.

A assimetria negativa indica que a média dos dados é menor que a mediana e, portanto, que a distribuição dos dados é assimétrica à esquerda. Já a assimetria positiva indica que a média dos dados é maior que a mediana e, portanto, que a distribuição dos dados é assimétrica à direita.

Tabela 5: Coeficiente de assimetria para as variáveis quantitativas de usairpollution

knitr::kable(apply(usairpollution[,-1],2,skewness),col.names = "Coeficiente de assimetria")

	Coeficiente de assimetria
SO2	1.5841126
temp	0.8229757
manu	3.4846033
popul	2.9412580
wind	0.0026751
precip	-0.6925181
predays	-0.5500923

De acordo com a Tabela 5, temos que:

• A distribuição de SO2 é assimétrica à direita;
• A distribuição de temp é assimétrica à direita;
• A distribuição de manu é assimétrica à direita;
• A distribuição de popul é assimétrica à direita;
• A distribuição de wind é assimétrica à direita;
• A distribuição de precip é assimétrica à esquerda;
• A distribuição de predays é assimétrica à esquerda.

Os gráficos abaixo são relativos às densidades dessas variáveis. Através deles é possível visualizar a informação contida nos coeficientes de assimetria.

par(mfrow=c(4,2))

ggplot(usairpollution, aes(x=SO2)) + 
  geom_density()+
scale_x_continuous(breaks = seq(0,120,by=10), limits= c(0,120))+
labs(title = "Densidade da variável SO2 \n",
      x = "\nTeor de SO2 em µg/m³ \n",
      y = "Densidade\n")+
  theme_classic()

ggplot(usairpollution, aes(x=temp)) + 
  geom_density()+
scale_x_continuous(breaks = seq(0,120,by=10), limits= c(0,120))+
labs(title = "Densidade da variável temp \n",
     x = "\nTemperatura média anual (°F)",
      y = "Densidade\n")+
  theme_classic()

ggplot(usairpollution, aes(x=manu)) + 
  geom_density()+
scale_x_continuous(breaks = seq(0,3400,by=200), limits= c(0,3400))+
labs(title = "Densidade da variável manu \n",
     x = "\n Nº de empresas manufatureiras que empregam 20 ou mais trabalhadores",
      y = "Densidade\n")+
  theme_classic()
  
ggplot(usairpollution, aes(x=popul)) + 
  geom_density()+
scale_x_continuous(breaks = seq(0,3400,by=200), limits= c(0,3400))+
labs(title = "Densidade da variável popul \n",
     x = "\n Tamanho da população (Censo de 1970) em milhares",
      y = "Densidade\n")+
  theme_classic()
    
ggplot(usairpollution, aes(x=wind)) + 
  geom_density()+
scale_x_continuous(breaks = seq(0,120,by=10), limits= c(0,120))+
labs(title = "Densidade da variável wind \n",
     x = "\n Velocidade média anual do vento em milhas por hora",
      y = "Densidade\n")+
  theme_classic()


ggplot(usairpollution, aes(x=precip)) + 
  geom_density()+
scale_x_continuous(breaks = seq(0,120,by=10), limits= c(0,120))+
labs(title = "Densidade da variável precip \n",
     x = "\n Precipitação média anual em polegadas",
      y = "Densidade\n")+
  theme_classic()

ggplot(usairpollution, aes(x=predays)) + 
  geom_density()+
scale_x_continuous(breaks = seq(0,200,by=10), limits= c(0,200))+
labs(title = "Densidade da variável predays \n",
     x = "\n Nº médio de dias com precipitação por ano",
      y = "Densidade\n")+
  theme_classic()

2.4 Boxplots

Os gráficos abaixo são relativos aos boxplots das variáveis quantitativas do dataset usairpollution. Através deles é possível visualizar a distribuição dos dados, assimetria, quartis e outliers.

par(mfrow=c(4,2))
cores <- brewer.pal(n = 7, name = "Set1")


a<-ggplot(data = usairpollution, aes (y = SO2))+
geom_boxplot(fill=cores[1], color="black") +
labs(title = "Distribuição da variável SO2 \n",
       y = " Teor de SO2 em µg/m³ \n")+
scale_y_continuous(breaks = seq(0,130,by=10), limits= c(0,130))+
scale_x_continuous(labels=NULL)+
theme_classic()
a

b<- ggplot(data = usairpollution, aes (y = temp))+
geom_boxplot(fill=cores[2], color="black") +
labs(title = "Distribuição da variável temp \n",
       y = "Temperatura média anual (°F) \n")+
scale_y_continuous(breaks = seq(0,130,by=10), limits= c(0,130))+
scale_x_continuous(labels=NULL)+
theme_classic()
b
c<- ggplot(data = usairpollution, aes (y = manu))+
geom_boxplot(fill=cores[3], color="black") +
labs(title = "Distribuição da variável manu \n",
       y = " Nº de empresas manufatureiras que empregam 20 ou mais trabalhadores \n")+
scale_y_continuous(breaks = seq(0,3400,by=200), limits= c(0,3400))+
scale_x_continuous(labels=NULL)+
theme_classic()
c

d<- ggplot(data = usairpollution, aes (y = popul))+
geom_boxplot(fill=cores[4], color="black") +
labs(title = "Distribuição da variável popul \n",
       y = " Tamanho da população (Censo de 1970) em milhares \n")+
scale_y_continuous(breaks = seq(0,3400,by=200), limits= c(0,3400))+
scale_x_continuous(labels=NULL)+
theme_classic()
d

e<-ggplot(data = usairpollution, aes (y = wind))+
geom_boxplot(fill=cores[5], color="black") +
labs(title = "Distribuição da variável wind \n",
       y = " Velocidade média anual do vento em milhas por hora \n")+
scale_y_continuous(breaks = seq(0,50,by=2), limits= c(0,50))+
scale_x_continuous(labels=NULL)+
theme_classic()
e


f<-ggplot(data = usairpollution, aes (y = precip))+
geom_boxplot(fill=cores[6], color="black") +
labs(title = "Distribuição da variável precip \n",
       y = " Precipitação média anual em polegadas \n")+
scale_y_continuous(breaks = seq(0,130,by=10), limits= c(0,130))+
scale_x_continuous(labels=NULL)+
theme_classic()
f

g<-ggplot(data = usairpollution, aes (y = predays))+
geom_boxplot(fill=cores[7], color="black") +
labs(title = "Distribuição da variável predays \n",
       y = " Nº médio de dias com precipitação por ano \n")+
scale_y_continuous(breaks = seq(0,180,by=10), limits= c(0,180))+
scale_x_continuous(labels=NULL)+
theme_classic()
g

2.5 Outliers

Para a variável SO2, os outliers correspondem às medições de dióxido de enxofre nas cidades de Chicago, Philadelphia e Providence.

Tabela 6: Outliers para a variável SO2

out<- ggplot_build(a)[["data"]][[1]][["outliers"]]
knitr::kable(usairpollution[which(usairpollution$SO2==out[[1]][1]| usairpollution$SO2==out[[1]][2] |usairpollution$SO2==out[[1]][3]),c(1,2)],row.names=FALSE)

City	SO2
Chicago	110
Philadelphia	69
Providence	94

Para a variável temp, o outlier corresponde a temperatura média anual(°F) na cidade de Miami.

Tabela 7: Outliers para a variável temp

out<- ggplot_build(b)[["data"]][[1]][["outliers"]]
knitr::kable(usairpollution[which(usairpollution$temp==out[[1]][1]),c(1,3)],row.names=FALSE)

City	temp
Miami	75.5

Para a variável manu, os outliers correspondem ao nº de empresas manufatureiras que empregam 20 ou mais trabalhadores nas cidades de Chicago, Cleveland, Detroit e Philadelphia.

Tabela 8: Outliers para a variável manu

out<- ggplot_build(c)[["data"]][[1]][["outliers"]]
knitr::kable(usairpollution[which(usairpollution$manu==out[[1]][1]|usairpollution$manu==out[[1]][2]|usairpollution$manu==out[[1]][3]|usairpollution$manu==out[[1]][4]),c(1,4)],row.names=FALSE)

City	manu
Chicago	3344
Cleveland	1007
Detroit	1064
Philadelphia	1692

Para a variável popul, os outliers correspondem ao tamanho da população (milhare) nas cidades de Chicago, Detroit e Philadelphia.

Tabela 9: Outliers para a variável popul

out<- ggplot_build(d)[["data"]][[1]][["outliers"]]
knitr::kable(usairpollution[which(usairpollution$popul==out[[1]][1]|usairpollution$popul==out[[1]][2]|usairpollution$popul==out[[1]][3]),c(1,5)],row.names=FALSE)

City	popul
Chicago	3369
Detroit	1513
Philadelphia	1950

A variável wind não possui outliers.

Para a variável precip, os outliers correspondem a precipitação média anual em polegadas nas cidades de Albuquerque e Phoenix.

Tabela 10: Outliers para a variável precip

out<- ggplot_build(f)[["data"]][[1]][["outliers"]]
knitr::kable(usairpollution[which(usairpollution$precip==out[[1]][1]|usairpollution$precip==out[[1]][2]),c(1,7)],row.names=FALSE)

City	precip
Albuquerque	7.77
Phoenix	7.05

Por fim, para a variável predays os outliers correspondem ao número médio de dias com precipitação por ano nas cidades de Albuquerque, Phoenix e San Francisco.

Tabela 11: Outliers para a variável predays

out<- ggplot_build(g)[["data"]][[1]][["outliers"]]
knitr::kable(usairpollution[which(usairpollution$predays==out[[1]][1]|usairpollution$predays==out[[1]][2] | usairpollution$predays==out[[1]][3]),c(1,8)],row.names=FALSE)

City	predays
Albuquerque	58
Buffalo	166
Phoenix	36

2.6 Gráficos de dispersão

Foram feitos gráficos de dispersão entre a variável SO2 e as demais a fim de verificar possíveis relações, especialmente as lineares. Entretanto, ao observar os gráficos não foi possível identificar visualmente quaisquer relações.

par(mfrow=c(3,2))
cores2 <- brewer.pal(n = 6, name = "Dark2")
ggplot(usairpollution, aes(x=temp, y=SO2)) +
geom_point(color=cores2[1])+
labs(title = "Gráfico de dispersão entre as variáveis SO2 e temp",
      y = " Teor de SO2 em µg/m³ \n",
      x = "\nTemperatura média anual (°F)")+
theme_classic()

ggplot(usairpollution, aes(x=manu, y=SO2)) +
geom_point(color=cores2[2])+
labs(title = "Gráfico de dispersão entre as variáveis SO2 e manu",
      y = " Teor de SO2 em µg/m³ \n",
      x = "\n Nº de empresas manufatureiras que empregam 20 ou mais trabalhadores")+
theme_classic()

ggplot(usairpollution, aes(x=popul, y=SO2)) +
geom_point(color=cores2[3])+
labs(title = "Gráfico de dispersão entre as variáveis SO2 e popul",
      y = " Teor de SO2 em µg/m³ \n",
      x = "\n Tamanho da população (Censo de 1970) em milhares")+
theme_classic()

ggplot(usairpollution, aes(x=wind, y=SO2)) +
geom_point(color=cores2[4])+
labs(title = "Gráfico de dispersão entre as variáveis SO2 e wind",
      y = " Teor de SO2 em µg/m³ \n",
      x = "\n Velocidade média anual do vento em milhas por hora")+
theme_classic()

ggplot(usairpollution, aes(x=precip, y=SO2)) +
geom_point(color=cores2[5])+
labs(title = "Gráfico de dispersão entre as variáveis SO2 e precip",
      y = " Teor de SO2 em µg/m³ \n",
      x = "\n Precipitação média anual em polegadas")+
theme_classic()

ggplot(usairpollution, aes(x=predays, y=SO2)) +
geom_point(color=cores2[6])+
labs(title = "Gráfico de dispersão entre as variáveis SO2 e predays",
      y = " Teor de SO2 em µg/m³ \n",
      x = "\n Nº médio de dias com precipitação por ano")+
theme_classic()

2.7 Correlograma

O correlograma abaixo mostra os coeficientes de correlação de Pearson para todas as variáveis quantitativas do dataset.

corr <- cor(usairpollution[,-1]) #Matriz de correlação
#p.mat <- cor_pmat(usairpollution[,-1]) # Matriz de p-valores das correlações
ggcorrplot(corr,outline.color = "white",type = "lower",lab = TRUE)

• A variável de interesse SO2 possui maior correlação com a variável manu, 0.64.

• As variáveis manu e popul possuem correlação 0.96, o que significa que essas variáveis possuem uma relação linear positiva: quanto maior a população, maior o número de empresas manufatureiras.

2.8 Distribuição geográfica da variável SO2

Para produzir o mapa de calor para a variável teor de dióxido de enxofre (SO2), foi utilizado um dataset que contém populações de cidades dos E.U.A de 1790 a 2010 e suas respectivas latitudes e longitudes, aqui nomeado popul_1790.2010.

O primeiro passo da construção do mapa de calor consistiu em fazer um merge entre os datasets popul_1790.2010 e usairpollution pela variável chave City, a fim de obter as latitudes e longitudes. O resultado foi um dataset com 130 observações, visto que existem cidades de mesmo nome mas que pertencem a estados diferentes.

Com o objetivo de eliminar múltiplas entradas para a mesma cidade, foi feita uma comparação entre variáveis popul e X1970, ambas referentes à população dessas cidades de acordo com o Censo de 1970. Foram mantidas entradas únicas para cada cidade, de forma que o quadrado da diferença entre popul e X1970 fosse mínima.

popul_1790.2010 <- read.csv("E:/Multivariada/Atividade 01/1790-2010_MASTER.txt")
x <- merge(usairpollution,popul_1790.2010,by="City")
x$popul <- x$popul*1000

cities <- as.character(unique(x$City))
y <- NULL
for(i in 1:length(cities)){
  
  keep <- which(x$City==cities[i])[which.min((x$popul[which(x$City==cities[i])] - x$X1970[which(x$City==cities[i])])^2)]
  y <- rbind(y,x[keep,])
}

Com os dados das latitudes e longitudes em mãos foi possível construir o mapa abaixo, que contém os teores de \(SO_2\) para as 41 cidades do dataset usairpollution. Cada ponto representa uma cidade e, quanto maior o tamanho do ponto ou mais próximo de amarelo a sua cor, maior o teor de \(SO_2\) observado. As cidades com as maiores medições de \(SO_2\) são Chicago, Providence, Philadelphia, Cleveland e Pittsburgh.

ggplot(y, aes(LON, LAT,label=City)) + borders("state") + 
  geom_point(aes(size = SO2, fill=SO2),colour="black", shape=21) + 
  coord_quickmap() +
  geom_label_repel(aes(label=ifelse(SO2>60,as.character(City),'')), 
                   size = 2,
                   box.padding   = 0.35, 
                   point.padding = 0.5,
                   segment.color = 'grey50')+
  labs(title = "Distribuição geográfica da variável SO2\n",
      y = " Latitude\n",
      x = "\nLongitude")+
  theme(plot.title = element_text(hjust = 0.5))+
  scale_fill_viridis(option = "C")+
  theme_classic()

3 Análise de Componentes Principais (PCA)

As componentes principais são novas variáveis que são construídas como combinações lineares das variáveis iniciais.

Essas combinações são feitas de forma que essas novas variáveis não sejam correlacionadas e a maioria das informações contidas nas variáveis iniciais seja armazenada nos primeiros componentes.

Então, a ideia é que dados k-dimensionais forneçam k componentes principais. O método busca colocar o máximo de informações possíveis nas primeiras componentes, para que, se quisermos reduzir a dimensionalidade de um conjunto de dados, possamos focar a análise nos primeiros componentes sem sofrer uma grande penalidade em termos de perda de informação.

As variáveis do conjunto de dados devem ser padronizadas previamente à análise. No caso desse estudo, a variável de interesse, SO2, foi retirada da análise de componentes principais, pois não é nosso interesse agregá-la às demais.

3.1 Vetores de Carga

Os vetores de carga contêm os coeficientes da combinação linear que gerou cada componente principal para as variáveis utilizadas. A primeira componente, por exemplo, pode ser descrita como

\(0.32temp-0.61manu-0.57popul-0.35wind+0.04precip-0.24predays\)

Tabela 12: Vetores de carga

res.pca <- prcomp(usairpollution[,-c(1,2)], scale = TRUE, center = TRUE)
knitr::kable(res.pca$rotation)

	PC1	PC2	PC3	PC4	PC5	PC6
temp	0.3296461	-0.1275974	0.6716861	-0.3064573	-0.5580564	0.1361878
manu	-0.6115424	-0.1680577	0.2728863	0.1368408	0.1020421	0.7029705
popul	-0.5778220	-0.2224533	0.3503741	0.0724813	-0.0780655	-0.6946413
wind	-0.3538388	0.1307915	-0.2972533	-0.8694258	-0.1132669	0.0245250
precip	0.0408070	0.6228578	0.5045629	-0.1711483	0.5681834	-0.0606222
predays	-0.2379159	0.7077653	-0.0930885	0.3113069	-0.5800039	0.0219606

3.2 Scree Plot

Podemos observar que a primeira componente explica cerca de 36.6% da variância total dos dados, a segunda 25% e a terceira, 23.2%.

fviz_eig(res.pca,addlabels = TRUE, xlab = "Dimensões", ylab = "% da variância explicada", main = "Scree Plot" )

Utilizando as 4 primeiras componentes principais, a proporção da variância explicada (PVE) acumulada é de 97.52%.

Tabela 13: Resumo PCA

summary(res.pca)

## Importance of components:
##                          PC1   PC2    PC3    PC4    PC5     PC6
## Standard deviation     1.482 1.225 1.1810 0.8719 0.3385 0.18560
## Proportion of Variance 0.366 0.250 0.2324 0.1267 0.0191 0.00574
## Cumulative Proportion  0.366 0.616 0.8485 0.9752 0.9943 1.00000

3.3 Qualidade da representação das variáveis (cos2)

A qualidade da representação de uma variável é medida pelo quadrado das cargas ou coordenadas, chamada de cos2.

O gráfico abaixo mostra o total de cos2 para as duas primeiras componentes, ou seja, o valor para cada variável corresponde a soma de cos2 para as duas primeiras componentes.

fviz_pca_var(res.pca, col.var = "cos2",
gradient.cols = viridis_pal()(12),
repel = TRUE # Avoid text overlapping
) + labs(title =" Qualidade da representação (cos2) - Variáveis")

As variáveis mais próximas da circunferência são predays, precip, manu e popul, ou seja sua representação no mapa de fatores é boa, então são mais importantes para as primeiras componentes. Já as variáveis wind e temp estão mais próximas da origem, o que significa não foram bem representadas pelas componentes 1 e 2. Podemos observar que as variáveis popul, manu, wind e predays são variáveis positivamente correlacionadas, pois estão do mesmo lado da circunferência.

3.4 Contribuição das variáveis nas componentes (contrib)

As variáveis que mais contribuem na construção da primeira e segunda componentes (Dim.1 e Dim.2) são as mais importantes para explicar a variabilidade no conjunto de dados. Essas variáveis são:

• manu
• popul
• precip
• predays

var <- get_pca_var(res.pca)

corrplot(var$contrib, is.corr=FALSE,tl.col = 'black',col.lim = c(0, 80))

O gráfico abaixo contém as contribuições das variáveis na primeira componente (%). A linha tracejada em vermelho representa a contribuição média esperada (%).

Se todas as variáveis contribuíssem de forma igual para a construção das componentes, o valor da contribuição média esperada (%) seria \(\frac{1}{6} \approx 16.67\%\).

Portanto, se uma determinada variável tem contribuição acima de \(16.67\%\), esta pode ser considerada importante na contribuição da componente.

Na primeira componente, as variáveis manu e popul ultrapassam a contribuição média esperada (%), logo, são importantes em sua construção.

fviz_contrib(res.pca, choice = "var", axes = 1) + labs(title ="Contribuição das variáveis na Primeira Componente", y = "Contribuição em %")

Já na componente 2, as variáveis predays e precip são importantes para sua construção.

fviz_contrib(res.pca, choice = "var", axes = 2) + labs(title ="Contribuição das variáveis na Segunda Componente", y = "Contribuição em %")

Se considerarmos as componentes 1 e 2 simultaneamente, as variáveis predays, manu e popul são as mais importantes.

fviz_contrib(res.pca, choice = "var", axes = 1:2) + labs(title ="Contribuição das variáveis na Primeira e Segunda Componentes", y = "Contribuição em %")

3.5 Qualidade da representação dos indivíduos (cos2)

O gráfico abaixo contém a qualidade da representação(cos2) dos indivíduos. Os individuos similares ficam mais próximos entre si. Quanto mais distante do centro estiver o indivíduo, melhor sua representação nas componentes 1 e 2: destacam-se as cidades de Philadelphia, Chicago, Phoenix e Albuquerque.

us <- data.frame(usairpollution, row.names = 1)
res.pca <- prcomp(us[,-2], scale = TRUE, center = TRUE)
ind = get_pca_ind(res.pca)
fviz_pca_ind(res.pca, repel = TRUE,
pointsize = "cos2",
pointshape = 21,
fill = "yellow"
) + theme_minimal() + labs(title =" Qualidade da representação (cos2) - Indivíduos")

O gráfico abaixo corrobora com o que foi observado anteriormente: as cidades melhor representadas nas componentes 1 e 2 são Chicago, San Francisco, Albuquerque, Philadelphia, Denver e Phoenix.

fviz_cos2(res.pca, choice = "ind", axes=1:2) + labs(title =" Qualidade da representação (cos2) - Indivíduos - Componentes 1 e 2", y = "Cos2 - Qualidade da representação")

3.6 Contribuição dos índividuos nas componentes (contrib)

O gráfico abaixo contém as contribuições dos indivíduos nas duas primeiras componente (%). A linha tracejada em vermelho representa a contribuição média esperada (%).

fviz_contrib(res.pca, choice = "ind", axes = 1:2) +labs(title =" Contribuição dos Indivíduos nas componentes 1 e 2", y = "Contrib - Contribuição (%)")

Se todas os indivíduos contribuíssem de forma igual para a construção das componentes, o valor da contribuição média esperada (%) seria \(\frac{1}{41} \approx 2.44\%\).

Portanto, se um determinad indivíduo tem contribuição acima de \(2.44\%\), este pode ser considerada importante na contribuição da componente.

Para as duas primeiras componente, os indivíduos Chicago, Phoenix, Philadelphia, Albuquerque, San Francisco, Denver, Cleveland, e Salt Lake City ultrapassam a contribuição média esperada (%), logo, são importantes em sua construção.

4 Escalonamento Multidimensional (MDS)

O MDS retorna uma solução ótima para representar os dados em um espaço de dimensão inferior, onde o número de dimensões k é pré-especificado pela pessoa fazendo a análise. Por exemplo, a escolha de \(k = 2\) otimiza a localização dos objetos para um gráfico bidimensional de dispersão.

O algoritmo MDS toma como dado de entrada a matriz de dissimilaridade, representando as distâncias entre pares de objetos. A ideia é encontrar uma representação dos dados em menor dimensão e que preserve ao máximo as distâncias entre pontos.

O MDS será aplicado nas variáveis quantitativas do dataset usairpollution, com exceção da variável SO2. Previamente a utilização do algoritmo, os dados foram padronizados para a métrica [0,1]: \(X' = \frac{x-\min(x)}{\max(x)-min(x)}\). Foi utilizada a distância euclidiana para cálculo de dissimilaridades.

4.1 Escalonamento Multidimensional Clássico

O MDS Clássico preserva a métrica de distância original, entre pontos, da melhor forma possível. Isto é, as distâncias ajustadas no mapa MDS e as distâncias originais estão na mesma métrica. O MDS Clássico pertence à chamada categoria de escalonamento multidimensional métrico. É também conhecido como análise de coordenadas principais e é adequado para dados quantitativos.

Para demonstrar que é possível reconstruir a matriz de distância original utilizando as coordenadas geradas a partir do MDS, calcula-se a diferença máxima absoluta entre a matriz de distâncias original e a reconstruída.

O gráfico abaixo mostra as diferenças máximas observadas entre a matriz de distância original e a reconstruída a a partir das coordenadas MDS ao variar o valor de k entre 2 e 7.

As diferenças próximas de zero para \(k=\{5,6,7\}\) expressam que, a partir de \(k=5\) coordenadas, a matriz de distância original pode ser recuperada quase que totalmente.

P<- data.frame(usairpollution[,-2], row.names = 1)
P <- apply(P,2, function(x) (x-min(x))/(max(x)-min(x)))

valor <- NULL
for(i in 2:7){
  valor<-rbind(valor,c(i,max(abs(dist(P) - dist(cmdscale(dist(P),k=i, eig=TRUE)$points)))))
}

df<-data.frame(k = valor[,1], `dif_max` = valor[,2])

ggplot(data=df, aes(x=k, y=dif_max, group=1)) +
  geom_line()+
  geom_point()+
  theme_bw()+
  labs(title="Diferenças máximas absolutas entre matriz de distância original e matriz \nreconstruída a partir das coordenadas MDS, para k entre 2 e 7\n", y="Diferença máxima absoluta\n",x="\nk")+
  ylim(0,1.25)

Aplicando a função cmdscale (stats), com \(k=5\), abaixo temos o gráfico de dispersão das observações para as coordenadas 1 e 2. Pode ser observado que as cidades com teores de \(SO_2\) mais similares ficam mais agrupadas, enquanto que as que possuem teores mais discrepantes ficam mais afastadas.

mds <- P %>%
  dist() %>%          
  cmdscale(k=5, eig=TRUE) 


mds_label<- as.data.frame(mds$points)
colnames(mds_label)<- paste("Coord",1:5)
mds_label$City <- usairpollution[,1]
mds_label$SO2 <- usairpollution[,2]

ggplot(mds_label, aes(x=`Coord 1`, y = `Coord 2`)) +
  geom_point(aes(size=SO2), colour = "#32CD32")+
  labs(title= "Gráfico de dispersão - MDS Clássico com k=5", x= "Coord 1", y="Coord 2")+
  theme(plot.title = element_text(hjust = 0.5))+
  geom_label_repel(aes(label=City), 
                   size = 2,
                   box.padding   = 0.35, 
                   point.padding = 0.5,
                   segment.color = 'grey50')+
  theme_bw()

5 t-Distributed Stochastic Neighbor Embedding (t-SNE)

É uma técnica de visualização em 1D, 2D ou 3D de datasets de altas dimensões. No gráfico abaixo, podemos perceber que algumas cidades com valor de SO2 baixos ficaram concentradas na parte superior do gráfico. Comparando os resultados gráficos entre MDS e t-SNE, nota-se que as observações estão mais uniformemente distribuídas no t-SNE, enquanto que no MDS estão mais agrupadas.

tic("tsne")
toc(log = TRUE, quiet = TRUE)

tsne_fit <- Rtsne(P, perplexity=13)
tsne_coord <- data.frame(tsne_fit$Y)
tsne_coord$City <- usairpollution[,1]
tsne_coord$SO2 <- usairpollution[,2]

ggplot(data=tsne_coord, aes(x=X1, y=X2))+ 
  geom_point(aes(size=SO2), colour = "#32CD32") + 
  geom_label_repel(aes(label=City), 
                   size = 2,
                   box.padding   = 0.35, 
                   point.padding = 0.5,
                   segment.color = 'grey50')+ 
  labs(x="Coord 1\n",y="\nCoord 2", title="t-SNE para o dataset usairpollution") +
  theme_bw() 

6 Conclusão

O dataset usairpollution contém 8 variáveis e 41 observações relativas a cidades estadunidenses. A variável de interesse é SO2, o teor de dióxido de enxofre em microgramas por metro cúbico, uma medida de poluição.

Na Análise Descritiva, vimos estatísticas resumo, assimetria e densidades, boxplots, outliers, gráficos de dispersão, correlograma e distribuição geográfica para as 7 variáveis quantitativas: SO2, temp, manu, popul, wind, precip e predays.

Em seguida, na Análise de Componentes Principais, vimos que com 4 componentes conseguimos explicar cerca de 97.52% da variabilidade total do conjunto de dados.

No Escalonamento Multidimensional, vimos que a partir de \(k=5\) coordenadas, a matriz de distância original pode ser recuperada quase que totalmente.

E, por fim, vimos no t-SNE que podemos representar dados de altas dimensionalidades em 2 dimensões.

Como próximos passos, poderíamos, por exemplo, fazer uma regressão linear entre SO2 e as seis componentes principais calculadas ou então uma clusterização com os resultados do MDS.

Todas as análises realizadas abriram diversas possibilidades para explorar melhor a relação entre SO2 e as demais variáveis.