# Multivariate multiple linear regression model

I have 2 response variables (Y1, Y2) and some independent variables. I need to predict both Y1 and Y2 using the same set of predictors.In other words I need to fit the model:
$Y=XB+E$
where:
$Y$ has $n$ rows and $q$ columns,
$X$ has $n$ rows and $p$ columns,
$B$ has $p$ rows and $q$ columns,
$E$ has $n$ rows and $q$ columns.
In this case, $q=2$.

One way to do this in R is by calling the lm() function:

m1=lm(cbind(Y1,Y2)~x1+x2+x3)


This is equivalent to fitting two univariate multiple regression models, one for Y1, and the other for Y2:

m1.1=lm(Y1~x1+x2+x3)
m1.2=lm(Y2~x1+x2+x3)


How can I fit a multivariate multiple linear regression model for both Y1 and Y2 which permits correlation between Y1 and Y2?
In this type of problem you have two covariance matrices involved, one that is $n \times n$ over the rows of $Y$, and another that is $q \times q$ over the columns of $Y$. If the $q \times q$ covariance matrix is the identity matrix, then you're right, your solution for $B$ is equivalent to running separate regressions.