O objetivo principal é a visualização de conjunto multivariados. Nesta tarefa, os gráficos bivariados são extermamente importantes e costumam ser a primeira alternativa para a obtençao de padrão gráfico dos dados. A exibição da relação entre duas variáveis dependerá das medidas de cada variável, ou seja, se ela são nominais, ordinais, discretas ou contínuas. Há tipos de gráficos e cuidados adequados a cada uma dessas situações. Na visualização bivariada, podemos resslatar as técnicas baseadas em pontos, como o caso do scatter plot matrix e as técnicas de redução de dimensionalidade que buscam projetar em um espaço 2D os pontos de um espaço multidimensional. Essas técnicas são extensamente exploradas em multivariada e podemos citar: análise de componentes principais, análise fatorial, escalonamento multidimentsional. Essa última uma técnica bastante interessante em sua concepção de manter a ordenação dos pontos. Vamos aqui discutir as possibilidades e os problemas na construção de gráficos de dispersão e de scatter plot matrix. É sempre importante insepcionar a matriz de correlações dos dados nesse contexto de gráficos bivariados. Nesse sentido estaremos discutindo algumas possibilidades na visualização das matrizes de correlações que facilitarão bastante a exploração multivariada.
Apresentamos a seguir comandos e detalhes de gráficos bivariados usando o bastante conhecido conjunto de dados `iris, que apresenta 4 variáveis contínuas relacionadas com medidas morfológicas das flores iris e uma categórica que tem como níveis três de suas espécies.
Tipo bastante conhecido de gráfico, mas que eventualmente pode resumir demais os dados.
# Box-plots - quantitativa vs. categórica
variaveis <- names(iris[-5])
par(mfrow = c(2, 2))
for(i in 1: length(variaveis)) {
with(iris, {
dados <- eval(parse(text = variaveis[i]))
boxplot(dados ~ Species, data = iris, main = variaveis[i],
cex.axis = 0.85)
})
}
par(mfrow = c(1, 1))
O violin plot é um método para visualizar dados numéricos. É semelhante a um box-plot, agregando entretanto uma suavização por núcleo da densidade, rebatida girado em cada lado. Enquanto um box-plot mostra apenas as estatísticas resumo, como mediana e intervalos interquartílicos, o violin plot apresenta as distribuição completa dos dados. assim, o box-plot é conveniente para comparar estatísticas, mas não permite ver variações nos dados. Para distribuições multimodais (aquelas com vários picos), isso pode ser particularmente limitante. Por outro lado, o violin plot é um híbrido ddo box plot com a suavização por núcleo-estimador, podendo mostrar os picos dos dados
library(vioplot)
nomes <- levels(iris$Species)
with(iris,{
#for(i in 1:3) assign(paste0("x",i), Petal.Width[Species == nomes[i]])
for(i in 1:3) assign(paste("x",i, sep=""),
Petal.Width[Species == nomes[i]])
vioplot(x1, x2, x3, names = nomes, col = "lightblue", horizontal = TRUE)
# col = "gold"
title("Violin Plots de Petal.Width")
}) # end with
O diagrama de dispersão (scatter plot) é a forma mais simples de visualização gráfica de duas variáveis quantitativas, com cada variável representada em um dos eixos. Por exemplo, com o conjunto de dados iris, plotamos Petal.Length vs. Petal.Width para as 150 flores do conjunto de dados.
Uso do formato fórmula (Petal.Width ~ Petal.Length) como argumento, já que dessa maneira pode-se declarar o conjunto de dados com o argumento data = iris
plot(Petal.Width ~ Petal.Length, data = iris, main="Conjunto de Dados - Íris",
xlab = "Comprimento pétala (cm)", ylab = "Largura pétala (cm)")
Scatter plot das variáveis Petal.Widthe Petal.Length, com o fator Species representado por carateres distintos (argumento pch)
# Scatter plot com o fator 'Species' - Pétalas (por caracter)
plot(Petal.Width ~ Petal.Length, data = iris, xlab = "Comprimento (cm)",
pch = c(23, 24, 25)[unclass(iris$Species)], ylab = "Largura (cm)",
main = "Conjunto de Dados Íris - Pétalas")
# legenda
legend("bottomright", legend = c("Setosa", "Versicolor", "Virginica"),
pch = c(23, 24, 25), bty = "n")
Scatter plot das variáveis Sepal.Widthe Sepal.Length, com o fator Species representado por cores distintas (argumento bg) e uso de caractere que pode ser preenchido (pch = 21)
# Scatter plot com fator 'Species' - Sépalas (por cor)
plot(Sepal.Width ~ Sepal.Length, data = iris, xlab = "Comprimento (cm)",
bg = c("red","green3","blue")[unclass(iris$Species)], pch = 21,
ylab = "Largura (cm)", main = "Conjunto de Dados Íris - Sépalas")
# legenda
legend("bottomright", legend = c("Setosa", "Versicolor", "Virginica"),
pch = 20, col = c("red", "green3", "blue"), cex = 0.9, bty = "n")
Tenta-se uma visualização das 5 variáveis, nos eixos utilizam-se a largura e o comprimento e aos invés de pontos usam-se letras para cada espécie (´Ss,VvéRr, com letras maiúsculas para as medidas das Sépalas e minúsculas, para as sépalas. São usadas as funçõesmatplot()ematpoints()`.
# Scatter plot com matplot - Todos comprimentos e larguras
iS <- iris$Species == "setosa"
iV <- iris$Species == "versicolor"
iG <- iris$Species == "virginica"
op <- par(bg = "bisque")
matplot(c(1, 8), c(0, 4.5), type = "n",
xlab = "Comprimento (cm)", ylab = "Largura (cm)",
main = "Dimensões de Pétalas e Sépalas de Flores de Íris")
matpoints(iris[iS,c(1,3)], iris[iS,c(2,4)], pch = "sS", col = c(2,4))
matpoints(iris[iV,c(1,3)], iris[iV,c(2,4)], pch = "vV", col = c(2,4))
matpoints(iris[iG,c(1,3)], iris[iG,c(2,4)], pch = "rR", col = c(2,4))
legend(1, 4, c(" Pétalas Setosa", " Sépalas Setosa",
"Pétalas Versicolor", "Sépalas Versicolor",
" Pétalas Virginica", " Sépalas Virginica"), cex=0.9,
pch = "sSvVrR", col = rep(c(2,4), 3), bty = "n")
par(op)
Uso do comando matplot() de array para criação do plot, com uma camada para cada espécie.
# Scatter plot com array e matplot
# Criação da array
nome.var <- colnames(iris)[-5]
nome.esp <- as.character(unique(iris[,"Species"]))
# criação array (uma matriz 50x3 para cada espécie)
iris.S <- array(NA, dim = c(50, 4, 3),
dimnames = list(NULL, nome.var, nome.esp))
# alocação dos valores no array, por camada
for(i in 1:3) iris.S[, , i] <- data.matrix(iris[1:50 + 50*(i-1), -5])
# 1a. observação, 2a. camada
iris.S[1, , 2]
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 7.0 3.2 4.7 1.4
# 1a. observação, todas as camadas
iris.S[1, , ]
## setosa versicolor virginica
## Sepal.Length 5.1 7.0 6.3
## Sepal.Width 3.5 3.2 3.3
## Petal.Length 1.4 4.7 6.0
## Petal.Width 0.2 1.4 2.5
# scatter plot - bg apenas na região do gráfico
matplot(iris.S[, "Petal.Length", ], iris.S[, "Petal.Width", ], type = "n",
sub = paste(c("S: ", "C: ", "V: "), dimnames(iris.S)[[3]],
sep = "", collapse= ", "),
xlab = "Petal.Length", ylab = "Petal.Width",
main = "Conjunto de Dados Iris")
rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4],
col = "bisque")
matpoints(iris.S[, "Petal.Length", ], iris.S[, "Petal.Width", ], pch = "SCV",
col = rainbow(3, start = 0.8, end = 0.1))
ggplot2:Necessária a instalação dos pacotes ggplot2e gridExtra
# Scatter plot com pacote ggplot2
library(ggplot2)
library(gridExtra)
Objeto gráfico (graphical object - grob) sp1 contém os pares ordenados das variáveis relacionados com a fórmula Petal.Width ~ Petal.Length, com representação simples por pontos (geo_point).
# Plot com pontos default
sp1 <- ggplot(iris, aes(x = Petal.Length, y = Petal.Width))
sp1 + geom_point() +
xlab("Comprimento (cm)") + ylab("Largura (cm)") +
ggtitle("Conjunto de Dados Íris - Pétalas")
O grob sp2 é criado a partir do objeto sp1e contém os pontos coloridos por default, de acordo com os níveis do fator SPecies.
# cor p/ nível do fator Species
sp2 <- sp1 + geom_point(aes(color = factor(Species))) +
xlab("Comprimento (cm)") + ylab("Largura (cm)") +
ggtitle("Conjunto de Dados Íris - Pétalas")
sp2
O grob sp3é criado a partir do objetos sp1e contém os pontos coloridos conforme o usuário, de acordo com os níveis do fator Species.
# Cores conforme usuário
sp3 <- sp1 + geom_point(aes(color=factor(Species))) +
scale_color_manual(values = c("orange", "purple", "gray")) +
xlab("Comprimento (cm)") + ylab("Largura (cm)") +
ggtitle("Conjunto de Dados Íris - Pétalas")
sp3
O grob sp4é criado a partir do objetos sp1e contém os pontos com formato de caractere conforme o usuário, de acordo com os níveis do fator Species.
# forma e tamanho dos ponto conforme nível fator
sp4 <- sp1 + geom_point(aes(shape = factor(Species))) +
xlab("Comprimento (cm)") + ylab("Largura (cm)") +
ggtitle("Conjunto de Dados Íris - Pétalas")
sp4
Gráfico com cores default e formato dos pontos conforme usuário.
# Formato pontos cfe. usuário
sp1 + geom_point(aes(shape = factor(Species), color = factor(Species))) +
scale_shape_manual(values = c(0, 2, 3)) +
theme(legend.position = c(0,1), legend.justification = c(0,1)) +
theme(legend.title = element_blank()) +
xlab("Comprimento (cm)") + ylab("Largura (cm)") +
ggtitle("Conjunto de Dados Íris - Pétalas")
Uso da função ggplot() para criação do gráfico
# Scatterplot comprimentos + espécies
ggplot(iris, aes(x = Sepal.Length, y = Petal.Length, color = Species)) +
geom_point()+
xlab("Comprimento da sépala (cm)") + ylab("Comprimento da pétala (cm)") +
ggtitle("Conjunto de Dados Íris") +
scale_color_discrete(name ="Espécies") +
theme(legend.position = c(1, 0), legend.justification = c(1, 0))
Criação de gráfico múltiplo dos grobs sp2, sp3e sp4, com a função grid.arrange().
# painel com os gráficos
grid.arrange(sp2, sp3, sp4, nrow=1)
Uso da função rug()para apresentar dispersão dos valores da variável nas margens do gráfico.
plot(Petal.Width ~ Petal.Length, data = iris, main = "Conjunto de Dados Iris",
xlab = "Petal.Length", ylab = "Petal.Width")
with(iris, {
rug(jitter(Petal.Length, amount = 0.05), side = 1, ticksize = 0.02,
col = "light blue")
rug(jitter(Petal.Width, amount = 0.05), side = 2, ticksize = 0.02,
col = "light blue")
}) # end with
Construção de áreas para 3 gráficos por meio da função layout() e construção sequenciada dos três gráficos: gráfico de dispersão de Petal.Lengthe Petal.Width, histograma da variável Petal.Lengthe box-plot da variável Petal,Width.
# Dispersão com marginais (histograma e boxplot)
layout(matrix(c(2,0,1,3), nrow = 2, byrow = T), respect = T,
widths = c(2, 1), heights = c(1, 2))
xlim <- range(iris$Petal.Length) * 1.1
plot(Petal.Width ~ Petal.Length, data = iris, type = "n", xlim = xlim,
cex.lab = 0.9, xlab = "Petal.Length", ylab = "Petal.Width")
points(Petal.Width ~ Petal.Length, data = iris, cex = 0.6,
pch = c("S", "C", "V")[unclass(Species)],
col = rainbow(3, start = 0.8, end = 0.1)[unclass(Species)])
hist(iris$Petal.Length, freq = F, ylab = "Densidade", main = NULL, xlab = NULL)
boxplot(iris$Petal.Width)
ggpubr e ggExtra:O pacote ggExtra permite o uso da função ggMarginal() que constrói vários tipos de gráficos univariados.
# Dispersão com marginais - pacote ggpubr
library(ggpubr) # comando ggscatter
library(ggExtra) # comando ggMarginal
O pacote ggpubroferece a função ggscatter(), que constrói o objetos gráficos (grobs) de diagramas de dispersão, com inúmeros argumentos que permitem flexibidade na construção dos gráficos. O objetos p tem características relacionados com cores, tamnaho e trasparência das cores.
p <- ggscatter(iris, x = "Petal.Length", y = "Petal.Width",
color = "Species", palette = "jco",
size = 3, # size: tamanho dos pontos
alpha = 0.6) # alpha: transparência cores
Marginal plot construído com a função ggMarginal(), com gráficos de densidade em amabas as margens, obtidas com argumento type = "density".
ggMarginal(p, type = "density")
Margens do gráficos com Box-plots, obtidos com o argumento type = "boxplot", da função ggMarginal().
# marginais: box-plot
ggMarginal(p, type = "boxplot")
Margens do gráficos com histogramas, obtidos com o argumento type = "histogram", da função ggMarginal().
# marginais: histograma
ggMarginal(p, type = "histogram")
Gráficos de dispersão combinados com gráficos marginais, usando os pacotes ggpubr e cowplot com construção de gráfico por meio de graphical objects (grobs ). Nos casos abaixo usa-se o pacote cowplotporque o pacote ggextranão lida com vários grupos. Os gráficos apresentam um aparência melhor.
library(ggpubr)
library(cowplot) # ggExtra não lida com vários grupos
Gráfico de dispersão das variáveis Petal.Lengthe Petal.Width, combinados com o histograma de Petal.Length na posição 3 e o box-plot de Petal.Width, na posição 4.
# marginais plots - x: histograma e y: box-plot
# gráfico paralelo ao eixo x - histograma
xplot <- gghistogram(iris, x = "Petal.Length", bins = 15, ylab = "",
ggtheme = theme_bw())
# gráfico paralelo ao eixo y
yplot <- ggboxplot(iris, y = "Petal.Width", alpha = 0.5,
xlab = "", ggtheme = theme_bw())
# limpeza dos gráficos
p.b <- p + rremove("legend")
xplot <- xplot + clean_theme() + rremove("legend")
yplot <- yplot + clean_theme() + rremove("legend")
# montagen gráfico com pacote cowplot
plot_grid(xplot, NULL, p.b, yplot, ncol = 2, align = "hv",
rel_widths = c(2, 1), rel_heights = c(1, 2))
Gráfico de dispersão das variáveis Petal.Lengthe Petal.Width, combinados com gráficos suavizados de densidade por grupo coloridos. Gráficos de suavização de densidade facilitam a visualização dos dados e auxiliam muito na identificão de padrão dos dados.
# Scatter plot e density marginal plots coloridos por grupos
# scatter plot por grupo
sp1 <- ggscatter(iris, x = "Petal.Length", y = "Petal.Width", color = "Species",
palette = "jco", size = 3, alpha = 0.6) + border()
# densidades marginais por grupos
xplot1 <- ggdensity(iris, "Petal.Length", fill = "Species",
palette = "jco")
yplot1 <- ggdensity(iris, "Petal.Width", fill = "Species",
palette = "jco") + rotate()
# limpeza dos gráficos
sp1 <- sp1 + rremove("legend")
yplot1 <- yplot1 + clean_theme() + rremove("legend")
xplot1 <- xplot1 + clean_theme() + rremove("legend")
# montagem do gráfico com pacote cowplot
plot_grid(xplot1, NULL, sp1, yplot1, ncol = 2, align = "hv",
rel_widths = c(2, 1), rel_heights = c(1, 2))
Gráfico de dispersão das variáveis Petal.Lengthe Petal.Width, combinados com box-plots coloridos por grupo.
# Scatter plot e box-plots por grupos
# scatter plot por grupo + ggtheme
sp2 <- ggscatter(iris, x = "Petal.Length", y = "Petal.Width", color= "Species",
palette = "jco", size = 3, alpha = 0.6, ggtheme = theme_bw())
# box-plots marginais
xplot2 <- ggboxplot(iris, x = "Species", y = "Petal.Length", color = "Species",
fill = "Species", palette = "jco", alpha = 0.5,
ggtheme = theme_bw()) + rotate()
yplot2 <- ggboxplot(iris, x = "Species", y = "Petal.Width", color = "Species",
fill = "Species", palette = "jco", alpha = 0.5,
ggtheme = theme_bw())
# limpeza dos gráficos
sp2 <- sp2 + rremove("legend")
yplot2 <- yplot2 + clean_theme() + rremove("legend")
xplot2 <- xplot2 + clean_theme() + rremove("legend")
# montagem do gráfico com pacote cowplot
plot_grid(xplot2, NULL, sp2, yplot2, ncol = 2, align = "hv",
rel_widths = c(2, 1), rel_heights = c(1, 2))
Dependendo da situação a reta de regressão ajuda a entender o padrão dos dados, no caso indicando a presença de clusters nos dados.
# Scatter plot com linha de regressão
# reta de regressão y = Petal.Width, x = Petal.Length
petal.lm <- lm(Petal.Width ~ Petal.Length, data = iris)
summary(petal.lm)
##
## Call:
## lm(formula = Petal.Width ~ Petal.Length, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.56515 -0.12358 -0.01898 0.13288 0.64272
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.363076 0.039762 -9.131 4.7e-16 ***
## Petal.Length 0.415755 0.009582 43.387 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2065 on 148 degrees of freedom
## Multiple R-squared: 0.9271, Adjusted R-squared: 0.9266
## F-statistic: 1882 on 1 and 148 DF, p-value: < 2.2e-16
# coeficiente de determinação da regressão
cor(iris$Petal.Length, iris$Petal.Width)^2
## [1] 0.9271098
# scatter plot com reta de regressão
r2 <- summary(petal.lm)$r.squared
# scatter plot e reta de regressão.
plot(Petal.Width ~ Petal.Length, data = iris,
main = "Conjunto de Dados Íris - Pétalas", xlab = "Comprimento (cm)",
ylab = "Largura (cm)")
# reta de regressão
abline(petal.lm, lty = 2, col = "blue")
# R2 informado no plot com uso de anotação matemática e variável
text(x = 1, y = 2.25, bquote(r^2 == .(round(r2,2), 3)), pos = 4, cex = 0.90)
O pacote ggplot2amplia as possibilidades de visualização da relação entre as duas variáveis.
library(ggplot2)
O default do ggplot2é adicionar linha de regressão com banda de confiança de 95%.
# scatter plot e regressão c/ banda de confiança
ggplot(iris, aes(x = Petal.Length, y = Petal.Width)) +
geom_point(shape = 1) + # círculos vazios
geom_smooth(method = lm) # adiciona linha de regressão
## `geom_smooth()` using formula 'y ~ x'
# default: banda de confianção de 95%
O argumento se = FALSE retira a banda de confiança do gráfico.
# scatter plot e regressão s/ banda de confiança
ggplot(iris, aes(x = Petal.Length, y = Petal.Width)) +
geom_point(shape = 1) + # círculos vazios
geom_smooth(method = lm, # adiciona linha de regressão
se = FALSE) # sem banda de confiança
## `geom_smooth()` using formula 'y ~ x'
A função geom_smooth() adiciona linha de regressão suavizada. Seu argumento methodtem como valores possíveis "lm" (a reta de regressão), "glm", "gam", "loess" ou uma função de regressão qualquer. A função geom_smooth()sem explicitação de arumento, desenhará uma linha loess com banda de confiança, entre os dados. LOWESS (Locally Weighted Scatterplot Smoothing), algumas vezes denominado LOESS (locally weighted smoothing), é uma ferramenta popular usada em análise de regressão que cria uma linha de regressão suavizada em um scatter plot. Ela ajuda a visualizar a relação entre as variáveis.
# scatter plot e loess c/ banda de confiança
ggplot(iris, aes(x = Petal.Length, y = Petal.Width)) +
geom_point(shape = 1) + # círculos vazios
geom_smooth() # adiciona curva suavizada loess
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
# com banda de confiança
O argumento se = FALSE retira a banda de confiança do gráfico.
# scatter plot e loess s/ banda de confiança
ggplot(iris, aes(x = Petal.Length, y = Petal.Width)) +
geom_point(shape = 1) + # círculos vazios
geom_smooth(method = "loess", # adiciona curva suavizada loess
se = FALSE) # sem banda de confiança
## `geom_smooth()` using formula 'y ~ x'
No ggplot2os códigos dos formatos dos pontos são dados pela tabela abaixo:
Assim, pode-se customizar o formato dos pontos por nível do fator Species, optando-se por: círculo vazio (1), triângulo (2) e cruz (3).
# ajusta pontos por Species, com formatos diferentes
ggplot(iris, aes(x = Petal.Length, y = Petal.Width, shape = Species)) +
geom_point() +
scale_shape_manual(values = c(1, 2, 3)) # círculo vazio, triângulo e
# cruz
Aqui, o argumento shape = 1da função geom_point() plotará círculos vazios nos pontos dos dados, os quais são coloridos de acordo com o fator Species. A função aes(group = 1) imporá uma única reta de regressão pelos dados.
# scatter plot com cores por Species
ggplot(iris, aes(x = Petal.Length, y = Petal.Width, color = Species)) +
geom_point(shape = 1) + # círculos vazios
geom_smooth(method = lm , # adiciona reta de regressão
aes(group = 1), # única reta
se = FALSE) # sem banda de confiança
## `geom_smooth()` using formula 'y ~ x'
O default é uma reta de regressão por nível da variável Speciesjá que foi imposto no comando aes()geral. Foi usado um pallete um pouco mais escuro que o normal.
# scatter plot com cores e linhas de regressão distintos p/ grupo
ggplot(iris, aes(x = Petal.Length, y = Petal.Width, color = Species)) +
geom_point(shape = 1) + # círculos vazios
scale_colour_hue(l = 50) + # paleta um pouco mais escura que a normal
geom_smooth(method = lm) # adiciona retas de regressão
## `geom_smooth()` using formula 'y ~ x'
# com banda de confiança
O argumento se = FALSE retira as bandas de confiança do gráfico.
# scatter plot com cores diferentes e linhas de regressão s/ ic
ggplot(iris, aes(x = Petal.Length, y = Petal.Width, color = Species)) +
geom_point(shape = 1) +
scale_colour_hue(l = 50) + # paleta um pouco mais escura que a normal
geom_smooth(method = lm, # adiciona retas de regressão
se = FALSE) # sem banda de confiança
## `geom_smooth()` using formula 'y ~ x'
O argumento fullrange = TRUE estende as retas de regressão pela área do gráfico, além do domínio dos dados.
# Estende linhas de regressão além do domínio dos dados
ggplot(iris, aes(x = Petal.Length, y = Petal.Width, color = Species)) +
geom_point(shape = 1) + # círculos vazios
scale_colour_hue(l = 50) + # paleta um pouco mais escura que a normal
geom_smooth(method = lm, # adiciona retas de regressão
se = FALSE, # sem banda de confiança
fullrange = TRUE) # estende as retas de regressão
## `geom_smooth()` using formula 'y ~ x'
A função facet_wrap()apresenta o diagrama de dispersão das variáveis, com linha de regressão por cada nível do fator Species.
# painéis por Species c/ linhas de regressão
ggplot(iris, aes(x = Petal.Length, y = Petal.Width, group = Species)) +
geom_point(shape = 1, size = 1.5) + # tamanho dos pontos
stat_smooth(method = lm, se = FALSE) + # retas regressão s/ ic
facet_wrap( ~ Species) # painéis por Species
## `geom_smooth()` using formula 'y ~ x'
Uma questão importante ao trabalharmos com diagramas de dispersão é quando há muitos empates ou dados com valores muito próximos. Para exemplificar esse efeito vamos simular um conjunto de 20 observações crescendo linearmente, com um ruído aditivo
# Conjunto de dados com 20 observações com sobreposição
set.seed(666)
# geração dados crescente com ruído (n = 20)
dados <- data.frame(cond = rep(c("A", "B"), each = 10),
xvar = 1:20 + rnorm(20, sd = 3),
yvar = 1:20 + rnorm(20, sd = 3))
O diagrama de dispersão com linha de regressão estimada por loess será obtida pelos comandos abaixo:
# scatter plot com loess dos dados originais
ggplot(dados, aes(x = xvar, y = yvar)) +
geom_point(shape = 19) + # círculos sólidos
geom_smooth(se = FALSE) # suavização por loess
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
A perturbação afetou a estimação da linha de regressão, principalmente na região com menos informação (menor quantidade de pontos). Agora, arredondamos os valores do conjunto de dados para o múltiplo de 5 mais próximo, forçando empates.
# arredondando xvar and yvar para o 5 mais próximo
dados$xarr <- round(dados$xvar/5)*5
dados$yarr <- round(dados$yvar/5)*5
O argumento àlpha` controla a opacidade da cor dos pontos, ou seja, com mais pontos a cor ficará menos opaca devido à sobreposição obtida.
# pontos parcialmente transparente, c/ opacidade 1/4
ggplot(dados, aes(x = xarr, y = yarr)) +
geom_point(shape = 19, # círculos sólidos
alpha = 1/4) + # opacidade 1/4
geom_smooth(method = "loess", se = FALSE)
## `geom_smooth()` using formula 'y ~ x'
Percebe-se a região com maior concentração de pontos, auxiliado pelo valor escolhido de alpha = 0.25. Uma outra possibilidade de gráfico é introduzir uma perturbação na posição do ponto por meio da função osition_jitter() aplicada ao argumento position.
# Perturbação dos pontos
ggplot(dados, aes(x = xarr, y = yarr)) +
geom_point(shape = 1, # Use hollow circles
position = position_jitter(width = 1, height = 0.5)) +
geom_smooth(method = "loess", se = FALSE)
## `geom_smooth()` using formula 'y ~ x'
Pode-se perceber um pouco mais as concentrações de pontos nas sub-regiões do gráfico.
ggpubr:library(ggpubr)
Uso da função ggscatter(), com reta de regressão (add = "reg.line"), intervalo de confiança (conf.int = TRUE) e cores da linha de regressão e fundo. A função stat_cor() adiciona o valor da correlação de p-valor de teste de significância.
# plot básico
ggscatter(iris, x = "Petal.Length", y = "Petal.Width",
add = "reg.line", # Adiciona reta de regressão
conf.int = TRUE, # Adiciona intervalo de confiança
add.params = list(color = "blue",
fill = "lightgray")
) +
stat_cor(method = "pearson", label.x = 2.5, label.y = 2) # Adiciona correlação
## `geom_smooth()` using formula 'y ~ x'
O formato dos pontos do pacote ggpubrsão dados pelos códigos abaixo:
# formatos de pontos
show_point_shapes()
Optou-se mudar o formato dos pontos para losângulos (shape = 18).
# Muda formato do ponto
ggscatter(iris, x = "Petal.Length", y = "Petal.Width", shape = 18,
add = "reg.line", # Adiciona reta de regressão
conf.int = TRUE, # Adiciona intervalo de confiança
add.params = list(color = "blue",
fill = "lightgray")
)
## `geom_smooth()` using formula 'y ~ x'
Gráfico com formato, cores, correlações e retas de regressão distintas, por nível do fator Species. OfatorSpeciesestá atribuído aos argumentoscoloreshapee usa-se a funçãoaes(color = Species)como arumneto da funçãostat_cor()`.
# gráfico colorido por grupos
ggscatter(iris, x = "Petal.Length", y = "Petal.Width",
add = "reg.line",
conf.int = TRUE,
color = "Species", palette = "jco", # Cor por Species
shape = "Species" # Formato pontos por Species
) +
stat_cor(aes(color = Species), label.x = 3) # correlação por Species
## `geom_smooth()` using formula 'y ~ x'
As retas de regressão são estendidas com o argumento fullrange = TRUE, com a adição de desnsidade marginal com o argumento rug = TRUE.
# Estendendo reta de regressão e adição de densidade marginal (rug)
ggscatter(iris, x = "Petal.Length", y = "Petal.Width",
add = "reg.line",
color = "Species", palette = "jco",
shape = "Species",
fullrange = TRUE, # Estendendo a reta de regressão
rug = TRUE # Adicionando rug marginal
)+
stat_cor(aes(color = Species), label.x = 3)
## `geom_smooth()` using formula 'y ~ x'
Traçar elipses de probabilidade pode ser uma maneira útil de visualizar distribuições normais bivariadas. Uma elipse de probabilidade representa uma curva de probabilidade constante dentro do qual se encontra uma certa porcentagem da distribuição. A largura e a orientação das elipses de probabilidade podem fornecer informações sobre a correlação entre as variáveis de interesse. É usada a função ggscatterdo pacote ggpubr.
library(ggpubr)
Elipses de 95% de confiança são plotadas com os argumentos ellipse = TRUEe ellipse.level = 0.95.
# elipses de concentração
ggscatter(iris, x = "Petal.Length", y = "Petal.Width",
color = "Species", palette = "jco",
shape = "Species",
ellipse = TRUE, ellipse.level = 0.95)
O polígono convexo é o polígono cujos lados conectam O menor conjunto convexo de pontos que contém todas as observações. Ele é obtido com o argumento ellipse.type = "convex".
# polígono convexo
ggscatter(iris, x = "Petal.Length", y = "Petal.Width",
color = "Species", palette = "jco",
shape = "Species",
ellipse = TRUE, ellipse.type = "convex")
Os argumentos mean.point = TRUEe star.plot = TRUE plotam o centróide dos dados, conectando-o a cada ponto.
# Adiciona média dos pontos por grupos e estrelas
ggscatter(iris, x = "Petal.Length", y = "Petal.Width",
color = "Species", palette = "jco",
shape = "Species",
ellipse = TRUE,
mean.point = TRUE,
star.plot = TRUE)
Há técnicas de estimação não-paramétrica de funções que são bastante utilizadas na visualização de funções de densidade de probabilidade a partir dos dados. Essas técnicas permitem identificar padrões na densidade que podem não ser percebidos com as ferramentas habituais, como o histograma. Dentre outras metodologias de estimação não paramétrica, destaca-se o núcleo estimador da função de densidade. Seja uma função de densidade de probabilidades \(f_X\), desconhecida, \(x\) um ponto qualquer do contradomínio da variável aleatória \(X\) e \(x_1, x_2, \dots, x_n\) uma amostra aleatória de população com densidade \(f\). A estimativa por núcleo da densidade no ponto \(x\) do contradomínio de \(X\) é dada pela expressão a seguir:
\[ \hat{f}_X(x) = \frac{1}{nh} \sum_{i=1}^nK \left( \frac{x - x_i}{h} \right) \]
onde, \(K\) é uma função de densidade qualquer denominada núcleo e \(h\) é o parâmetro de suavização, largura de banda ou janela. O núcleo estimador pode ser encarado com uma média ponderada pela densidade, de maneira que os vizinhos mais próximos do ponto \(x\) influencial mais a densidade do que os pontos mais distantes. Em geral, usa-se a função de densidade normal como o núcleo \(K\), embora não haja restrição para usar qualquer outra função de densiade. Por sua vez, a janela \(h\) controla a influência dos vizinhos: valores grande da janela tendem a considerar a influência de vizinhos mais distantes, enquanto os valores pequenos da janela dão menos peso aos vizinhos mais distante. A escolha do valor da janela é crucial no uso da técnica e há várias metodologias para buscar seu valor ótimo. Valores pequenos da janela deixam passar mais ruído na estimação e valores grandes suavizam demais aestimação da densidade. Sejam os dados abaixo relativos a uma amostra aleatória de uma população com função de densidade de probabilidade \(f\).
# densidade univariada - funcionamento do núcleo
# observações
amostra<- c(0, 1, 1.1, 1.5, 1.9, 2.8, 2.9, 3.5)
n <- length(amostra)
Para considerar o contradomínio da variável aleatória, alarga-se um pouco a amplitude observada doss dados,
xgrid <- seq(from = min(amostra) -1, to = max(amostra) + 1, by = 0.01)
Com a função sapply(), calcula-se facilmente, as estimativas por núcleo da função de densidade de probabilidade desconhecida \(f\). Usa-se o valor da janela \(h = 0,4\).
# janela
h <- 0.4
# aplicação do núcleo no grid
# contribuição núcleo de cada valor observado a cada ponto (matriz 551 x 8)
contr <-sapply(amostra, function(a) dnorm((xgrid - a)/h)/(n * h))
# gráfico da estimação e das contribuições
plot(xgrid, rowSums(contr), lwd = 2, xlab = "x",
ylab = expression(hat(f)(x)), type = "l",
main = "Estimação Densidade Univariada\nNúcleo Estimador" )
# observações
#rug(amostra, lwd = 2)
A estimação da função de densidade, baseada nos 8 valores oservados e numa janela \(h = 0,4\), apresenta multimodalidade. Poderíamos entender a estimativa da função de densidade como a soma, em cada ponto \(x\), do contradomínio, calibradas por \(1/nh\).
# aplicação do núcleo no grid
# contribuição núcleo de cada valor observado a cada ponto (matriz 551 x 8)
contr <-sapply(amostra, function(a) dnorm((xgrid - a)/h)/(n * h))
# gráfico da estimação e das contribuições
plot(xgrid, rowSums(contr), lwd = 2, xlab = "x",
ylab = expression(hat(f)(x)), type = "l",
main = "Estimação Densidade Univariada\nNúcleo Estimador" )
# observações
rug(amostra, lwd = 2)
# plot das contribuições
out <- apply(contr, 2, function(b) lines(xgrid, b, lty = 2))
Podemos entender então a estimativa de densidade, como a soma de oito densidades com desvio padrão \(h\), cada uma participando com peso \(1/n\). Ou seja janela maiores, produzirão densidades mais achatadas em torno de cada ponto e produzirão estimativas mais suaves. Por exemplo, no caso em que \(h = 0,7\):
# janela mais larga - maior suavização
h <- 0.7
# contribuição núcleo de cada valor observado a cada ponto (matriz 551 x 8)
contr <-sapply(amostra, function(a) dnorm((xgrid - a)/h)/(n * h))
# gráfico da estimação e das contribuições
plot(xgrid, rowSums(contr), lwd = 2, xlab = "x",
ylab = expression(hat(f)(x)), type = "l",
main = "Estimação Densidade Univariada\nNúcleo Estimador" )
# observações
rug(amostra, lwd = 2)
# plot das contribuições
out <- apply(contr, 2, function(b) lines(xgrid, b, lty = 2))
Por outrolado, se usarmos janelas menores as densidades serão concentradas em torno de cada ponto observado e aestimativa da densidade terá menos suavidade, pois terá deixado passar mais rúido. Seja por exempelo a estimativa por núcleo usando uma janela \(h=0.2\):
# janela mais larga - maior suavização
h <- 0.2
# contribuição núcleo de cada valor observado a cada ponto (matriz 551 x 8)
contr <-sapply(amostra, function(a) dnorm((xgrid - a)/h)/(n * h))
# gráfico da estimação e das contribuições
plot(xgrid, rowSums(contr), lwd = 2, xlab = "x",
ylab = expression(hat(f)(x)), type = "l",
main = "Estimação Densidade Univariada\nNúcleo Estimador" )
# observações
rug(amostra, lwd = 2)
# plot das contribuições
out <- apply(contr, 2, function(b) lines(xgrid, b, lty = 2))
Há vários métodos para escolha da janela ótima na estimação da função de densidade de uma amostra de dados. Afunção kde(), do pacote ksfaz essa escolha por default ou aceita os valores desejado de janela.
# Estimação densidade
library(ks)
Obtemos, com os comandos abaixo, a estimativa da variável Petal.Length, do conjunto de dados iris, que sabemos apresenta clusters significativos e bem discriminados:
# Densidade univariada - Petal.Length
fhat <- kde(x = iris[, 3])
# densidade suavizada por núcleo estimador
plot(fhat, cont = 50, col.cont = "blue", cont.lwd = 2,
ylab = expression(hat(f)(x)),
xlab = "Comprimento pétala (cm)")
Gráfico da densidade de Petal.Length apresenta bimodalidade, indicando graficament os intervalos de concentração de densidade. Percebe-se também uma oscilação na densidade nas regiões nas quais sabemos há misturas de espécies diferentes. Repetimos o procedimento para a variável Petal.Widthe observamos comportamento similar.
# Densidade univariada - Petal.Width
fhat <- kde(x = iris[, 4])
plot(fhat, cont = 50, col.cont = "blue", cont.lwd = 2,
ylab = expression(hat(f)(x)),
xlab = "Largura pétala (cm)")
Já conhecemos o grafico de dispersão das variáveis Petal.Length e Petal.Width. Uma visualização interessabte seria de sua densidade bivariada. Os gráficos das curvas de densidade é dado pela mesma função kde(), só que agora usando uma matriz dde dados contendo seuus valores observdos.
# Densidade bivariada - Petal.Length e Petal.WWidth
# contour plot
fhat <- kde(x = iris[, 3:4])
# contour plot
plot(fhat, display = "filled.contour2", cont = seq(10, 90, by = 10))
O gráfico das curvas de densidade permite que visualizemos as regiões de concentrações de densidade, indicando um terceiro grupo (com uma concemtração mais tênue de densidade) que não apareceu nos gráficos univariados anteriores dessas variáveis. Pode-se obter a visão em perspectiva da superfície da densidade dos dados usando o argumento display = "persp.
# Densidade bivariada - Petal.Length e Petal.WWidth
# perspectiva
plot(fhat, display = "persp", thin = 3, border = 1, col = "white")
O pacote MASS tem a função kde2d() que constrói a matriz com o grid de pontos e os valores de densidade necessária para usar como argumento da função `contour().
library(MASS)
library(colorspace)
As estimativas densidade bivariada por núcleo necessárias para traçar as curvas de mesma densidade são:
# estimativas densidade bivariada por núcleo - contour plot
plot(Petal.Width ~ Petal.Length, data = iris, cex = 0.75,
xlab = "Comprimento pétala (cm)",
ylab = "Largura pétala (cm)", main = "Estimação de Densidade Bivariada")
# estimação densidade - package: MASS
k <- kde2d(iris[, 3], iris[, 4])
cnt <- contourLines(k$x, k$y, k$z)
n <- length(cnt)
# definição de cores - package: colorspace
cores <- rev(sequential_hcl(n))
# curvas de nível - desenho e cores
for( i in seq_len(n) ) lines(cnt[[i]], col=cores[i])
O pacote hexbin permite a construção de mapa de calor bivariado de muita qualidade e de uma maneira bem simplificada, com a função hexbinplot()
library(hexbin)
# estimativa densidade codificada por cor
hexbinplot(Petal.Width ~ Petal.Length, iris, aspect = 1)
É possível obtermos a densidade trivariada em objeto do tipo rglque facilita bastante a visualização da superfície da função de densidade trivariada. Instalem o pacote misc3de executem os comandos abaixo.
library(misc3d)
Executem os comandos fora do R Markdown e visualizem e analisem amplamente a densidade trivariada.
# Densidade trivariada - Objeto rgl
fhat <- kde(x = iris[,2:4])
plot(fhat, drawpoints = TRUE)
lattice:Usaremos a função xyplot () do pacote lattice para estudar a relação entre duas variáveis quantitativas, em conjunto ou não com uma ou mais fatores.
O conjunto de dados gapminderFive.txt, com 1.704 observações de 6 variáveis, referentes à expectativa de vida e renda per capita de 141 países, de 1952 a 2007. As variáveis desse conjunto de dados são:
country: identificação dos 141 países do estudo.
year: ano das observações, a partir de 1952, medidas em intervalos de 5 anos.
pop: população total, em unidades.
continent: continente, com níveis "Asia", "Europe", "Africa", "Americas"e "Oceania".
lifeExp: expectativa de vida ao nascer (caso os padrões de mortalidade permaneçam os mesmos)
gdpPercap: renda per capita, em US$.
O carregamento dos dados e verificação do conjunto de dados importado
# carregamento com read.delim()
gapminder <- read.delim("..//dados//gapminderDataFiveYear.txt")
# dimensões do conjunto de dados
dim(gapminder)
## [1] 1704 6
# estrutura do conjunto de dados
str(gapminder)
## 'data.frame': 1704 obs. of 6 variables:
## $ country : chr "Afghanistan" "Afghanistan" "Afghanistan" "Afghanistan" ...
## $ year : int 1952 1957 1962 1967 1972 1977 1982 1987 1992 1997 ...
## $ pop : num 8425333 9240934 10267083 11537966 13079460 ...
## $ continent: chr "Asia" "Asia" "Asia" "Asia" ...
## $ lifeExp : num 28.8 30.3 32 34 36.1 ...
## $ gdpPercap: num 779 821 853 836 740 ...
Devido à pequena quantidade de informações serão retiradas todas as observações referentes à Oceania.
gap5Dados <- droplevels(subset(gapminder, continent != "Oceania"))
dim(gap5Dados)
## [1] 1680 6
# valores das variáveis categóricas
unique(gap5Dados$continent)
## [1] "Asia" "Europe" "Africa" "Americas"
length(unique(gap5Dados$country))
## [1] 140
unique(gap5Dados$year)
## [1] 1952 1957 1962 1967 1972 1977 1982 1987 1992 1997 2002 2007
Vamos conhecer algumas funções gráfica do pacote lattice, principalmente xyplot().
library(lattice)
Para construção do diagrama de dispersão das variáveis contínuas do conjuntode dados gap5Dados, com a função xyplot(). Pode-se inserir linhas de grade na figura, por meio do argumento grid = TRUE ou linhas horizontais (grid = "h") e linhas verticias (grid = "v").
# sem guias
xyplot(lifeExp ~ gdpPercap, gap5Dados)
# com guias
xyplot(lifeExp ~ gdpPercap, gap5Dados, grid = TRUE)
Aparentemente, é necessário a log-transformação do eixo horizontal, isto é, para a variável de renda per capita (gdpPercap). Pode-se transformar os dados diretamente na fórmula, mas não é ideal, devido às marcas do eixo, que podem não ser fáceis de serem lidas. O ideal é comandar a transformação log por meio do argumento scale que, além de trsnaformar os dados, apresenta no eixo os valores da escala original.
# ---- Log-transformação do eixo x
# log-transformação pela fórmula
xyplot(lifeExp ~ log10(gdpPercap), gap5Dados, grid = TRUE)
# log-transformação - argumento scale
xyplot(lifeExp ~ gdpPercap, gap5Dados, grid = TRUE,
scales = list(x = list(log = 10)))
Uma outra estratégia para rotular os eixos é ajustar o argumento equispaced.log = FALSE. Há outras maneiras disponíveis no pacote latticeExtra. I wish the tick marks and grid lined up, but I’m trying not to let the perfect be the enemy of the good. Let’s press on.
# Rótulo de eixo de variável log-transformada
xyplot(lifeExp ~ gdpPercap, gap5Dados, grid = TRUE,
scales = list(x = list(log = 10, equispaced.log = FALSE)))
O argumento typepode ser usado para acrescentar elementos no gráfico, com, por exemplo, type = c("p", "r") que acrescenta uma linha de regressão no gráfico
# pontos (default)
xyplot(lifeExp ~ gdpPercap, gap5Dados, grid = TRUE,
scales = list(x = list(log = 10, equispaced.log = FALSE)),
type = "p")
# pontos e a reta de regressão
xyplot(lifeExp ~ gdpPercap, gap5Dados,
grid = TRUE,
scales = list(x = list(log = 10, equispaced.log = FALSE)),
type = c("p", "r"))
No primeiro gráfico desenha-se a linha de regressão com cor mais laranja (col.line = "darkorange") e mais grossa (lwd = 3). E, no gráfico da direita, sobrepõe-se uma linha suavizada de regressão por meio do comando type = c("p", "smooth"), capturando uma mudança sutil na direção da tendência dos dados.
# mudança da cor da reta de regressão
xyplot(lifeExp ~ gdpPercap, gap5Dados, grid = TRUE,
scales = list(x = list(log = 10, equispaced.log = FALSE)),
type = c("p", "r"), col.line = "darkorange", lwd = 3)
# sobreposição de curva suavizada
xyplot(lifeExp ~ gdpPercap, gap5Dados, grid = TRUE,
scales = list(x = list(log = 10, equispaced.log = FALSE)),
type = c("p", "smooth"), col.line = "darkorange", lwd = 3)
O argumento group permite a especificação de variável de agrupamento dos pontos do diagrama de dispersão. Os dados dos vários grupos serão sobrepostos e distinguíveis visualmente. Uma legenda pode ser obtida pelo argumento auto.key = TRUE, embora haja outras alternativas para customizar a legenda.
# agrupamento por continentes
xyplot(lifeExp ~ gdpPercap, gap5Dados, grid = TRUE,
scales = list(x = list(log = 10, equispaced.log = FALSE)),
group = continent)
## legenda
xyplot(lifeExp ~ gdpPercap, gap5Dados, grid = TRUE,
scales = list(x = list(log = 10, equispaced.log = FALSE)),
group = continent, auto.key = TRUE)
O argumento group pode ser usado em conjuntp com o argumento type e modificada a formatação da legenda.
# curva suavizada por grupo (groups + type = "smooth")
xyplot(lifeExp ~ gdpPercap, gap5Dados, grid = TRUE,
scales = list(x = list(log = 10, equispaced.log = FALSE)),
group = continent, auto.key = TRUE,
type = c("p", "smooth"), lwd = 4)
# alterando a legenda
xyplot(lifeExp ~ gdpPercap, gap5Dados, grid = TRUE,
scales = list(x = list(log = 10, equispaced.log = FALSE)),
group = continent,
auto.key = list(columns = nlevels(as.factor(gap5Dados$continent))),
type = c("p", "smooth"), lwd = 4)
A mesma função xyplot() permite também que os diagramas de dispersão dos grupos sejam desenhados em diferentes painéis, ao invés de serem sobrepostos, especificando o fator de condicionamento \(z\) na fórmula do gráfico como em y ~ x | z. Esse tipo de gráfico é denominado condicionamento multipainel, Algumas a variável de condicionamento é a mesma de agrupamento, para se manter o esquema cromático constante e visível.
# multipainéis por continente - mesma cor
xyplot(lifeExp ~ gdpPercap | continent, gap5Dados, grid = TRUE,
scales = list(x = list(log = 10, equispaced.log = FALSE)))
# variável de condicionamento e de agrupamento é o mesmo
xyplot(lifeExp ~ gdpPercap | continent, gap5Dados, grid = TRUE,
group = continent,
scales = list(x = list(log = 10, equispaced.log = FALSE)))
Pode-se combinar o condicionamento, agrupamento e o argumento type.
# condicionamento + type "r" ou "smooth"
xyplot(lifeExp ~ gdpPercap | continent, gap5Dados, grid = TRUE,
scales = list(x = list(log = 10, equispaced.log = FALSE)),
type = c("p", "smooth"), col.line = "darkorange", lwd = 4)
# condicionamento + agrupamento + type "r" ou "smooth"
xyplot(lifeExp ~ gdpPercap | continent, gap5Dados, grid = TRUE,
group = continent,
scales = list(x = list(log = 10, equispaced.log = FALSE)),
type = c("p", "smooth"), lwd = 4)
Em gráficos univariados pode-se combater o overplotting com jitter ou estimando a curva de densidade. Em diagramas de dispersão, a solução é regular a transparência dos pontos com o argumento alpha. No gráfico da direita usa-se alpha = 0.5. Para se interpretar essa característica do gráfico, imagie que se a trasnparência for igual a \(1/k\) então a área torna-se opaca quando são sobrepostos \(k\) pontos. ( then an area becomes opaque when \(k\) points overlap. (Nem todos os dispositivos gráficos suportam o controle da transparência). No gráfico da direita, utiliza-se a função panel.smoothScatter() que estima a densidade bivariada por núcleo.
# Multi-painés com ajuste da transparência dos pontos
xyplot(lifeExp ~ gdpPercap | continent, gap5Dados, grid = TRUE,
scales = list(x = list(log = 10, equispaced.log = FALSE)),
type = c("p", "smooth"), lwd = 4, alpha = 1/2)
# Multipainéis com transparência e suavização da densidade bivariada
xyplot(lifeExp ~ gdpPercap | continent, gap5Dados, grid = TRUE,
scales = list(x = list(log = 10, equispaced.log = FALSE)),
panel = panel.smoothScatter)
A hexagonalização hexagonal binning é uma maneira interessante de lidar com overplotting. O plano é dividido em uma rede hexagonal e os hexágonos individuais são sombreados de acordo com quantidade de pontos que que pertencem a ele. Isso requer o pacote hexbin e com a função hexbinplo() como um substituto para a função xyplot().
library(hexbin)
hexbinplot(lifeExp ~ gdpPercap, gap5Dados,
scales = list(x = list(log = 10, equispaced.log = FALSE)),
aspect = 1, bins=50)
Divirtam-se bastante!!