A Tutorial On Restricted Maximum Likelihood Estimation In .

A Tutorial on Restricted Maximum Likelihood Estimationin Linear Regression and Linear Mixed-Effects ModelXiuming Zhangzhangxiuming@u.nus.eduA*STAR-NUS Clinical Imaging Research CenterOctober 12, 2015SummaryThis tutorial derives in detail an estimation procedure—restricted maximum likelihood (ReML) [Patterson and Thompson, 1971] [Harville, 1974]—that is able to produceunbiased estimates for variance components of an linear model. We first introduce theconcept of bias in variance components by maximum likelihood (ML) estimation insimple linear regression and then discuss a post hoc correction. Next, we apply ReMLto the same model and compare the ReML estimate with the ML estimate followed bypost hoc correction. Finally, we explain the linear mixed-effects (LME) model for longitudinal analysis [Bernal-Rusiel et al., 2013] and demonstrate how to obtain unbiasedestimators of the parameters with ReML.1 Linear RegressionFamiliarity with basic linear regression facilitates the understanding of more complexlinear models. We Therefore, start with this and introduce the concept of bias inestimating variance components.1.1 The ModelThe simplest linear regression is of the form y β0 β1 x , where y is named response(or dependent variable or prediction), x is named regressor (or explanatory variable,independent variable), β’s are regression coefficients, and is called residual (or1

error), which is assumed to distribute as a zero-mean Gaussian with an unknownvariance, i.e., N (0, σ 2 ).When we have more than one regressor (a.k.a. multiple linear regression1 ), themodel comes in its matrix formy Xβ ,(1)where y is the response vector, X is the design matrix with each its row specifyingunder what design or conditions the corresponding response is observed (hence thename), β is the vector of regression coefficients, and is the residual vector distributingas a zero-mean multivariable Gaussian with a diagonal covariance matrix N (0, σ 2 IN ),where IN is the N N identity matrix. Thereforey N (Xβ, σ 2 IN ),(2)meaning that linear combination Xβ explains (or predicts) response y with uncertaintycharacterized by a variance of σ 2 .As seen, this assumption on the covariance matrix mandates that (i) residualsamong responses are independent, and (ii) residuals of all the responses have the samevariance σ 2 . This assumption makes parameter estimations straightforward, but meanwhile imposes some limitations on the model. For example, the model cannot handleproperly intercorrelated responses, such as the longitudinal measurements of one individual. This motivates the usage of linear mixed-effects (LME) model in analyzinglongitudinal data [Bernal-Rusiel et al., 2013].1.2 Parameter EstimationUnder the model assumptions, we aim to estimate the unknown parameters (β and σ 2 )from the data available (X and y). Maximum likelihood (ML) estimation is the mostcommon estimator. We maximize the log-likelihood w.r.t. β and σ 2L(β, σ 2 y, X) N1Nlog 2π log σ 2 2 (y Xβ)T (y Xβ)222σand obtainβ̂ (X T X) 1 X T y1σ̂ 2 (y X β̂)T (y X β̂),N(3)(4)where N is the number of responses. β̂ is simply the ordinary least squares (OLS)estimator of β, and we compute σ̂ 2 with the value of β̂. As we will prove in the1A single-response special case of general linear model, which itself is a special case of generalizedlinear model with identity link and normally distributed responses.2

following section, estimator σ̂ 2 is biased downwards as compared with real value σ 2 ,because we neglect the loss of degree of freedom (DoF) for estimating β.1.3 Estimation Bias in Variance ComponentThe bias of an estimator refers to the difference between this estimator’s expectation(here E{σ̂ 2 }) and the true value (here σ 2 ). To facilitate computing E{σ̂ 2 }, we definematrix A X(X T X) 1 X T , so Ay X β̂ ŷ. In fact, A is an orthogonal projectionthat satisfies A2 A (idempotent) and A A (Hermitian or self-adjoint or simplysymmetric for real matrices). We can then verify that (IN A)T (IN A) IN A.E{σ̂ 2 } 1E{(y X β̂)T (y X β̂)}N1E{(y Ay)T (y Ay)}N1E{y T (I A)T (I A)y}N1E{y T (I A)y}N 1E{y T y} E{y T Ay}NTheorem 1.1If y N (0, IN ), and A is an orthogonal projection, then y T Ay χ2 (k) withk rank(A).Proof. If A is idempotent, its eigenvalues satisfy λ2 λ. If A is also Hermitian,its eigenvalues are reala . Hence, in eigendecomposition A QΛQ 1 QΛQT b ,Λ is a diagonal matrix containing either 0 or 1.(i) When A is full-rank, Λ IN , A has to be IN , and y T Ay y T y. Then theresult follows immediately from the definition of chi-square distribution.(ii) When A is of rank k N ,TTTy Ay y QΛQ y WT"#Ik 000W WkT Wk ,where W QT y, and Wk denotes the vector containing the first k elements ofW . Since Q is orthogonal, W N (0, IN ), so Wk N (0, Ik ). The result againfollows directly from the definition of chi-square distribution.abSee http://proofwiki.org/wiki/Hermitian Matrix has Real Eigenvalues.Eigenvectors of a real symmetric matrix are orthogonal, i.e., Q 1 QT .3

Recall Equation (2). By Theorem 1.1, y Xβ T y Xβ χ2 (N ).σσHence,E{y T y} N σ 2 (Xβ)T (Xβ).Similarly, y Xβσ T Ay Xβσ χ2 (k),where k rank(A). Hence,E{y T Ay} kσ 2 (Xβ)T (Xβ).Substituting E{y T y} and E{y T Ay} givesE{σ̂ 2 } N k 2σ σ2,N(5)where k rank(A) is just the number of columns (regressors) in X. Therefore, ourestimation of the variance component is biased downwards. This bias is especiallysevere when we have many regressors (a large k), in which case we need to correctthis bias by simply multiplying a factor of N/(N k). Hence, the corrected, unbiasedestimator becomes1(y X β̂)T (y X β̂)N k1 (y X(X T X) 1 X T y)T (y X(X T X) 1 X T y),N k2 σ̂unbiased(6)which is classically used in linear regression [Verbeke and Molenberghs, 2009].2 Restricted Maximum LikelihoodIn simple problems where solutions to variance components are closed-form (like linear regression above), we can remove the bias post hoc by multiplying a correctionfactor. However, for complex problems where closed-form solutions do not exist, weneed to resort to a more general method to obtain a bias-free estimation for variancecomponents. Restricted maximum likelihood (ReML) [Patterson and Thompson,1971] [Harville, 1974] is one such method.2.1 The TheoryGenerally, estimation bias in variance components originates from the DoF loss inestimating mean components. If we estimated variance components with true mean4

component values, the estimation would be unbiased. The intuition behind ReML is tomaximize a modified likelihood that is free of mean components instead of the originallikelihood as in ML.Consider a general linear regression modely Xβ ,where y is still an N -vector of responses, X is still an N k design matrix, but residual is no longer assumed to distribute as N (0, σ 2 IN ), but rather N (0, H(θ)), where H(θ)is a general covariance matrix parametrized by θ. For simplicity, H(θ) is often writtenas just H. Previously, we have been referring to θ as “variance components.”If vector a is orthogonal to all columns of X, i.e., aT X 0, then aT y is known as anerror contrast.We can find at mostN k such vectors that are linearly independent2 .hiDefine A a1 a2 . . . aN k . It follows that AT X 0 and E{AT y} 0. S IN X(X T X) 1 X T is a candidate for A, as SX 0. Furthermore, it can be shownAAT S and AT A IN .The error contrast vectorw AT y AT (Xβ ) AT N (0, AT HA)is free of β. [Patterson and Thompson, 1971] has proven that in the absence of information on β, no information about θ is lost when inference is based on w rather than ony. We can now directly estimate θ by maximizing a “restricted” log-likelihood functionLw (θ AT y). This bypasses estimating β first and can Therefore, produce unbiasedestimates for θ.Once H(θ) is known, the generalized least squares (GLS) solution to β minimizing squared Mahalanobis length of the residual (Y Xβ)T H 1 (Y Xβ) is justβ̂ (X T H 1 X) 1 X T H 1 y.(7)We now derive a convenient expression for Lw (θ AT y) [Harville, 1974].Lw (θ AT y) log fw (AT y θ)TZ log fw (A y θ)fβ̂ (β̂ β, θ) dβ̂Z log fw (AT y θ) fβ̂ (GT y β, θ) dβZ log fw (AT y θ)fβ̂ (GT y β, θ) dβZ log fw,β̂ (AT y, GT y β, θ) dβ2(β̂, β exchangeable here)Imagine a plane spanned by k 2 linearly independent vectors in three-dimensional space (N 3). Wecan find at most N k 1 vector (passing the origin) orthogonal to the plane.5

Z logfy hiT y β, θ dβA GZ1hi logfy (y β, θ) dβ. det A G Interlude 2.1hiWe express det A G in terms of X. hiT hi 12 det A G det A GA Ghi det"#! 12AT G A T GGT A G T G 1 1 det AT A 2 det GT G GT A(AT A) 1 AT G 2 11 (det I) 2 det GT G GT AI 1 AT G 2 1 det GT G GT SG 2 1 det X T X 2We continue deriving1hiLw (θ A y) log det A G Z 12ZTT log det X Xfy (y β, θ) dβfy (y β, θ) dβ 1T 1p log det X Xexp (y Xβ) H (y Xβ) dβ2(2π)N det H Z 2111 T NT 1 log det X X (2π) 2 (det H) 2 exp (y Xβ) H (y Xβ) dβ2T 12Z1Interlude 2.2We can decompose (y Xβ)T H 1 (y Xβ) into(y X β̂)T H 1 (y X β̂) (β β̂)T (X T H 1 X)(β β̂)with Equation (7).We resumeTTLw (θ A y) log det X X log det X T X 21(2π) 12(2π) 2 (det H) 2 N2 12(det H)N16 1T 1exp (y Xβ) H (y Xβ) dβ2 1exp (y X β̂)T H 1 (y X β̂)2Z

1exp (β β̂)T (X T H 1 X)(β β̂) dβ2 2111 NT T 1 log det X X (2π) 2 (det H) 2 exp (y X β̂) H (y X β̂)2Zk1(2π) 2 (det X T H 1 X) 2(a Gaussian integral) 111 log(2π)det X T X 2 (det H) 2 (det X T H 1 X) 2 1T 1exp (y X β̂) H (y X β̂)21111 (N k) log(2π) log det X T X log det H log det X T H 1 X22221T 1(8) (y X β̂) H (y X β̂),2 12 (N k)whereβ̂ (X T H 1 X) 1 X T H 1 y.(9)With this convenient expression, we can maximize the restricted log-likelihoodLw (θ AT y) w.r.t. variance components θ to obtain an unbiased estimate for the covariance matrix H(θ̂) and the corresponding regression coefficient estimates β̂. NewtonRaphson method is usually employed. For more computational details, see [Lindstromand Bates, 1988].2.2 Applied to Simplest Linear RegressionWe have seen the estimation bias in θ by ML from Equation (5). In the simplestform of linear regression where we assume H σ 2 IN , estimation σ̂ 2 is closed-form(Equation (4)), allowing us to correct the bias simply with a multiplicative factor. Inthis section, we verify that, in the simplest form of linear regression, the ReML methodproduces exactly the same solutions as ML method followed by the post hoc correction.Setdd111Lw (σ 2 AT y) log det H log det X T H 1 X (y X β̂)T H 1 (y X β̂)22dσdσ222d11 (N k) log σ 2 2 (y X β̂)T (y X β̂)dσ 2 22σ 0We obtain exactly the same result as Equation (6), produced by post hoc correction.It is worth noticing that in this simplest linear regression case, the mean estimateβ̂ is independent of the variance component θ (Equation (3)). This implies althoughthe ML and ReML estimates of σ̂ 2 are different, the estimates of β̂ are the same. Thisis no longer true for more complex regression models, such as the linear mixed-effectsmodel, as to be seen in the next section. Thus, for those complex models, we have a7