Aula Prática 2 (GET00126) - Testes para comparação de duas populações normais multivariadas

# librarias
library(nortest)
library(corrplot)
library(car)
library(ellipse)
library(mvnormalTest)
library(ACSWR)
library(Hotelling)
library(psych)
library(biotools)

Os dados a seguir contém informações sobre duas variáveis medidas em 45 aves fêmeas e 45 aves machos da especie hook-billed kites (veja a segunda questão 2 da primeira aula prática), a saber, comprimento da cauda \((x_1)\) e comprimento das asas \((x_2)\). Examine a normalidade bivariada dos dados, e identifique valores discrepantes.

birds_data_female<- read.table("~/Documents/Cursos_GET_UFF/20241GET00126A1/Datasets/birds_female.dat")
colnames(birds_data_female) <- c("tail_length","wing_length")

birds_data_male<- read.table("~/Documents/Cursos_GET_UFF/20241GET00126A1/Datasets/birds_male.dat")
colnames(birds_data_male) <- c("tail_length","wing_length")

group <- factor(c(rep("F",45), rep("M",45)))

birds_data_complete <- data.frame(rbind(birds_data_female,birds_data_male))
birds_data_complete$group <- group
birds_data_complete

pairs(birds_data_female, pch = 21, bg = "coral")

pairs(birds_data_male,pch = 21, bg = "purple")

pairs(birds_data_complete[,c(1,2)], pch = 21, bg=c("coral","purple")[birds_data_complete$group])

Aves Fêmeas

par(mfrow=c(2,2))
nvar  <- ncol(birds_data_female)
nomes <- names(birds_data_female)
for (i in 1:ncol(birds_data_female)) {
  hist(birds_data_female[,i], main = paste("Histograma de",nomes[i]), xlab = paste(nomes[i]))
  plot(density(birds_data_female[,i]), main = paste("Densidade de",nomes[i]))
}

Aves Machos

par(mfrow=c(2,2))
nvar  <- ncol(birds_data_male)
nomes <- names(birds_data_male)
for (i in 1:ncol(birds_data_male)) {
  hist(birds_data_male[,i], main = paste("Histograma de",nomes[i]), xlab = paste(nomes[i]))
  plot(density(birds_data_male[,i]), main = paste("Densidade de",nomes[i]))
}

mvnormalTest::mardia(birds_data_female)

$mv.test

$uv.shapiro
            W      p-value UV.Normality
tail_length 0.9698 0.2857  Yes         
wing_length 0.9813 0.6726  Yes         

Teste \(T^2\) de Hotelling

X1 <- birds_data_female
X2 <- birds_data_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 (pooled estimador)
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)
# Transformando T2 pela F
fcalc = (n1+n2-1-p)/((n1+n2-2)*p)*t2calc
fcalc

         [,1]
[1,] 1.800573

Teste M de Box

boxM(data = birds_data_complete[,-3],grouping = birds_data_complete[,3])

    Box's M-test for Homogeneity of Covariance Matrices

data:  birds_data_complete[, -3]
Chi-Sq (approx.) = 27.536, df = 3, p-value = 4.546e-06

1. Teste \(H_0: \boldsymbol{\mu}_1 - \boldsymbol{\mu}_2 = 0\). Qual a sua decisão?
   
2. Qual o p-valor do teste? Lembre que \(F_{p,p(n_1+n_2-2),\alpha} = \frac{(n_1+n_2-1-p)}{(n_1+n_2-2)p}T^2_{p,n_1+n_2-2,\alpha}\)
3. Usando que, \(T^2 = \frac{n_1n_2}{n_1+n_2}(\bar{\mathbf X}_1 - \bar{\mathbf X}_2)^\top\mathbf S_{pl}^{-1}(\bar{\mathbf X}_1 - \bar{\mathbf X}_2) \approx_{\mathcal D} \chi^2_p\) sob \(H_0\), para \(n_i - p\), \(i=1,2\) são grandes, realize um teste assintótico para testar \(H_0: \boldsymbol{\mu}_1 - \boldsymbol{\mu}_2 = 0\).
4. A suposição de igualdade das matrizes de covariâncias das fêmeas e dos machos é razoável?
5. Determine a elipse de 95% de confiança para \(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2\) e os intervalos de confiança simultâneos para os componentes de \(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2\).
6. Removendo o(s) valor(es) discrepantes, refaça (a)-(d).
7. Consulte sobre as funções hotelling.stat() e hotelling.test() do pacote Hotelling.