1 Objetivos de aprendizaje

Formular y ajustar un modelo lineal multivariado (MLM): \(\mathbf{Y} = \mathbf{X}B + E\).
Realizar pruebas globales (Wilks) y pruebas por predictor (Type I vs Type II).
Interpretar coeficientes por respuesta y pruebas univariadas auxiliares.
Evaluar supuestos mediante análisis de residuales y normalidad multivariada (Mardia).
Comunicar hallazgos con una discusión integrada.

2 Paquetes y datos

Teoría: car provee MANOVA/ANOVA multivariada para objetos mlm y MVN implementa pruebas de normalidad multivariada coherentes con el supuesto \(E\sim \mathcal{MN}\).

library(car)   # Anova (Type II/III) y métodos para 'mlm'
library(MVN)   # Normalidad multivariada (Mardia, HZ, Royston)
data(iris)     # Conjunto clásico de medidas en cm
str(iris)

## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

3 Especificación del modelo

Modelo:

Respuestas (p=2): \(Y_1=\texttt{Sepal.Length}\), \(Y_2=\texttt{Sepal.Width}\)
Predictores: \(X_1=\texttt{Petal.Length}\), \(X_2=\texttt{Petal.Width}\)

``` \(\begin{bmatrix} \text{Sepal.Length} & \text{Sepal.Width} \end{bmatrix} ~ \sim ~ \text{Petal.Length} + \text{Petal.Width}\) ```

Tarea A (antes de correr): prediga el signo esperado de cada coeficiente, justificando con argumentos biológicos o de correlación simple.

Predicción de signos (anote aquí):
- Efecto de Petal.Length sobre Sepal.Length: ____
- Efecto de Petal.Width sobre Sepal.Length: ____
- Efecto de Petal.Length sobre Sepal.Width: ____
- Efecto de Petal.Width sobre Sepal.Width: ____

3.1 Ajuste del modelo multivariado

fit1 <- lm(cbind(Sepal.Length, Sepal.Width) ~ Petal.Length + Petal.Width, data = iris)
coef(fit1)

##              Sepal.Length Sepal.Width
## (Intercept)     4.1905824   3.5870492
## Petal.Length    0.5417772  -0.2571378
## Petal.Width    -0.3195506   0.3640421

Guía: La matriz de coeficientes \(\hat B\) reporta, por columna, la ecuación de cada respuesta. Interprete magnitud y signo controlando por el otro predictor.

Anote su interpretación breve (2–3 líneas por respuesta):
- Ecuación Sepal.Length:
- Ecuación Sepal.Width:

4 Pruebas globales multivariadas

4.1 Type I (secuencial) con `manova`

summary(manova(fit1), test = "Wilks")

##               Df   Wilks approx F num Df den Df    Pr(>F)    
## Petal.Length   1 0.11668   552.65      2    146 < 2.2e-16 ***
## Petal.Width    1 0.85265    12.62      2    146 8.836e-06 ***
## Residuals    147                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Teoría: Wilks \(\Lambda\) compara variación explicada vs residual; valores pequeños ⇒ mayor efecto conjunto. Type I depende del orden de los términos en la fórmula.

Preguntas rápidas:
1) ¿Rechaza \(H_0: B_{\text{Petal.Length}}=0\)? ¿Y para Petal.Width?
2) ¿Por qué los resultados de Type I pueden cambiar si reordenas los predictores?

4.2 Type II (marginal) con `car::Anova`

Anova(fit1, test = "Wilks")  # por defecto, Type II

## 
## Type II MANOVA Tests: Wilks test statistic
##              Df test stat approx F num Df den Df    Pr(>F)    
## Petal.Length  1   0.43886   93.342      2    146 < 2.2e-16 ***
## Petal.Width   1   0.85265   12.616      2    146 8.836e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Teoría: Type II prueba cada término ajustando por el resto (sin interacciones no incluidas). Suele ser preferible cuando el orden no debe influir.

Compara Type I vs Type II: ¿Qué cambia en los tamaños de efecto (\(\Lambda\), F)? ¿Qué predictor muestra mayor efecto conjunto ajustado?

5 Pruebas univariadas auxiliares (por respuesta)

sum_SL <- summary(lm(Sepal.Length ~ Petal.Length + Petal.Width, data = iris))
sum_SW <- summary(lm(Sepal.Width  ~ Petal.Length + Petal.Width, data = iris))
sum_SL; sum_SW

## 
## Call:
## lm(formula = Sepal.Length ~ Petal.Length + Petal.Width, data = iris)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.18534 -0.29838 -0.02763  0.28925  1.02320 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.19058    0.09705  43.181  < 2e-16 ***
## Petal.Length  0.54178    0.06928   7.820 9.41e-13 ***
## Petal.Width  -0.31955    0.16045  -1.992   0.0483 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4031 on 147 degrees of freedom
## Multiple R-squared:  0.7663, Adjusted R-squared:  0.7631 
## F-statistic:   241 on 2 and 147 DF,  p-value: < 2.2e-16

## 
## Call:
## lm(formula = Sepal.Width ~ Petal.Length + Petal.Width, data = iris)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.06198 -0.23389  0.01982  0.20580  1.13488 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.58705    0.09373  38.272  < 2e-16 ***
## Petal.Length -0.25714    0.06691  -3.843  0.00018 ***
## Petal.Width   0.36404    0.15496   2.349  0.02014 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3893 on 147 degrees of freedom
## Multiple R-squared:  0.2131, Adjusted R-squared:  0.2024 
## F-statistic:  19.9 on 2 and 147 DF,  p-value: 2.238e-08

Guía: Útiles para dirección/magnitud por respuesta, pero no sustituyen a la prueba multivariada. Relacione las R² y los signos con su predicción inicial (Tarea A).

Reflexión: ¿Por qué Sepal.Length tiene R² mucho mayor que Sepal.Width? Proponga una hipótesis.

6 Diagnóstico de supuestos

6.1 Gráficos de residuales

res1 <- residuals(fit1)       # n x 2
op <- par(mfrow = c(2,2))
qqnorm(res1[,1]); qqline(res1[,1], lwd = 2)
qqnorm(res1[,2]); qqline(res1[,2], lwd = 2)
plot(fitted(fit1)[,1], res1[,1], xlab="Ajustados Sepal.Length", ylab="Residuales")
plot(fitted(fit1)[,2], res1[,2], xlab="Ajustados Sepal.Width",  ylab="Residuales")

par(op)

Qué buscar:
- QQ-plots ~ recta ⇒ normalidad marginal razonable.
- Res vs Ajustados ~ nube sin forma ⇒ linealidad y varianza constante.

Checklist:
- [ ] Sin curvatura sistemática
- [ ] Sin embudo pronunciado (heterocedasticidad)
- [ ] Pocos atípicos extremos

6.2 Normalidad multivariada (Mardia)

MVN::mvn(as.data.frame(res1), mvnTest = "mardia", multivariatePlot = "qq")

## $multivariateNormality
##              Test         Statistic           p value Result
## 1 Mardia Skewness 0.908428482782596 0.923348514485471    YES
## 2 Mardia Kurtosis 0.255111395155304 0.798637063234119    YES
## 3             MVN              <NA>              <NA>    YES
## 
## $univariateNormality
##               Test     Variable Statistic   p value Normality
## 1 Anderson-Darling Sepal.Length    0.2631    0.6964    YES   
## 2 Anderson-Darling Sepal.Width     0.6786    0.0749    YES   
## 
## $Descriptives
##                n         Mean   Std.Dev      Median       Min      Max
## Sepal.Length 150 1.341285e-17 0.4003412 -0.02762794 -1.185344 1.023195
## Sepal.Width  150 1.501042e-17 0.3866448  0.01981651 -1.061983 1.134881
##                    25th      75th       Skew   Kurtosis
## Sepal.Length -0.2983839 0.2892468 0.05922405 -0.3047638
## Sepal.Width  -0.2338942 0.2058010 0.08376038  0.7534951

Interpretación: p-valores grandes ⇒ sin evidencia contra normalidad multivariada. Si fallara, considere transformaciones o términos no lineales/interacciones.

7 Discusión integrada (esqueleto de informe)

Hallazgo 1 (Efecto conjunto): …
Hallazgo 2 (Por respuesta): …
Supuestos (diagnóstico): …
Limitaciones: correlación entre predictores, posible multicolinealidad, escala de medidas.
Extensiones: agregar interacción Petal.Length:Petal.Width, o incluir Species como factor (MANOVA con covariables).

8 Extensión opcional (para crédito extra)

Calcule la distancia de Mahalanobis de los residuales para detectar observaciones potencialmente influyentes:
\(MD^2_i = e_i^\top \hat\Sigma^{-1} e_i\).

Sigma_hat <- cov(res1)
md2 <- mahalanobis(res1, center = c(0,0), cov = Sigma_hat)
cut <- qchisq(0.975, df = ncol(res1))
which(md2 > cut)
hist(md2, breaks = 15, main = "MD² residuales", xlab = "MD²")

8.1 Rúbrica breve (auto-evaluación)

Interpreta \(\hat B\) correctamente por respuesta.
Distingue Type I vs Type II y justifica su elección.
Integra resultados multivariados y univariados sin contradicciones lógicas.
Evalúa supuestos con evidencia gráfica y pruebas.
Redacta una discusión clara y breve (8–12 líneas).

Laboratorio: Regresión Lineal Multivariada (Cap. 7) — iris

Guía de trabajo

2025-08-21

1 Objetivos de aprendizaje

2 Paquetes y datos

3 Especificación del modelo

3.1 Ajuste del modelo multivariado

4 Pruebas globales multivariadas

4.1 Type I (secuencial) con `manova`

4.2 Type II (marginal) con `car::Anova`

5 Pruebas univariadas auxiliares (por respuesta)

6 Diagnóstico de supuestos

6.1 Gráficos de residuales

6.2 Normalidad multivariada (Mardia)

7 Discusión integrada (esqueleto de informe)

8 Extensión opcional (para crédito extra)

8.1 Rúbrica breve (auto-evaluación)

Laboratorio: Regresión Lineal Multivariada (Cap. 7) — iris

Guía de trabajo

2025-08-21

1 Objetivos de aprendizaje

2 Paquetes y datos

3 Especificación del modelo

3.1 Ajuste del modelo multivariado

4 Pruebas globales multivariadas

4.1 Type I (secuencial) con manova

4.2 Type II (marginal) con car::Anova

5 Pruebas univariadas auxiliares (por respuesta)

6 Diagnóstico de supuestos

6.1 Gráficos de residuales

6.2 Normalidad multivariada (Mardia)

7 Discusión integrada (esqueleto de informe)

8 Extensión opcional (para crédito extra)

8.1 Rúbrica breve (auto-evaluación)

4.1 Type I (secuencial) con `manova`

4.2 Type II (marginal) con `car::Anova`