La revue canadienne de statistique Multivariate flattening for better predictions

The Canadian Journal of StatisticsVol.28,No.2,2000,Pages???-???La revue canadienne de statistiqueMultivariateflattening for better pre-dictionsMartin BILODEAUUniversit´e de Montr´e alKey words and phrases:Multivariateflattening,multivariate regression,prediction, risk.Mathematics subject classification codes(1991):primary62J07;secondary62H12. ABSTRACTMultivariate regression with p responses as opposed to p multiple regressions is getting increasingly more attention,especially in the context of prediction.Multivariateflat-tening methods are investigated as a way to obtain improved predictions over ordinary least-squares.With respect to sum of squares of prediction error,or SPE risk,an unbi-ased estimate of the risk is derived for two recent prediction methods,OPT and GCV, proposed by Breiman&Friedman(1997).Expressions for the exact SPE risk of OPT and GCV are derived generally for p≥1and evaluated numerically for p=1.R´ESUM´ELes m´e thodes de pr´e vision de p variables utilisant la r´e gression multivari´e e,plutˆo t que p r´e gression multiples,sont de plus en plus utilis´e es.Le r´e tr´e cissement multivari´e est´e tudi´e comme m´e thode de pr´e vision pour am´e liorer la m´e thode des moindres carr´e s.On obtient l’estimateur sans biais du risque SPE de deux m´e thodes r´e centes de pr´e vision,la m´e thode OPT et celle du GCV de Breiman et Friedman(1997).Une expression pour le risque exact de OPT et GCV est obtenue lorsque p≥1et une´e valuation num´e rique en est faite quand p=1.1.INTRODUCTIONThe problem of predicting several,possibly correlated,responses from the same set of predictors is becoming increasingly important.Applications by Breiman& Friedman(1997)(denoted hereafter by BF)include prediction of changes in the valuations of stocks in60industry groups by using over100econometric variables as predictors.Or,in chemometrics the prediction of6output characteristics of the polymers produced as predicted by22predictor variables.The model is the common multivariate regression model with p responses and1k predictors on n independent observations,Y=XB+E,where the error terms are assumed to be normally distributed,E∼N n p(0,I⊗Σ),Σ>0(positive definite),and X is an n×k matrix of rank k.The design matrix X is usually random in which case we assume that it is independently distributedof the error terms.Let|=E.=denote mutual independence so that X| The ordinary least-squares estimate of B isˆB OLS=(X X)−1X Y=(ˆb1,...,ˆb p),whereˆb i=(X X)−1X y i.Prediction of a“future”observation of p responses,y, at a given value of the predictor,x,assuming the same modely=B x+eis still valid,based on OLS is done in the obvious wayˆy OLS=ˆB OLS x=(ˆb 1x,...,ˆb p x) .We assume that“future”is independent of the“past”;i.e.,(x,e)|=(X,E).So,prediction of y1,...,y p is done independently via each of the p multiple regres-sions.This paper looks into improved predictions,˜y,by multivariate“flattening”techniques˜y =ˆy OLS A,x ˜B=x ˆB OLS A,which means prediction of the i th response is accomplished via a linear combina-tion,through the i th column of A:p×p,of all the p response predictions.To evaluate the performance of a prediction we evaluate an expected loss func-tion(risk).Consider the prediction of y by˜y=˜B x for any estimate˜B of B with respect to SPER(sum of squares of prediction error when the independent variable is random)motivated by Stein(1960)E(y−˜y) Σ−1(y−˜y)=p+E tr(˜B−B)Σ−1(˜B−B) Ω,whereΩ=E xx .SPER is obtained on subtracting p from the aboveR SPER(˜B)=E tr(˜B−B)Σ−1(˜B−B) Ω.A closely related risk function denoted SPE by Dempster et al.(1976)isR SPE(˜B)=E tr(˜B−B)Σ−1(˜B−B) X X.The SPE risk was investigated by van der Merwe&Zidek(1980)and Bilodeau& Kariya(1989).It is the most convenient risk for mathematical reasons.A canonical form of this problem which simplify notations and derivations is as follows.LetV=Y [I−X(X X)−1X]Y2be the residuals sum of squares.Then,given X ,ˆB OLS ∼N k p (B ,(X X )−1⊗Σ),V ∼W p (n −k,Σ),ˆBOLS |=V .The problem in its canonical form will be to estimate B ,with respect to SPE riskR SPE (˜B)=E tr (˜B −B )Σ−1(˜B −B ) C −1,based on Z and V whose conditional distribution given C :k ×k isZ ∼N k p (B ,C ⊗Σ),V∼W p (m,Σ),Z|=V .Results in terms of the canonical problem can be transposed to the original formu-lation with the correspondanceZ ↔ˆBOLS ,C ↔(X X )−1,m ↔(n −k ),V ↔V ,k ↔k,p ↔p.Van der Merwe &Zidek (1980)proposed FICYREGˆB FICYREG =Z I −r (Z C −1Z )−1V,where r =(k −p −1)/(m +p +1)and establisedR SPE (ˆBFICYREG )<R SPE (Z )=pk,or,FICYREG improves uniformly over OLS.Bilodeau &Kariya (1989)producedmany other improved predictions and showed the Efron-Morris (EM)estimateˆBEM =ˆB FICYREG −s V /tr(Z C −1Z ),s =(p −1)/(m +p +1),improves uniformly over FICYREG.Attention to this problem was revived recently by BF,in a paper commentedby 26discussants,who proposed two multivariate “flattening”techniques,ˆBGCV and ˆBCV based,respectively,on generalized cross-validation and ordinary cross-validation.A very extensive simulation convincingly reports notable risk reduction over other commonly used prediction methods.In the idealized situation,assuming B ,Σ,and Ω,are known they considered as well the optimization problemmin A(y −A ˆyOLS ) (y −A ˆy OLS )and obtainedA =[I +r (B ΩB )−1Σ]−1,for a certain scalar r .This leads for p ≤k to the sample versionˆBOPT =Z I +r (Z C −1Z )−1V−1.32.MAIN RESULTSIn this section we investigate conditions under which ˆBOPT and ˆB GCV (2)lead to improved predictions over OLS.Both ˆBOPT and ˆB GCV have an expansion as r ↓0whose two first terms give ˆBFICYREG so from the result of van der Merwe &Zidek (1980)we expect there exists an interval 0<r <r 0such that ˆBOPT and ˆBGCV are uniformly better than OLS.The main tool is that of unbiased risk estimation developed by Stein (1981).Using a result of Golub &van Loan (1989)there exists a nonsingular W :p ×p such thatWVW=diag(λi )W (Z C −1Z )W=diag(γi ),(1)and thusW −1(Z C −1Z )−1VW =diag(x i ),where x i =λi /γi >0,i =1,...,p .Proposition 1.[van der Merwe &Zidek 1980]The SPE risk difference between FICYREG and OLS is unbiasedly estimated byˆR(ˆB FICYREG )−pk =x i ·r [r (m +p +1)−2(k −p −1)].This implies that for0<r <2(k −p −1)/(m +p +1),k >p +1,FICYREG has a SPE risk uniformly less than that of OLS.Proposition 2.The SPE risk difference between OPT and OLS is unbiasedlyestimated byˆR(ˆB OPT )−pk =2r x i (1+rx i )−1· (1+rx i )−1+2r x i (1+rx i )−2+2k (1+rx i )−1−2kp −2r (1+rx i )−2· (1+rx i )−1+2r (1+rx i )−1 2+2r (1+rx i )−2−2r (1+rx i )−3+(m −p −1) x −1i 1−(1+rx i )−1 2.The choice of the flattening matrix A can also be made with the GCV criteriaproposed by Craven &Wahba (1979).BF equation (3.9)showed this approximation reduces toˆBGCV =Z r 2(1−r )−1I + 2−(1−r )−1 I +(Z C −1Z )−1V −1 −1· −r I + I +(Z C −1Z )−1V−1 .(2)4The prescription for r is r =k/n ,n =k +m .For GCV the unbiased risk estimation is contained in the following proposition.Leth (x )=−r +(1+x )−1g (x )=r 2(1−r )−1+ 2−(1−r )−1(1+x )−1.Proposition 3.The SPE risk difference between GCV and OLS is unbiasedly estimated byˆR(ˆB GCV )−pk =−2 2−(1−r )−1x i (1+x i )−2h (x i )g (x i )−2+ (1+x i )−1g (x i )−1·x i (1+x i )−1h (x i )g (x i )−1+2kh (x i )g (x i )−1−2kp +2x i (1+x i )−2g (x i )−1+2 x i (1+x i )−1· (1+x i )−1g (x i )−1+2 2−(1−r )−1(1+x i )−1h (x i )g (x i )−2· (1+x i )−1h (x i )g (x i )−1− (1+x i )−1g (x i )−1·(1+x i )−1h (x i )g (x i )−1+ (1+x i )−2h (x i )2g (x i )−3−(1+x i )−2h (x i )g (x i )−2+2(1+x i )−1·(1+x i )−1g (x i )−1−2 (1+x i )−1· (1+x i )−1h (x i )g (x i )−2−2 (1+x i )−2h (x i )g (x i )−2+2(1+x i )−2g (x i )−1+(m −p −1) x −1i1−h (x i )g (x i )−1 2.Propositions 2and 3have important consequences regarding conditional (givenC )risk,as well as unconditional risk;i.e.,they apply when the predictor is consid-ered fixed or random.Suppose,under a condition on r ,that ˆR(ˆB GCV )−pk ≤0,for all x 1,...,x p >0.Then,averaging over the conditional distribution of the x i ’s given C ,R (ˆB GCV |C )−pk ≤0,for all positive definite matrix C .Averaging further over the distribution of C ,whatever it may be,R (ˆB GCV )−pk ≤0,for all possible parameters B ,Σ,and all distributions on the positive definitematrix C .For example,when the rows,x i ,of X are i.i.d.N k (0,Ω),then C−1has a W k (n,Ω)distribution.The determination of an interval for r such that the unbiased estimate of the risk difference between OPT (or GCV)with OLS is always (for all x i >0)non-positive is more difficult due to the complexity of the expressions.However,on the domain x 1=···=x p ≡x >0such an interval exists.Consideration of this domain leads to simplifications of estimated risks difference in Propositions 2and 3,to which Baranchik’s (1970)technique then applies.The hope is that estimated risks difference will remain nonpositive outside this particular domain.5This can be checked for a given value of r by plotting the function,representing the estimated risks difference,on several rays coming out of the origin.The domain x1=···=x p is only one such ray,which makes a45degrees angle with all the axes. Proposition4.Let k>p+1and definer OPT=2(k−p−1)/(m+p+1),r GCV=1+(m+p+1)2(k−p−1)+(m+p+1)2(k+14)(k−2)/[2(k−4)],where x denotes the smallest integer greater than or equal to x,known as the “ceiling”function.A plot of this minimum value of m versus k in Figure1shows that m need not be very large compared to k.Suchfinite sample properties of GCV,even for p=1,seem novel.6mkFigure1:Minimum value of m as a function of k≥5when p=1so that r=k/n ofˆB GCV belongs to the specified interval in Proposition4.In general,for given values of p,k,and m,a plot ofˆR(˜B)−pk versus x1,...,x p would show that˜B is minimax if the surface is always nonpositive.Figure2shows the surface plot of GCV for p=2and k=4.In the worst scenario m=1and r=r GCV,GCV is seen to be minimax.A similar result can be observed for OPT which is also minimax.The GCV,however,may not be minimax for the prescription r=k/n for small values of m=n−k as shown in Figure3.The surface plot is positive for small values of x1and x2which indicate GCV may not be minimax.Figure2:Unbiased estimate of the risk difference between GCV and OLS,ˆR(ˆB GCV)−pk,for p=2,k=4,m=1,and r=r GCV.7Figure 3:Unbiased estimate of the risk difference between GCV and OLS,ˆR(ˆB GCV )−pk ,for p =2,k =4,m =1,and r =k/n .A surface plot of the function ˆR(˜B )−pk for large values of p is not possible.However,on a given ray emerging from the origin in the positive quadrant and determined by a unit vector u ,one can plot this function at points t u ,on an interval 0≤t ≤a .For GCV,Figure 4shows that the prescription r =r GCV leads to a minimax prediction procedure;this may not be the case,however,for r =k/n when n is too small.12345t -300-250-200-150-100-50U12345t -400-300-200-100100200UFigure 4:Representation of U =ˆR(ˆB GCV )−pk on 10rays chosen at random in the positive quadrant.The values are p =20,k =50,and m =50.On the left r =r GCV =0.348,whereas on the right r =k/n =1/2.To get some insights into the new GCV method we made a comparison of exact8SPE risks when p=1.We use the notations Z=z and(Z C−1Z)−1V=x.The comparison is made between the positive part procedures(let[t]+=max(t,0))ˆb FICYREG+=1−k−2m+2−1z,ˆb GCV+=[h(x)/g(x)]+z,where in h(·)and g(·)we set r=k/n.When p=1the positive parts FI-CYREG+and GCV+are,respectively,better than FICYREG and GCV follows from Baranchik(1970).Note thatˆb OPT does not require a positive part correc-tion.The exact conditional risks are evaluated using the unbiased risk estimates,ˆR≡ˆR(x),and averaging over the conditional distribution of x.It is well known1/x follows a canonical F(k,m;d)conditional distribution[Johnson&Kotz(1970),p. 191]with noncentrality parameter given by c=b C−1b/σ2.Thus,the conditional risk is given byR(·|C)=∞i=0e−c/2(c/2)i2,m(1+y)(n+2i)/2dy.(3)Evaluation of unconditional risk requires averaging R in(3)over the distribution of c.For example,this can be done if C−1∼W k(n,Ω)in which case the distribution is simply c∼(b Ωb/σ2)χ2n.The unconditional risk is then a function of the signal-to-noise ratio d=b Ωb/σ2,R(·)=∞i=0d in/2+i)B−1n2,m(1+y)(n+2i)/2dy.(4)Examination of Figure5reveals for p=1,k=5,and m=40,that GCV+is better than FICYREG+and OPT for small signal-to-noise ratios,whereas OPT is the preferred method for large signal-to-noise ratios.90.51 1.52d12345Risk Figure 5:Comparison of exact risks as a function of signal-to-noise ratio d =b Ωb /σ2.Joined line is for OPT,dashed line corresponds to FICYREG+,whereas dot-dashed line represents GCV+.The values are p =1,k =5,and m =40.The evaluation of the exact risk of FICYREG,OPT,or GCV,when p >1would use the same idea.The exact joint conditional distribution of the eigenvalues x 1,...,x p given in Muirhead (1982,Theorem 10.4.2)is in terms of hypergeometric functions.Although it is known that this joint distribution depends on B and Σonly through the eigenvalues of the signal-to-noise matrix (B C −1B )Σ−1the evaluation of exact conditional risk still remains a difficult computational problem:R (·|C )=y 1>···>y p >0ˆR (y −11,...,y −1p )·etr−12n ;12Θ,F (I +F )−1·Γp1Γp12(m ) πp2/22pp i =1y (k −p −1)/2iso that the flattening matrix A is positive semidefinite.Risk of these estimates are minimum at B =0and they are said to shrink towards 0.Multivariate flattening yields substantial risk reduction especially when the estimate shrinks to a point close to the “true”B .We have (conditionally and,thus,unconditionally)for a general estimate ˜B=˜B (Z ,V )R SPE ˜B (Z ,V );B −B (0),Σ =R SPE B (0)+˜B (Z −B (0),V );B ,Σ .Thus,ˆB (0)FICYREG+=B (0)+ Z −B (0) W (0) diag (1−rx (0)i )+ W (0) −1,ˆB (0)GCV+=B (0)+ Z −B (0) W (0) diag h (x (0)i )/g (x (0)i ) + W (0) −1,ˆB (0)OPT=B (0)+ Z −B (0) I +r Z −B (0) C −1 Z −B (0) −1V −1,whereW (0) −1 Z −B (0) C −1 Z −B (0) −1VW (0) =diag x (0)i,should have good risk performance around B =B (0).ACKNOWLEDGEMENTSThis research was supported in part by a grant from the Natural Sciences and Engineering Research Council of Canada.APPENDIXLemma 1.[Bilodeau &Kariya 1989]Let ˜B =Z +G (Z ,V )be an arbitrary estimate,where G =(g ij ),Z =(z ij ),and V =(v ij ).Assume E g 2ij <∞,E |∂g ij /∂z kl |<∞,E (∂g ij /∂v kl )2<∞,and the conditions in Theorem 2.1of Haff(1976)are satisfied.The unique unbiased estimate of the risk difference R (˜B)−pk is ˆR(˜B )−pk =2tr A +2tr TC −1+(m −p −1)tr GV −1G C −1,with A =(a ij ),T =(t ij )anda ij=p k =1∂g ik /∂z jk t ij =pk =1p l =1[g ik ·∂g jl /∂v kl +∂g ik /∂v kl ·g jl ](1+δkl )/2.Proof of Proposition 2:DefineG (Z ,V )=Z I +r (Z C −1Z )−1V−1−Z 11andU =I +r (Z C −1Z )−1V .Then the derivatives in Lemma 1are obtained as ∂G /∂z jk =Z ∂∂z jkU ·U −1+e j e k U −1−e j e k =−r ZU −1∂∂z jk (Z C −1Z )·(Z C −1Z )−1VU −1+e j e k U −1−e j e k =r ZU −1(Z C −1Z )−1 Z C −1e j e k+e k e j C −1Z ·(Z C −1Z )−1VU −1+e j e k U −1−e j e k .One then extracts ∂g ik /∂z jk from ∂G /∂z jk and easily checks thatA =r tr (Z C −1Z )−1VU −1 ·ZU −1(Z C −1Z )−1Z C −1+r ZU −1(Z C −1Z )−1U −1V (Z C −1Z )−1Z C −1+tr U −1 ·I −p I .Diagonalizing with the matrix W (1)the trace of A reduces totr A =r x i (1+rx i )−1· (1+rx i )−1+r x i (1+rx i )−2+k (1+rx i )−1−kp.Now,we evaluate T with∂G /∂v kl =Z ∂U −1/∂v kl =−ZU −1∂U /∂v kl ·U −1=−r ZU −1(Z C −1Z )−1∂V /∂v kl ·U −1=−r ZU −1(Z C −1Z )−1{e k e l +e l e k }1Then after very lengthy but straightforward algebra one can check thatˆR(ˆB OPT)−pk=px(m+p+1)w1(x)2−2(k−p−1)w1(x)+2(p+1)pw 1(x)[1+w1(x)]x2.(5) Since w1(x)>0,w 1(x)<0,and lim x↓0w1(x)=r,thenˆR(ˆB OPT)−pk<0if 0<r<2(k−p−1)/(m+p+1),k>p+1.For GCV the equation(5)is still validbut withw2(x)=x−11−h(x)g(x)−1.Since w2(x)>0,w 2(x)<0,and lim x↓0w2(x)=r/(r−1)2,we getˆR(ˆB GCV)−pk< 0if0<r/(r−1)2<2(k−p−1)/(m+p+1).Solving this quadratic equation for r gives the result.REFERENCESBaranchik,A.J.(1970).A family of minimax estimators of the mean of a multivariate normal distribution.Ann.Math.Statist.,41,642–645.Bilodeau,M.,&Kariya,T.(1989).Minimax estimators in the normal MANOVA model.J.Multivariate Anal.,28,260–270.Breiman,L.,&Friedman,J.H.(1997).Predicting multivariate responses in multiple linear regression.J.R.Statist.Soc.B,59,3–54.Craven,P.,&Wahba,G.(1979).Smoothing noisy data with spline functions:estimat-ing the correct degree of smoothing by the method of generalized cross-validation.Numer.Math.,31,317–403.Dempster,A.P.,Schatzoff,M.,&Wermuth,M.(1977).A simulation study of alterna-tives to ordinary least-squares.J.Amer.Statist.Assoc.,72,77–106.Golub,H.G.,&van Loan,C.F.(1984).Matrix Computations2nd ed..John Hopkins University Press,Baltimore.Haff,L.R.(1976).An identity for the Wishart distribution with applications.J.Mul-tivariate Anal.,9,531–544.Johnson,N.L.,&Kotz,S.(1970).Distributions in Statistics:Continuous Univariate Distributions-2.Wiley.Muirhead,R.J.(1982).Aspects of Multivariate Analysis.Wiley.Stein,C.(1960).Multiple regression.Contrib.Prob.Statist.(Hotelling volume), 424–443.Stein,C.(1981).Estimation of the mean of a multivariate normal distribution.Ann.Statist.,9,1135–1151.van der Merwe,A.,&Zidek,J.V.(1980).Multivariate regression analysis and canonical variates.Canad.J.Statist.,8,27–39.Received10December1998D´e partement de math´e matiques et de statistique Accepted17June1999Universit´e de Montr´e alC.P.6128,succursale Centre-villeMontr´e al(Qu´e bec)Canada H3C3J7e-mail:bilodeau@dms.umontreal.ca13。