library(tidyr) #Allows for us to manipulate the data structure
library(data.table) #Allows for us to manipulate the data structure
library(lme4) #Allows us to fit mixed effects models
library(lattice) #for plotting random effects
library(nlme)6 Modeling the Covariance
Recall (from notes): Using the gls() function to fit a linear model using generalized least squares. The syntax for the function gls() in nlme package is
gls(model, data, correlation, weights, subset, method, na.action, control, verbose)
- Description
model: A two-sided linear formula object describing the modeldata: Optional dataframecorrelation: See description for details.- The
corStructobject describes the within-group correlation structure. - Note we use
corCompSymmwhich has the following syntax: corCompSymm(value, form, fixed) compound symmetry structure corresponding to a constant correlation.
- The
weights: See description for details.- Describes the variance function.
- Note we use
varClasseswhich defines standard classes of variance function structures and use “varIdent” which describes constant variance(s), that is generally used to allow different variances.
Model 1: \(E(y_{ij} \mid X_{ij}) = \beta_0 + \beta_{trt}trt + \beta_{T_1}T_1 + \beta_{T_4}T_4 + \beta_{T_6}T_6 + \beta_{trt:T_1}trt*T_1 + \beta_{trt:T_4}trt*T_4 + \beta_{trt:T_6}trt*T_6\).
This model includes time as a categorical covariate and also includes an interaction term to account and adjust for any possible discrepancies between the treatment and control.
# unstructured covariance
mod_1_unstruc <- gls(measurements ~ factor(group) + factor(level) + factor(group):factor(level),
data = long_TLC, corr = corSymm(form = ~1 | id), weights = varIdent(form = ~1 |
as.factor(level)), method = "REML")
summary(mod_1_unstruc)Generalized least squares fit by REML
Model: measurements ~ factor(group) + factor(level) + factor(group):factor(level)
Data: long_TLC
AIC BIC logLik
2452.076 2523.559 -1208.038
Correlation Structure: General
Formula: ~1 | id
Parameter estimate(s):
Correlation:
1 2 3
2 0.571
3 0.570 0.775
4 0.577 0.582 0.581
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | as.factor(level)
Parameter estimates:
lead0 lead1 lead4 lead6
1.000000 1.325877 1.370471 1.524816
Coefficients:
Value Std.Error t-value p-value
(Intercept) 26.272 0.7102892 36.98775 0.0000
factor(group)Treatment 0.268 1.0045006 0.26680 0.7898
factor(level)lead1 -1.612 0.7919151 -2.03557 0.0425
factor(level)lead4 -2.202 0.8149218 -2.70210 0.0072
factor(level)lead6 -2.626 0.8885198 -2.95548 0.0033
factor(group)Treatment:factor(level)lead1 -11.406 1.1199371 -10.18450 0.0000
factor(group)Treatment:factor(level)lead4 -8.824 1.1524735 -7.65658 0.0000
factor(group)Treatment:factor(level)lead6 -3.152 1.2565568 -2.50844 0.0125
Correlation:
(Intr) fct()T fct()1 fct()4 fct()6
factor(group)Treatment -0.707
factor(level)lead1 -0.218 0.154
factor(level)lead4 -0.191 0.135 0.680
factor(level)lead6 -0.096 0.068 0.386 0.385
factor(group)Treatment:factor(level)lead1 0.154 -0.218 -0.707 -0.481 -0.273
factor(group)Treatment:factor(level)lead4 0.135 -0.191 -0.481 -0.707 -0.272
factor(group)Treatment:factor(level)lead6 0.068 -0.096 -0.273 -0.272 -0.707
f()T:()1 f()T:()4
factor(group)Treatment
factor(level)lead1
factor(level)lead4
factor(level)lead6
factor(group)Treatment:factor(level)lead1
factor(group)Treatment:factor(level)lead4 0.680
factor(group)Treatment:factor(level)lead6 0.386 0.385
Standardized residuals:
Min Q1 Med Q3 Max
-2.1756516 -0.6849999 -0.1515556 0.5294195 5.6327729
Residual standard error: 5.022503
Degrees of freedom: 400 total; 392 residual# exchangeable (constant variances)
mod_1_exch_const <- gls(measurements ~ factor(group) + factor(level) + factor(group):factor(level),
data = long_TLC, corr = corCompSymm(form = ~factor(level) | id), weights = varIdent(form = ~1),
method = "REML")
summary(mod_1_exch_const)Generalized least squares fit by REML
Model: measurements ~ factor(group) + factor(level) + factor(group):factor(level)
Data: long_TLC
AIC BIC logLik
2480.621 2520.334 -1230.311
Correlation Structure: Compound symmetry
Formula: ~factor(level) | id
Parameter estimate(s):
Rho
0.5954401
Coefficients:
Value Std.Error t-value p-value
(Intercept) 26.272 0.9370175 28.037898 0.0000
factor(group)Treatment 0.268 1.3251428 0.202242 0.8398
factor(level)lead1 -1.612 0.8428574 -1.912542 0.0565
factor(level)lead4 -2.202 0.8428574 -2.612541 0.0093
factor(level)lead6 -2.626 0.8428574 -3.115592 0.0020
factor(group)Treatment:factor(level)lead1 -11.406 1.1919804 -9.568950 0.0000
factor(group)Treatment:factor(level)lead4 -8.824 1.1919804 -7.402807 0.0000
factor(group)Treatment:factor(level)lead6 -3.152 1.1919804 -2.644339 0.0085
Correlation:
(Intr) fct()T fct()1 fct()4 fct()6
factor(group)Treatment -0.707
factor(level)lead1 -0.450 0.318
factor(level)lead4 -0.450 0.318 0.500
factor(level)lead6 -0.450 0.318 0.500 0.500
factor(group)Treatment:factor(level)lead1 0.318 -0.450 -0.707 -0.354 -0.354
factor(group)Treatment:factor(level)lead4 0.318 -0.450 -0.354 -0.707 -0.354
factor(group)Treatment:factor(level)lead6 0.318 -0.450 -0.354 -0.354 -0.707
f()T:()1 f()T:()4
factor(group)Treatment
factor(level)lead1
factor(level)lead4
factor(level)lead6
factor(group)Treatment:factor(level)lead1
factor(group)Treatment:factor(level)lead4 0.500
factor(group)Treatment:factor(level)lead6 0.500 0.500
Standardized residuals:
Min Q1 Med Q3 Max
-2.5147478 -0.6973588 -0.1498707 0.5542798 6.5106947
Residual standard error: 6.625714
Degrees of freedom: 400 total; 392 residual# exchangeable (heterogeneous variances)
mod_1_exch_heter <- gls(measurements ~ factor(group) + factor(level) + factor(group):factor(level),
data = long_TLC, corr = corCompSymm(form = ~factor(level) | id), weights = varIdent(form = ~1 |
factor(level)), method = "REML")
summary(mod_1_exch_heter)Generalized least squares fit by REML
Model: measurements ~ factor(group) + factor(level) + factor(group):factor(level)
Data: long_TLC
AIC BIC logLik
2459.96 2511.587 -1216.98
Correlation Structure: Compound symmetry
Formula: ~factor(level) | id
Parameter estimate(s):
Rho
0.6102708
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | factor(level)
Parameter estimates:
lead0 lead1 lead4 lead6
1.000000 1.279672 1.323197 1.519220
Coefficients:
Value Std.Error t-value p-value
(Intercept) 26.272 0.7237951 36.29757 0.0000
factor(group)Treatment 0.268 1.0236008 0.26182 0.7936
factor(level)lead1 -1.612 0.7506796 -2.14739 0.0324
factor(level)lead4 -2.202 0.7713883 -2.85459 0.0045
factor(level)lead6 -2.626 0.8726932 -3.00908 0.0028
factor(group)Treatment:factor(level)lead1 -11.406 1.0616213 -10.74394 0.0000
factor(group)Treatment:factor(level)lead4 -8.824 1.0909078 -8.08868 0.0000
factor(group)Treatment:factor(level)lead6 -3.152 1.2341746 -2.55393 0.0110
Correlation:
(Intr) fct()T fct()1 fct()4 fct()6
factor(group)Treatment -0.707
factor(level)lead1 -0.211 0.149
factor(level)lead4 -0.181 0.128 0.402
factor(level)lead6 -0.060 0.043 0.383 0.383
factor(group)Treatment:factor(level)lead1 0.149 -0.211 -0.707 -0.285 -0.270
factor(group)Treatment:factor(level)lead4 0.128 -0.181 -0.285 -0.707 -0.271
factor(group)Treatment:factor(level)lead6 0.043 -0.060 -0.270 -0.271 -0.707
f()T:()1 f()T:()4
factor(group)Treatment
factor(level)lead1
factor(level)lead4
factor(level)lead6
factor(group)Treatment:factor(level)lead1
factor(group)Treatment:factor(level)lead4 0.402
factor(group)Treatment:factor(level)lead6 0.383 0.383
Standardized residuals:
Min Q1 Med Q3 Max
-2.1429199 -0.6927789 -0.1528885 0.5263121 5.5480303
Residual standard error: 5.118004
Degrees of freedom: 400 total; 392 residual