How to do multivariate regression in R? I need to preform a multivariate normal regression in R. The question is: 

Let $Y_1$, $Y_2$, and $Y_3$ follows multivariate normal distribution. What is 
  
  
*
  
*the conditional of $Y_3$ given $Y_1$ and $Y_2$
  
*the conditional of $Y_2$ given $Y_1$
From these two, derive:
  
*the joint distribution of $Y_3$ and $Y_2$ given $Y_1$. 

Now suppose you have a sample of size $n$ from the multivariate normal distribution. Do the two regressions in (1) and (2). How can I combine them to get (3), the regression of $Y_3$ and $Y_2$ on $Y_1$?
library(mvtnorm)
mu  <- c(1,2,3)
Sig <- matrix(c(4,2,1,2,4,-1,1,-1,4), nrow=3, ncol=3) 
Y   <- rmvnorm(20, mean=mu, sigma=Sig) #generate multivariate normal distribution
y3  <- lm(Y[,3]~Y[,1] + Y[,2]) 
y2  <- lm(Y[,2]~Y[,1])

 A: To follow up on the answer above, if you are interested in the multivariate regression of $(y_2, y_3)$ on $y_1$, we can show using properties of the multivariate normal distribution that this is equivalent to the univariate regressions of $y_2$ on $y_1$ and of $y_3$ on $y_1$. 
Let $Y = (y_1, y_2, y_3)^T \sim N\left((\mu_1, \mu_2, \mu_3)^T,\left(\begin{array}{ccc}\sigma_1^2 & \sigma_{12} & \sigma_{13} \\ \sigma_{12} & \sigma_2^2 & \sigma_{23} \\ \sigma_{13} & \sigma_{23} & \sigma_3^2 \end{array} \right) \right)$. 
By properties of the multivariate normal distribution, the joint regression of $y_2, y_3$ on $y_1$ is $E((y_2,y_3)^T|y_1) = (\mu_2,\mu_3)^T +  (\sigma_{12}, \sigma_{13})^T \frac{y_1 - \mu_1}{\sigma_1^2}$. After a little additional algebra, we can show that this reduces to a 2x1 vector of the univariate regressions of $y_3$ and $y_2$ respectively on $y_1$. 
Therefore as shown above by the previous answer, it is sufficient to regress each of $y_3$ and $y_2$ on $y_1$ individually to get the joint regression. 
