Os passos contidos nesse tutorial foram retirados da seguinte página:
[https://rstudio-pubs-static.s3.amazonaws.com/574220_2ddd7eb0d872432886602576d6f9801a.html]
Os dados consistem de 50 unidades amostrais de três espécies (setosa, virginica, versicolor) de íris (uma espécie de planta). Temos então um total de 150 unidades amostrais.
head(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
summary(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100
## 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300
## Median :5.800 Median :3.000 Median :4.350 Median :1.300
## Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
## 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
## Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
## Species
## setosa :50
## versicolor:50
## virginica :50
##
##
##
Se forem observar a função head informa os primeiros valores da tabela de dados. Percebe-se que temos 4 variáveis que foram analisadas:
Vamos explorar e conhecer um pouco mais sobre esses dados?
table(iris$Species)
##
## setosa versicolor virginica
## 50 50 50
Podemos verificar que temos 50 amostras por espécie da planta. Isso foi apresentado no início desse estudo, mas agora com a função “table” podemos comprovar na prática.
Vamos investigar um pouco mais?? Inicialmente com a fução “summary” obtivemos a estatística descritiva de cada variável independentemente de cada espécie. Agora vamos obter os dados de maneira mais detalhada.
tapply(X = iris$Sepal.Length, INDEX = list(iris$Species), FUN = mean)
## setosa versicolor virginica
## 5.006 5.936 6.588
Utilizamos acima a função “tapply”. Vamos agora proceder o mesmo cálculo, porém com outra função. Veja como fica:
aggregate(Sepal.Length ~ Species, data = iris, mean)
## Species Sepal.Length
## 1 setosa 5.006
## 2 versicolor 5.936
## 3 virginica 6.588
Vamos agora utilizar esse mesmo comando para o cálculo da média das outras variáveis:
aggregate(Sepal.Length ~ Species, data = iris, mean)
## Species Sepal.Length
## 1 setosa 5.006
## 2 versicolor 5.936
## 3 virginica 6.588
aggregate(Sepal.Width ~ Species, data = iris, mean)
## Species Sepal.Width
## 1 setosa 3.428
## 2 versicolor 2.770
## 3 virginica 2.974
aggregate(Petal.Length ~ Species, data = iris, mean)
## Species Petal.Length
## 1 setosa 1.462
## 2 versicolor 4.260
## 3 virginica 5.552
Vamos agora utilizar dessa mesma função “aggregate” para calcular o desvio padrão:
aggregate(Sepal.Length ~ Species, data = iris, sd)
## Species Sepal.Length
## 1 setosa 0.3524897
## 2 versicolor 0.5161711
## 3 virginica 0.6358796
aggregate(Sepal.Width ~ Species, data = iris, sd)
## Species Sepal.Width
## 1 setosa 0.3790644
## 2 versicolor 0.3137983
## 3 virginica 0.3224966
aggregate(Petal.Length ~ Species, data = iris, sd)
## Species Petal.Length
## 1 setosa 0.1736640
## 2 versicolor 0.4699110
## 3 virginica 0.5518947
Interessante né? Mas vamos tentar otimizar e colocar esses resultados de médias e desvio padrão numa tabela só:
medias <- matrix(NA, ncol = 3, nrow = 4)
colnames(medias) <- unique(iris$Species)
rownames(medias) <- names(iris)[-5]
for (i in 1:4){medias[i,] <- tapply(iris[,i], iris$Species, mean)}
medias
## setosa versicolor virginica
## Sepal.Length 5.006 5.936 6.588
## Sepal.Width 3.428 2.770 2.974
## Petal.Length 1.462 4.260 5.552
## Petal.Width 0.246 1.326 2.026
Agora que aprendemos alguns conceitos básicos de manipulação no R, vamos falar da Estatística descritiva.
Vamos falar um pouco sobre:
| Parâmetro | Descrição | Função no R |
|---|---|---|
| Média | Média aritmética | mean() |
| Mediana | valor central | median() |
| Moda | Valor mais frequente | sort(table(), decreasing = TRUE)[1] |
| Desvio padrão | variação em torno da média | sd() |
| Quantis | pontos de corte dividindo uma distribuição de probabilidade | quantile() |
A média é um dos fundamentos para obtenção de medidas de tendência central.
\(\overline{x}=\frac{\sum_{i=1}^{n}\times _{i}}{n}\)
vars <- iris[, -5]
apply(vars, 2, mean)
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 5.843333 3.057333 3.758000 1.199333
Vamos ao cálculo da mediana. Como dito anteriormente, é o valor central.
apply(vars, 2, median)
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 5.80 3.00 4.35 1.30
Calculando agora a moda da minha amostra:
moda <- sort(table(iris$Petal.Length), decreasing = TRUE)
moda[1]
## 1.4
## 13
Agora que vimos alguns cálculos básicos de medidas de tendência central é necessário verificar quanto desses valores estão espalhados em torno da média. Isso chama-se na linguagem estatística de medidas de dispersão.
A variância mede o quanto os dados estão dispersos em torno da média, ou seja, do valor esperado.
Monitorar a variância é essencial para as indústrias de manufatura e qualidade porque a redução da variância dos processos aumenta a precisão e diminui o número de defeitos dos produtos manufaturados.
\(s^{2}=\sum_{i=1}^{n}\frac{(x_{i}-x)^{2}}{n-1}\)
apply(vars, 2, var)
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.6856935 0.1899794 3.1162779 0.5810063
O desvio padrão é uma medida que indica a dispersão dos dados dentro de uma amostra com relação à média.
\(s=\sqrt{\sum_{i=1}^{n}\frac{(x_{i}-x)^{2}}{n-1}}\)
sd01 <- apply(vars, 2, sd)
sd01
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.8280661 0.4358663 1.7652982 0.7622377
Podemos avaliar a correlação entre as variáveis.
Tipos de Correlação
cor(vars)
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 1.0000000 -0.1175698 0.8717538 0.8179411
## Sepal.Width -0.1175698 1.0000000 -0.4284401 -0.3661259
## Petal.Length 0.8717538 -0.4284401 1.0000000 0.9628654
## Petal.Width 0.8179411 -0.3661259 0.9628654 1.0000000
Vamos agora falar dos métodos gráficos:
| Tipo de gráfico | Descrição |
|---|---|
| barras | cada barra representa a frequência de observações de um grupo |
| histograma | cada barra representa a frequência de observalçoes em um dado conjunto de valores |
| densidade | estimativa da distribuição estatística baseada nos dados |
| quantil-quantil | comparação dos dados em relação a uma distribuição normal |
| box-plot | representação visual da média, quartis, simetria e outliers |
| dispersão | display gráfico de variáveis contínuas nos eixos x e y |
O R em seu pacote base possui funções para execução dos gráficos acima. Contudo, pode ser instalado o pacote “ggplot2” que é uma ferramenta poderosa para ambiente gráfico.
Na tabela abaixo segue um compilado das funções do Rbase e do pacote ggplot2 para obtenção dos gráficos apresentados anteriormente.
| Tipo de gráfico | Rbase | ggplot2 |
|---|---|---|
| barras | barplot() | geom_bar() |
| histograma | histogram() | geom_histogram() |
| densidade | plot(density()) | geom_density() |
| quantil-quantil | qqnorm() | geom_qq() |
| box-plot | boxplot() | geom_boxplot() |
| dispersão | plot() | geom_point() |
barplot(table(iris$Species))
barplot(table(iris$Sepal.Width))
par(mfrow=c(2, 2))
hist(iris$Sepal.Length)
hist(iris$Sepal.Width)
hist(iris$Petal.Length)
hist(iris$Petal.Length)
par(mfrow=c(1, 2))
hist(iris$Sepal.Width)
hist(iris$Sepal.Width, breaks = 4)
A curva de densidade mostra a probabilidade de observar determinado valor. Em comparação ao histograma, no eixo y, ao invés de termos a frequência, temos a densidade probabilística.
par(mfrow=c(1, 2))
hist(iris$Sepal.Width)
hist(iris$Sepal.Width, freq = FALSE)
Vamos utilizar a função density.
par(mfrow=c(1, 2))
# plot da curva de densidade
plot(density(iris$Sepal.Width))
# plot da curva de densidade sobre o histograma de densidade
hist(iris$Sepal.Width, freq = FALSE)
lines(density(iris$Sepal.Width), col="blue")
O boxplot nos fornece uma análise visual da posição, dispersão, simetria, caudas e valores discrepantes (outliers) do conjunto de dados.
Posição: Em relação à posição dos dados, observa-se a linha central do retângulo (a mediana ou segundo quartil).
Dispersão: A dispersão dos dados pode ser representada pelo intervalo interquatilico que é a diferença entre o terceiro quartil e o primeiro quartil (tamanho da caixa), ou ainda pela amplitude que é calculada da seguinte maneira: valor máximo – valor mínimo. Embora a amplitude seja de fácil entendimento, o intervalo interqualitico é uma estatística mais robusta para medir variabilidade uma vez que não sofre influência de outliers.
Simetria: Um conjunto de dados que tem uma distribuição simétrica, terá a linha da mediana no centro do retângulo. Quando a linha da mediana está próxima ao primeiro quartil, os dados são assimétricos positivos e quando a posição da linha da mediana é próxima ao terceiro quartil, os dados são assimétricos negativos. Vale ressaltar que a mediana é a medida de tendência central mais indicada quando os dados possuem distribuição assimétrica, uma vez que a média aritmética é influenciada pelos valores extremos.
Caudas: As linhas que vão do retângulo até aos outliers podem fornecer o comprimento das caudas da distribuição.
Outliers: Já os outliers indicam possíveis valores discrepantes. No boxplot, as observações são consideradas outliers quando estão abaixo ou acima do limite de detecção de outliers. O limite de detecção de outliers é construído utilizando o intervalo interquartílico, dado pela distância entre o primeiro e o terceiro quartil. Sendo assim, os limites inferior e superior de detecção de outlier são dados por:
Limite Inferior = Primeiro Quartil – 1,5 * (Terceiro Quartil – Primeiro Quartil)
Limite Superior = Terceiro Quartil + 1,5 * (Terceiro Quartil – Primeiro Quartil)
[Oliveira, Bruno. “O que é Desvio Padrão e Erro Padrão? - Blog da Oper”. Oper, 14 de agosto de 2019, https://operdata.com.br/blog/desvio-padrao-e-erro-padrao/.]
Vamos elaborar os box-plots de nossas variáveis:
boxplot(iris$Sepal.Length)
boxplot(iris$Sepal.Width)
boxplot(iris$Petal.Length)
boxplot(iris$Petal.Width)
Vamos rodar por espécie:
boxplot(Sepal.Length ~ Species, data = iris)
boxplot(Sepal.Width ~ Species, data = iris)
boxplot(Petal.Length ~ Species, data = iris)
boxplot(Petal.Width ~ Species, data = iris)
Podemos perceber a existência de alguns outliers. Vamos usar a própria função boxplot para indentificá-los.
boxplot(iris$Sepal.Width)
my_boxplot <- boxplot(iris$Sepal.Width, plot=FALSE)
my_boxplot
## $stats
## [,1]
## [1,] 2.2
## [2,] 2.8
## [3,] 3.0
## [4,] 3.3
## [5,] 4.0
##
## $n
## [1] 150
##
## $conf
## [,1]
## [1,] 2.935497
## [2,] 3.064503
##
## $out
## [1] 4.4 4.1 4.2 2.0
##
## $group
## [1] 1 1 1 1
##
## $names
## [1] "1"
# o objeto é uma lista e os valores outliers estão guardados no elemento $out da lista
outliers <- my_boxplot$out
#qual a posicao dos outliers
which(iris$Sepal.Width %in% outliers)
## [1] 16 33 34 61
# vamos usar a posicao para indexar o objeto
iris[which(iris$Sepal.Width %in% outliers), c("Sepal.Width", "Species")]
## Sepal.Width Species
## 16 4.4 setosa
## 33 4.1 setosa
## 34 4.2 setosa
## 61 2.0 versicolor
No caso anterior consideramos outliers em relação à distribuição da variável para todas as espécies juntas. É razoável assumir que cada espécie tenha um padrão morfométrico distinto de modo que poderíamos identificar outliers de maneira específica por espécie.
boxplot(Sepal.Width ~ Species, data = iris)
my_boxplot2 <- boxplot(Sepal.Width ~ Species, data=iris, plot=FALSE)
my_boxplot2
## $stats
## [,1] [,2] [,3]
## [1,] 2.9 2.0 2.2
## [2,] 3.2 2.5 2.8
## [3,] 3.4 2.8 3.0
## [4,] 3.7 3.0 3.2
## [5,] 4.4 3.4 3.8
##
## $n
## [1] 50 50 50
##
## $conf
## [,1] [,2] [,3]
## [1,] 3.288277 2.688277 2.910622
## [2,] 3.511723 2.911723 3.089378
##
## $out
## [1] 2.3
##
## $group
## [1] 1
##
## $names
## [1] "setosa" "versicolor" "virginica"
# o objeto é uma lista e os valores outliers estão guardados no elemento $out da lista
outliers2 <- my_boxplot2$out
# neste caso, queremos apenas os outliers da especie setosa
# vamos usar a posicao para indexar o objeto
iris[iris$Sepal.Width %in% outliers2 &
iris$Species == "setosa",
c("Sepal.Width", "Species")]
## Sepal.Width Species
## 42 2.3 setosa
A análise da distribuição dos dados é fundamental. Os testes paramétricos por exemplo só podem ser aplicados em dados que apresentem uma distribuição normal. Observe a tabela abaixo com algumas distribuições:
| Distribuição | Tipo | E (X) | \(\sigma^{2}(X)\) | Uso | Exemplo |
|---|---|---|---|---|---|
| normal | contínua | \(\sigma^{2}\) | A distribuição normal é simétrica em torno da média o que implica que e média, a mediana e a moda são todas coincidentes. | distribuição de tamanho | |
| binomial | discreta | \(np\) | \(np(1-p)\) | número de sucessos em \(n\) tentativas | presença ou ausência de espécies |
| Poisson | discreta | \(\lambda\) | \(\lambda\) | modelar o número de ocorrências de um evento por um certo período de tempo ou por um certo volume ou por uma certa área. | número de chegas em um estabelecimento bancário. |
| Log-normal | contínua | \(log(\mu)\) | \(log(\sigma^{2})\) | curva assimétrica | Distribuição de abundância de espécies |
Vamos verificar a normalidade através dos seguintes gráficos:
par(mfrow = c(1,3))
qqnorm(iris$Sepal.Length[iris$Species == "setosa"],
main = "setosa")
qqline(iris$Sepal.Length[iris$Species == "setosa"])
qqnorm(iris$Sepal.Length[iris$Species == "versicolor"],
main = "versicolor")
qqline(iris$Sepal.Length[iris$Species == "versicolor"])
qqnorm(iris$Sepal.Length[iris$Species == "virginica"],
main = "virginica")
qqline(iris$Sepal.Length[iris$Species == "virginica"])
Esse gráfico mostra as correelações anteriormente calculadas. Podemos observar uma forte correlação entre as variáveis Petalwidth e Petallenght.
pairs(vars)
library("GGally")
## Loading required package: ggplot2
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
ggpairs(vars)