Overview of this module Overview of this module Course 02429 Analysis of correlated data: Mixed Linear Models 1 Different view on the random effects approach Example: Activity of rats 2 Gaussian model of spatial correlation Example: Activity of rats 3 Other spatial correlation structures 4 Diagram of analysis 5 The semi-variogram 6 Analysing the time structure by polynomial regression Module 12: Repeated measures II, advanced methods Per Bruun Brockhoff DTU Compute Building 324 - room 220 Technical University of Denmark 2800 Lyngby – Denmark e-mail: [email protected] Per Bruun Brockhoff ([email protected]) Mixed Linear Models, Module 12 Fall 2014 1 / 24 Per Bruun Brockhoff ([email protected]) Aim of this module Mixed Linear Models, Module 12 Fall 2014 2 / 24 Aim of this module Aim of this module Remember the rats data set? Summary of experiment: 3 treatments: 1, 2, 3 (concentration) 10 cages per treatment 10 contiguous months The response is activity (log(count) of intersections of light beam during 57 hours) Extend the toolbox for dealing with repeated measures Introduce models where the covariance structure is directly specified 10.5 Improve our ability to handle long series 10.0 log(count) See how to specify these models in R 9.5 9.0 8.5 1 2 3 4 5 6 7 8 9 10 Month Per Bruun Brockhoff ([email protected]) Mixed Linear Models, Module 12 Fall 2014 3 / 24 Per Bruun Brockhoff ([email protected]) Mixed Linear Models, Module 12 Fall 2014 4 / 24 Different view on the random effects approach Different view on the random effects approach Different view on the random effects approach Activity of rats analyzed via compound symmetry model Any mixed model can be expressed as: The model is the same as the random effects model, but specified directly y ∼ N (Xβ, ZGZ0 + R), The total covariance of all observations are described by V = ZGZ0 + R The ZGZ0 part is specified through the random effects of the model The R part has so far been σ 2 I, but in this module we will put some structure into R For instance the structure known from the random effects model   0 σ2 cov(yi1 , yi2 ) =  individual 2 σindividual + σ2 , if individuali1 6= individuali2 and i1 6= i2 , if individuali1 = individuali2 and i1 6= i2 , if i1 = i2 This structure is known as compound symmetry Per Bruun Brockhoff ([email protected]) Mixed Linear Models, Module 12 Fall 2014 6 / 24 lnc ∼ N (µ, V), where µi = µ + α(treatmi ) + β(monthi ) + γ(treatmi , monthi ), , if cagei1 6= cagei2 and i1 6= i2  0 σd2 , if cagei1 = cagei2 and i1 6= i2 Vi1 ,i2 =  2 2 σd + σ , if i1 = i2 lme(lnc ~ month + treatm + month:treatm, random = ~1 | cage, data = rats) OR directly into the R-matrix: gls(lnc~month+treatm+month:treatm, correlation=corCompSymm(form=~1|cage),data=rats) Per Bruun Brockhoff ([email protected]) 7 / 24 Example: Activity of rats , if individual i1 6= individual i2 and i1 6= i2 , if individual i1 = individual i2 and i1 6= i2 ν2 + τ2 , if i1 = i2 The entire model is: lnc µi ∼ = Vi1 ,i2 = N (µ, V), where µ + α(treatmi ) + β(monthi ) + γ(treatmi , monthi ), and    0 , if cagei1 6= cagei2 and i1 6= i2  −(monthi1 −monthi2 )2 2 2 , if cagei1 = cagei2 and i1 6= i2 ν + τ exp ρ2    2 ν + τ 2 + σ2 , if i1 = i2 This model is implemented by: lme(lnc~month+treatm+month:treatm, random=~1|cage, correlation=corGaus(form=~as.numeric(month)|cage,nugget=T), data=rats) ν2 ν2 + 0.5τ2 0 Covariance = Fall 2014 Rats data via spatial Gaussian correlation model Covariance structures depending on “how far” observations are apart are known as spatial The following covariance structure has been proposed for repeated measurements  Vi1 ,i2 Mixed Linear Models, Module 12 Gaussian model of spatial correlation Gaussian model of spatial correlation and Implemented in R traditionally by a random effect: Gaussian model of spatial correlation 0    −(ti1 −ti2 )2 ν 2 + τ 2 exp 2 ρ    2 ν + τ 2 + σ2 Example: Activity of rats 0 Per Bruun Brockhoff ([email protected]) Distance t − t 0.83ρ i1 12 i2 Mixed Linear Models, Module Fall 2014 9 / 24 Per Bruun Brockhoff ([email protected]) Mixed Linear Models, Module 12 Fall 2014 10 / 24 Gaussian model of spatial correlation Example: Activity of rats Gaussian model of spatial correlation Example: Activity of rats Reduction from spatial Gaussian to random effects? The relevant part of the R output is: Random effects: Formula: ~1 | cage (Intercept) Residual √ ˆ 2 + τˆ2 ) StdDev: 0.1404056 (= νˆ) 0.2171559 (= σ Correlation Structure: Gaussian spatial correlation Formula: ~as.numeric(month) | cage Parameter estimate(s): range nugget 2.3863954 (= ρˆ2 ) 0.2186744 (= σ ˆ 2 /(ˆ σ 2 + τˆ2 )) Number of Observations: 300 Number of Groups: 30 The spatial Gaussian model (A) is an extension of the random effects model (B) Use the restricted/residual likelihood ratio test (A) 2 GA→B = 2`(B) re − 2`re , where GA→B ∼ χ2 For the rats data we get: Model (A) Spatial Gaussian (B) Random effects Notice the R parametrization of the variance parameters 2`re -105.3 8.6 G–value GA→B = 113.9 df 2 P–value PA→B < 0.0001 The spatial Gaussian structure is clearly a better description of the covariance structure in the rats data set Per Bruun Brockhoff ([email protected]) Mixed Linear Models, Module 12 Fall 2014 11 / 24 Per Bruun Brockhoff ([email protected]) Other spatial correlation structures Mixed Linear Models, Module 12 Fall 2014 12 / 24 Diagram of analysis Other spatial correlation structures Diagram of analysis R has several build–in correlation structures. A few examples are: Write in R Name corGaus Gaussian corExp corAR1 corSymm exponential autoregressive(1) unstructured Correlation term Identify "individuals" −(ti1 −ti2 )2 } ρ2 −|ti1 −ti2 | 2 τ exp{ } ρ τ 2 ρ|i1 −i2 | τ 2 exp{ Select covariance structure from knowledge about the experiment guided by information criteria τi21 ,i2 Unfortunately is can be very difficult to choose — especially for “short” individual series General advice: Keep it simple: Numerical problems often occur with (too) complicated structures Graphical methods: Especially for “long” series the (semi)–variogram is useful Information criteria: AIC or BIC = “2` + penalty(#par)” can be used as guideline Try to cross–validate your main conclusion(s) by one of the “simple” methods Per Bruun Brockhoff ([email protected]) Mixed Linear Models, Module 12 Fall 2014 14 / 24 Covariance parameters are tested by likelihood ratio test The green arrow is often omitted by the argument that a non–significant simplification of the mean structure should not change the covariance structure much Select fixed effects Select covariance structure Change model Change model Test covariance parameters Test fixed effects Interpret results Per Bruun Brockhoff ([email protected]) Mixed Linear Models, Module 12 Fall 2014 16 / 24 The semi-variogram The semi-variogram The semi-variogram The semi-variogram print(plot(Variogram(model2, form =˜as.numeric(month)|cage, data = rats))) ● Plotting of the empirical covariation (of residuals) versus "time" Plot of σ 2 + γ(u), where ● ● γ(u) = τ 2 1 − λ(u) λ(u) = exp{−u2 /ρ2 } ● ● ● 0.5 ● 2 3 ν2 ● 1 τ2 0.0 σ2 0 Covariance 4 Semivariogram 1.0 2 4 6 8 Distance 0 1 Per Bruun Brockhoff ([email protected]) 2 3 4 Mixed Linear Models, Module 12 Fall 2014 18 / 24 Per Bruun Brockhoff ([email protected]) Mixed Linear Models, Module 12 Fall 2014 19 / 24 Time difference The semi-variogram The semi-variogram The semi-variogram The semi-variogram model3 <- lme(lnc month + treatm + month:treatm, random = 1 | cage, correlation = corExp(form =˜as.numeric(month) | cage, nugget = T), data = rats) print(plot(Variogram(model3, form =˜as.numeric(month)|cage, data = rats))) 1.0 model4 <- lme(lnc month + treatm + month:treatm, random = 1 | cage, correlation = corAR1(form = as.numeric(month) | cage), data = rats) print(plot(Variogram(model4, form =˜as.numeric(month)|cage, data = rats))) ● ● 1.0 ● ● ● 0.8 ● ● ● Semivariogram Semivariogram 0.8 ● 0.6 ● 0.4 0.2 ● ● ● ● 0.6 ● 0.4 0.2 ● ● 0.0 2 2 4 6 8 Mixed Linear Models, Module 12 6 8 Distance Distance Per Bruun Brockhoff ([email protected]) 4 Fall 2014 20 / 24 Per Bruun Brockhoff ([email protected]) Mixed Linear Models, Module 12 Fall 2014 21 / 24 Analysing the time structure by polynomial regression Analysing the time structure by polynomial regression Analysing the time structure by polynomial regression Overview of this module 1 Do the factor based analysis as shown above. 2 Explorative plotting of individual and treatment average regression lines/curves. 3 Potentially make a ”high degree” decomposition based on the simple ”split-plot” repeated measures model using lmer and lmerTest. 4 Check if a linear or quadratic regression model could be used as an alternative to the factor based model 5 IF a regression approach seems to capture what is going on, then try to fit the random coefficient model as an alternativ to the correlation structure used from above - chose the best one at the end. 6 A possibility is that a factor based model is needed for the main effect of time, whereas a quantitative model would fit the interaction effect. Per Bruun Brockhoff ([email protected]) Mixed Linear Models, Module 12 Fall 2014 23 / 24 1 Different view on the random effects approach Example: Activity of rats 2 Gaussian model of spatial correlation Example: Activity of rats 3 Other spatial correlation structures 4 Diagram of analysis 5 The semi-variogram 6 Analysing the time structure by polynomial regression Per Bruun Brockhoff ([email protected]) Mixed Linear Models, Module 12 Fall 2014 24 / 24