Régression multiple multivariée dans R

J'ai 2 variables dépendantes (DV) dont chacune des notes peut être influencée par l'ensemble des 7 variables indépendantes (IV). Les DV sont continus, alors que l'ensemble des IV consiste en un mélange de variables codées continues et binaires. (Dans le code ci-dessous, les variables continues sont écrites en majuscules et les variables binaires en minuscules.)

Le but de l’étude est de découvrir comment ces variables sont influencées par les variables des solutions intraveineuses. J'ai proposé le modèle de régression multiple multivariée (MMR) suivant:

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

Pour interpréter les résultats, j'appelle deux déclarations:

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

Les sorties des deux appels sont collées ci-dessous et sont considérablement différentes. Quelqu'un peut-il s'il vous plaît expliquer quelle déclaration parmi les deux devrait être choisie pour résumer correctement les résultats du ROR, et pourquoi? Toute suggestion serait grandement appréciée.

Sortie utilisant summary(manova(my.model))statement:

> 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

Sortie utilisant Manova(my.model)statement:

> 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 

En bref, c'est parce que la base-R manova(lm())utilise des comparaisons de modèle séquentiel pour que l' on appelle somme de type I de places, alors que car« s Manova()par défaut utilise des comparaisons de modèles pour la somme de type II de carrés.

Je suppose que vous connaissez bien l'approche de comparaison de modèles à l'analyse de la variance ou à la régression. Cette approche définit ces tests en comparant un modèle restreint (correspondant à une hypothèse nulle) à un modèle sans restriction (correspondant à l'hypothèse alternative). Si vous n'êtes pas familier avec cette idée, je vous recommande l'excellent "Conception d'expériences et analyse de données" de Maxwell & Delaney (2004).

Pour le type SS de type I, le modèle restreint dans une analyse de régression pour votre premier prédicteur cest le modèle nul qui utilise uniquement le terme absolu:, lm(Y ~ 1)où, Ydans votre cas, correspond à la variable multidimensionnelle DV définie par cbind(A, B). Le modèle sans restriction ajoute ensuite un prédicteur c, c'est-à-dire lm(Y ~ c + 1).

Pour les SS de type II, le modèle sans restriction utilisé dans une analyse de régression pour votre premier prédicteur cest le modèle complet, qui inclut tous les prédicteurs à l'exception de leurs interactions, c'est-à-dire lm(Y ~ c + d + e + f + g + H + I). Le modèle restreint supprime le prédicteur cdu modèle sans restriction, c'est-à-dire lm(Y ~ d + e + f + g + H + I).

Étant donné que les deux fonctions reposent sur des comparaisons de modèles différentes, elles conduisent à des résultats différents. Il est difficile de répondre à la question à laquelle on préfère - cela dépend vraiment de vos hypothèses.

Ce qui suit suppose que vous sachiez comment les statistiques de test multivariées, telles que Pillai-Bartlett Trace, sont calculées en fonction du modèle null, du modèle complet et de la paire de modèles restreints sans restriction. Par souci de brièveté, je ne considère que les prédicteurs cet H, et ne teste que 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                                           

A titre de comparaison, le résultat de carla Manova()fonction de s utilisant SS type 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 ***

Maintenant, vérifiez manuellement les deux résultats. Construisez d'abord la matrice de conception et comparez-la à la matrice de conception de R.$X$

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

Définissez maintenant la projection orthogonale pour le modèle complet ( , en utilisant tous les prédicteurs). Cela nous donne la matrice .$P_{f} = X (X'X)^{-1} X'$$W = Y' (I-P_{f}) Y$

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

Restreintes et modèles sans restriction pour le type SS I plus leurs projections et , ce qui conduit à la matrice .$P_{rI}$$P_{uI}$$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

Restreints et modèles sans restriction pour le type SS II plus leurs projections et , conduisant à la matrice .$P_{rI}$$P_{uII}$$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-trace Bartlett pour les deux types de SS: trace 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

Notez que les calculs pour les projections orthogonales imitent la formule mathématique, mais sont une mauvaise idée numériquement. On devrait vraiment utiliser des décompositions QR ou SVD en combinaison avec crossprod().

Eh bien, je n’ai toujours pas assez de points pour commenter la réponse précédente et c’est pourquoi je l’écris en tant que réponse séparée, alors veuillez me pardonner. (Si possible, veuillez me repasser les 50 points de repère;)

Voici donc les 2cents: Les tests d’erreurs des types I, II et III sont essentiellement des variations dues au déséquilibre des données. (Defn Unbalanced: ne pas avoir le même nombre d'observations dans chacune des strates). Si les données sont équilibrées, le test d'erreur des types I, II et III donne exactement les mêmes résultats.

Alors que se passe-t-il lorsque les données sont déséquilibrées?

Considérons un modèle qui comprend deux facteurs A et B; il y a donc deux effets principaux et une interaction, AB. SS (A, B, AB) indique le modèle complet SS (A, B) indique le modèle sans interaction. SS (B, AB) indique le modèle qui ne tient pas compte des effets du facteur A, etc.

Cette notation a maintenant un sens. Gardez cela à l'esprit.

SS(AB | A, B) = SS(A, B, AB) - SS(A, B)

SS(A | B, AB) = SS(A, B, AB) - SS(B, AB)

SS(B | A, AB) = SS(A, B, AB) - SS(A, AB)

SS(A | B)     = SS(A, B) - SS(B)

SS(B | A)     = SS(A, B) - SS(A)

Type I, aussi appelé somme "séquentielle" de carrés:

1) SS(A) for factor A.

2) SS(B | A) for factor B.

3) SS(AB | B, A) for interaction AB.

Nous estimons donc l’effet principal de A d’abord, l’effet de B donné A, puis l’interaction AB donnée A et B (C’est là qu’il s’agit de données déséquilibrées, les différences s’inscrivent. puis interaction dans une "séquence")

Type II:

1) SS(A | B) for factor A.

2) SS(B | A) for factor B.

Le type II teste la signification de l'effet principal de A après B et B après A. Pourquoi n'y a-t-il pas de SS (AB | B, A)? La mise en garde est que la méthode de type II ne peut être utilisée que lorsque nous avons déjà testé que l'interaction est non significative. Etant donné qu'il n'y a pas d'interaction (SS (AB | B, A) n'est pas significatif), le test de type II a un meilleur pouvoir que celui de type III

Type III:

1) SS(A | B, AB) for factor A.

2) SS(B | A, AB) for factor B.

Nous avons donc testé l'interaction pendant le type II et l'interaction était significative. Nous devons maintenant utiliser le type III car il prend en compte le terme d'interaction.

Comme @caracal l'a déjà dit, lorsque les données sont équilibrées, les facteurs sont orthogonaux et les types I, II et III donnent tous les mêmes résultats. J'espère que ça aide !

Divulgation: La plupart de ce n'est pas mon propre travail. J'ai trouvé cette excellente page liée et j'ai eu l'impression de la réduire davantage pour la simplifier.