Mixed Models – Part 2

1. Mixed Models – Part 2

2. Learning Outcomes – Part 2 • The participant will learn: • In what circumstances LS residuals provide direct information on how the model for the fixed effects is misspecified • How and why to transform a GLS model into a LS model • The different types of residuals in mixed models and their uses • SAS diagnostics for fixed effects and covariance parameters in mixed models • How PRISM plots provide information about the covariance structure • General strategies for building mixed models • How penalised spline regression can be used as a nonparametric fitting method

3. Department of Statistics TEXAS A&M UNIVERSITY Choosing the models for the mean and the covariance structure

4. Three approaches to modeling the covariance • Random effects covariance structures:- random intercepts and/or random slopes • Covariance pattern models • Unstructured covariance

5. Verbeke, G. & Molenbergs, G. (2000), Linear Mixed Models for Longitudinal Data, Springer, New York

6. 5.4.1 Model-Building Guidelines “Fitting linear mixed models implies that an appropriate mean structure as well as a covariance structure--i.e., both for random effects and error components--needs to be specified. These steps are not independent of each other. First, unless robust inference is used, an appropriate covariance model is essential to obtain valid inferences for the parameters in the mean structure, which is usually of primary interest. … Once an appropriate mean structure Xiβ for E(Yi) has been selected, we use the ordinary least squares (OLS) method to estimate β, and we hereby ignore the fact that not all measurements are independent. It follows from the theory of generalized estimating equations (GEE) that this OLS estimator is consistent for β (Liang and Zeger, 1986). A helpful tool for deciding which time-varying covariates should be included in the model is a plot of the OLS residual profiles versus time. When this plot shows constant variability over time, we assume stationarity and we do not include random effects other than intercepts. In cases where the variability varies over time and where there is still some remaining systematic structure in the residual profiles (i.e., where the between-subject variability is large in comparison to the overall variation), the following guidelines can be used to select one or more random effects additional to the random intercepts.” Dmitrienko, A., Molenberghs, G. Chuang-Stein, C. & Offen, W. (2005), Analysis of Clinical Trials Using SAS – A Practical Guide, SAS, Cary NC

7. Department of Statistics TEXAS A&M UNIVERSITY A cautionary example

8. Interest Rate Example Tryfos, P. (1998), Methods for Business Analysis and Forecasting: Text & Cases, Wiley, New York

9. Interest Rate Example Conclusions:

10. Interest Rate Example Next Steps: Add a quadratic term in Loans Closed to the model Or ?

11. Interest Rate Example

12. Interest Rate Example: Model with AR(1) errors Generalized least squares fit by maximum likelihood Model: InterestRate ~ LoansClosed + I(LoansClosed^2) + VacancyIndex Correlation Structure: AR(1) Formula: ~Month Parameter estimate(s): Phi 0.9568896 Coefficients: Value Std.Error t-value p-value (Intercept) 7.163999 0.4505311 15.901231 0.0000 LoansClosed -0.004723 0.0044814 -1.053932 0.3086 I(LoansClosed^2) 0.000009 0.0000285 0.298965 0.7691 VacancyIndex -0.079368 0.1349623 -0.588074 0.5652 What happened to the quadratic term in LoansClosed that is evident in the LS residual plots?

13. Interest Rate Example: Model with AR(1) errors Generalized least squares fit by maximum likelihood Model: InterestRate ~ LoansClosed + VacancyIndex Correlation Structure: AR(1) Formula: ~Month Parameter estimate(s): Phi 0.9572093 Coefficients: Value Std.Error t-value p-value (Intercept) 7.122990 0.4182065 17.032232 0.0000 LoansClosed -0.003432 0.0011940 -2.874452 0.0110 VacancyIndex -0.076340 0.1307842 -0.583710 0.5676

14. Interest Rate Example after transforming a model with AR(1) errors into a model with iid errors Call: lm(formula = ystar ~ xstar - 1) Coefficients: Estimate Std. Error t value Pr(>|t|) xstar(Intercept) 7.122990 0.418207 17.032 1.12e-11 xstarLoansClosed -0.003432 0.001194 -2.874 0.011 xstarVacancyIndex -0.076340 0.130784 -0.584 0.568 Generalized least squares fit by maximum likelihood Coefficients: Value Std.Error t-value p-value (Intercept) 7.122990 0.4182065 17.032232 0.0000 LoansClosed -0.003432 0.0011940 -2.874452 0.0110 VacancyIndex -0.076340 0.1307842 -0.583710 0.5676

15. Interest Rate Example after transforming a model with AR(1) errors into a model with iid errors

16. Interpreting residual plots Li, K.C. and Duan, N. (1989), Regression analysis under link violation, Annals of Statistics, 17, 1009-1052.

17. Interpreting residual plots Cook, R.D. and Weisberg, S. (1999), Graphs in statistical analysis: is the medium the message? The American Statistician, 53(1), 29-37.

18. Interpreting residual plots

19. Using generalized least squares when the errors are AR(1)

20. Transforming a model with AR(1) errors into a model with iid errors

21. A general approach to transforming GLS into LS

22. Coming in 2008

23. Lessons learnt from the cautionary example

24. Department of Statistics TEXAS A&M UNIVERSITY Random effects covariance structures

25. Random effects covariance structures: Orthodontic Growth Data Investigators at the University of North Carolina Dental School followed the growth of 27 children (16 males, 11 females) from age 8 top age 14. Every two years they measured the distance between the pituitary and pterygomaxillary fissure, two points that are easily identified on x-ray exposures of the side of the head. These data are reported in Pothoff and Roy (1964) and analyzed in detail in Pinheiro and Bates (2000) and elsewhere. Pinheiro, J.C. & Bates, D.M. (2000). Mixed Effects Models in S and S-PLUS, Springer, New York. Pothoff, R.F. & Roy, S.N. (1964). A generalized multivariate analysis of variance model especially useful for growth curve problems, Biometrika 51, 313-326.

26. Random effects covariance structures: Orthodontic Growth Data distance age Subject Sex 1 26.0 8 M01 Male 2 25.0 10 M01 Male 3 29.0 12 M01 Male 4 31.0 14 M01 Male . 61 22.0 8 M16 Male 62 21.5 10 M16 Male 63 23.5 12 M16 Male 64 25.0 14 M16 Male 65 21.0 8 F01 Female 66 20.0 10 F01 Female 67 21.5 12 F01 Female 68 23.0 14 F01 Female . 105 24.5 8 F11 Female 106 25.0 10 F11 Female 107 28.0 12 F11 Female 108 28.0 14 F11 Female

27. Random effects covariance structures: Orthodontic Growth Data

29. Random effects covariance structures: Orthodontic Growth Data – Plots of Y= Distance Correlations T8 T10 T12 T14 T8 1.000 0.626 0.711 0.600 T10 0.626 1.000 0.635 0.759 T12 0.711 0.635 1.000 0.795 T14 0.600 0.759 0.795 1.000 The correlations are relatively constant across time. This suggests that we fit ______ as the covariance structure.

30. Random intercepts model Compound Symmetry = Random intercepts

31. Example Covariance Matrices Compound Symmetry = Random intercepts

32. Random intercepts covariance structure: Orthodontic Growth Data – R output Linear mixed-effects model fit by REML Data: Orthodont AIC BIC logLik 445.7572 461.6236 -216.8786 Random effects: Formula: ~1 | Subject (Intercept) Residual StdDev: 1.816214 1.386382 Fixed effects: distance ~ age * Sex Value Std.Error DF t-value p-value (Intercept) 16.340625 0.9813122 79 16.651810 0.0000 age 0.784375 0.0775011 79 10.120823 0.0000 SexFemale 1.032102 1.5374208 25 0.671321 0.5082 age:SexFemale -0.304830 0.1214209 79 -2.510520 0.0141

33. Random intercepts covariance structures: Orthodontic Growth Data – SAS output

34. Random intercepts covariance structure: Orthodontic Growth Data

35. Random intercepts covariance structure: Orthodontic Growth Data: Fitted lines for each subject Examples of “shrinkage” in random intercepts

36. Likelihood ratio tests For models with the same fixed effects, the REML log-likelihoods can be used to compare nested models for the covariance. A formal statistical test is obtained by taking twice the difference in the respective maximized log-likelihoods and comparing it to a chi-squared distribution with degrees of freedom equal to the difference between the number of covariance parameters in the full and reduced models. This test is called the likelihood ratio test. Fitzmaurice, G.M., Laird, N.N. & Ware, J.H. (2004). Applied Longitudinal Analysis. Wiley, New York (Chapter 4)

37. Random intercepts versus unstructured covariance: Orthodontic Growth Data Conclusion:

38. Testing on the boundary of the parameter space In general, when testing a null hypothesis that is on the boundary of the parameter space (e.g., the variance of a random effect equals zero), the usual null distribution for the likelihood ratio test is no longer a chi-squared distribution with degrees of freedom equal to the difference between the number of parameters in the full and reduced models; instead, the null distribution is a mixture of chi-squared distributions. For example, when comparing two nested models, one with q correlated random effects, the other with q+1 correlated random effects, the null distribution of the likelihood ratio test is a 50:50 mixture of chi-squared distributions with q and q+1 degrees of freedom. Fitzmaurice, G.M., Laird, N.N. & Ware, J.H. (2004). Applied Longitudinal Analysis. Wiley, New York (Chapter 8)

39. Random intercepts covariance structure: Orthodontic Growth Data – R output Linear mixed-effects model fit by REML Fixed effects: distance ~ age * Sex Value Std.Error DF t-value p-value (Intercept) 16.340625 0.9813122 79 16.651810 0.0000 age 0.784375 0.0775011 79 10.120823 0.0000 SexFemale 1.032102 1.5374208 25 0.671321 0.5082 age:SexFemale -0.304830 0.1214209 79 -2.510520 0.0141 lm(formula = distance ~ age * Sex, data = Orthodont) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 16.3406 1.4162 11.538 < 2e-16 age 0.7844 0.1262 6.217 1.07e-08 SexFemale 1.0321 2.2188 0.465 0.643 age:SexFemale -0.3048 0.1977 -1.542 0.126 Model df AIC BIC logLik Test L.Ratio p-value RI 1 6 445.7572 461.6236 -216.8786 lm0 2 5 493.5591 506.7811 -241.7796 1 vs 2 49.80187 <.0001 Critical value (alpha=0.001) is 9.55 Conclusion:

40. Comparing models using AIC & BIC “Often it is of interest to compare non-nested models for the covariance. To compare non-nested models, an … approach is the Akaike Information Criterion (AIC). According to the AIC, given a set of competing models for the covariance, one should select the model that minimizes AIC = -2l + 2c where l is the maximized REML log- likelihood and c is the number of covariance parameters. Fitzmaurice, G.M., Laird, N.M. & Ware, J.H. (2004). Applied Longitudinal Analysis. Wiley, New York (page 176)

41. Comparing models using AIC & BIC

42. Comments on AIC and BIC “BIC extracts a very large penalty for the estimation of each additional covariance parameter. In general we do not recommend the use of BIC for covariance model selection as it entails a high risk of selecting a model that is too simple or parsimonious for the data at hand.” Fitzmaurice, G.M., Laird, N.M. & Ware, J.H. (2004). Applied Longitudinal Analysis. Wiley, New York (page 177)

43. Residuals in mixed models

45. Scaled or Cholesky residuals

47. Random intercepts covariance structure: Orthodontic Growth Data – SAS output Conclusions:

48. Random intercepts covariance structure: Orthodontic Growth Data – SAS output Conclusions:

49. Random intercepts covariance structure: Orthodontic Growth Data – Marginal & conditional residuals Conclusions:

50. Random intercepts covariance structure: Orthodontic Growth Data – Mixed model diagnostics in SAS Conclusions: