1
$\begingroup$

I'm trying to run a multivariate multiple regression in R, i.e. including multiple predictors and multiple outcome variables in the same linear regression model. Does anybody know how to pull out the coefficients and p-values for the relationship between each predictor and outcome pair in a multivariate multiple regression? I cannot seem to work out how to do that (I've been trying!).

Let me explain with an example dataset, if it helps:

#Create example dataset
df <- data.frame(pid=factor(501), y1=numeric(501), y2=numeric(501), 
y3=numeric(501), y4=numeric(501), x1=factor(501), x2=factor(501), 
x3=factor(501), x4=numeric(501), x5=numeric(501))
df$pid <- seq(1,501, by=1)
df$y1 <- seq(1,101, by=0.2)
df$y2 <- seq(401,201, by=-0.4)
df$y3 <- sqrt(rnorm(501, 7, 0.5))^3
df$x1 <- c(rep(c("sad","happy"), each=250), "sad")
df$x2 <- c(rep(c("human","vehicle","animal"), each=167))
df$x3 <- c(rep(seq(1,10, by=0.1), each=5), seq(1,46, by=1))
df$x4 <- rnorm(501, 3, .24)
df$x5 <- sqrt(rnorm(501, 23, 3.5))    

I then create the model using this:

#Specify the regression model
model <- lm(cbind(y1, y2, y3) ~ x1 + x2 + x3 + x4 + x5, data=df)

I can't simply use summary(lm) as doing so runs separate regressions without accounting for familywise error, nor does it account for the dependent variables possibly being correlated.

Reiterating my question. Does anybody know how to pull output so I can work out the coefficient and p-values but doing so in the same model? For example, I want to work out the coefficients and p-values of:

x1, x2, x3, x4 and x5 on y1
x1 and x2 on y2
x1, x2, x3, x4 and x5 on y3
... etc etc

I tried the car package:

modelanova <- car::Anova(model)
summary(modelanova)

However, I couldn't get it to break down to a particular outcome variable, it'd only produce coefficients overall (as if a composite outcome variable had been created)

Any ideas would be wonderful. I know I could run several univariate multiple regressions but I am particularly interested in running a single multivariate multiple regression.

$\endgroup$
2
$\begingroup$

The p-values are correct - they are the same as when running each of your $Y$s as a separate regression. Check out posts that have been made regarding the purpose of multivariate regression for more details, like this one, or this one, or even this one, but here's a good reference with examples:

https://socserv.socsci.mcmaster.ca/jfox/Books/Companion/appendix/Appendix-Multivariate-Linear-Models.pdf

Consider the the setup for linear regression, where $Y$ is $n\times m$ and $X$ is $n \times k$:

$Y=X\beta+\varepsilon$

$\beta=(X'X)^{-1}X'Y$

$\varepsilon=Y-X\beta=(I-X(X'X)^{-1}X')Y$

$\hat{\sigma}^2=\frac{\varepsilon'\varepsilon}{n-k}$

Notice that $\hat{\sigma}^2$ is $m \times m$, so its a scalar when $m=1$. To get standard errors of the betas in this case, the formula is:

$\hat{\sigma}_\beta=diagonal(\sqrt{\hat{\sigma}^2(X'X)^{-1}})$

But this doesn't work when $m>1$. To get the same standard errors as running $m$ univariate regressions, you could take the kronecker product like this:

$\hat{\sigma}_\beta=diagonal(\sqrt{\hat{\sigma}^2\otimes(X'X)^{-1}})$

But the advantage of doing multivariate regression is to be able to do combination and joint tests of the betas, like to test whether $\beta_{j,y_1}$ and $\beta_{j,y_2}$ are different from each other. For that, you'll need to do Multivariate Analysis of Variance (MANOVA) tests which require the extra information in $\varepsilon'\varepsilon$ and $Y'Y$ to perform.

$\endgroup$
2
+25
$\begingroup$

If you do not want a model that treats each outcome separately, maybe canonical correlation analysis would suit your situation better? Basically, it builds linear combinations of $Y$ and linear combinations of $X$ such that the correlation between these variates would be maximized. One can then test the significance of correlation associated with each variate, or induce sparsity to select the "most important" variables with LASSO, for example.

R implementations are available in base (well, stats, cancor) or several packages (CCA, or mixOmics for sparse CCA).

$\endgroup$
  • $\begingroup$ Thank you for suggesting this as an alternative approach. I am still hopeful to use multivariate multiple regression (not just for this problem in particular, but to determine if my suggestion is possible). $\endgroup$ – SimonsSchus May 15 '17 at 13:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.