PUNTO 3
Una empresa farmaceutica ha realizado un estudio con n = 30 empleados que han consumido un nuevo suplemento durante 6 meses. Para cada sujeto se registraron las siguientes variables: Universidad del Rosario 2 de 3 a) X1: Resistencia f´ısica (minutos). b) X2: Tiempo de reacci´on (ms). c) X3: Memoria a corto plazo (puntuacion)
Cargar datos
setwd("C:/Users/prestamour/Downloads")
library(readxl)## Warning: package 'readxl' was built under R version 4.4.3library(MVN)## Warning: package 'MVN' was built under R version 4.4.3farmaceutica = read_excel("farmaceutica.xlsx")a)
# ---- (a) Cargar datos y calcular media, n y matriz de covarianza ----
# Cambia la ruta según dónde tengas el archivo
datos <- read_excel("farmaceutica.xlsx")
n <- nrow(datos) # tamaño de muestra
media <- colMeans(datos) # vector de medias
S <- cov(datos) # matriz de covarianza
cat("n =", n, "\n")## n = 30cat("Media muestral:\n")## Media muestral:print(media)## Resistencia Tiempo_Reaccion Memoria
## 35.65105 249.68097 80.45005cat("Matriz de covarianza S:\n")## Matriz de covarianza S:print(S)## Resistencia Tiempo_Reaccion Memoria
## Resistencia 13.277140 6.323688 9.47725
## Tiempo_Reaccion 6.323688 22.634125 10.48137
## Memoria 9.477250 10.481371 28.21304b)
# ---- (b) Verificar simetría y valores propios ----
cat("¿Es es simetrica? :", all(S == t(S)), "\n")## ¿Es es simetrica? : TRUEvalprop <- eigen(S)$values
cat("Valores propios de S:\n")## Valores propios de S:print(valprop)## [1] 40.944273 14.656308 8.523729c)
# ---- (c) Prueba de normalidad multivariada (Mardia) ----
# Ejecutar Mardia sin especificar mvnTest
mardia <- mvn(data = datos)
# Mostrar todo el resultado
print(mardia)## $multivariate_normality
## Test Statistic p.value Method MVN
## 1 Henze-Zirkler 0.626 0.42 asymptotic ✓ Normal
##
## $univariate_normality
## Test Variable Statistic p.value Normality
## 1 Anderson-Darling Resistencia 0.241 0.753 ✓ Normal
## 2 Anderson-Darling Tiempo_Reaccion 0.402 0.337 ✓ Normal
## 3 Anderson-Darling Memoria 0.425 0.297 ✓ Normal
##
## $descriptives
## Variable n Mean Std.Dev Median Min Max 25th 75th
## 1 Resistencia 30 35.651 3.644 35.688 27.512 44.067 33.371 38.038
## 2 Tiempo_Reaccion 30 249.681 4.758 250.441 240.030 259.326 246.247 253.278
## 3 Memoria 30 80.450 5.312 81.821 68.369 89.909 76.794 84.476
## Skew Kurtosis
## 1 -0.068 3.158
## 2 -0.235 2.387
## 3 -0.424 2.446
##
## $data
## # A tibble: 30 × 3
## Resistencia Tiempo_Reaccion Memoria
## <dbl> <dbl> <dbl>
## 1 38.5 251. 82.0
## 2 33.5 252. 81.6
## 3 32.6 240. 76.4
## 4 34.5 246. 83.9
## 5 34.0 246. 83.7
## 6 33.3 240. 73.8
## 7 36.5 247. 78.1
## 8 38.1 256. 84.7
## 9 38.3 255. 79.8
## 10 40.5 255. 75.9
## # ℹ 20 more rows
##
## $subset
## NULL
##
## $outlierMethod
## [1] "none"
##
## attr(,"class")
## [1] "mvn"d)
# ---- (d) Prueba T^2 de Hotelling ----
mu0 <- c(35, 250, 80) # vector bajo H0
d <- as.matrix(media - mu0) # diferencia
T2 <- n * t(d) %*% solve(S) %*% d # estadístico T^2
p <- ncol(datos)
Fstat <- ((n - p) / (p * (n - 1))) * T2 # conversión a F
gl1 <- p
gl2 <- n - p
pvalue <- 1 - pf(Fstat, df1 = gl1, df2 = gl2)
cat("T^2 =", T2, "\n")## T^2 = 1.610505cat("F =", Fstat, " con gl =", gl1, ",", gl2, "\n")## F = 0.4998119 con gl = 3 , 27cat("p-valor =", pvalue, "\n")## p-valor = 0.6855578# Valor crítico para comparación
Fcrit <- qf(1 - 0.05, df1 = gl1, df2 = gl2)
T2crit <- (p * (n - 1) / (n - p)) * Fcrit
cat("Valor critico F =", Fcrit, "\n")## Valor critico F = 2.960351cat("Valor critico T^2 =", T2crit, "\n")## Valor critico T^2 = 9.53891