Multivariado regresión múltiple en R

Question

Tengo 2 variables dependientes (DVs), cada una de cuyas partituras pueden ser influenciados por el conjunto de las 7 variables independientes (IVs). DVs son continuos, mientras que el conjunto de IVs se compone de una mezcla de continuo y con codificación binaria de las variables. (En el código a continuación las variables continuas se han escrito en letras mayúsculas y variables binarias en letras minúsculas.)

El objetivo del estudio es descubrir cómo estos DVs son influenciados por IVs variables. Me propuso el siguiente multivariante de regresión múltiple (MMR) modelo:

my.model <- lm(cbind(A, B) ~ c + d + e + f + g + H + I)

Para interpretar los resultados que llamar a dos declaraciones:

1. summary(manova(my.model))
2. Manova(my.model)

Los resultados de ambas llamadas se pegan a continuación y son significativamente diferentes. Puede alguien por favor explicar que la declaración de entre los dos deben escogerse adecuadamente resumir los resultados de la vacuna MMR, y por qué? Cualquier sugerencia será muy apreciada.

Salida usando summary(manova(my.model)) declaración:

> summary(manova(my.model))
           Df   Pillai approx F num Df den Df    Pr(>F)    
c           1 0.105295   5.8255      2     99  0.004057 ** 
d           1 0.085131   4.6061      2     99  0.012225 *  
e           1 0.007886   0.3935      2     99  0.675773    
f           1 0.036121   1.8550      2     99  0.161854    
g           1 0.002103   0.1043      2     99  0.901049    
H           1 0.228766  14.6828      2     99 2.605e-06 ***
I           1 0.011752   0.5887      2     99  0.556999    
Residuals 100                                              
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Salida usando Manova(my.model) declaración:

> library(car)
> Manova(my.model)

Type II MANOVA Tests: Pillai test statistic
  Df test stat approx F num Df den Df    Pr(>F)    
c  1  0.030928   1.5798      2     99   0.21117    
d  1  0.079422   4.2706      2     99   0.01663 *  
e  1  0.003067   0.1523      2     99   0.85893    
f  1  0.029812   1.5210      2     99   0.22355    
g  1  0.004331   0.2153      2     99   0.80668    
H  1  0.229303  14.7276      2     99 2.516e-06 ***
I  1  0.011752   0.5887      2     99   0.55700    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Answer 1

Brevemente dicho, esto es debido a que la base-R manova(lm()) utiliza secuencial modelo de comparaciones para los llamados de Tipo I de la suma de los cuadrados, mientras que car's Manova() por defecto utiliza comparaciones del modelo para el Tipo II de la suma de los cuadrados.

Supongo que usted está familiarizado con el modelo de comparación de enfoque ANOVA o análisis de regresión. Este enfoque define estas pruebas mediante la comparación de un modelo restringido (correspondiente a una hipótesis nula) a un modelo sin restricciones (correspondiente a la hipótesis alternativa). Si usted no está familiarizado con esta idea, recomiendo Maxwell Y Delaney excelente "Diseño de experimentos y análisis de datos" (2004).

Para el tipo I SS, el modelo restringido en un análisis de regresión para su primer predictor c es el null-modelo que utiliza solamente el absoluto plazo: lm(Y ~ 1), donde Y sería en su caso, el multivariante DV definido por cbind(A, B). El modelo sin restricciones, a continuación, agrega predictor c, es decir, lm(Y ~ c + 1).

Para el tipo II SS, el modelo sin restricciones en un análisis de regresión para su primer predictor c es el modelo completo que incluye todos los predictores, excepto para sus interacciones, es decir, lm(Y ~ c + d + e + f + g + H + I). El modelo restringido quita predictor c de la modelo sin restricciones, es decir, lm(Y ~ d + e + f + g + H + I).

Puesto que ambas funciones se basan en diferentes comparaciones del modelo, que conducen a diferentes resultados. La pregunta que uno es preferible es difícil de responder, realmente depende de la validez de sus hipótesis.

Lo que sigue se supone que usted está familiarizado con la forma multivariante de la estadística de prueba como la Pillai-Bartlett Traza se calcula con base en el null-modelo, el modelo completo, y el par de restringido no restringido de modelos. Por razones de brevedad, sólo he considerar predictores c y H, y sólo prueba para c.

> N <- 100                             # generate some data: number of subjects
> c <- rbinom(N, 1, 0.2)               # dichotomous predictor c
> H <- rnorm(N, -10, 2)                # metric predictor H
> A <- -1.4*c + 0.6*H + rnorm(N, 0, 3) # DV A
> B <-  1.4*c - 0.6*H + rnorm(N, 0, 3) # DV B
> Y <- cbind(A, B)                     # DV matrix
> my.model <- lm(Y ~ c + H)            # the multivariate model
> summary(manova(my.model))            # from base-R: SS type I
          Df  Pillai approx F num Df den Df  Pr(>F)    
c          1 0.06835   3.5213      2     96 0.03344 *  
H          1 0.32664  23.2842      2     96 5.7e-09 ***
Residuals 97                                           

Para la comparación, el resultado de car's Manova() función utilizando SS tipo II.

> library(car)                           # for Manova()
> Manova(my.model, type="II")
Type II MANOVA Tests: Pillai test statistic
  Df test stat approx F num Df den Df  Pr(>F)    
c  1   0.05904   3.0119      2     96 0.05387 .  
H  1   0.32664  23.2842      2     96 5.7e-09 ***

Ahora manualmente verificar los resultados. Construir el diseño de la matriz de $X$ primero y comparar a R del diseño de la matriz.

> X  <- cbind(1, c, H)
> XR <- model.matrix(~ c + H)
> all.equal(X, XR, check.attributes=FALSE)
[1] TRUE

Ahora definir la proyección ortogonal para el modelo completo ($P_{f} = X (X X)^{-1} X'$, utilizando todos los predictores). Esto nos da la matriz $W = Y' (I-P_{f})$Y.

> Pf  <- X %*% solve(t(X) %*% X) %*% t(X)
> Id  <- diag(N)
> WW  <- t(Y) %*% (Id - Pf) %*% Y

Con y sin restricciones de los modelos para SS tipo I, además de sus proyecciones de $P_{rI}$ y $P_{uI}$, que conduce a la matriz $B_{I} = Y' (P_{uI} - P_{PrI})$Y.

> XrI <- X[ , 1]
> PrI <- XrI %*% solve(t(XrI) %*% XrI) %*% t(XrI)
> XuI <- X[ , c(1, 2)]
> PuI <- XuI %*% solve(t(XuI) %*% XuI) %*% t(XuI)
> Bi  <- t(Y) %*% (PuI - PrI) %*% Y

Con y sin restricciones de los modelos de SS de tipo II, además de sus proyecciones de $P_{rI}$ y $P_{uII}$, que conduce a la matriz $B_{II} = Y' (P_{uII} - P_{PrII})$Y.

> XrII <- X[ , -2]
> PrII <- XrII %*% solve(t(XrII) %*% XrII) %*% t(XrII)
> PuII <- Pf
> Bii  <- t(Y) %*% (PuII - PrII) %*% Y

Pillai-Bartlett de seguimiento para ambos tipos de SS: seguimiento de $(B + W)^{-1} B$.

> (PBTi  <- sum(diag(solve(Bi  + WW) %*% Bi)))   # SS type I
[1] 0.0683467

> (PBTii <- sum(diag(solve(Bii + WW) %*% Bii)))  # SS type II
[1] 0.05904288

Tenga en cuenta que los cálculos de las proyecciones ortogonales de imitar la fórmula matemática, sino que son una mala idea numéricamente. Uno realmente debe usar QR-descomposición o de enfermedad vesicular porcina en combinación con crossprod() lugar.