Análise dos dados morfológicos dos BlueJays

# Carregar librarias
library(tidyverse)
library(grid)
library(Stat2Data)
library(magrittr)
library(broom)
library(gridExtra)
library(corrplot)
library(car)
library(nortest)
library(ellipse)
library(mvnormalTest)

O conjunto de dados a seguir consiste de 123 observações sobre as seguintes 8 variáveis:

• BirdID: rótulo para identificação da ave
• KnowSex: sexo da ave codificado como F ou M
• BillDepth: Espessura do bico medida na narina (em mm)
• BillWidth: Largura do bico (em mm)
• BillLength: Comprimento do bico (em mm)
• Head: Distância da ponta do bico até a parte de trás da cabeça (em mm)
• Mass: Massa corporal (em gramas)
• Skull: Distância da base do bico até a parte de trás do crânio (em mm)

Carregando os dados

data("BlueJays")
blue_jays_data<- BlueJays[,-9]
head(blue_jays_data)

## Boxplots para as aves fêmeas

female <- blue_jays_data[blue_jays_data$KnownSex == "F",]
female <- female[,-c(1,2)]
boxplot(female)

Densidades para as aves fêmeas

female <- blue_jays_data[blue_jays_data$KnownSex == "F",]
female <- female[,-c(1,2)]
nomes_female <- colnames(female)
par(mfrow=c(4,2))
for (i in 1:ncol(female)) {
  plot(density(female[,i]), main = paste("Densidade de",nomes_female[i]))
  polygon(density(female[,i]), col = "coral")
}

## Boxplots para as aves machos

male <- blue_jays_data[blue_jays_data$KnownSex == "M",]
male <- male[,-c(1,2)]
boxplot(male)

Densidades para as aves machos

male <- blue_jays_data[blue_jays_data$KnownSex == "M",]
male <- male[,-c(1,2)]
nomes_male <- colnames(male)
par(mfrow=c(4,2))
for (i in 1:ncol(male)) {
  plot(density(female[,i]), main = paste("Densidade de",nomes_male[i]))
  polygon(density(female[,i]), col = "purple")
}

Diagrama de dispersão

blue_jays_data_numeric <- blue_jays_data[c(-1,-2)]
pairs(blue_jays_data_numeric, 
      pch=21,
      bg = c("coral","purple")[as.factor(blue_jays_data$KnownSex)])

Gráfico de Correlação

corrplot(cor(blue_jays_data_numeric), method = "square" , order = "hclust", tl.col='black', tl.cex=.75) 

Testando a igualdade dos grupos

Nesta seção iremos realizar um teste de hipóteses para comparar se os dois grupos são iguais em relação ao seu vetor médios. A suposição de normalidade multivariada não parece ser satisfeita. Contudo, por razões didáticas faremos o teste sob essa suposição. Com efeito, iremos testar:

\[\begin{align}H_0:\boldsymbol{\mu}_1 = \boldsymbol{\mu}_2\\ H_1:\boldsymbol{\mu}_1 \neq \boldsymbol{\mu}_2\end{align},\] em que \(\boldsymbol{\mu}_1\) e \(\boldsymbol{\mu}_2\) são os vetores de medias de [BillDepth, BillWidth, BillLength, Mass,Skull] para as fêmeas e os machos, respectivamente.

X1 <- female
X2 <- male
#vetor de médias por sexo
X1_barra <- apply(X1,2,mean)
X2_barra <- apply(X2,2,mean)
# matrizes de covariâncias por sexo
S1 <- cov(X1)
S2 <- cov(X2)
# estimativa da matriz de covariâncias (estimador comum)
p  <- ncol(X1) 
n1 <- nrow(X1)
n2 <- nrow(X2)
W1 <- (n1-1)*S1
W2 <- (n2-1)*S2
# matriz de covariâncias pooled 
Spl <- (W1 + W2)/(n1+n2-2) 
invSpl <- solve(Spl)
# Estatística do teste T2 de Hotelling
t2calc <- ((n1*n2)/(n1+n2))*t(X1_barra - X2_barra)%*%invSpl%*%(X1_barra-X2_barra)
fcalc <- (n1+n2-1-p)/((n1+n2-2)*p)*t2calc
pvalor <- pf(fcalc,df1=p,df2 = n1+n2-p-1,lower.tail = FALSE)
resultados<-list(t2=as.numeric(t2calc),
     f=as.numeric(fcalc),
     pvalor=as.numeric(pvalor))
resultados

$t2
[1] 122.4434

$f
[1] 19.56395

$pvalor
[1] 1.137787e-15

Pelos resultados do teste, rejeitamos \(H_0\), portanto há diferença significativa entre as aves fêmeas e machos.