R中的多元多元回歸

我有 2 個因變量 (DV)，每個變量的分數都可能受到 7 個自變量 (IV) 的影響。DV 是連續的，而 IV 集由連續和二進制編碼變量的混合組成。（在下面的代碼中，連續變量用大寫字母書寫，二進制變量用小寫字母書寫。）

該研究的目的是揭示這些 DVs 如何受 IVs 變量的影響。我提出了以下多元多元回歸 (MMR) 模型：

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

為了解釋結果，我稱之為兩個陳述：

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

兩個調用的輸出都粘貼在下面，並且有很大不同。有人可以解釋一下應該選擇兩者中的哪一個來正確總結 MMR 的結果，為什麼？任何建議將不勝感激。

使用summary(manova(my.model))語句輸出：

> 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

使用Manova(my.model)語句輸出：

> 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 

簡而言之，這是因為 base-Rmanova(lm())對所謂的 I 類平方和使用順序模型比較，而默認情況下，car對Manova()II 類平方和使用模型比較。

我假設您熟悉方差分析或回歸分析的模型比較方法。這種方法通過將受限模型（對應於零假設）與非受限模型（對應於備擇假設）進行比較來定義這些檢驗。如果你不熟悉這個想法，我推薦 Maxwell & Delaney 的優秀“設計實驗和分析數據”（2004 年）。

對於 I 型 SS，您的第一個預測變量的回歸分析中的受限模型c是僅使用絕對項的空模型：lm(Y ~ 1)，Y在您的情況下，將是由 定義的多元 DV cbind(A, B)。不受限制的模型然後添加預測器c，即lm(Y ~ c + 1)。

對於 II 型 SS，您的第一個預測變量的回歸分析中的無限制模型c是完整模型，其中包括除了它們的交互作用之外的所有預測變量，即lm(Y ~ c + d + e + f + g + H + I). 受限模型從非受限模型中移除預測變量c，即lm(Y ~ d + e + f + g + H + I)。

由於這兩個函數都依賴於不同的模型比較，因此它們會導致不同的結果。哪個更可取的問題很難回答——這實際上取決於你的假設。

接下來的內容假設您熟悉如何基於空模型、完整模型和一對受限-非受限模型計算多變量檢驗統計量（如 Pillai-Bartlett Trace）。為簡潔起見，我只考慮預測變量c和H，並且只測試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 

為了比較，使用 SS 類型 IIcar的函數的結果。Manova()

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 ***

現在手動驗證這兩個結果。構建設計矩陣首先與 R 的設計矩陣進行比較。

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

現在定義完整模型的正交投影（，使用所有預測變量）。這給了我們矩陣.

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

SS 類型 I 的受限和非受限模型及其預測和，導致矩陣.

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

SS II 型的受限和非受限模型及其預測和，導致矩陣.

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

兩種類型 SS 的 Pillai-Bartlett 跡線：.

(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

請注意，正交投影的計算模擬了數學公式，但在數值上是一個壞主意。應該真正結合使用 QR 分解或 SVD crossprod()。

引用自：https://stats.stackexchange.com/questions/11127