Inciso (a)

Se cargan los datos desde un archivo .DAT. Se quedan solo tres variables de interés: Fuel, Repair y Capital.

library(readr)
# Cargar datos de T6-10.DAT
Datos<-read.table("T6-10.DAT",header = TRUE, 
                   col.names = c("Fuel", "Repair", "Capital","Combustible"))
datos<- Datos[, -4]

Se generan Q-Q plots para comparar cada variable con la distribución normal.
Si los puntos siguen la línea diagonal indican distribución normal.
Desviaciones grandes indican colas pesadas o presencia de outliers.

#Inciso a) Q-Q plots y diagramas de dispersión

par(mfrow=c(1,3))
qqnorm(datos$Repair, main="Q-Q Repair"); qqline(datos$Repair)
qqnorm(datos$Capital, main="Q-Q Capital"); qqline(datos$Capital)
qqnorm(datos$Fuel, main="Q-Q Fuel"); qqline(datos$Fuel)

Se dibujan diagramas de dispersión entre cada par de variables.Sirven para detectar relaciones lineales y puntos atípicos. Así como la distancias de mahalanobis para visualizar los outliers

pairs(datos, main="Scatterplots")

#Calcular la matriz de covarianzas y la media
Sigma <- cov(datos)
mu <- colMeans(datos)
#Calcular la distancia de Mahalanobis para cada observación
md2 <- mahalanobis(datos, center = mu, cov = Sigma)
p <- ncol(datos)
cut <- qchisq(0.95, df = p)  
which(md2 > cut)
## [1]  8 20 40 46 55
# Visualizacion
hist(md2, breaks = 20, main = "Distancias de Mahalanobis", xlab = "MD²")

Se eliminan las observaciones 8,20,40,46 y 55, previamente identificadas como outliers.
Y se repite el análisis gráfico, pero sin los outliers.

# Q-Q y scatterplots sin los outliers (obs 8,20,40,46 y 55) ------
datos_sin<- datos[-c(8,20,40,46,55),]

par(mfrow=c(1,3))
qqnorm(datos_sin$Repair, main="Q-Q Repair sin outliers"); qqline(datos_sin$Repair)
qqnorm(datos_sin$Capital, main="Q-Q Capital sin outliers"); qqline(datos_sin$Capital)
qqnorm(datos_sin$Fuel, main="Q-Q Fuel sin outliers"); qqline(datos_sin$Fuel)

pairs(datos_sin, main="Scatterplots sin outliers")

Los Q-Q plots se ajustan mejor a la recta indicando mayor evidencia de normalidad.

Se calculan los valores críticos (basados en la distribución F) para construir las elipses de confianza al 95%.

# Parámetros
n <- nrow(datos_sin)
p <- 2
f.inv <- qf(1-0.05, p, n-p)
ratio <- (p*(n-1))/(n*(n-p))
factor <- sqrt(ratio*f.inv)

# Librería ellipse
library(ellipse)
## 
## Adjuntando el paquete: 'ellipse'
## The following object is masked from 'package:graphics':
## 
##     pairs
#Repair vs Capital
cov_RC <- cov(cbind(datos_sin$Repair, datos_sin$Capital))
medias_RC <- colMeans(cbind(datos_sin$Repair, datos_sin$Capital))

plot(ellipse(cov_RC, centre=medias_RC, t=factor),
     type="l", col="lightblue", lwd=2,
     xlab="Repair", ylab="Capital",
     main="Elipse 95% Repair vs Capital")
points(medias_RC[1], medias_RC[2], col="red", pch=19)

# Repair vs Fuel
cov_RF <- cov(cbind(datos_sin$Repair, datos_sin$Fuel))
medias_RF <- colMeans(cbind(datos_sin$Repair, datos_sin$Fuel))
plot(ellipse(cov_RF, centre=medias_RF, t=factor),
     type="l", col="lightblue", lwd=2,
     xlab="Repair", ylab="Fuel",
     main="Elipse 95% Repair vs Fuel")
points(medias_RF[1], medias_RF[2], col="red", pch=19)

# Capital vs Fuel 
cov_CF <- cov(cbind(datos_sin$Capital, datos_sin$Fuel))
medias_CF <- colMeans(cbind(datos_sin$Capital, datos_sin$Fuel))

plot(ellipse(cov_CF, centre=medias_CF, t=factor),
     type="l", col="lightblue", lwd=2,
     xlab="Capital", ylab="Fuel",
     main="Elipse 95% Capital vs Fuel")
points(medias_CF[1], medias_CF[2], col="red", pch=19)

Se grafica la elipse de confianza bivariada para el par (Repair, Capital). La elipse delimita la región donde se espera que esté la media conjunta al 95%.
El punto rojo marca el vector de medias. Esto mismo se repite después para (Repair vs Fuel) y (Capital vs Fuel).

Inciso (b)

Se calculan medias, matriz de covarianza S, número de observaciones y variables, y el nivel de significancia

#Inciso b) Intervalos de confianza 
library(MVN)
medias <- colMeans(datos_sin)
S <- cov(datos_sin)
n <- nrow(datos_sin)
p <- ncol(datos_sin)
alfa <- 0.05

Intervalos de confianza Bonferron

Se construyen intervalos de confianza Bonferroni (95%) para cada media:

# Bonferroni 
t_crit <- qt(1-0.05/(2*p), df=n-1)
Bonf <- matrix(NA, nrow=p, ncol=2)
for (i in 1:p){
  error <- t_crit*sqrt(S[i,i]/n)
  Bonf[i,] <- c(medias[i]-error, medias[i]+error)
}
rownames(Bonf) <- names(medias)
Bonf
##              [,1]      [,2]
## Fuel     9.786834 11.615430
## Repair   7.103967  9.917543
## Capital 10.294745 13.810161

Usan la t de Student ajustada y son más amplios para protegen contra el error de considerar varias comparaciones simultáneas.

Intervalos de confianza Hotelling T^2

Se construyen los intervalos simultáneos T^2 de Hotelling al 95%

#  Hotelling T^2
f_crit <- qf(1-0.05, p, n-p)
factor <- (p*(n-1))/(n-p)
T2 <- matrix(NA, nrow=p, ncol=2)
for (i in 1:p){
  error <- sqrt(factor*f_crit*S[i,i]/n)
  T2[i,] <- c(medias[i]-error, medias[i]+error)
}
rownames(T2) <- names(medias)
T2
##             [,1]     [,2]
## Fuel    9.610731 11.79153
## Repair  6.833005 10.18850
## Capital 9.956193 14.14871

Basados en la distribución F y son más eficientes y estrechos que Bonferroni cuando hay correlación entre variables.

En resumen general
En los Q-Q plots iniciales se notan desviaciones y por eso mismo segun el inciso se eliminan las observaciones 8, 20, 40, 46 y 55, asi los datos se ajustan mejor a la normalidad multivariada.
Las elipses de confianza muestran gráficamente la variabilidad conjunta de pares de variables y los intervalos Bonferroni son más conservadores (más largos).
Por ultimo los intervalos Hotelling T² son más cortos y eficientes, pero menos conservadores.