Workshop análises básicas no R

Raquel Divieso

2023-07-10

Agora que já temos uma base sobre manipulação de dados no R, vamos ver algumas análises de dados, normalmente utilizadas em estudos ecológicos.

  • Análises de variância;
  • Correlações
  • Regressões
library(datasets)
library(dplyr)
library(tidyr)
library(multcompView)
library(viridis)

Iris é um conjunto de dados morfológicos de três espécies de plantas do gênero Iris disponível no R muito acessível para testes e treinamentos.

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

Para um conjunto de dados com várias espécies (ou outro tipo de grupos) a função agragate facilita o cáculo de diferentes estatísticas

sep_mean<-aggregate(iris$Sepal.Length, list(iris$Species), FUN=mean) # Média
sep_mean
##      Group.1     x
## 1     setosa 5.006
## 2 versicolor 5.936
## 3  virginica 6.588
sep_var<-aggregate(iris$Sepal.Length, list(iris$Species), FUN=var) # Variância
sep_var
##      Group.1         x
## 1     setosa 0.1242490
## 2 versicolor 0.2664327
## 3  virginica 0.4043429
sep_sd<-aggregate(iris$Sepal.Length, list(iris$Species), FUN=sd) #Desvio padrão 
sep_sd
##      Group.1         x
## 1     setosa 0.3524897
## 2 versicolor 0.5161711
## 3  virginica 0.6358796

Visualização dos dados: Um plot violino básico

Primeiramente vamos visitar o site: https://www.data-to-viz.com/ Aqui temos uma coleção de tipos de gráficos e os códigos para fazê-los no R.

# install.packages("ggplot2")
# install.packages("gridExtra")
library(ggplot2)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine

Também podemos usar plots violinos para representar duas variáveis. Vamos usar novamente os dados das plantas Iris para visualizar:

p <- ggplot(iris, aes(x=Species, y=Sepal.Length, fill=Species)) + 
  geom_violin()

aes: função que determina a estética do seu plot fill: argumento de preenchimento do plot - os grupos de interesse.

p2 <- ggplot(iris, aes(x=Species, y=Petal.Length, fill=Species)) + 
  geom_violin()
grid.arrange(p, p2, ncol=2)

Podemos melhorar o nosso gráfico:

p <- ggplot(iris, aes(x=Species, y=Sepal.Length, fill=Species)) + 
  geom_violin()+
  xlab("")+
  ylab("comprimento")+
  annotate("text", x=1, y =7.7, label = "Sépala", size = 5)+ 
  theme(legend.position = "none")
  
p2 <- ggplot(iris, aes(x=Species, y=Petal.Length, fill=Species)) + 
  geom_violin()+
  xlab("")+
  ylab("")+
  annotate("text", x=1, y =6.55, label = "Pétala", size = 5)+ 
  theme(legend.position = "none")
pdf("Flores.pdf")
grid.arrange(p, p2, ncol=2)
dev.off()
## png 
##   2

Testes de hipóteses

1. One-way Anova

ANOVA one-way (unidirecional): Quando há uma única variável independente categórica (também conhecida como fator) e uma única variável dependente contínua. Uma análise de variância procura testar variações significativas nas médias da variável dependente (y) ao longo dos níveis da variável independente (x).

Vamos testar se diferentes dietas em galinhas tem efeito no peso dos animais:

dados_ch<-ChickWeight
head(dados_ch)
##   weight Time Chick Diet
## 1     42    0     1    1
## 2     51    2     1    1
## 3     59    4     1    1
## 4     64    6     1    1
## 5     76    8     1    1
## 6     93   10     1    1
dim(dados_ch)
## [1] 578   4
str(dados_ch)
## Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame':   578 obs. of  4 variables:
##  $ weight: num  42 51 59 64 76 93 106 125 149 171 ...
##  $ Time  : num  0 2 4 6 8 10 12 14 16 18 ...
##  $ Chick : Ord.factor w/ 50 levels "18"<"16"<"15"<..: 15 15 15 15 15 15 15 15 15 15 ...
##  $ Diet  : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
##  - attr(*, "formula")=Class 'formula'  language weight ~ Time | Chick
##   .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv> 
##  - attr(*, "outer")=Class 'formula'  language ~Diet
##   .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv> 
##  - attr(*, "labels")=List of 2
##   ..$ x: chr "Time"
##   ..$ y: chr "Body weight"
##  - attr(*, "units")=List of 2
##   ..$ x: chr "(days)"
##   ..$ y: chr "(gm)"
dados_ch$peso_idade<-dados_ch$weight/dados_ch$Time #criando uma nova coluna de dados
hist(dados_ch$peso_idade) #chacando a distribução dos dados criados

model <- aov(dados_ch$peso_idade ~ dados_ch$Diet) # a função aov ajusta  uma análise de variância.

Encontramos um erro, o modelo não pode ser criado porque a variável y possui valores não-finitos. Vamos checar:

table(is.finite(dados_ch$peso_idade))
## 
## FALSE  TRUE 
##    50   528

50 valores não podem ser analisados como estão. Uma solução é a transformação de dados. Outra é a remoção destas amostras, que é o que vamos fazer

#Vamos transformar todos os dados em NA
dados_ch$peso_idade[!is.finite(dados_ch$peso_idade)] <- NA
#E remover os NAs de todo o conjunto de dados
dados_ch<-na.omit(dados_ch)
# Testando se funcionou:
table(is.finite(dados_ch$peso_idade))
## 
## TRUE 
##  528

E agora sim podemos criar nosso modelo usando a função aov:

model <- aov(dados_ch$peso_idade ~ dados_ch$Diet)

# E ver os resultados usando a função summary:
summary(model)
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## dados_ch$Diet   3    756   252.0   11.35 3.19e-07 ***
## Residuals     524  11634    22.2                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

O ajuste dos dados ao modelo indica o quanto o tipo da dieta explica sobre a variabilidade no ganho de peso. Olhando os resultados acima e considerando os seguintes resultados:

Df: Os graus de liberdade do modelo Sum Sq: soma dos quadrados - é a variância que o modelo não explicou. Mean Sq: media ao quadrado - é a variância explicada por cada componente . F-value: estatística F - é a medida usada para comparar as médias ao quadrados dentro e entre os grupos. Pr(>F): Tvalor de p da estatística F = significância estatística Residuals: Desvio relativo da média do grupo.

Nós podemos concluir, portanto que existem diferenças significativas entre os tipos de tratamento.

Comparações múltiplas pareadas entre as médias de grupo:

Um p-value significativo indica a diferença entre os grupos, mas ainda precisamos definir estes grupos. As comparações pareadas são usadas para identificarmos as diferenças nas médias entre certos pares de grupos e se estas são estatisticamente significativas.

Uma das formas de realizarmos estas comparações par-a-par é utilizando o teste de Tukey HSD (Tukey Honest Significant Differences) com a função TukeyHSD

tuk<-TukeyHSD(model, conf.level=.95)

Apenas a diferença entre os grupos 1-3 e 1-4 são significativamente diferentes.

Outra forma de testarmos as diferenças entre as médias dos grupos é utilizando um teste T

pairwise.t.test(dados_ch$peso_idade, dados_ch$Diet,
                 p.adjust.method = "BH")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  dados_ch$peso_idade and dados_ch$Diet 
## 
##   1       2     3    
## 2 0.036   -     -    
## 3 8.3e-06 0.036 -    
## 4 8.3e-06 0.036 0.945
## 
## P value adjustment method: BH

A utilização de cada um destes métodos (e de outros) depende muito dos seus dados, principalmente da distribuição. Nesta página temos um resumo sobre as diferentes casos para aplicação dos testes: https://statplace.com.br/blog/qual-o-melhor-teste-para-a-comparacao-de-medias/

Plotando os resultados da ANOVA

Para plotar este resultados vamos um pouquinho mais longe. A função multcompLetters4 (do pacote…) cria letras de acordo com as categorias de similaridades baseadas nos resultados da ANOVA e do teste Tukey.

cat <- multcompLetters4(model, tuk)
print(cat)
## $`dados_ch$Diet`
##    3    4    2    1 
##  "a"  "a" "ab"  "b"

Agora vamos criar uma tabela com os fatores

Tk <- group_by(dados_ch, Diet) %>%
  summarise(mean=mean(peso_idade), quant = quantile(peso_idade, probs = 0.75)) %>%
  arrange(desc(mean))

E incluir o display de letras na tabela criada

cld <- as.data.frame.list(cat$`dados_ch$Diet`)
Tk$cat <- cld$Letters
ggplot(dados_ch, aes(Diet, peso_idade)) + 
  geom_boxplot()+
  labs(x="Tipo de dieta", y="Peso por idade") +
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_text(data = Tk, aes(x = Diet, y = quant, label = cat), size = 4, vjust=-1, hjust =-1)

Testando os resíduos

Homogeneidade das variâncias

plot(model, 1)

Normalidade

plot(model, 2)

Levene’s Test (center = median)

#leveneTest(ados_ch$peso_idade ~ dados_ch$Diet)

2. Anova two-way

Um teste de ANOVA bivariada é usada para avaliar o efeito simultâneo de duas variáveis preditoras categóricas (x) em uma variável resposta (y)

library(carData)
data(OBrienKaiserLong)

Tratamentos: control, A, B. gender: F, M phase: pre, post, fup.

table(OBrienKaiserLong$treatment, OBrienKaiserLong$phase)
##          
##           pre post fup
##   control  25   25  25
##   A        20   20  20
##   B        35   35  35

A construção do seu modelo depende das suas hipóteses!

mod_aov1 <- aov(score ~ treatment + phase , data = OBrienKaiserLong)
summary(mod_aov1)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## treatment     2  186.7   93.37   24.39 2.38e-10 ***
## phase         2  167.5   83.75   21.87 1.94e-09 ***
## Residuals   235  899.8    3.83                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod_aov2 <- aov(score ~ treatment * phase , data = OBrienKaiserLong)
summary(mod_aov2)
##                  Df Sum Sq Mean Sq F value   Pr(>F)    
## treatment         2  186.7   93.37  26.216 5.48e-11 ***
## phase             2  167.5   83.75  23.514 5.07e-10 ***
## treatment:phase   4   77.0   19.25   5.405 0.000354 ***
## Residuals       231  822.8    3.56                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod_aov3 <- aov(score ~ phase +treatment*gender , data = OBrienKaiserLong)
summary(mod_aov3)
##                   Df Sum Sq Mean Sq F value   Pr(>F)    
## phase              2  167.5   83.75   27.32 2.22e-11 ***
## treatment          2  186.7   93.37   30.46 1.80e-12 ***
## gender             1   58.3   58.29   19.01 1.95e-05 ***
## treatment:gender   2  130.2   65.12   21.24 3.38e-09 ***
## Residuals        232  711.2    3.07                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(OBrienKaiserLong) +
  aes(x = treatment, y = score, fill = phase) +
  geom_boxplot()+
  scale_fill_brewer(viridis(3))+
  guides(fill=guide_legend(title="Phases"))

3. Correlação e Regressão

Enquanto uma correlação mede a força ou grau de relação entre duas variáveis. A regressão é uma equação (modelo) que descreve essa relação com valores específicos.

Testando a regra de Bergman em Aves.

A regra de Bergman identifica um padrão de diminuição do tamanho corporal entre as espécies com a diminuição da latitude. Espécies com latitudes mais próximas de zero teriam menor tamanho de corpo.

Dados utilizados: <AVONET: morphological, ecological and geographical data for all birds>[https://onlinelibrary.wiley.com/doi/full/10.1111/ele.13898]

aves<-read.csv("C:/Users/raque/OneDrive - ufpr.br/Área de Trabalho/Curso_R_basico_CAEB/Aves_dados.csv")

Vamos criar um conjunto de dados somente com os traits que vamos utilizar.

dat<-select(aves, Specie, order, Mass, Centroid.Latitude)
dat<-na.omit(dat) #e remover todos os NAs

Vamos checar se existe algum nome duplicado na mossa coluna com os nomes das espécies:

table(duplicated(dat$Specie)) 
## 
## FALSE 
##  9855

Conferindo a dimensão dos nossos dados:

dim(dat)
## [1] 9855    4

Temos 9855 especies para realizar as análises

Teste de Correlação

Para o teste vamos utilizar a função cor_res:

cor_res<-cor.test(abs(dat$Centroid.Latitude), log(dat$Mass))
cor_res
## 
##  Pearson's product-moment correlation
## 
## data:  abs(dat$Centroid.Latitude) and log(dat$Mass)
## t = 17.11, df = 9853, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1506242 0.1889728
## sample estimates:
##       cor 
## 0.1698628

Legal! Encontramos uma correlação entre as variáveis. Agora, uma regreção vai nos mostrar a direção desta correlação (para sabermos se nossa hipótese pode ser aceita) e principalmente o quanto um modelo linear explica o padrão dos dados.

3. Regressão linear

A função lm calcula uma equação (um modelo linear) com base nos dados e, então com base no ajuste dos dados, consegue demonstrar o quanto uma variável (x) pode explicar da outra (y).

O ~ é utilizado para separar as variáveis y e x, nesta ordem. No nosso caso, imaginamos que a latitude (correspondente as variações climaticas e ambientais) terá uma influência no tamanho das espécies. Então latitude=X, Massa=Y

reg_res<-lm(log(dat$Mass)~abs(dat$Centroid.Latitude))
summary(reg_res)
## 
## Call:
## lm(formula = log(dat$Mass) ~ abs(dat$Centroid.Latitude))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4593 -1.1857 -0.3099  0.9317  7.9517 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                3.606679   0.023519  153.35   <2e-16 ***
## abs(dat$Centroid.Latitude) 0.017756   0.001038   17.11   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.553 on 9853 degrees of freedom
## Multiple R-squared:  0.02885,    Adjusted R-squared:  0.02875 
## F-statistic: 292.7 on 1 and 9853 DF,  p-value: < 2.2e-16

Apesar de observarmos uma correlação significativa, algumas estatísticas nos mostram que o poder de explicação da latitude no tamanho dos organismos é muito baixo (R²=2.8%). O valor do estimate também é baixo (0.017756) mas positivo, mostrando um aumento na massa, com o aumento da latitude central da espécie.

Plotando
plot(abs(dat$Centroid.Latitude), log(dat$Mass), xlab= "Latitude central absoluta", ylab = "Massa corporal (log)", pch=20)
abline(reg_res, col="red", lwd= 2)

Para fins exploratórios, vamos fazer um plot único com todas as retas de regressão para cada ordem.

ggplot(dat,aes(abs(Centroid.Latitude), log(Mass), colour=factor(order))) +
  geom_point() +
  geom_smooth(method='lm')+ 
  theme(legend.position = "none") #este comando remove a legenda
## Warning in qt((1 - level)/2, df): NaNs produzidos

## Warning in qt((1 - level)/2, df): NaNs produzidos

## Warning in qt((1 - level)/2, df): NaNs produzidos
## Warning in max(ids, na.rm = TRUE): nenhum argumento não faltante para max;
## retornando -Inf

## Warning in max(ids, na.rm = TRUE): nenhum argumento não faltante para max;
## retornando -Inf

## Warning in max(ids, na.rm = TRUE): nenhum argumento não faltante para max;
## retornando -Inf

Existe uma varição entre as retas calculadas pelo modelo linear para cada ordem. Vamos então usar um loop para calcular as estatísticas (p-value, slope e R²) para cada ordem com um número mínimo de espécies (50 espécies)

Primeiro vamos filtrar as ordens com mais de 50 espécies:

orders <- data.frame(table(dat$order)) %>% 
         filter(Freq>= 50) %>% as.list()
length(orders$Var1)
## [1] 19

Criando uma tabela para guardar os resultados:

tab_res<-data.frame(matrix(ncol = 4, nrow = 19))
colnames(tab_res) <- c("Order", "p_value", "slope", "R2")
rownames(tab_res) <- orders$Var1

tab_res
##                   Order p_value slope R2
## Accipitriformes      NA      NA    NA NA
## Anseriformes         NA      NA    NA NA
## Apodiformes          NA      NA    NA NA
## Bucerotiformes       NA      NA    NA NA
## Caprimulgiformes     NA      NA    NA NA
## Charadriiformes      NA      NA    NA NA
## Columbiformes        NA      NA    NA NA
## Coraciiformes        NA      NA    NA NA
## Cuculiformes         NA      NA    NA NA
## Falconiformes        NA      NA    NA NA
## Galliformes          NA      NA    NA NA
## Gruiformes           NA      NA    NA NA
## Passeriformes        NA      NA    NA NA
## Pelecaniformes       NA      NA    NA NA
## Piciformes           NA      NA    NA NA
## Procellariiformes    NA      NA    NA NA
## Psittaciformes       NA      NA    NA NA
## Strigiformes         NA      NA    NA NA
## Suliformes           NA      NA    NA NA

Em loop simples uma variável pode ser substituida por uma sequencia de (no nosso caso) números, a cada rodada o i é substituído pelo número seguinte.

Note que nós poderíamos fazer o mesmo código abaixo 19 vezes, substituindo o nome das ordens e assim preencher a tabela com os resultados:

dat_o<-subset(aves, order=="Accipitriformes") #esta ordem é a primeira da lista
reg<-lm(log(dat_o$Mass)~abs(dat_o$Centroid.Latitude))
res<-summary(reg)
# cada uma dos seguintes códigos vai salvar o objeto seguinte ao <- na respectiva coluna ($Order, $p_value, $slope e $R2) e linha correspondente ordem (linha 1, neste caso).
tab_res$Order[1]<- "Accipitriformes" 
tab_res$p_value[1] <- res$coefficients[2,4] 
tab_res$slope[1] <- res$coefficients[2,1]
tab_res$R2[1] <- res$r.squared*100 

res$coefficients
##                               Estimate  Std. Error   t value      Pr(>|t|)
## (Intercept)                  6.3920230 0.099990768 63.926132 3.878778e-153
## abs(dat_o$Centroid.Latitude) 0.0130809 0.004212123  3.105535  2.126697e-03

A parte res\(coefficients corresponde ao local onde está o p-value dentro do resultado da regressão (linha 2, coluna 4). res\)coefficients[2,1] é onde está o valor do estimate (slope) dentro do resultado da regressão (linha 2, coluna 1) res$r.squared é onde está o valor do R² (slope) dentro do resultado da regressão. Podemos multiplicar o valor por 100 e obter a % de explicação.

Veja que só a primeira linha foi preenchida:

tab_res
##                             Order     p_value     slope       R2
## Accipitriformes   Accipitriformes 0.002126697 0.0130809 3.847822
## Anseriformes                 <NA>          NA        NA       NA
## Apodiformes                  <NA>          NA        NA       NA
## Bucerotiformes               <NA>          NA        NA       NA
## Caprimulgiformes             <NA>          NA        NA       NA
## Charadriiformes              <NA>          NA        NA       NA
## Columbiformes                <NA>          NA        NA       NA
## Coraciiformes                <NA>          NA        NA       NA
## Cuculiformes                 <NA>          NA        NA       NA
## Falconiformes                <NA>          NA        NA       NA
## Galliformes                  <NA>          NA        NA       NA
## Gruiformes                   <NA>          NA        NA       NA
## Passeriformes                <NA>          NA        NA       NA
## Pelecaniformes               <NA>          NA        NA       NA
## Piciformes                   <NA>          NA        NA       NA
## Procellariiformes            <NA>          NA        NA       NA
## Psittaciformes               <NA>          NA        NA       NA
## Strigiformes                 <NA>          NA        NA       NA
## Suliformes                   <NA>          NA        NA       NA

Como temos 19 ordens para refazer a mesma análise sendo que a única mudança será o nome da ordem - e a posição dos resultados na tabela que podem corresponder ao mesmo número de objetos na lista de ordens. Nestes casos é muito simples construir um loop, substintuindo o i pelos valores de 1 a 19 a cada rodada:

for (i in 1:19){
  ord<-as.character(orders$Var1[i])
  dat_o<-subset(aves, order==ord)
  reg<-lm(log(dat_o$Mass)~abs(dat_o$Centroid.Latitude))
  res<-summary(reg)
  tab_res$Order[i]<- ord
  tab_res$p_value[i] <- res$coefficients[2,4] 
  tab_res$slope[i] <- res$coefficients[2,1]
  tab_res$R2[i] <- res$r.squared*100}

Agora podemos ver e comparar as estatpisticas para cada ordem:

tab_res
##                               Order      p_value         slope           R2
## Accipitriformes     Accipitriformes 0.0021266973  0.0130808961  3.847822218
## Anseriformes           Anseriformes 0.0008585930  0.0097858620  6.810821840
## Apodiformes             Apodiformes 0.1285208290  0.0070326704  0.526671773
## Bucerotiformes       Bucerotiformes 0.9028599144  0.0030349813  0.024217482
## Caprimulgiformes   Caprimulgiformes 0.5834112248 -0.0030746377  0.287383811
## Charadriiformes     Charadriiformes 0.0369858106  0.0049205511  1.189853556
## Columbiformes         Columbiformes 0.9235173862 -0.0004018295  0.003046816
## Coraciiformes         Coraciiformes 0.4393613121  0.0058355076  0.401847535
## Cuculiformes           Cuculiformes 0.3226731165 -0.0069997871  0.724553229
## Falconiformes         Falconiformes 0.0319743560  0.0137397474  7.206017836
## Galliformes             Galliformes 0.3893321824  0.0031043661  0.261966098
## Gruiformes               Gruiformes 0.0004893399  0.0255380624  7.661439986
## Passeriformes         Passeriformes 0.5074520380 -0.0005285277  0.007471743
## Pelecaniformes       Pelecaniformes 0.0772727872  0.0133255654  2.970705069
## Piciformes               Piciformes 0.0026419488  0.0107875495  2.193874312
## Procellariiformes Procellariiformes 0.0013405839  0.0231958641  7.929340084
## Psittaciformes       Psittaciformes 0.0409523232  0.0097201439  1.191687920
## Strigiformes           Strigiformes 0.0083820099  0.0130137354  3.440629939
## Suliformes               Suliformes 0.0236785302  0.0076625715 10.013556434

Exercício

Escolha um conjunto de dados (no pacote datasets) ou utilize seus próprios dados para testar uma hipótese de sua escolha.