---
title: "Aula TP 5 -- Resolução FT 3"
author: "João Malato e Tiago A. Marques"
date: "Outubro de 2019"
output:
prettydoc::html_pretty:
theme: cayman
highlight: github
toc: true
toc_depth: 3
---
## Disclaimer
Este documento fornece algumas propostas de resolução aos exercícios da folha de trabalho TP3.pdf, resolvida durante a aula prática da quinta semana (14 a 18 de Outubro de 2019). É de notar que pode haver (quase) infinitas maneiras de solucionar cada um dos exercícios. Se esta resolução não for *exactamente igual* à que foi apresentada na aula, ou à que cada um de vós conseguiu criar, não tem mal. O importante aqui é que sejam capazes de compreender o que fizeram e entender o que é apresentado neste documento para que possam mais à frente reproduzi-lo de forma independente.
Neste documento assumimos que os ficheiros de dados se encontram na mesma directoria do .Rmd. Caso não seja o caso, será necessário primeiro identificar essa directoria.
Além do HTML em si, vamos disponibilizar o .Rmd também para poderem ver como se faz este relatório mais bonitinho do que a versão sem sal que costumam implementar nas aulas. Nota: precisam instalar o package "prettydoc" para isso acontecer!
Feito este preliminar, passemos à resolução.
# Exercício 1.
**Importe os dados `DataTP3densidadescoruja.csv`. Os dados são relativos à densidade de coruja das torres em 4 habitats distintos: montado, zona agrícola, matos, eucaliptal. As densidades expressam o número de indivíduos por $\text{km}^2$.**
**Obtenha a média, a mediana, os quartis, o desvio-padrão, a variância e a distância inter-quartil da densidade para os vários habitats. Caracterize os vários conjuntos de dados (habitats).**
### Importar os dados
> **Pro tip:** Dar um nome curto e facilmente identificável aos dados. Ninguém consegue escrever "*DataTP3densidadescorujas*" rapidamente e sempre que for necessário aceder à data frame. Usem sempre nomes concisos e simples de perceber.
```{r}
# dados separados por ';', logo é possível usar a função read.csv2()
dcoruja <- read.csv2("DataTP3densidadescorujas.csv")
# ou então usa-se read.csv() mas é necessário escrever o tipo de separação no código e os decimais
# dcoruja <- read.csv("DataTP3densidadescorujas.csv", sep = ";", dec = ",")
```
> **Pro tip:** Quando importamos um conjunto de dados novos, vale sempre a pena dar uma vista de olhos quanto à sua estrutura e conjunto de valores (que é basicamente o que estes exercícios pretendem ensinar).
```{r}
# para ver as primeiras linhas da data frame
head(dcoruja)
# para ver a estrutura de cada uma das colunas
str(dcoruja)
# para estudar de um modo geral como se comportam os dados em cada um dos habitats
summary(dcoruja)
```
Através destas funções (pode-se também usar a função `tail()` para ver as últimas linhas da data frame) fica-se com uma visão global sobre o tipo de dados com que estamos a trabalhar e de como poderemos lidar com eles.
### Média
Usa-se a função `mean()` mas como são várias colunas pode-se usar a função `apply()` para optimizar o código.
```{r}
apply(dcoruja, 2, mean)
# para cada parâmetro da função lê-se:
# 1. X usar os dados dcoruja,
# 2. MARGIN aplicar a função FUN a cada coluna (número 2, se fosse o número 1 era linhas),
# 3. FUN função que se pretende aplicar nas linhas/colunas dos dados
```
O uso desta função evita a escrita repetida e algo cansativa (muito *copy-paste*) da função para a média.
```{r, eval=F}
# Como seria chato escrever para se obter a mesmíssima coisa:
mean(dcoruja$montado)
mean(dcoruja$zona.agricola)
mean(dcoruja$matos)
mean(dcoruja$eucaliptal)
```
### Mediana
Obtém-se a mediana ao aplicar a função `median()`.
```{r}
apply(dcoruja, 2, median)
```
### Quartis
Cada um dos quartis calcula-se através da função `quantile()` onde se pode seleccionar o quartil de probabilidade pretendido através do parâmetro `probs` da função (ex.: `quantile(dados$coluna, probs = c(0.25, 0.75))` devolve os quantis de probabilidade 25% e de 75%).
```{r}
# para obter todos os quartis de cada coluna
apply(dcoruja, 2, quantile)
```
Ou obtendo os quantis pedidos para a distância inter-quartil. Para este caso temos de definir o que faz a função `FUN`, no último parâmetro da `apply()`).
```{r}
# opção mais avançada para obter APENAS os quartis de 25% e 75%
apply(dcoruja, 2, function(x) quantile(x, probs = c(0.25, 0.75)))
```
### Desvio-padrão
O desvio padrão calcula-se através da função `sd()`.
```{r}
apply(dcoruja, 2, sd)
```
### Distância interquartil
Pode-se calcular a distância interquartil através da subtracção do quartil de probabilidade de 75% com o quartil de 25% (ver como calcular quartis de probabilidade específicos em cima com a função `quantile()`).
```{r}
# função IQR para interquartile range
apply(dcoruja, 2, quantile, probs = 0.75) - apply(dcoruja, 2, quantile, probs = 0.25)
```
Podemos também aplicar a função `IQR()`, que calcula essa diferença de forma automática.
```{r}
# função IQR para interquartile range
apply(dcoruja, 2, IQR)
```
Poderíamos até definir a nossa propria função para fazer o mesmo que a `IQR()` (util se não soubéssemos que tal função existia). Vamos então criar a função `DIQ()`.
```{r}
DIQ <- function(x) {
diff(quantile(x, prob = c(0.25, 0.75)))
}
```
Usando a função para obter a distância inter-quartil.
```{r}
apply(X = dcoruja, MARGIN = 2, FUN = DIQ)
```
***
# Exercício 2.
**Obtenha os histogramas e gráficos de caixas e bigodes (boxplots) para os vários habitats. Formate os gráficos indicando as respectivas legendas dos eixos e os títulos dos gráficos. O que é que estes gráficos permitem acrescentar às particularidades dos dados descritas em 1.**
### Histogramas
Um histograma padrão para cada coluna pode ser facilmente feito através da função `hist()`. Como mostra o exemplo em baixo, para os Montados.
```{r, fig.align='center', fig.width=8}
hist(dcoruja$montado)
```
Mas como é pedido, é necessário manipular os histogramas de modo a torná-los mais legíveis (alterar eixos, legendas e títulos). Tudo isso pode ser feito através do uso dos parâmetros da função `hist()`.
> **Pro tip:** A grande maioria das funções de visualização (`plot()`, `hist()`, `barplot()`, `boxplot()`, etc) tem uma estrutura de ajuste semelhantes. Aprendam bem quais os parâmetros que mais usam e quais os que mais gostam, para poderem personalizar os gráficos ao vosso bel-prazer sempre que vos for pedido.
```{r, fig.width=8, fig.height=7, fig.align='center'}
# definir uma sequência de quatro cores para se poder usar
colour <- c("#004589", "#0080ff", "#89c4ff", "#d8ebff")
par(mfrow=c(2,2)) # modificar a estrutura dos gráficos para ter 2 linhas e 2 colunas (4 gráficos)
hist(dcoruja$montado, # conjunto de dados
xlab = "Densidade", # label do eixo dos x
ylab = "Frequência", # label do eixo dos y
main = "Montado", # título do histograma
col = colour[1], # cor do histograma
las = 1) # alteração da posição do eixo dos y
hist(dcoruja$zona.agricola, xlab = "Densidade", ylab = "Frequência", main = "Zona Agrícola", col = colour[2], las = 1)
hist(dcoruja$matos, xlab = "Densidade", ylab = "Frequência", main = "Matos", col = colour[3], las = 1)
hist(dcoruja$eucaliptal, xlab = "Densidade", ylab = "Frequência", main = "Eucaliptal", col = colour[4], las = 1)
```
Na realidade podiamos fazer isto tudo de uma vez, com um código um pouco mais avançado, através da aplicação de um ciclo for (*for loop*).
```{r, fig.width=8, fig.height=7, fig.align='center'}
# Alternativa
par(mfrow = c(2,2))
for(i in seq_along(colnames(dcoruja))) {
hist(dcoruja[, i],
xlab = "Densidade",
ylab = "Frequência",
main = colnames(dcoruja)[i],
col = colour[i],
las = 1)
box(lwd = 2) # proposta para alterar a 'line width' (espessura) da caixa do histograma
}
```
### Boxplots
Podemos obter um boxplot dos dados através da função `boxplot()`. Para esta função também se pode seleccionar cada coluna, mas se as colunas estiverem bem difinidas basta indicar-se a data frame que os boxplots surgem categorizados para cada habitat.
```{r, fig.align='center', fig.width=8, fig.height=7}
boxplot(dcoruja)
```
Sendo também necessário alterar a estrutura dos gráficos podemos usar opções semelhantes às que se usou anteriormente nos histogramas.
```{r, fig.align='center', fig.width=8, fig.height=7}
boxplot(dcoruja,
ylim = c(0, 25), # amplitude do eixo dos y
lty = 1, # tipo de linha que os quartis apresentam
las = 1,
col = colour,
names = c("Montado", "Zona Agrícola", "Matos", "Eucaliptal"),
main = "Boxplots para densidades dos habitats",
pch = 19) # tipo de ponto dos outliers
```
**Não se esqueçam:** Os gráficos aqui criados servem para algo. Façam o favor de *olhar* para eles e interpretar algo. O gráfico mais bonito do mundo não nos serve de nada se não conseguirmos obter informação que esclareça as nossas questões. Interpretem e escrevam!
Neste caso, poderiamos dizer que a maior densidade média era encontrada em eucaliptal e montado, seguido de matos, sendo as densidades menores em zonas agricolas. Era ainda de realçar que a variabilidade do Eucaliptal era muito maior que no montado. assim sendo, se por exemplo nos dissessem que queriam maximizar a probabilidade de ir a um local e ver mais do que 6 corujas, qual seria o local a visitar? Isso mesmo, um montado!
Nos matos há duas observações discrepantes, que são muito superiores à média dos outros matos. Poder-se-ia dar o caso de estes matos ficarem perto de um eucaliptal?
### Extra
Na realidade, o ficheiro de dados não estava organizado tal como é recomendado que um ficheiro de excel para ser analisado seja organizado, isto é, tendo cada variável apenas numa coluna -- aqui tinhamos a mesma variável, `densidade`, em 4 colunas diferentes. Mais informação sobre como organizar folhas de cálculo:
Neste caso, para termos o conjunto de dados organizado como é recomendado teriamos de criar uma coluna para a `zona` (com os nomes dos quatro habitats distintos) e outra coluna para os valores de densidades associados a cada observação.
```{r}
# funções usadas
# rep() - repete algo
# names() - devolve o nome das várias colunas
# nrow() - devolve o número de linhas da data frama
newdcoruja <- data.frame(zona = rep(names(dcoruja), each = nrow(dcoruja)),
densidade = c(dcoruja[, 1], dcoruja[, 2], dcoruja[, 3], dcoruja[, 4]))
```
Podemos agora ver o que são os novos dados
```{r}
newdcoruja
```
e por fim usar este novo conjunto de dados para fazer os boxplot, por exemplo.
```{r, fig.align='center', fig.width=8, fig.height=7}
boxplot(newdcoruja$densidade ~ newdcoruja$zona,
col = colour[c(4, 3, 1, 2)],
ylim = c(0, 25),
ylab = "Densidade",
xlab = "Habitats")
```
***
# Exercício 3.
**Importe os dados `DataTP3coleopteros.csv`. Caracterize os dados das diferentes áreas geográficas no que se refere aos indicadores de localização, aos indicadores de dispersão, à distribuição dos dados, a observações discrepantes e a observações influentes. Use um conjunto de estatísticas descritivas e métodos gráficos.**
### Importar os dados
```{r}
dcole <- read.csv2("DataTP3coleopteros.csv")
```
### Análise exploratória
```{r}
head(dcole)
str(dcole)
summary(dcole)
```
### Indicadores de localização
Os indicadores de localização podem ser:
- média
- mediana
```{r}
apply(dcole, 2, mean)
apply(dcole, 2, median)
```
Possíveis questões:
- Como varia a média/mediana entre os três países? Todas diferem? Quais diferem mais?
- Qual o que apresenta os maiores valores? Qual o que apresenta os menores?
- Será que os valores das médias são consistentes com os valores das observações? (*para isto é necessário ir estudar a dispersão dos dados.*)
### Localizadores de dispersão
Os indicadores de dispersão aqui estudados são:
- variância
- desvio padrão
```{r}
apply(dcole, 2, var)
apply(dcole, 2, sd)
```
### Distribuição dos dados, Observações discrepantes e Observações influentes
```{r, fig.align='center', fig.width=8, fig.height=7}
coleoptero.colours <- c("#bcddc7", "#b6d051", "#75be6f")
boxplot(dcole, ylim = c(0, 250), las = 1, col = coleoptero.colours, lty = 1)
```
Portugal é o país que apresenta maior dispersão dos dados, com observações a variarem entre valores de aproximadamente `r round(range(dcole$Portugal)[1], 2)` e `r round(range(dcole$Portugal)[2], 2)`, enquanto que os dados de Marrocos e Espanha são menos dispersos. Marrocos é o único país que apresenta observações discrepantes com uma observação (`r max(dcole$Marrocos)`) muito acima das restantes.
>**Nota**: Estes valores discrepantes poderão ser potenciais outliers -- mas *outlier* é uma palavra muito forte que nem sequer gostamos de usar, pois (1) ou se está a falar de erros; ou (2) de valores que não são consistentes com uma distribuição Normal/Gaussiana, mas porque é que uma variável aleatória ecológica haveria de ser Normal?
Espanha e Marrocos são os países que mais diferem em termos de registo de dados.
```{r include=FALSE}
sdsem <- sd(dcole$Marrocos[-c(5, 10)])
sdcom <- sd(dcole$Marrocos)
```
A presença de observações discrepantes faz variar significativamente o desvio padrão de Marrocos. Retirando os dois pontos extremos.
|Com observações discrepantes|Sem observações discrepantes|
|------------|------------|
|`r sdcom`|`r sdsem`|
A pergunta inicial era traiçoeira. Era perguntado quais eram as **observações discrepantes** e as **influentes**. As discrepantes podemos identificar, como fizémos acima, mas a observações influentes dependem do objectivo. Só se poderia dizer quais observações são influentes depois de se ter feito uma análise, e uma vez que essa análise não foi feita, aqui apenas se podia dizer que as *observações discrepantes tinham o potencial para ser influentes*.
***
# Exercicio 4.
**A partir dos dados `DataTP3coleopteros.csv` relativos a Portugal, aplique as seguintes transformações: raiz quadrada, logarítmo, arco seno e arco seno raiz quadrada. Verifique os efeitos destas transformações através de estatísticas descritivas e de métodos gráficos.**
Sendo que vamos trabalhar só com a coluna referente a Portugal podemos usar essa coluna como variável.
```{r}
colePT <- dcole$Portugal
```
### Raíz quadrada
No R calcula-se a raíz quadrada através da função `sqrt()`.
```{r}
sqrt(colePT)
```
Para estudar visualmente os efeitos da aplicação da raíz quadrada podemos realizar alguns boxplots aos dados originais e compará-los com os dados transformados.
```{r, fig.align='center', fig.width=8, fig.height=7}
# juntar os dois conjuntos de dados através da função c(olumn)bind()
sqrtPT <- cbind(Original = colePT,
Transformado = sqrt(colePT))
# outra alternativa seria criar uma data frame
# sqrtPT <- data.frame(Original = colePT, Transformado = sqrt(colePT))
boxplot(sqrtPT, ylab = "Densidade", main = "Original vs. Raíz Quadrada", col = coleoptero.colours[1], las = 1)
# Para ver os dados transformados numa escala mais ajustada
# boxplot(sqrtPT[, 2], ylab = expression(sqrt(Densidade)), main = "Dados transformados", col = coleoptero.colours[1], las = 1, ylim = c(2,14))
```
É possível ver que a transformação dos dados provocou uma grande diminuição dos valores e o aparecimento de uma observação baixa discrepante (um candidato a outlier).
Podemos estudar a média e o desvio padrão dos dados originais e compará-los com os dados transformados.
```{r}
# juntar os dois conjuntos de dados através da função r(ow)bind()
sqrtPT2 <- rbind(mean = apply(sqrtPT, 2, mean),
sd = apply(sqrtPT, 2, sd))
# para ver a tabela criada
sqrtPT2
# outra alternativa seria criar uma data frame
# sqrtPT2 <- data.frame(apply(sqrtPT, 2, mean),
# apply(sqrtPT, 2, sd),
# row.names = c("Original", "Transformado"))
# colnames(sqrtPT2) <- c("Média", "Desvio Padrão")
```
### Logaritmo
```{r, fig.align='center', fig.width=8, fig.height=7}
# juntar os dois conjuntos de dados através da função cbind()
logPT <- cbind(Original = colePT,
Transformado = log(colePT))
# outra alternativa seria criar uma data frame
# logPT <- data.frame(Original = colePT, Transformado = log(colePT))
boxplot(logPT, ylab = "Densidade", main = "Original vs. Logaritmo", col = coleoptero.colours[1], las = 1)
# Para ver os dados transformados numa escala mais ajustada
# boxplot(logPT[, 2], ylab = "log(Densidade)", main = "Dados transformados", col = coleoptero.colours[1], las = 1, ylim = c(2,6))
```
```{r}
# juntar os dois conjuntos de dados através da função rbind()
logPT2 <- rbind(mean = apply(logPT, 2, mean),
sd = apply(logPT, 2, sd))
# para ver a tabela criada
logPT2
# outra alternativa seria criar uma data frame
# logPT2 <- data.frame(apply(logPT, 2, mean),
# apply(logPT, 2, sd),
# row.names = c("Original", "Transformado"))
# colnames(logPT2) <- c("Média", "Desvio Padrão")
```
### Arco seno
Para se poder trabalhar com a função arco seno (`asin()`) é necessário transformar os dados para que fiquem sob a forma de frequências (entre 0 e 1). Para tal dividimos os valores do objecto `colePT` pelo seu valor máximo (encontramos o máximo através da função `max()`).
Apesar de não ser pedido, podemos também calcular também o seno. Esta transformaçã faz-se através da função `sin()`, ficando assim com duas transformações trigonométricas para comparar.
```{r, fig.align='center', fig.width=8, fig.height=7}
# Qual o valor máximo registado em Portugal?
max(colePT)
# Dividir pelo máximo para obter as frequências
colePT.freq <- colePT/max(colePT)
# criar uma tabela com os dados originais e alterados pelas funções trigonométricas
trigPT <- cbind(Original = colePT,
Seno = sin(colePT.freq),
ArcoSeno = asin(colePT.freq))
# criar dois boxplots:
par(mfrow = c(1, 2)) # dividir o espaço dos gráficos em 1 linha com 2 colunas (dois gráficos lado a lado)
# gráfico para comparar variação das transformações
boxplot(trigPT, ylab = "Densidades", main = "Comparação", col = coleoptero.colours[1], las = 1, lty = 1)
# só alterações trigonométricas (colunas 2 e 3 da tabela trigPT)
boxplot(trigPT[, c(2, 3)], ylim = c(0, 2), main = "Alterações trigonométricas", ylab = "Densidades", col = coleoptero.colours[1], las = 1, lty = 1)
```
> **Pro tip:** Às vezes a escala para um valor que se pretende apresentar não é a ideal para outro valor também apresentado (exemplo em esquerda acima, relativamente às transformações). Nessa altura pode ser necessário arranjar soluções que mostrem melhor os valores obtidos. Para este exemplo podemos simplesmente criar um gráfico (direita) complementar com os valores das transformações trigonométricas.
```{r}
# juntar os dois conjuntos de dados através da função r(ow)bind()
trigPT2 <- rbind(mean = apply(trigPT, 2, mean),
sd = apply(trigPT, 2, sd))
# para ver a tabela criada (arredondada com apenas 3 casas decimais)
round(trigPT2, 3)
```
### Arco seno de raíz quadrada
Nesta situação podemos usar uma função *dentro* de uma função, mais especificamente, calcular a raíz quadrada dos dados observados de Portugal e calcular o arco seno desse valor (`asin(sqrt(dados$coluna))`).
Não esquecer que temos de usar os valores das frequências dos dados (`colePT.freq`) para calcular o arco seno.
```{r, fig.align='center', fig.width=8, fig.height=7}
# podemos calcular funções dentro de funções (perguntem-me sobre o package 'magrittr' mais logo)
asin(sqrt(colePT.freq))
# criando uma tabela
trigsqrtPT <- cbind(Frequencias = colePT.freq,
AsinSqrt = asin(sqrt(colePT.freq)))
# boxplot
boxplot(trigsqrtPT, ylab = "Densidades", main = "Comparação", col = coleoptero.colours[1], las = 1, lty = 1, ylim = c(0,2))
```
E como se comportam as médias e desvios padrão?
```{r}
# juntar as médias e desvios padrão dos dados
trigsqrtPT2 <- rbind(mean = apply(trigsqrtPT, 2, mean),
sd = apply(trigsqrtPT, 2, sd))
# para ver a tabela criada arredondada
round(trigsqrtPT2, 3)
```
***
# Exercício 5.
**Pretende utilizar um procedimento de teste que exige que os dados sejam normais e homocedásticos. Avalia se os dados `DataTP3pesos.csv` cumprem esses requisitos, através de métodos gráficos e de procedimentos de teste.**
```{r}
# ler os dados
pesos <- read.csv2("DataTP3pesos.csv")
```
### Análise exploratória
```{r}
head(pesos)
str(pesos)
summary(pesos)
```
> **Nota:** À medida que vão ficando mais à vontade com o R podem ir extendendo esta rápida análise exploratória inicial. Assim garantem que começam logo a trabalhar dados com o maior conhecimento possível.
### Testar Normalidade
Queremos inferir sobre a normalidade dos dados. Para tal podemos aplicar testes estatísticos para testar uma hipótese nula ($H_0$) para essa distribuição normal (Gaussiana) dos dados. Neste caso, a nossa hipótese nula e consequente hipótese alternativa ($H_1$) são:
$$ H_0: \text{As amostras provêm de uma população com distribuição Gaussiana}\\vs.\\H_1: \text{As amostras não provêm de uma população com distribuição Gaussiana}\ ,$$
Alguns testes de normalidade usados são:
- teste de normalidade Shapiro-Wilk -- `shapiro.test()`
- testes de Kolgomorov-Smirnov -- `ks.test()`
Para estas análises vamos trabalhar com o nível de significância de 5%, isto é, $\alpha = 0.05$.
```{r}
# Shapiro-Wilk normality test
shapiro.test(pesos$dieta.A)
shapiro.test(pesos$dieta.B)
```
Os dados da `dieta.A` apresentam um p-value ('valor-p' em português) menor que o $\alpha$ escolhido ($\text{p-value} < 0.05$) logo temos evidência estatística suficiente para rejeitar a hiópese nula (ou $H_0$). Podemos assim assumir que esta amostra não tem uma distribuição Normal.
Já os dados da `dieta.B` apresentam um p-value maior que o nível de significancia ($\text{p-value} = 0.07944$), logo não temos informação estatística suficiente para rejeitar $H_0$, não existindo assim evidências que a amostra não tenha uma distribuição normal.
> **Pro tip:** Num teste de hipóteses fazemos sempre inferência sobre a hipótese nula $H_0$. Ou seja, dependendo do teste de hióteses, se $\text{p-value} < \alpha$ podemos afirmar que existe evidência estatística suficiente para rejeitar $H_0$, o que por consequinte nos leva a poder considerar a hipótese alternativa proposta, $H_1$. Caso $\text{p-value} \geq \alpha$, já não podemos dizer o mesmo. Neste caso dizemos que não temos evidências para rejeitar $H_0$ sendo essa a que assumimos ser a hipótese mais correcta.
Podemos inferir algo sobre a distribuição dos dados vendo como se comportam graficamente.
```{r, fig.align='center', fig.width=8, fig.height=5}
# Representação gráfica dos dois tipos de dietas
par(mfrow=c(1,2))
hist(pesos$dieta.A, main = "Dieta A", ylab = "Frequência", xlab = "Pesos", col = "cornflowerblue", ylim = c(0, 7), xlim = c(6, 16))
hist(pesos$dieta.B, main = "Dieta B", ylab = "Frequência", xlab = "Pesos", col = "darksalmon", ylim = c(0, 3), xlim = c(0, 20))
```
Estes dados apresentam gráficos que dificilmente serão de uma distribuição normal, dado a existência de classes sem valores e um declive tão abrupto.
Na realidade, o que fizemos acima ignora um facto fundamental sobre o pressuposto da Gaussianidade. É que esse pressuposto *nunca* é sobre os dados em si, mas sobre os residuos de um modelo - são sobre esses que existe um pressuposto de gaussianidade, ou não.
Neste caso, para testar a Gaussianidade, deveriamos fazer um *teste unico*, sobre a totalidade dos resíduos. Assim sendo, o teste correcto seria não sobre os valores originais, mas **sobre os valores originais subtraidos da sua média** (`dados - mean(dados)`), resultando assim nos resíduos associados a um teste de t para duas amostras.
```{r}
residuos <- c(pesos$dieta.A - mean(pesos$dieta.A),
pesos$dieta.B - mean(pesos$dieta.B))
shapiro.test(residuos)
```
Podemos comprovar a mesma conclusão com o teste de Kolmogorov-Smirnov:
```{r, warning=F}
ks.test(residuos, "pnorm")
```
Através dos dois testes vamos rejeitar a hipótese nula para os níveis usuais de significância ($\text{testes de Shapiro e K.S.: p-value} < 0.05$).
### Testar Homocedastidade
Queremos estudar a homocedasticidade dos dados (*homoskedasticity* em ingês). Isso corresponde a testar a igualdade de variâncias das variáveis. Assim, as nossas hipóteses são:
$$ H_0: \text{As amostras provêm de populações com variâncias iguais (homocedasticidade)}\\vs.\\H_1: \text{As amostras provêm de populações com variâncias diferentes (heterocedasticidade)}\ ,$$
Podemos aplicar o teste específico para homogeneidade de variâncias, o teste de Bartlett (através da função `bartlett.test()`).
Considerando o nível de significância ainda 5% ($\alpha = 0.05$),
```{r}
# Barlett test of homogeneity of variances
bartlett.test(x = c(pesos$dieta.A, pesos$dieta.B), # vector com os dados
g = rep(c("A", "B"), each = 12)) # como agrupamos os elementos
```
Como o valor-p é maior que o nível de significância ($\text{p-value} = 0.08502$) não existe evidência estatística que nos leve à rejeição de $H_0$, isto é, não podemos rejeitar a hipótese de que variâncias seja igual.
***
# Exercício 6.
**Aplique as transformações que achar convenientes aos dados analisados na alínea anterior e verifique o cumprimento dos pressupostos de normalidade e homogeneidade de variâncias, através de testes e de gráficos. O que seria melhor efectuar nesta situação?**
A transformação de um conjunto de dados pode servir para obter dados que respeitem determinados pressupostos, mas também nos permite estudar e inferir sobre determinadas características estatísticas que possam não ser observáveis ou estar "disfarçadas" na sua forma original.
Existem várias transformações possíveis, estas aqui demonstradas são das mais usuais.
Para os testes de hipóteses realizados vamos assumir um nível de significância de 10% ($\alpha = 0.10$).
### Transformação raíz quadrada
Implementando a transformação e vendo o seu efeito.
```{r}
sqrtpesos <- sqrt(pesos)
summary(sqrtpesos)
```
Refazendo o teste sobre a Gaussianidade dos residuos, claculamos os resíduos com a raiz quadrada
```{r}
res_sqrtpesos <- c(sqrtpesos$dieta.A - mean(sqrtpesos$dieta.A),
sqrtpesos$dieta.B - mean(sqrtpesos$dieta.B))
```
```{r}
# Shapiro-Wilk test para as duas dietas
shapiro.test(res_sqrtpesos)
```
É possível ver que esta transformação não alterou o resultado do teste da normalidade.
```{r, fig.align='center', fig.width=8, fig.height=5}
# Representação gráfica
hist(res_sqrtpesos, main = "", ylim = c(0, 12), ylab = "Frequência", xlab = expression(sqrt(Pesos)), col = "lightcoral", las = 1)
```
### Transformação logarítmica
```{r}
logpesos <- log(pesos)
summary(logpesos)
res_logpesos <- c(logpesos$dieta.A - mean(logpesos$dieta.A),
logpesos$dieta.B - mean(logpesos$dieta.B))
```
Testando a normalidade dos dados:
$$ H_0: \text{A população de onde veio a amostra tem uma distribuição Gaussiana.}$$
```{r}
# Shapiro-Wilk test para as duas dietas
shapiro.test(res_logpesos)
```
Com o resultado do teste de normalidade, não rejeitamos a hipótese nula $H_0$ de que a amostra provém de uma distribuição Gaussiana.
```{r, fig.align='center', fig.width=8, fig.height=5}
# Representação gráfica
hist(res_logpesos, main = "", ylab = "Frequência", xlab = expression(log(Pesos)), col = "lightcoral", las = 1, ylim = c(0, 10), xlim = c(-1, 1))
```
### Transformação de standadização
A standardização faz-se através da função `scale()`, que vai calcular a equação
$$Z = \frac{(x-\bar{x})}{\sigma_x}\ ,$$
onde $\bar x$ e $\sigma_x$ representam respectivamente a média e o desvio padrão dos dados $x$.
```{r}
# aplicação da função scale()
scale(pesos)
# tratar os dados standardizados como uma data frame (função as.data.frame())
stpesos <- as.data.frame(scale(pesos))
summary(stpesos)
res_stpesos <- c(stpesos$dieta.A - mean(stpesos$dieta.A),
stpesos$dieta.B - mean(stpesos$dieta.B))
```
```{r}
# Shapiro-Wilk Normality Test para as dietas
shapiro.test(res_stpesos)
```
Esta transformação não alterou o resultado do teste da normalidade que temos vindo a obter. A amostra da dieta.A ($\text{p-value} < 0.05$) continua a garantir a rejeição da hipótese nula para normalidade enquanto a amostra da dieta.B ($\text{p-value} = 0.07944$) continua sem evidências que possam afirmar que não é normal.
```{r, fig.align='center', fig.width=9, fig.height=5}
# Representação gráfica
hist(res_stpesos, main = "", ylab = "Frequência", xlab = "Pesos", col = "lightcoral", las = 1, ylim = c(0, 12))
box()
```
Por fim, revejam os resultados, os testes de hipóteses e realizem as interpretações que acharem adequadas.
> **Nota:** Depois de todas estas transformações, apenas a logaritmização alterou as conclusões referentes aos dados obtidas no exercício anterior. Mas todo o processo não foi em vão! Às vezes acontece passarmos muito tempo às voltas com os dados, e faz parte da aprendizagem tentar todas as ferramentas à nossa disposição para garantir que os dados se comportam de uma maneira ou de outra.
> **Pro tip:** Enquanto biólogos e conhecedores da biologia e da relação dos dados com o mundo real -- mais que a maioria de quem faz este tipo de análises --, podem e devem aplicar a vossa experiência do curso nas análises e tentar ir mais além da mera apresentação dos resultados. Perguntem-se sobre o porquê dos resultados obtidos serem como são e quais as suas implicações biológicas. E por favor, comentem o vosso código!
Qualquer questão relativamente a esta resolução poderá ser discutida com qualquer um de nós.