Aula Prática 1 (GET00126) - Normal Multivariada

# librarias
library(nortest)
library(corrplot)
library(car)
library(ellipse)
library(mvnormalTest)
library(ACSWR)

1. O objeto air_polution contém 42 observações de 7 variáveis sobre a poluição de ar registradas às 12:00 PM em Los Angeles em diferentes dias. Com os dados sobre a poluição do ar examine a suposição de normalidade bivariada para o par de variáveis \(X_5\) = NO2 e \(X_6\) = O3.

1. Examine a normalidade univariada de \(X_5\) e \(X_6\).
2. Calcule as formas quadráticas \(d_i^2 = (\mathbf{x}_i - \mathbf{\bar{x}})^\top \mathbf{S}^{-1}(\mathbf{x}_i - \mathbf{\bar{x}})\), \(i = 1,\ldots, 42\), \(\mathbf{x}_j^\top = [x_{i5},x_{i6}]\) e construa um Q-Q plot apropriado para a forma quadrática \(d^2\). A suposição de normalidade bivariada para os dados é razoável?
3. Determine a proporção de observações de \(\mathbf{x}_i^\top = [x_{i5},x_{i6}]\), \(i = 1,\ldots, 42\) dentro do gráfico de contorno de uma normal bivariada com uma probabilidade de 50%.

air_pollution <- read.table("~/Documents/Cursos_GET_UFF/20241GET00126A1/Datasets/air_pollution.dat")
colnames(air_pollution) <- c("wind","solar_radiation","CO","NO","NO2","O3","HC")
head(air_pollution)

air_pollution_56 <- air_pollution[,c(5,6)]
pairs(air_pollution_56, pch = 19)

boxplot(air_pollution_56, pch = 19)

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

par(mfrow=c(1,2))
for (j in 1:nvar)
{
  qqPlot(scale(air_pollution_56[,j]),dist="norm",mean=0,sd=1,col.lines=1,grid="FALSE",
         xlab="quantil da N(0,1)",ylab=nomes[j],cex=1.2,pch=20,id=FALSE)
}

# Shapiro Wilks Test
shapiro.test(air_pollution_56[,1])

    Shapiro-Wilk normality test

data:  air_pollution_56[, 1]
W = 0.93004, p-value = 0.013

shapiro.test(air_pollution_56[,2])

    Shapiro-Wilk normality test

data:  air_pollution_56[, 2]
W = 0.89603, p-value = 0.001098

n <- nrow(air_pollution_56)
nvar <- ncol(air_pollution_56)
vmu <- colMeans(air_pollution_56)
mcov <- cov(air_pollution_56)
# forma quadrática
d2 <- n*mahalanobis(air_pollution_56,center=vmu,cov=mcov)
qqPlot(d2,dist="chisq",df=nvar,col.lines=1,grid="FALSE",
       xlab="quantil da distribuição qui-quadrado",
       ylab="d2",
       cex=1.2, 
       pch = 20,
       id = FALSE)

Mardia_Test<-mardia(air_pollution_56)
Mardia_Test$mv.test

rho <- cor(air_pollution_56)
var <- diag(mcov)
i=1;j=2
plot(ellipse::ellipse(x = rho[i,j],
                      scale = c(sqrt(var[i]),sqrt(var[j])),
                      centre = c(vmu[i],vmu[j]),
                      level = 0.5
                      ),
     type = "l",
     xlab = nomes[i],
     ylab = nomes[j])
lines(vmu[i],vmu[j], pch = 8, col = "red", type = "p")

points(air_pollution_56[,i],air_pollution_56[,j],
       pch=20,
       cex=1.2,
       cex.lab=1.2,
       cex.axis=1.2,
       col="darkgreen")

2. Um ecologista de vida selvagem, mediu \(x_1\): comprimento do cauda (tail_length) em milimetros e \(x_2\): comprimento das asas (wing_length) em milimetros em uma amostra de 45 aves fêmeas da especie hook-billed kites (Gavião-caracoleiro em português). As observações estão armazenadas no objeto birds_data.

1. Examine a normalidade univariada das variáveis \(X_1\): tail_length e \(X_2\) = wing_length.
2. Examine a normalidade bivariada de \([X_1,X_2]^\top\). É um modelo viável para \([X_1,X_2]^\top\)?
3. Suponha que é sabido que \(\mu_1\) = 190 mm e \(\mu_2\) = 275 mm para machos da espécie hook-billed kites. São esses valores razoáveis para o comprimento médio da cauda e o comprimento médio das assas para aves fêmeas? Explique.

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

3. O objeto psy_profile_data consiste de 130 observações geradas por pontuações (scores) de um teste psicológico administrado em adolescentes peruanos (15, 16 e 17 anos de idade). Para cada adolescente o sexo [gender (masc. = 1 , fem. = 2)] e estado socioeconômico [socio (baixo = 1, médio = 2)] também foi registrado. As pontuações foram acumuladas em 5 subescalas rotuladas como: independência (indep), suporte (supp), benevolência (benev), conformidade (conform), e liderança (leader).

1. Examine a normalidade univariada de cada uma das seguintes variáveis: indep, supp, benev, conform, e leader.
2. Usando todas as 5 variáveis, verifique a suposição de normalidade multivariada.
3. Construa uma elipse de 99% de confiança para cada par de variáveis.
4. Assumindo normalidade multivariada. Determine os comprimentos e direções dos eixos do ellipsóide de 95% de confiança para \(\mathbf{\mu}\).

psy_profile_data <- read.table("~/Documents/Cursos_GET_UFF/20241GET00126A1/Datasets/psy_profile.dat")
colnames(psy_profile_data) <- c("indep","supp","benev","conform","leader","gender","socio")
head(psy_profile_data)

4. Um estudo psicológico formado por 4 testes foi conduzido em 32 mulheres e 32 homens. As pontuações de cada teste para cada indivíduo foram coletadas em um vetor [pict, formas,brinqu,vocab]. As componentes são descritas a seguir:

• pict: score obtido no teste de inconsistência pictoriais,
• formas: score obtido no teste de placas de formas,
• brinqu: score obtido no teste de reconhecimento de brinquedos,
• vocab: score obtido no teste de vocabulário.

Um dos interesses do estudo é saber se o vetor de médias das pontuações dos testes é o mesmo para homens e mulheres. Para isso conduza um teste de hipóteses apropriado. Não esqueça de checar todas as suposições.

data(mfp)
homens = mfp[,1:4]
mulheres = mfp[,5:8]
colnames(homens) <- c("pict", "formas", "brinqu", "vocab")
colnames(mulheres) <- c("pict", "formas", "brinqu", "vocab")
psy_tests <- rbind(homens,mulheres)
sexo <- rbind(matrix("M",nrow = 32),matrix("F",nrow=32))
psy_tests$sexo <- sexo
head(psy_tests)

n1  <- nrow(homens); n2 <- nrow(mulheres); p <- ncol(mulheres)
m1  <- matrix(colMeans(homens),ncol = 1); m2 <- matrix(colMeans(mulheres),ncol = 1) 
S1  <- var(homens); S2 <- var(mulheres)
Spl <- (1/(n1+n2-2))*((n1-1)*S1 + (n2-1)*S2); invSpl <- solve(Spl)
t2_calc <- ((n1*n2)/(n1+n2))*t(m1 - m2)%*%invSpl%*%(m1-m2)
f_calc  <- (n1+n2-1-p)/((n1+n2-2)*p)*t2calc
p_valor <- pf(fcalc,df1=p,df2 = n1+n2-p-1,lower.tail = FALSE)
data.frame(t2_calc,f_calc,p_valor)