[R] Multivariate multiple regression

[Iain Pardoe]

Thanks Henric.

If I may I'd like to go a little further ...  For example, Johnson and
Wichern's example 7.10:

> ex7.10 <-
+   data.frame(y1 = c(141.5, 168.9, 154.8, 146.5, 172.8, 160.1, 108.5),
+              y2 = c(301.8, 396.1, 328.2, 307.4, 362.4, 369.5, 229.1),
+              z1 = c(123.5, 146.1, 133.9, 128.5, 151.5, 136.2, 92),
+              z2 = c(2.108, 9.213, 1.905, .815, 1.061, 8.603, 1.125))
> attach(ex7.10)
> f.mlm <- lm(cbind(y1,y2)~z1+z2)
> y.hat <- c(1, 130, 7.5) %*% coef(f.mlm)
> round(y.hat, 2)
         y1     y2
[1,] 151.84 349.63
> qf.z <- t(c(1, 130, 7.5)) %*%
+   solve(t(cbind(1,z1,z2)) %*% cbind(1,z1,z2)) %*%
+   c(1, 130, 7.5)
> round(qf.z, 5)
        [,1]
[1,] 0.36995
> n.sigma.hat <- SSD(f.mlm)$SSD
> round(n.sigma.hat, 2)
     y1    y2
y1 5.80  5.22
y2 5.22 12.57
> F.quant <- qf(.95,2,3)
> round(F.quant, 2)
[1] 9.55

This gives me all the information I need to calculate a 95% confidence
ellipse for y=(y1,y2) at (z1,z2)=(130,7.5) using JW's equation (7-46):

(y-y.hat) %*% ((n-r-1) * solve(n.sigma.hat)) %*% t(y-y.hat)
<= qf.z * (m*(n-r-1)/(n-r-m)) * F.quant

But, what if instead I'd like to sample (y1,y2) values from this
distribution?  I can sample from an F(m,n-r-m) distribution easily
enough, but then how do I transform this to a single point in (y1,y2)
space?

Any ideas would be gratefully appreciated.  Thanks.