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.
# 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).
## montado zona.agricola matos eucaliptal
## 1 13.2 2.1 2.1 1.4
## 2 11.1 0.1 1.8 16.7
## 3 9.8 0.6 1.9 12.2
## 4 10.7 0.5 2.5 23.6
## 5 7.5 0.2 2.2 2.9
## 6 15.2 0.5 19.2 18.7
## 'data.frame': 11 obs. of 4 variables:
## $ montado : num 13.2 11.1 9.8 10.7 7.5 15.2 10.2 11.5 16.1 9.8 ...
## $ zona.agricola: num 2.1 0.1 0.6 0.5 0.2 0.5 1.8 1.1 1.4 2.1 ...
## $ matos : num 2.1 1.8 1.9 2.5 2.2 19.2 2.4 1.7 2.7 0.9 ...
## $ eucaliptal : num 1.4 16.7 12.2 23.6 2.9 18.7 0.2 10.2 19.1 3.2 ...
## montado zona.agricola matos eucaliptal
## Min. : 7.50 Min. :0.1000 Min. : 0.90 Min. : 0.20
## 1st Qu.:10.00 1st Qu.:0.3500 1st Qu.: 1.85 1st Qu.: 2.70
## Median :10.90 Median :0.6000 Median : 2.20 Median :10.20
## Mean :11.45 Mean :0.9545 Mean : 5.50 Mean :10.06
## 3rd Qu.:12.35 3rd Qu.:1.6000 3rd Qu.: 2.60 3rd Qu.:17.70
## Max. :16.10 Max. :2.1000 Max. :23.10 Max. :23.60
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.
## montado zona.agricola matos eucaliptal
## 11.4545455 0.9545455 5.5000000 10.0636364
# 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.
Mediana
Obtém-se a mediana ao aplicar a função median()
.
## montado zona.agricola matos eucaliptal
## 10.9 0.6 2.2 10.2
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%).
## montado zona.agricola matos eucaliptal
## 0% 7.50 0.10 0.90 0.2
## 25% 10.00 0.35 1.85 2.7
## 50% 10.90 0.60 2.20 10.2
## 75% 12.35 1.60 2.60 17.7
## 100% 16.10 2.10 23.10 23.6
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()
).
# 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)))
## montado zona.agricola matos eucaliptal
## 25% 10.00 0.35 1.85 2.7
## 75% 12.35 1.60 2.60 17.7
Desvio-padrão
O desvio padrão calcula-se através da função sd()
.
## montado zona.agricola matos eucaliptal
## 2.496944 0.782769 7.801282 8.462183
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()
).
# função IQR para interquartile range
apply(dcoruja, 2, quantile, probs = 0.75) - apply(dcoruja, 2, quantile, probs = 0.25)
## montado zona.agricola matos eucaliptal
## 2.35 1.25 0.75 15.00
Podemos também aplicar a função IQR()
, que calcula essa diferença de forma automática.
## montado zona.agricola matos eucaliptal
## 2.35 1.25 0.75 15.00
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()
.
Usando a função para obter a distância inter-quartil.
## montado zona.agricola matos eucaliptal
## 2.35 1.25 0.75 15.00
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.
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.
# 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).
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.
Sendo também necessário alterar a estrutura dos gráficos podemos usar opções semelhantes às que se usou anteriormente nos histogramas.
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.
# 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
## zona densidade
## 1 montado 13.2
## 2 montado 11.1
## 3 montado 9.8
## 4 montado 10.7
## 5 montado 7.5
## 6 montado 15.2
## 7 montado 10.2
## 8 montado 11.5
## 9 montado 16.1
## 10 montado 9.8
## 11 montado 10.9
## 12 zona.agricola 2.1
## 13 zona.agricola 0.1
## 14 zona.agricola 0.6
## 15 zona.agricola 0.5
## 16 zona.agricola 0.2
## 17 zona.agricola 0.5
## 18 zona.agricola 1.8
## 19 zona.agricola 1.1
## 20 zona.agricola 1.4
## 21 zona.agricola 2.1
## 22 zona.agricola 0.1
## 23 matos 2.1
## 24 matos 1.8
## 25 matos 1.9
## 26 matos 2.5
## 27 matos 2.2
## 28 matos 19.2
## 29 matos 2.4
## 30 matos 1.7
## 31 matos 2.7
## 32 matos 0.9
## 33 matos 23.1
## 34 eucaliptal 1.4
## 35 eucaliptal 16.7
## 36 eucaliptal 12.2
## 37 eucaliptal 23.6
## 38 eucaliptal 2.9
## 39 eucaliptal 18.7
## 40 eucaliptal 0.2
## 41 eucaliptal 10.2
## 42 eucaliptal 19.1
## 43 eucaliptal 3.2
## 44 eucaliptal 2.5
e por fim usar este novo conjunto de dados para fazer os boxplot, por exemplo.
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
Análise exploratória
## Portugal Espanha Marrocos
## 1 146.62775 187.1135 18.07499
## 2 136.67867 194.4733 21.80148
## 3 93.69134 193.7110 26.39281
## 4 57.54435 185.1346 20.22935
## 5 184.87033 175.7776 228.20000
## 6 180.22817 159.7605 18.33291
## 'data.frame': 11 obs. of 3 variables:
## $ Portugal: num 146.6 136.7 93.7 57.5 184.9 ...
## $ Espanha : num 187 194 194 185 176 ...
## $ Marrocos: num 18.1 21.8 26.4 20.2 228.2 ...
## Portugal Espanha Marrocos
## Min. : 10.64 Min. :150.3 Min. : 1.234
## 1st Qu.: 92.91 1st Qu.:167.2 1st Qu.: 19.281
## Median :124.79 Median :185.1 Median : 26.393
## Mean :120.94 Mean :178.0 Mean : 41.204
## 3rd Qu.:163.43 3rd Qu.:190.4 3rd Qu.: 28.805
## Max. :193.21 Max. :199.9 Max. :228.200
Indicadores de localização
Os indicadores de localização podem ser:
- média
- mediana
## Portugal Espanha Marrocos
## 120.93767 177.95293 41.20379
## Portugal Espanha Marrocos
## 124.78882 185.13458 26.39281
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
## Portugal Espanha Marrocos
## 3162.7172 308.3572 3921.9483
## Portugal Espanha Marrocos
## 56.23804 17.56010 62.62546
Distribuição dos dados, Observações discrepantes e Observações influentes
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 10.64 e 193.21, 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 (228.2) 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.
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 |
---|---|
62.6254608 | 5.6175171 |
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.
Raíz quadrada
No R calcula-se a raíz quadrada através da função sqrt()
.
## [1] 12.108995 11.690965 9.679429 7.585799 13.596703 13.424909 11.170892
## [8] 13.899852 10.484098 9.598289 3.261227
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.
# 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.
# 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
## Original Transformado
## mean 120.93767 10.591014
## sd 56.23804 3.105623
Logaritmo
# 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))
# 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
## Original Transformado
## mean 120.93767 4.5990405
## sd 56.23804 0.8245156
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.
## [1] 193.2059
# 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.
# 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)
## Original Seno ArcoSeno
## mean 120.938 0.565 0.760
## sd 56.238 0.240 0.447
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.
# podemos calcular funções dentro de funções (perguntem-me sobre o package 'magrittr' mais logo)
asin(sqrt(colePT.freq))
## [1] 1.0575599 0.9992871 0.7703260 0.5772801 1.3615633 1.3086311 0.9334368
## [8] 1.5707963 0.8545259 0.7622239 0.2368309
# 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?
# 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)
## Frequencias AsinSqrt
## mean 0.626 0.948
## sd 0.291 0.377
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.
Análise exploratória
## dieta.A dieta.B
## 1 8.2 3.4
## 2 9.2 5.2
## 3 7.5 7.8
## 4 7.4 9.3
## 5 7.2 12.5
## 6 13.2 3.6
## 'data.frame': 12 obs. of 2 variables:
## $ dieta.A: num 8.2 9.2 7.5 7.4 7.2 13.2 14.1 7.1 8.1 7.9 ...
## $ dieta.B: num 3.4 5.2 7.8 9.3 12.5 3.6 4.7 9.8 16.4 3.8 ...
## dieta.A dieta.B
## Min. : 7.100 Min. : 3.400
## 1st Qu.: 7.475 1st Qu.: 4.100
## Median : 7.750 Median : 6.200
## Mean : 8.750 Mean : 7.325
## 3rd Qu.: 8.450 3rd Qu.: 9.425
## Max. :14.100 Max. :16.400
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\).
##
## Shapiro-Wilk normality test
##
## data: pesos$dieta.A
## W = 0.67165, p-value = 0.0004555
##
## Shapiro-Wilk normality test
##
## data: pesos$dieta.B
## W = 0.87666, p-value = 0.07944
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.
# 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.
residuos <- c(pesos$dieta.A - mean(pesos$dieta.A),
pesos$dieta.B - mean(pesos$dieta.B))
shapiro.test(residuos)
##
## Shapiro-Wilk normality test
##
## data: residuos
## W = 0.88157, p-value = 0.008938
Podemos comprovar a mesma conclusão com o teste de Kolmogorov-Smirnov:
##
## One-sample Kolmogorov-Smirnov test
##
## data: residuos
## D = 0.37493, p-value = 0.002348
## alternative hypothesis: two-sided
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\)),
# 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
##
## Bartlett test of homogeneity of variances
##
## data: c(pesos$dieta.A, pesos$dieta.B) and rep(c("A", "B"), each = 12)
## Bartlett's K-squared = 2.9661, df = 1, p-value = 0.08502
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.
## dieta.A dieta.B
## Min. :2.665 Min. :1.844
## 1st Qu.:2.734 1st Qu.:2.024
## Median :2.784 Median :2.482
## Mean :2.937 Mean :2.619
## 3rd Qu.:2.906 3rd Qu.:3.070
## Max. :3.755 Max. :4.050
Refazendo o teste sobre a Gaussianidade dos residuos, claculamos os resíduos com a raiz quadrada
res_sqrtpesos <- c(sqrtpesos$dieta.A - mean(sqrtpesos$dieta.A),
sqrtpesos$dieta.B - mean(sqrtpesos$dieta.B))
##
## Shapiro-Wilk normality test
##
## data: res_sqrtpesos
## W = 0.92937, p-value = 0.09435
É possível ver que esta transformação não alterou o resultado do teste da normalidade.
Transformação logarítmica
## dieta.A dieta.B
## Min. :1.960 Min. :1.224
## 1st Qu.:2.012 1st Qu.:1.410
## Median :2.048 Median :1.811
## Mean :2.142 Mean :1.861
## 3rd Qu.:2.133 3rd Qu.:2.243
## Max. :2.646 Max. :2.797
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.}\]
##
## Shapiro-Wilk normality test
##
## data: res_logpesos
## W = 0.95865, p-value = 0.4119
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.
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\).
## dieta.A dieta.B
## [1,] -0.2326631 -0.96419198
## [2,] 0.1903607 -0.52201477
## [3,] -0.5287797 0.11668565
## [4,] -0.5710820 0.48516667
## [5,] -0.6556868 1.27125949
## [6,] 1.8824556 -0.91506118
## [7,] 2.2631770 -0.64484177
## [8,] -0.6979892 0.60799367
## [9,] -0.2749654 2.22931012
## [10,] -0.3595702 -0.86593038
## [11,] -0.5287797 -0.76766878
## [12,] -0.4864773 -0.03070675
## attr(,"scaled:center")
## dieta.A dieta.B
## 8.750 7.325
## attr(,"scaled:scale")
## dieta.A dieta.B
## 2.363934 4.070766
# tratar os dados standardizados como uma data frame (função as.data.frame())
stpesos <- as.data.frame(scale(pesos))
summary(stpesos)
## dieta.A dieta.B
## Min. :-0.6980 Min. :-0.9642
## 1st Qu.:-0.5394 1st Qu.:-0.7922
## Median :-0.4230 Median :-0.2764
## Mean : 0.0000 Mean : 0.0000
## 3rd Qu.:-0.1269 3rd Qu.: 0.5159
## Max. : 2.2632 Max. : 2.2293
##
## Shapiro-Wilk normality test
##
## data: res_stpesos
## W = 0.80585, p-value = 0.0003677
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.
# 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.