- p. 1/??Statistics 203: Introduction to Regressionand Analysis of VarianceMixed EffectsJonathan Taylor- p. 2/??Today’s class■Mixed effects two-way ANOVA – what we didn’t cover lastclass.■Mixed effects models – random intercepts.■General form of a mixed effect model.■Details: MLE and REML.- p. 3/??Mixed effects model■In some studies, some factors can be thought of as fixed,others random.■For instance, we might have a study of the effect of astandard part of the brewing process on sodium levels in thebeer example.■Then, we might think of a model in which we have a fixedeffect for “brewing technique” and a random effect for beer.- p. 4/??Two-way mixed effects model■Yijk∼ µ··+ αi+ βj+ (αβ)ij+ εij, 1 ≤ i ≤ r, 1 ≤ j ≤ m, 1 ≤k ≤ n■εijk∼ N(0, σ2), 1 ≤ i ≤ r, 1 ≤ j ≤ m, 1 ≤ k ≤ n■αi∼ N(0, σ2α), 1 ≤ i ≤ r.■βj, 1 ≤ j ≤ m are constants.■(αβ)ij∼ N(0, (m − 1)σ2αβ/m), 1 ≤ j ≤ m, 1 ≤ i ≤ r.■Constraints:◆Pmj=1βj= 0◆Pri=1(αβ)ij= 0, 1 ≤ i ≤ r.◆Cov ((αβ)ij, (αβ)i0j0) = −σ2αβ/m■Cov(Yijk, Yi0j0k0) =δjj0σ2β+ δii0m−1mσ2αβ− (1 − δii0)1mσ2αβ+ δii0δkk0σ2- p. 5/??ANOVA tables: Two-way (mixed)SS df E(MS)SSA r − 1 σ2+ nmσ2αSSB m − 1 σ2+ nrPmj=1β2im−1+ nσ2αβSSAB (m − 1)(r − 1) σ2+ nσ2αβSSE =Pri=1Pmj=1Pnk=1(Yijk− Yij·)2(n − 1)rm σ2■To test H0: σ2α= 0 use SSA and SSE.■To test H0: β1= · · · = βm= 0 use SSB and SSAB.■To test H0: σ2αβuse SSAB and SSE.- p. 6/??Mixed linear models■Not every model is an ANOVA!■Suppose we study the effect of a blood pressure meant tolower blood pressure over time and we study r patients.■For each patient we record BP at regular intervals over aweek (every day, say).■Drug will have varying efficacy in the population.■ModelYij= β0+ γi+ β1Xij+ εij◆εij∼ N(0, σ2) i.i.d.◆γi∼ N(0, σ2γ) i.i.d. and independent of εij’s◆Xij= j, 1 ≤ j ≤ 7 in this example...- p. 7/??Random intercept model – balanced■Model is called a “random intercept” model.■It is “balanced” because every subject had the exact sameX’s.■E(Yij) = β0+ β1Xj.■Cov(Y ) = σ2Inr×nr+ σ2γJn0 . . . 00 Jn. . . 0............0 . . . . . . Jn= Ir⊗ V0where Jnis an n × n matrix with all 1’s;V0= σ2In+ σ2γJn.- p. 8/??Kronecker product■A useful bit of notation for multivariate statistics.■Given two matrices Am×nand Bp×q, the Kronecker productA ⊗ B is an mp × nq matrix with entries(A ⊗ B)(i−1)p+l,(j−1)q+k= Aij· Blk.■Properties:(A ⊗ B)t= At⊗ Bt(A ⊗ B)−1= A−1⊗ B−1(A ⊗ B)(X ⊗ Y ) = (AX) ⊗ (BY )rank(A ⊗ B) = rank(A)rank(B)Tr(A ⊗ B) = Tr(A)Tr(B)det(Aa×a⊗ Bb×b) = det(A)bdet Ba- p. 9/??Generalized least squaresIn a fixed effects model:Y = Xβ + ε■We know what OLS estimators of β are:bβOLS= (XtX)−1XtYif ε ∼ N(0, σ2) then they are MLE.■If ε ∼ N(0, D), D diagonalbβW LS= (XtD−1X)−1XtD−1Yare MLE estimators.■What if ε ∼ N(0, V ), V not diagonal? ThenbβGLS= (XtV−1X)−1XtV−1Yare MLE estimators. Proof is similar to OLS and WLS cases.- p. 10/??Fitting the random intercept model■Likelihood has parameters β0, β1, σ2, σ2γ.■For σ2, σ2γfixedbβ=(XtV−1X)(XtV−1Y ), V = Ir⊗ (σ2In+ σ2γJn)■ML estimates of β0, β1turn out to be identical to fixed effectsanalysisbβ1=Pnj=1XjY·j− nX · Y··Pnj=1X2j− nX2bβ0= Y··−bβ1XwhereY·j=1rrXi=1Yij.- p. 11/??Variance-covariance ofbβ■What about variance? In the generalized least squaressettingVar(bβ) = (XtV−1X)−1.■In the balanced random intercept modelVar(bβ0) =σ2γ+ σ2/nr+σ2X2rPnj=1(Xj− X)2Var(bβ1) =σ2rPnj=1(Xj− X)2Cov(bβ0,bβ1) = X · Var(bβ1).■Similar to fixed effects: only difference is σ2γin Var(β0).- p. 12/??Estimating σ2, σ2γ■If you go through the calculus, ignoring constraints σ2, σ2γ≥ 0the maximum likelihood estimator of σ2(˙· denotes a“pre”-MLE) is˙σ2=Pi,j(Yij−bYij)2r(n − 1)=SSEr(n − 1)and˙σ2γ=1nSSAm− bσ2whereSSE =Xi,j(Yij−bYi,j)2, SSA =Xi,j(Yi·− Y··)2■If ˙σ2γ< 0, then bσ2γ= 0 andbσ2=SSA + SSErn.- p. 13/??Unbalanced data■If not all the X’s are the same for each subject, or someobservations are missing, things are more complicated.■Likelihood has the same parameters only the variancecovariance matrix is no longer a Kronecker product and MLEestimates are not closed form anymore.■This is what computers are for!■Two different techniques: MLE and REML.- p. 14/??Residual maximum likelihood■Suppose that we have an n − p × n matrix K of rank n − psuch that K0X = 0.■ThenK0Y = K0(Xβ + ε) = Kε ∼ N(0, K0V K).■The distribution of K0Y doesn’t have any β’s in it.■We can use K0Y to estimate parameters in V (σ2, σ2γ)). Thistechnique is called REML – residual maximum likelihood.Also works if K is n × n, only then the inverses have to bethought of as generalized inverses.■Once V has been estimated bybVREM Lwe estimatebβREM L= (XtbV−1REM LX)−1XtbV−1REM LY.■We can approximate variance-covariance matrix asbβREM L∼ (XtbV−1REM LX)−1.- p. 15/??General form of a mixed effects model■The random intercept model can be generalized quite a bit:random slopes, more than one predictor, etc.■The random/mixed ANOVA models and random interceptmodel all have the formYn×1= Xn×pβp×1+ Zn×qγq×1+ εn×1where◆ε ∼ N(0, σ2I);◆γ ∼ N(0, D) for some covariance D: simplest model, D isdiagonal■In this modelY ∼ N(Xβ, ZDZ0+ σ2I).■Parameters to estimate: β and any parameters in D. If D isdiagonal, then there are q “variance” components