Multivariate multiple regression in R

Question

I have 2 dependent variables (DVs) each of whose score may be influenced by the set of 7 independent variables (IVs). DVs are continuous, while the set of IVs consists of a mix of continuous and binary coded variables. (In code below continuous variables are written in upper case letters and binary variables in lower case letters.)

The aim of the study is to uncover how these DVs are influenced by IVs variables. I proposed the following multivariate multiple regression (MMR) model:

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

To interpret the results I call two statements:

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

Outputs from both calls are pasted below and are significantly different. Can somebody please explain which statement among the two should be picked to properly summarize the results of MMR, and why? Any suggestion would be greatly appreciated.

Output using 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

Output using 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 

Answer 1

Briefly stated, this is because base-R's manova(lm()) uses sequential model comparisons for so-called Type I sum of squares, whereas car's Manova() by default uses model comparisons for Type II sum of squares.

I assume you're familiar with the model-comparison approach to ANOVA or regression analysis. This approach defines these tests by comparing a restricted model (corresponding to a null hypothesis) to an unrestricted model (corresponding to the alternative hypothesis). If you're not familiar with this idea, I recommend Maxwell & Delaney's excellent "Designing experiments and analyzing data" (2004).

For type I SS, the restricted model in a regression analysis for your first predictor c is the null-model which only uses the absolute term: lm(Y ~ 1), where Y in your case would be the multivariate DV defined by cbind(A, B). The unrestricted model then adds predictor c, i.e. lm(Y ~ c + 1).

For type II SS, the unrestricted model in a regression analysis for your first predictor c is the full model which includes all predictors except for their interactions, i.e., lm(Y ~ c + d + e + f + g + H + I). The restricted model removes predictor c from the unrestricted model, i.e., lm(Y ~ d + e + f + g + H + I).

Since both functions rely on different model comparisons, they lead to different results. The question which one is preferable is hard to answer - it really depends on your hypotheses.

What follows assumes you're familiar with how multivariate test statistics like the Pillai-Bartlett Trace are calculated based on the null-model, the full model, and the pair of restricted-unrestricted models. For brevity, I only consider predictors c and H, and only test for 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                                           

For comparison, the result from car's Manova() function using 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 ***

Now manually verify both results. Build the design matrix $X$ first and compare to R's design matrix.

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

Now define the orthogonal projection for the full model ($P_{f} = X (X'X)^{-1} X'$, using all predictors). This gives us the matrix $W = Y' (I-P_{f}) Y$.

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

Restricted and unrestricted models for SS type I plus their projections $P_{rI}$ and $P_{uI}$, leading to matrix $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

Restricted and unrestricted models for SS type II plus their projections $P_{rI}$ and $P_{uII}$, leading to matrix $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 trace for both types of SS: trace of $(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

Note that the calculations for the orthogonal projections mimic the mathematical formula, but are a bad idea numerically. One should really use QR-decompositions or SVD in combination with crossprod() instead.

Answer 2

Well, I still don't have enough points to comment on previous answer and thats why I am writing it as a separate answer, so please pardon me. (If possible please push me over the 50 rep points ;)

So here are the 2cents: Type I , II and III errors testing are essentially variations due to data being unbalanced. (Defn Unbalanced: Not having equal number of observations in each of the strata). If the data is balanced Type I , II and III error testing gives exact same results.

So what happens when the data is imbalanced?

Consider a model that includes two factors A and B; there are therefore two main effects, and an interaction, AB. SS(A, B, AB) indicates full model SS(A, B) indicates the model with no interaction. SS(B, AB) indicates the model that does not account for effects from factor A, and so on.

This notation now makes sense. Just keep it in mind.

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, also called "sequential" sum of squares:

1) SS(A) for factor A.

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

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

So we estimate main effect of A first them, effect of B given A, and then estimate interaction AB given A and B (This is where being imbalanced data, the differences kick in. As we estimate main effect first and then main of other and then interaction in a "sequence")

Type II:

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

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

Type II tests significance of main effect of A after B and B after A. Why is there no SS(AB | B, A) ? Caveat is that type II method can be used only when we have already tested for interaction to be insignificant. Given that there is no interaction (SS(AB | B, A) is insignificant) type II test has better power over type III

Type III:

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

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

So we tested for interaction during type II and interaction was significant. Now we need to use type III as it takes into account the interaction term.

As @caracal has said already, When data is balanced, the factors are orthogonal, and types I, II and III all give the same results. I hope this helps !

Disclosure: Most of it is not my own work. I found this excellent page linked and felt like boiling it down further to make it simpler.