Multivariate linear regression

● Interested in the relationship between random q-vector [Y1 , . . . , Yp]T  and random q-vector [X1 , . . . , Xq ]T

● Model

_ Population  version:    [Y1 , . . . , Yp ]T     |   X1 , . . . , Xq     ~   (BT [1, X1 , . . . , Xq ]T , σ 2 ),  where  B   = [βkj ](q+1)xp, i.e.,

*  E([Y1 , . . . , Yp ]T  | X1 , . . . , Xq ) = BT [1, X1 , . . . , Xq ]T

* cov ([Y1 , . . . , Yp]T  | X1 , . . . , Xq ) = 2 > 0, i.e., the conditional covariance of [Y1 , . . . , Yp]T  does not depend on X1 , . . . , Xq

v                x               B               E

n x p = n x (q + 1) (q + 1) x p + n x p

* v = [Yij ]nxp

*  Design matrix

x =  ┌'  1..     ┐'

'  1   Xn1     . . .    Xnq   'nx(q+1)

' rk(x) = q + 1 < p + q + 1 < n

*  E = [E1 . , . . . , En . ]T , where Ei(T)  is the ith row of E

*  Assume the independence across i , i.e.,

'  [Yi1, . . . , Yip , Xi1, . . . , Xiq ]T    [Y1 , . . . Yp , X1 , . . . , Xq ]T

' Ei .    (0p , 2)

●  Relationship with univariate linear regression

_  If 2 is diagonal, the multivariate model reduces to v.j  = xB.j  + E.j , j = 1, . . . , p

*  v.j : the jth column of v

*  B.j : the jth column of B

*  E.j  ~ (0n , σj(2)j Ⅰn )

'  σj(2)j :  (j, j)-entry of 2

●  Relationship with MANOVA

_ MANOVA models can be expressed as multivariate linear regression with carefully selected dummy (explanatory) variables.

Exercise: translate the following 1-way MANOVA model

vij  = ↓ + |i + Eij ,    i = 1, . . . , m,    j = 1, . . . , ni

into a multivariate linear regression model, where Eij   MVNp (0, 2) and     i |i  = 0.

● Least squares (LS) estimation (no need of (conditional) normality) _ LS = (xT x)_ 1 xTv

*  E(LS ) = B

_ LS = (n _ q _ 1)_ 1 (v _ xLS )T (v _ xLS ) = (n _ q _ 1)_ 1vT (Ⅰ _ 工)v * E( LS ) = 2

●  Maximum likelihood (ML) estimation (in need of (conditional) normality)

_ ML = (xT x)_ 1 xTv = LS

_  ML  = n_ 1 vT (Ⅰ _ 工)v = n_ 1 (n _ q _ 1) LS

*  Given x, n ML  ~ Wp (2, n _ q _ 1)

● Inference (in need of (conditional) normality) _ Inference on BT a, given a e 政q+1

*  Estimator M(T)L a

*  100(1 _ α)% confidence region for BT a

╱u e 政p  : (u _ M(T)L a)T L(_)S(1)(u _ M(T)L a) < aT (xT x)_ 1 a F1_α,p,n _p _q 、

_ Inference on v0 = BT x0 + E0  with a new observation vector x0 = [1, X01 , . . . , X0q]T  e 政q+1

*  Prediction vˆ 0  = BM(T)L x0

*  100(1 _ α)% prediction region for v0

╱u e 政p  : (u _ vˆ 0 )T L(_)S(1)(u _ vˆ 0 ) < (1 + x0(T)(xT x)_ 1 x0 } F1_α,p,n _p _q 、

_ Inference on aTv0  = aT BT x0    a e 政p  and a new observation vector x0  = [1, X01 , . . . , X0q]T  e 政q+1

*  Prediction aT vˆ 0  = aT BM(T)L x0

*  100(1 _ α)% Scheffé’s simultaneous prediction interval for aT v0

aTvˆ 0 土 入aT LS a(1 + x0(T)(xT x)_ 1 x0 } F1_α,p,n _p _q

install.packages(c ( !ellipse !))

options(digits  =  4)

tear  <-  c (

6.5 ,  6.2 ,  5.8 ,  6.5 ,  6.5 ,  6.9 ,  7.2 ,  6.9 ,  6.1 ,  6.3 ,

6.7 ,  6.6 ,  7.2 ,  7.1 ,  6.8 ,  7.1 ,  7.0 ,  7.2 ,  7.5 ,  7.6

)

gloss  <-  c (

9.5 ,  9.9 ,  9.6 ,  9.6 ,  9.2 ,  9.1 ,  10.0 ,  9.9 ,  9.5 ,  9.4 ,

9.1 ,  9.3 ,  8.3 ,  8.4 ,  8.5 ,  9.2 ,  8.8 ,  9.7 ,  10.1 ,  9.2

)

(plastic  <-  cbind(tear,  gloss))

rate  <-  factor(gl(2 , 10 ,length=nrow (plastic)),  labels=c ( "Low" ,  "High"))

# Model  fitting

fit  <-  lm(cbind(tear,  gloss)  ~  rate)

summary (fit)

# Prediction

(Obs_new  <-  data.frame(rate  =  factor(c ("High"),  levels  =  c ( "Low" ,  "High")))) (prediction  <-  t(predict(fit,  newdata  =  Obs_new)))

# Prediction  region

n  =  nrow (model.matrix(fit))