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.
Best Answer
I think you might do best to review multivariate regression (meaning no disrespect). There is a short tutorial at UVA's stats help page here: https://data.library.virginia.edu/getting-started-with-multivariate-multiple-regression/. They explain that multivariate regression is mostly the same as several univariate regressions, except that there are covariances between the betas for the different outcomes that need to be taken into account when testing the variables. In particular, they mention that the relevant plots are the same:
I will use their example to illustrate some data visualizations below (coded in
R). I'll start with a scatterplot matrix of the data.GENis binary, so I'll represent that with a different color and symbol. After fitting the model, if you try to runplot.lm()you'll get an error. However, it's easy enough to reproduce those plots manually. To plot a multiple regression model without interactions, you can pick a variable of interest and make a scatterplot with it and the response and draw the fitted function over it. Be sure to adjust the intercept by setting other variables at their means (see my answer to How to visualize a fitted multiple regression model?). You can also make scatterplots and qq-plots of the residuals (the latter lets you assess if their distribution is similar).