Designing Experiments and Analyzing Data: A Model Comparison Perspective (3rd edition) by Maxwell, Delaney, & Kelley
Information about the book is available at https://designingexperiments.com
R Code and Instructions to Accompany Chapter 15
These notes build on the prior chapters’ R Code and Instructions, where we discuss some operational aspects of R and detail the use of many functions and how to install packages. Here, we skip details that have already been discussed in the prior R Code and Instructions.
Load the AMCP package to use the data from Table 1, Chapter 15.
library(AMCP)
data(chapter_15_table_1)
chapter_15_table_1## Months30 Months36 Months42 Months48
## 1 108 96 110 122
## 2 103 117 127 133
## 3 96 107 106 107
## 4 84 85 92 99
## 5 118 125 125 116
## 6 110 107 96 91
## 7 129 128 123 128
## 8 90 84 101 113
## 9 84 104 100 88
## 10 96 100 103 105
## 11 105 114 105 112
## 12 113 117 132 130Margin means
# Means by row (individual here)
apply(chapter_15_table_1, MARGIN=1, FUN=mean)## [1] 109 120 104 90 121 101 127 97 94 101 109 123# Means by column (measurement occasion here)
apply(chapter_15_table_1, MARGIN=2, FUN=mean)## Months30 Months36 Months42 Months48
## 103 107 110 112First, we need to convert the structure of the data from “wide” or “multivariate” to “long” or “univariate.” We find the easiest way to do this is using the tidyr and dplyr packages.
library(tidyr)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union# Add an ID variable (by row number)
chapter_15_table_1$ID <- 1:nrow(chapter_15_table_1)
# Note that the "sep=-3" partitions between the 3rd and 2nd value from the end.
# The "-ID" here means to leave ID as it is (i.e., without "gathering" it)
data = chapter_15_table_1 %>%
gather(key=Unit.Month, value=Values, -ID) %>%
separate(col=Unit.Month, into=c("Unit", "Month"), sep=-3) %>%
arrange(ID, Month)
# Display the data
data## ID Unit Month Values
## 1 1 Months 30 108
## 2 1 Months 36 96
## 3 1 Months 42 110
## 4 1 Months 48 122
## 5 2 Months 30 103
## 6 2 Months 36 117
## 7 2 Months 42 127
## 8 2 Months 48 133
## 9 3 Months 30 96
## 10 3 Months 36 107
## 11 3 Months 42 106
## 12 3 Months 48 107
## 13 4 Months 30 84
## 14 4 Months 36 85
## 15 4 Months 42 92
## 16 4 Months 48 99
## 17 5 Months 30 118
## 18 5 Months 36 125
## 19 5 Months 42 125
## 20 5 Months 48 116
## 21 6 Months 30 110
## 22 6 Months 36 107
## 23 6 Months 42 96
## 24 6 Months 48 91
## 25 7 Months 30 129
## 26 7 Months 36 128
## 27 7 Months 42 123
## 28 7 Months 48 128
## 29 8 Months 30 90
## 30 8 Months 36 84
## 31 8 Months 42 101
## 32 8 Months 48 113
## 33 9 Months 30 84
## 34 9 Months 36 104
## 35 9 Months 42 100
## 36 9 Months 48 88
## 37 10 Months 30 96
## 38 10 Months 36 100
## 39 10 Months 42 103
## 40 10 Months 48 105
## 41 11 Months 30 105
## 42 11 Months 36 114
## 43 11 Months 42 105
## 44 11 Months 48 112
## 45 12 Months 30 113
## 46 12 Months 36 117
## 47 12 Months 42 132
## 48 12 Months 48 130#Note that we do not need "Unit" variable, so we drop it.
data$Unit <- NULL
# Examining the structure of the data.
str(data)## 'data.frame': 48 obs. of 3 variables:
## $ ID : int 1 1 1 1 2 2 2 2 3 3 ...
## $ Month : chr "30" "36" "42" "48" ...
## $ Values: int 108 96 110 122 103 117 127 133 96 107 ...# Examining the structure of the data we saw that "Month" is being treated as a character variable,
# yet we want it to be numeric:
data$Month <- as.numeric(data$Month)
str(data)## 'data.frame': 48 obs. of 3 variables:
## $ ID : int 1 1 1 1 2 2 2 2 3 3 ...
## $ Month : num 30 36 42 48 30 36 42 48 30 36 ...
## $ Values: int 108 96 110 122 103 117 127 133 96 107 ...There are multiple ways of producing trajectory plots. One convenient way is using the xyplot() function from the lattice. The lattice package provides a great deal of flexibility.
library(lattice)
# (each ID has it's own plot)
xyplot(Values ~ Month | ID, data=data, type="l", main="Trajectory Plot of McCarthy Scores", ylab="ID", xlab="Time in Months") 
Another way of plotting indivivdual trajectories is with ggplot2. The ggplot2 package provides a great deal of flexibility.
library(ggplot2)
Base.Plot <- ggplot(data = data, aes(x = Month, y = Values, group = ID))
Base.Plot + geom_line()
Yes another way to observe the trajectories is using vit() function from the MBESS package, which is easy to use but much less flexible than are the lattice or ggplot2 packages.
library(MBESS)
vit(id = "ID", occasion = "Month", score = "Values", Data = data)
#Or, without points
vit(id = "ID", occasion = "Month", score = "Values", pch="", Data = data)
A general linear model (GLM) appraoch to the data from Table 15.1 (compare this to Table 15.2 in MDK).
Model.GLM <- lm(Values ~ as.factor(Month) + as.factor(ID), data=data)
summary(Model.GLM)##
## Call:
## lm(formula = Values ~ as.factor(Month) + as.factor(ID), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.00 -4.25 0.00 4.25 14.00
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.040e+02 4.358e+00 23.862 < 2e-16 ***
## as.factor(Month)36 4.000e+00 3.183e+00 1.257 0.21769
## as.factor(Month)42 7.000e+00 3.183e+00 2.199 0.03498 *
## as.factor(Month)48 9.000e+00 3.183e+00 2.828 0.00791 **
## as.factor(ID)2 1.100e+01 5.513e+00 1.995 0.05433 .
## as.factor(ID)3 -5.000e+00 5.513e+00 -0.907 0.37102
## as.factor(ID)4 -1.900e+01 5.513e+00 -3.446 0.00157 **
## as.factor(ID)5 1.200e+01 5.513e+00 2.177 0.03676 *
## as.factor(ID)6 -8.000e+00 5.513e+00 -1.451 0.15619
## as.factor(ID)7 1.800e+01 5.513e+00 3.265 0.00255 **
## as.factor(ID)8 -1.200e+01 5.513e+00 -2.177 0.03676 *
## as.factor(ID)9 -1.500e+01 5.513e+00 -2.721 0.01031 *
## as.factor(ID)10 -8.000e+00 5.513e+00 -1.451 0.15619
## as.factor(ID)11 -6.527e-15 5.513e+00 0.000 1.00000
## as.factor(ID)12 1.400e+01 5.513e+00 2.539 0.01600 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.797 on 33 degrees of freedom
## Multiple R-squared: 0.7815, Adjusted R-squared: 0.6888
## F-statistic: 8.432 on 14 and 33 DF, p-value: 2.753e-07anova(Model.GLM)## Analysis of Variance Table
##
## Response: Values
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Month) 3 552 184.00 3.0269 0.04322 *
## as.factor(ID) 11 6624 602.18 9.9063 1.314e-07 ***
## Residuals 33 2006 60.79
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The lme4 package is a widely used package for mixed-effects models or hierarchical linear models, or what is often termed multilevel models, in R. Notice, below, that the fixed effects are provided (the intercept is included by default) and the random effects are included in parentheses, with an pipe symbol (i.e., “|”) used to denote the nesting strucutre. Here, the “(1|ID)” part, in words, means that the “intercept is nested within ID,” which means that there is a random intercept for each ID. We fit the model in two ways, one where Month is treated as a factor (directly below) and in the a future code-chunk, Month is treated continuous. To be clear, notice that the way SAS (as shown in the book) and R code the factor is different, as SAS takes the last ID whereas R takes the first, although the model parameterizations are done differently, they are equivalent. Compare this to Table 15.3 in MDK.
library(lme4)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following object is masked from 'package:tidyr':
##
## expandModel.Cat <- lmer(Values ~ as.factor(Month) + (1|ID), data=data)
summary(Model.Cat)## Linear mixed model fit by REML ['lmerMod']
## Formula: Values ~ as.factor(Month) + (1 | ID)
## Data: data
##
## REML criterion at convergence: 340.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.88627 -0.44709 -0.09063 0.53924 1.70501
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 135.35 11.634
## Residual 60.79 7.797
## Number of obs: 48, groups: ID, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 103.000 4.043 25.477
## as.factor(Month)36 4.000 3.183 1.257
## as.factor(Month)42 7.000 3.183 2.199
## as.factor(Month)48 9.000 3.183 2.828
##
## Correlation of Fixed Effects:
## (Intr) a.(M)3 a.(M)42
## as.fct(M)36 -0.394
## as.fct(M)42 -0.394 0.500
## as.fct(M)48 -0.394 0.500 0.500anova(Model.Cat)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## as.factor(Month) 3 552 184 3.0269Here we show selected outcomes from the model fit.
# Estimated coefficients (recall that was random) and slope (recall that was only fixed).
coef(Model.Cat)## $ID
## (Intercept) as.factor(Month)36 as.factor(Month)42 as.factor(Month)48
## 1 103.89905 4 7 9
## 2 113.78865 4 7 9
## 3 99.40378 4 7 9
## 4 86.81703 4 7 9
## 5 114.68770 4 7 9
## 6 96.70662 4 7 9
## 7 120.08203 4 7 9
## 8 93.11041 4 7 9
## 9 90.41324 4 7 9
## 10 96.70662 4 7 9
## 11 103.89905 4 7 9
## 12 116.48581 4 7 9
##
## attr(,"class")
## [1] "coef.mer"# Random effects (recall these are deviations from the fixed effect, of only intercept here).
ranef(Model.Cat)## $ID
## (Intercept)
## 1 0.8990539
## 2 10.7886474
## 3 -3.5962158
## 4 -16.1829711
## 5 11.6877013
## 6 -6.2933776
## 7 17.0820250
## 8 -9.8895934
## 9 -12.5867553
## 10 -6.2933776
## 11 0.8990539
## 12 13.4858092# Predicted values.
predict(Model.Cat)## 1 2 3 4 5 6 7
## 103.89905 107.89905 110.89905 112.89905 113.78865 117.78865 120.78865
## 8 9 10 11 12 13 14
## 122.78865 99.40378 103.40378 106.40378 108.40378 86.81703 90.81703
## 15 16 17 18 19 20 21
## 93.81703 95.81703 114.68770 118.68770 121.68770 123.68770 96.70662
## 22 23 24 25 26 27 28
## 100.70662 103.70662 105.70662 120.08203 124.08203 127.08203 129.08203
## 29 30 31 32 33 34 35
## 93.11041 97.11041 100.11041 102.11041 90.41324 94.41324 97.41324
## 36 37 38 39 40 41 42
## 99.41324 96.70662 100.70662 103.70662 105.70662 103.89905 107.89905
## 43 44 45 46 47 48
## 110.89905 112.89905 116.48581 120.48581 123.48581 125.48581# Outputs the data used in the model fit. This is useful for more complicated models with missing data.
model.frame(Model.Cat)## Values as.factor(Month) ID
## 1 108 30 1
## 2 96 36 1
## 3 110 42 1
## 4 122 48 1
## 5 103 30 2
## 6 117 36 2
## 7 127 42 2
## 8 133 48 2
## 9 96 30 3
## 10 107 36 3
## 11 106 42 3
## 12 107 48 3
## 13 84 30 4
## 14 85 36 4
## 15 92 42 4
## 16 99 48 4
## 17 118 30 5
## 18 125 36 5
## 19 125 42 5
## 20 116 48 5
## 21 110 30 6
## 22 107 36 6
## 23 96 42 6
## 24 91 48 6
## 25 129 30 7
## 26 128 36 7
## 27 123 42 7
## 28 128 48 7
## 29 90 30 8
## 30 84 36 8
## 31 101 42 8
## 32 113 48 8
## 33 84 30 9
## 34 104 36 9
## 35 100 42 9
## 36 88 48 9
## 37 96 30 10
## 38 100 36 10
## 39 103 42 10
## 40 105 48 10
## 41 105 30 11
## 42 114 36 11
## 43 105 42 11
## 44 112 48 11
## 45 113 30 12
## 46 117 36 12
## 47 132 42 12
## 48 130 48 12Here we combine the data used and the predicted values into Data.and.Predicted
Data.and.Predicted.Cat <- cbind(model.frame(Model.Cat), Predicted=predict(Model.Cat))
# Here we modify how lmer() has named month, due to how it is treated as a factor.
names(Data.and.Predicted.Cat)## [1] "Values" "as.factor(Month)" "ID"
## [4] "Predicted"names(Data.and.Predicted.Cat)[2] <- "Month"
# Here we turn month (back) into a numeric representation of what, in this model, was treated categorically.
Data.and.Predicted.Cat$Month <- as.numeric(Data.and.Predicted.Cat$Month)Now, like before, we use vit(). However, we can plot the predicted outcomes.
vit(id = "ID", occasion = "Month", score = "Predicted",
pch="", Data = Data.and.Predicted.Cat)
Random intercept only model with no effect of time.
Model.Int.Only <- lmer(Values ~ 1 + (1|ID), data=data)
summary(Model.Int.Only)## Linear mixed model fit by REML ['lmerMod']
## Formula: Values ~ 1 + (1 | ID)
## Data: data
##
## REML criterion at convergence: 361.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8488 -0.6912 0.1359 0.6088 1.7441
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 132.78 11.523
## Residual 71.06 8.429
## Number of obs: 48, groups: ID, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 108.000 3.542 30.49Data.and.Predicted.Int.Only <- cbind(Month=data$Month, model.frame(Model.Int.Only), Predicted=predict(Model.Int.Only))
vit(id = "ID", occasion = "Month", score = "Predicted", pch="", Data = Data.and.Predicted.Int.Only)
Here, for the above model, we can view the mean of the values and the prediction. Note that the predicted values for this model have a standard deviation of zero, as it is a means only model.
Data.and.Predicted.Int.Only %>%
group_by(ID) %>%
summarize(
Mean.Values = mean(Values),
S.Values = sd(Values),
Mean.Predicted = mean(Predicted),
S.Predicted = sd(Predicted))## # A tibble: 12 × 5
## ID Mean.Values S.Values Mean.Predicted S.Predicted
## <int> <dbl> <dbl> <dbl> <dbl>
## 1 1 109 10.645813 108.88200 0
## 2 2 120 13.114877 118.58404 0
## 3 3 104 5.354126 104.47199 0
## 4 4 90 6.976150 92.12394 0
## 5 5 121 4.690416 119.46604 0
## 6 6 101 8.981462 101.82598 0
## 7 7 127 2.708013 124.75806 0
## 8 8 97 12.780193 98.29797 0
## 9 9 94 9.521905 95.65196 0
## 10 10 101 3.915780 101.82598 0
## 11 11 109 4.690416 108.88200 0
## 12 12 123 9.416298 121.23005 0Compare the models using anova()
anova(Model.Int.Only, Model.Cat)## refitting model(s) with ML (instead of REML)## Data: data
## Models:
## Model.Int.Only: Values ~ 1 + (1 | ID)
## Model.Cat: Values ~ as.factor(Month) + (1 | ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Model.Int.Only 3 371.47 377.08 -182.73 365.47
## Model.Cat 6 368.71 379.94 -178.36 356.71 8.751 3 0.03279 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Treat Time Quantitatively.
We begin by treating Month as a factor for the fixed effects but as a continuous value for the random effect. Compare this to Table 15.3 in MDK (see especially the F-value from anova()).
Model.0 <- lmer(Values ~ as.factor(Month) + (1 + Month|ID), data=data)
summary(Model.0)## Linear mixed model fit by REML ['lmerMod']
## Formula: Values ~ as.factor(Month) + (1 + Month | ID)
## Data: data
##
## REML criterion at convergence: 336.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.81430 -0.47911 -0.05763 0.57397 1.53583
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 712.8192 26.6987
## Month 0.3892 0.6239 -0.90
## Residual 37.4364 6.1185
## Number of obs: 48, groups: ID, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 103.000 4.138 24.892
## as.factor(Month)36 4.000 2.722 1.470
## as.factor(Month)42 7.000 3.303 2.119
## as.factor(Month)48 9.000 4.092 2.199
##
## Correlation of Fixed Effects:
## (Intr) a.(M)3 a.(M)42
## as.fct(M)36 -0.421
## as.fct(M)42 -0.465 0.607
## as.fct(M)48 -0.471 0.595 0.749anova(Model.0)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## as.factor(Month) 3 200.07 66.689 1.7814Now, we plot the predicted scores with “Model.0,” which considers the fixed effect of month but a random slope for month.
Data.and.Predicted.0 <- cbind(model.frame(Model.0), Predicted=predict(Model.0))
vit(id = "ID", occasion = "Month", score = "Predicted", pch="", Data = Data.and.Predicted.0)
Here a model with fixed effects for linear, quadratic, and cubic effects are included, along with random effects for intercept and the linear effect.
Model.Trends <- lmer(Values ~ Month + I(Month^2) + I(Month^3) + (1+Month|ID), data=data)## Warning: Some predictor variables are on very different scales: consider
## rescalingsummary(Model.Trends)## Linear mixed model fit by REML ['lmerMod']
## Formula: Values ~ Month + I(Month^2) + I(Month^3) + (1 + Month | ID)
## Data: data
##
## REML criterion at convergence: 362.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.81429 -0.47911 -0.05762 0.57397 1.53581
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 712.7799 26.6979
## Month 0.3892 0.6238 -0.90
## Residual 37.4370 6.1186
## Number of obs: 48, groups: ID, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 6.800e+01 3.460e+02 0.197
## Month 1.583e+00 2.743e+01 0.058
## I(Month^2) -1.389e-02 7.135e-01 -0.019
## I(Month^3) 3.527e-13 6.095e-03 0.000
##
## Correlation of Fixed Effects:
## (Intr) Month I(M^2)
## Month -0.999
## I(Month^2) 0.997 -0.999
## I(Month^3) -0.994 0.998 -0.999
## fit warnings:
## Some predictor variables are on very different scales: consider rescalingNow we treat Month is treated as a quantitative variable in the straight-line change model formulation. Notice here that the covariance parameters (i.e., the variance of the random effects) is similar to that reported in the book. By default, the lmer() function uses a variation of maximum likelihood estimation termed “restricted information maximum likelihood.” Further, the way in which SAS and the lmer() function implement the estimation procedures are diffrent.
Model <- lmer(Values ~ Month + (1|ID), data=data)
summary(Model)## Linear mixed model fit by REML ['lmerMod']
## Formula: Values ~ Month + (1 | ID)
## Data: data
##
## REML criterion at convergence: 354.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.99786 -0.49141 -0.03432 0.54885 1.68963
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 136.13 11.668
## Residual 57.66 7.593
## Number of obs: 48, groups: ID, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 88.5000 7.2901 12.14
## Month 0.5000 0.1634 3.06
##
## Correlation of Fixed Effects:
## (Intr)
## Month -0.874Notice that the above output does not equal what is displayed in the book, as the intercept there is 103.5 but 88.5 here. Recall the interpretation of the intercept. This is due to the scaling of time. Here we will create a new time variable that simply subtracts the minimum value of time from all time values so that the intercept represents the mean at the initial timepoint.
data$Months.from.30 <- data$Month - 30
Model.30 <- lmer(Values ~ Months.from.30 + (1|ID), data=data)
summary(Model.30)## Linear mixed model fit by REML ['lmerMod']
## Formula: Values ~ Months.from.30 + (1 | ID)
## Data: data
##
## REML criterion at convergence: 354.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.99786 -0.49141 -0.03432 0.54885 1.68963
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 136.13 11.668
## Residual 57.66 7.593
## Number of obs: 48, groups: ID, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 103.5000 3.8350 26.99
## Months.from.30 0.5000 0.1634 3.06
##
## Correlation of Fixed Effects:
## (Intr)
## Mnths.fr.30 -0.383We can use “full information maximum likelihood” with REML=FALSE
Model.30_ML <- lmer(Values ~ Months.from.30 + (1|ID), REML=FALSE, data=data)
summary(Model.30_ML)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: Values ~ Months.from.30 + (1 | ID)
## Data: data
##
## AIC BIC logLik deviance df.resid
## 364.9 372.4 -178.5 356.9 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.03163 -0.49276 -0.03287 0.55799 1.70818
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 123.99 11.135
## Residual 56.06 7.487
## Number of obs: 48, groups: ID, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 103.5000 3.6881 28.063
## Months.from.30 0.5000 0.1611 3.104
##
## Correlation of Fixed Effects:
## (Intr)
## Mnths.fr.30 -0.393Here we show selected outcomes from the model fit using full information maximum likelihood.
# Estimated coefficients (recall that was random) and slope (recall that was only fixed).
coef(Model.30_ML)## $ID
## (Intercept) Months.from.30
## 1 104.39845 0.5
## 2 114.28140 0.5
## 3 99.90620 0.5
## 4 87.32790 0.5
## 5 115.17985 0.5
## 6 97.21085 0.5
## 7 120.57055 0.5
## 8 93.61705 0.5
## 9 90.92170 0.5
## 10 97.21085 0.5
## 11 104.39845 0.5
## 12 116.97675 0.5
##
## attr(,"class")
## [1] "coef.mer"# Random effects (recall these are deviations from the fixed effect, of only intercept here).
ranef(Model.30_ML)## $ID
## (Intercept)
## 1 0.8984501
## 2 10.7814010
## 3 -3.5938003
## 4 -16.1721015
## 5 11.6798511
## 6 -6.2891506
## 7 17.0705516
## 8 -9.8829509
## 9 -12.5783012
## 10 -6.2891506
## 11 0.8984501
## 12 13.4767513# Predicted values.
predict(Model.30_ML)## 1 2 3 4 5 6 7
## 104.39845 107.39845 110.39845 113.39845 114.28140 117.28140 120.28140
## 8 9 10 11 12 13 14
## 123.28140 99.90620 102.90620 105.90620 108.90620 87.32790 90.32790
## 15 16 17 18 19 20 21
## 93.32790 96.32790 115.17985 118.17985 121.17985 124.17985 97.21085
## 22 23 24 25 26 27 28
## 100.21085 103.21085 106.21085 120.57055 123.57055 126.57055 129.57055
## 29 30 31 32 33 34 35
## 93.61705 96.61705 99.61705 102.61705 90.92170 93.92170 96.92170
## 36 37 38 39 40 41 42
## 99.92170 97.21085 100.21085 103.21085 106.21085 104.39845 107.39845
## 43 44 45 46 47 48
## 110.39845 113.39845 116.97675 119.97675 122.97675 125.97675# Outputs the data used in the model fit. This is useful for more complicated models with missing data.
model.frame(Model.30_ML)## Values Months.from.30 ID
## 1 108 0 1
## 2 96 6 1
## 3 110 12 1
## 4 122 18 1
## 5 103 0 2
## 6 117 6 2
## 7 127 12 2
## 8 133 18 2
## 9 96 0 3
## 10 107 6 3
## 11 106 12 3
## 12 107 18 3
## 13 84 0 4
## 14 85 6 4
## 15 92 12 4
## 16 99 18 4
## 17 118 0 5
## 18 125 6 5
## 19 125 12 5
## 20 116 18 5
## 21 110 0 6
## 22 107 6 6
## 23 96 12 6
## 24 91 18 6
## 25 129 0 7
## 26 128 6 7
## 27 123 12 7
## 28 128 18 7
## 29 90 0 8
## 30 84 6 8
## 31 101 12 8
## 32 113 18 8
## 33 84 0 9
## 34 104 6 9
## 35 100 12 9
## 36 88 18 9
## 37 96 0 10
## 38 100 6 10
## 39 103 12 10
## 40 105 18 10
## 41 105 0 11
## 42 114 6 11
## 43 105 12 11
## 44 112 18 11
## 45 113 0 12
## 46 117 6 12
## 47 132 12 12
## 48 130 18 12Here we combine the data used and the predicted values into Data.and.Predicted
Data.and.Predicted <- cbind(model.frame(Model.30_ML), Predicted=predict(Model.30_ML))Now, like before, we use vit(). However, we can plot the predicted outcomes.
vit(id = "ID", occasion = "Months.from.30", score = "Predicted", pch="", Data = Data.and.Predicted, xlab="Number of Months From Baseline")
Now we fit a model in which there is a random intercept and a random linear trend (notice the parenthetical term within the lmer() function now has 1+Month). Note that this model treats Month as a numeric variable, whereas the output from Table 15.4 in MDK treats time as a categorical variable for the fixed effects but quantitatively for the random effects. Compare this output to Table 15.6 in MDK.
Model.Linear <- lmer(Values ~ Month + (1 + Month|ID), data=data)
summary(Model.Linear)## Linear mixed model fit by REML ['lmerMod']
## Formula: Values ~ Month + (1 + Month | ID)
## Data: data
##
## REML criterion at convergence: 349.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.79182 -0.43847 0.00418 0.64190 1.68326
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 735.5977 27.1219
## Month 0.4037 0.6354 -0.90
## Residual 34.8169 5.9006
## Number of obs: 48, groups: ID, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 88.5000 9.3028 9.513
## Month 0.5000 0.2231 2.241
##
## Correlation of Fixed Effects:
## (Intr)
## Month -0.925anova(Model.Linear)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Month 1 174.91 174.91 5.0237Here we combine the data used and the predicted values into Data.and.Predicted
Data.and.Predicted.Linear <- cbind(model.frame(Model.Linear), Predicted=predict(Model.Linear))As before, we use use vit(). There are now individual differences in both the intercepts and slopes.
vit(id = "ID", occasion = "Month", score = "Predicted", pch="", Data = Data.and.Predicted.Linear)
Now we fit a model in which there are fixed effects for the intercept and each month (i.e., month is treated categorically) and a random intercept and linear trend. Importantly, note that this model treats Month as a categorical variable for the fixed effect and as a numeric variable for the random effect. This model maps onto Table 15.4 in MDK.
Model.Categorical.Fix_Linear.Random <- lmer(Values ~ as.factor(Month) + (1 + Month|ID), data=data)
summary(Model.Categorical.Fix_Linear.Random)## Linear mixed model fit by REML ['lmerMod']
## Formula: Values ~ as.factor(Month) + (1 + Month | ID)
## Data: data
##
## REML criterion at convergence: 336.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.81430 -0.47911 -0.05763 0.57397 1.53583
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 712.8192 26.6987
## Month 0.3892 0.6239 -0.90
## Residual 37.4364 6.1185
## Number of obs: 48, groups: ID, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 103.000 4.138 24.892
## as.factor(Month)36 4.000 2.722 1.470
## as.factor(Month)42 7.000 3.303 2.119
## as.factor(Month)48 9.000 4.092 2.199
##
## Correlation of Fixed Effects:
## (Intr) a.(M)3 a.(M)42
## as.fct(M)36 -0.421
## as.fct(M)42 -0.465 0.607
## as.fct(M)48 -0.471 0.595 0.749anova(Model.Categorical.Fix_Linear.Random)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## as.factor(Month) 3 200.07 66.689 1.7814Here we create a Month variable that is treated as a factor (thus, using this variable make it unnecessary to use as.factor(Month)).
data2 <- cbind(data, Month_Factor=as.factor(data$Month))
# Using str() we can see the structure of each variable.
str(data2)## 'data.frame': 48 obs. of 5 variables:
## $ ID : int 1 1 1 1 2 2 2 2 3 3 ...
## $ Month : num 30 36 42 48 30 36 42 48 30 36 ...
## $ Values : int 108 96 110 122 103 117 127 133 96 107 ...
## $ Months.from.30: num 0 6 12 18 0 6 12 18 0 6 ...
## $ Month_Factor : Factor w/ 4 levels "30","36","42",..: 1 2 3 4 1 2 3 4 1 2 ...Notice that data2 is now being used, as it include a variable for Month treated categorically (for fixed effect) as well as quantitatively (for random effect). These results map onto the first part of Table 15.5 (with summary()) and then the second part (with anova()). Finally, confidence intervals are given using the `confint’ function.
Model.Categorical.Fix_Linear.Random <- lmer(Values ~ Month_Factor + (1 + Month|ID), data=data2)
summary(Model.Categorical.Fix_Linear.Random)## Linear mixed model fit by REML ['lmerMod']
## Formula: Values ~ Month_Factor + (1 + Month | ID)
## Data: data2
##
## REML criterion at convergence: 336.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.81430 -0.47911 -0.05763 0.57397 1.53583
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 712.8192 26.6987
## Month 0.3892 0.6239 -0.90
## Residual 37.4364 6.1185
## Number of obs: 48, groups: ID, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 103.000 4.138 24.892
## Month_Factor36 4.000 2.722 1.470
## Month_Factor42 7.000 3.303 2.119
## Month_Factor48 9.000 4.092 2.199
##
## Correlation of Fixed Effects:
## (Intr) Mn_F36 Mn_F42
## Mnth_Fctr36 -0.421
## Mnth_Fctr42 -0.465 0.607
## Mnth_Fctr48 -0.471 0.595 0.749anova(Model.Categorical.Fix_Linear.Random)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Month_Factor 3 200.07 66.689 1.7814CIs <- confint(Model.Categorical.Fix_Linear.Random)## Computing profile confidence intervals ...CIs## 2.5 % 97.5 %
## .sig01 25.5395248 25.694276
## .sig02 -0.9742166 -0.271004
## .sig03 0.1888073 1.091053
## .sigma 4.5233976 8.005324
## (Intercept) 94.6324000 111.367601
## Month_Factor36 -1.2524813 9.252481
## Month_Factor42 0.5367615 13.463239
## Month_Factor48 0.7050039 17.294997Here we combine the data used and the predicted values into Data.and.Predicted. This figure maps onto Figure 15.5.
Data.and.Predicted.Linear_CFLR <- cbind(model.frame(Model.Categorical.Fix_Linear.Random),
Predicted=predict(Model.Categorical.Fix_Linear.Random))
vit(id = "ID", occasion = "Month", score = "Predicted", pch="",
Data = Data.and.Predicted.Linear_CFLR)
Here we fit a straight-line change model with random intercept and slope. This maps onto the first part of Table 15.6 (with summary()) and then the second part (with anova())
Model.SLCM <- lmer(Values ~ Month + (1 + Month|ID), data=data)
summary(Model.SLCM)## Linear mixed model fit by REML ['lmerMod']
## Formula: Values ~ Month + (1 + Month | ID)
## Data: data
##
## REML criterion at convergence: 349.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.79182 -0.43847 0.00418 0.64190 1.68326
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 735.5977 27.1219
## Month 0.4037 0.6354 -0.90
## Residual 34.8169 5.9006
## Number of obs: 48, groups: ID, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 88.5000 9.3028 9.513
## Month 0.5000 0.2231 2.241
##
## Correlation of Fixed Effects:
## (Intr)
## Month -0.925anova(Model.SLCM)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Month 1 174.91 174.91 5.0237We can very explicetely create categorical variables for the months (to illustrate a more step-by-step approach).
data2$Month30 <- as.numeric(data2$Month==30)
data2$Month36 <- as.numeric(data2$Month==36)
data2$Month42 <- as.numeric(data2$Month==42)
data2$Month48 <- as.numeric(data2$Month==48)Compare the output below to the SAS output (which uses a different reference group) in Table 15.4 of MDK.
Model.Categorical.Fix_Linear.Random_Categorical <- lmer(Values ~ Month30 + Month36 + Month42 +
(1 + Month|ID), data=data2)
summary(Model.Categorical.Fix_Linear.Random_Categorical)## Linear mixed model fit by REML ['lmerMod']
## Formula: Values ~ Month30 + Month36 + Month42 + (1 + Month | ID)
## Data: data2
##
## REML criterion at convergence: 336.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.81430 -0.47911 -0.05762 0.57397 1.53582
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 712.8091 26.6985
## Month 0.3892 0.6238 -0.90
## Residual 37.4365 6.1185
## Number of obs: 48, groups: ID, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 112.000 4.231 26.470
## Month30 -9.000 4.092 -2.199
## Month36 -5.000 3.303 -1.514
## Month42 -2.000 2.722 -0.735
##
## Correlation of Fixed Effects:
## (Intr) Mnth30 Mnth36
## Month30 -0.506
## Month36 -0.493 0.749
## Month42 -0.434 0.595 0.607anova(Model.Categorical.Fix_Linear.Random_Categorical)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Month30 1 112.283 112.283 2.9993
## Month36 1 67.571 67.571 1.8049
## Month42 1 20.217 20.217 0.5400# Using Month continuously; random intercept and slope.
No.Random_Con <- lmer(Values ~ Month + (1+Month|ID), data=data2)
summary(No.Random_Con)## Linear mixed model fit by REML ['lmerMod']
## Formula: Values ~ Month + (1 + Month | ID)
## Data: data2
##
## REML criterion at convergence: 349.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.79182 -0.43847 0.00418 0.64190 1.68326
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 735.5977 27.1219
## Month 0.4037 0.6354 -0.90
## Residual 34.8169 5.9006
## Number of obs: 48, groups: ID, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 88.5000 9.3028 9.513
## Month 0.5000 0.2231 2.241
##
## Correlation of Fixed Effects:
## (Intr)
## Month -0.925anova(No.Random_Con)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Month 1 174.91 174.91 5.0237# Using the Months categorically; random intercept only.
No.Random_Cat <- lmer(Values ~ Month30 + Month36 + Month42 + (1|ID), data=data2)
summary(No.Random_Cat)## Linear mixed model fit by REML ['lmerMod']
## Formula: Values ~ Month30 + Month36 + Month42 + (1 | ID)
## Data: data2
##
## REML criterion at convergence: 340.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.88627 -0.44709 -0.09063 0.53924 1.70501
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 135.35 11.634
## Residual 60.79 7.797
## Number of obs: 48, groups: ID, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 112.000 4.043 27.703
## Month30 -9.000 3.183 -2.828
## Month36 -5.000 3.183 -1.571
## Month42 -2.000 3.183 -0.628
##
## Correlation of Fixed Effects:
## (Intr) Mnth30 Mnth36
## Month30 -0.394
## Month36 -0.394 0.500
## Month42 -0.394 0.500 0.500anova(No.Random_Cat)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## Month30 1 400 400 6.5803
## Month36 1 128 128 2.1057
## Month42 1 24 24 0.3948Generalized Least Squares
The nlme package predates the lmer package and hasdifferent capabilities. The “nl” part of nlme is “nonlinear.” Thus, models that cannot be written as a linear function (e.g., of time) can often be fit using nlme. However, there are other benefits and capabilities in nlme, such as the generalized least squares function, which is for fitting linear models using generalized least squares. As noted in the helpful for gls(): “the errors are allowed to be correlated and/or have unequal variances.” With some of the alternative models discussed in MDK, the gls() function can be used. Interestingly, as the gls() model is not exactly a mixed effects model, it is included in a mixed-effects modeling package (i.e., nlme but also SAS PROC MIXED) in much the same way such models are included in the mixed-effects modeling section part of MDK (recall we discuss how it can be considered a special case of such models).
# install.packages("nlme")
library(nlme)##
## Attaching package: 'nlme'## The following object is masked from 'package:lme4':
##
## lmList## The following object is masked from 'package:dplyr':
##
## collapseHere the variance is estimated differently for each month (which is considered here a factor) using the weights parameter and the varIdent() function. In words, the one-sided formula says that the variance is allowed to differ for each value of month. Thus, there is not homogeneity of variance across months assumed with the model below.
Model.GLS.1 <- gls(Values ~ Month_Factor,
weights = varIdent(form = ~ 1 | Month_Factor), data = data2)
summary(Model.GLS.1)## Generalized least squares fit by REML
## Model: Values ~ Month_Factor
## Data: data2
## AIC BIC logLik
## 382.9476 397.2211 -183.4738
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Month_Factor
## Parameter estimates:
## 30 36 42 48
## 1.0000000 1.0328266 0.9730408 1.0768354
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 103 3.958114 26.022493 0.0000
## Month_Factor36 4 5.690236 0.702959 0.4858
## Month_Factor42 7 5.522681 1.267500 0.2116
## Month_Factor48 9 5.816643 1.547284 0.1290
##
## Correlation:
## (Intr) Mn_F36 Mn_F42
## Month_Factor36 -0.696
## Month_Factor42 -0.717 0.499
## Month_Factor48 -0.680 0.473 0.488
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.6254850 -0.6933168 0.0000000 0.7119397 1.8962448
##
## Residual standard error: 13.71131
## Degrees of freedom: 48 total; 44 residualanova(Model.GLS.1)## Denom. DF: 44
## numDF F-value p-value
## (Intercept) 1 2865.6639 <.0001
## Month_Factor 3 0.9345 0.4321Other information is available within the object Model.GLS.1. The specific objects within Model.GLS.1 can be identified with names().
names(Model.GLS.1)## [1] "modelStruct" "dims" "contrasts" "coefficients"
## [5] "varBeta" "sigma" "apVar" "logLik"
## [9] "numIter" "groups" "call" "method"
## [13] "fitted" "residuals" "parAssign" "na.action"In addition to the heterogenious variances across the months, we can allow a correlational structure by including the correlation parameter in gls() information is available within the object Model.GLS.1. The specific objects within Model.GLS.1 can be identified with names(). Ths “1 | ID” part of the code says, in words, that “the errors are allowed to correlate.” Thus, combining weights = varIdent(form = ~ 1 | Month_Factor) and correlation=corSymm(form = ~ 1 | ID) in a gls() call leads to an unstructured covariance matrix for the errors. Compare these results to Table 15.10 in MDK. Note that the correlational structure part of the output corresponds to the “estimated correlation matrix for Subject 1” at the bottom of Table 15.10. Further, note that the anova() function produces output that corresponds to the “Type 3 Tests of Fixed Effects” in Table 15.10 (though the \(p\)-values differ because of some differences in operationalization of certain ideas in SAS PROC MIXED as compared to nlme).
Model.GLS <- gls(Values ~ Month_Factor,
weights = varIdent(form = ~ 1 | Month_Factor),
correlation=corSymm(form = ~ 1 | ID),
data = data2)
summary(Model.GLS)## Generalized least squares fit by REML
## Model: Values ~ Month_Factor
## Data: data2
## AIC BIC logLik
## 353.1128 378.0915 -162.5564
##
## Correlation Structure: General
## Formula: ~1 | ID
## Parameter estimate(s):
## Correlation:
## 1 2 3
## 2 0.795
## 3 0.696 0.760
## 4 0.599 0.466 0.853
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Month_Factor
## Parameter estimates:
## 30 36 42 48
## 1.000000 1.032823 0.973040 1.076839
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 103 3.958101 26.022581 0.0000
## Month_Factor36 4 2.579053 1.550957 0.1281
## Month_Factor42 7 3.045127 2.298755 0.0263
## Month_Factor48 9 3.692759 2.437202 0.0189
##
## Correlation:
## (Intr) Mn_F36 Mn_F42
## Month_Factor36 -0.275
## Month_Factor42 -0.419 0.530
## Month_Factor48 -0.381 0.087 0.797
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.625485e+00 -6.933197e-01 1.965979e-15 7.119438e-01 1.896251e+00
##
## Residual standard error: 13.71126
## Degrees of freedom: 48 total; 44 residualanova(Model.GLS)## Denom. DF: 44
## numDF F-value p-value
## (Intercept) 1 934.3057 <.0001
## Month_Factor 3 2.7379 0.0547There are auxilary functions that can be useful. One such function is the getVarCov() function, which extracts the covariance matrix for the errors. See Table 15.8 in MDK.
getVarCov(Model.GLS)## Marginal variance covariance matrix
## [,1] [,2] [,3] [,4]
## [1,] 188.00 154.360 127.36 121.180
## [2,] 154.36 200.540 143.63 97.452
## [3,] 127.36 143.630 178.00 168.090
## [4,] 121.18 97.452 168.09 218.000
## Standard Deviations: 13.711 14.161 13.342 14.765The cov2cor() function can be used for turning the extracted covariance matrix into the correlation matrix. See Table 15.9 in MDK.
cov2cor(getVarCov(Model.GLS))## Marginal variance covariance matrix
## [,1] [,2] [,3] [,4]
## [1,] 1.00000 0.79498 0.69623 0.59859
## [2,] 0.79498 1.00000 0.76023 0.46608
## [3,] 0.69623 0.76023 1.00000 0.85331
## [4,] 0.59859 0.46608 0.85331 1.00000
## Standard Deviations: 1 1 1 1As discussed in the book, there are yet other correlational structures for the errors that can be used. Here we show a structure that is called compound symmetry (equal covariances and equal variances). See Table 15.11 in MDK. Notice that the corCompSymm() function is used, which specifies that compound symmetry be used. Note too that the weights parameter is not needed as above (as the variances are the same for each month).
Model.GLS.CS <- gls(Values ~ Month_Factor,
correlation=corCompSymm(form = ~ 1 | ID),
data = data2)
summary(Model.GLS.CS)## Generalized least squares fit by REML
## Model: Values ~ Month_Factor
## Data: data2
## AIC BIC logLik
## 352.7563 363.4614 -170.3781
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | ID
## Parameter estimate(s):
## Rho
## 0.6900734
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 103 4.042858 25.477030 0.0000
## Month_Factor36 4 3.182972 1.256687 0.2155
## Month_Factor42 7 3.182972 2.199202 0.0332
## Month_Factor48 9 3.182972 2.827546 0.0070
##
## Correlation:
## (Intr) Mn_F36 Mn_F42
## Month_Factor36 -0.394
## Month_Factor42 -0.394 0.500
## Month_Factor48 -0.394 0.500 0.500
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.713690e+00 -6.604846e-01 2.029416e-15 7.140374e-01 1.856497e+00
##
## Residual standard error: 14.00487
## Degrees of freedom: 48 total; 44 residualanova(Model.GLS.CS)## Denom. DF: 44
## numDF F-value p-value
## (Intercept) 1 929.7391 <.0001
## Month_Factor 3 3.0269 0.0394getVarCov(Model.GLS.CS)## Marginal variance covariance matrix
## [,1] [,2] [,3] [,4]
## [1,] 196.14 135.35 135.35 135.35
## [2,] 135.35 196.14 135.35 135.35
## [3,] 135.35 135.35 196.14 135.35
## [4,] 135.35 135.35 135.35 196.14
## Standard Deviations: 14.005 14.005 14.005 14.005Another structure that may be of interest is an autoregressive structure, which can be implemented with the corAR1() function. Compare these results with Table 15.12. More generally, for other predefined correlational structures, see the corClasses() function from the nlme package.
Model.GLS.AR <- gls(Values ~ Month_Factor,
correlation=corAR1(form = ~ 1 | ID),
data = data2)
summary(Model.GLS.AR)## Generalized least squares fit by REML
## Model: Values ~ Month_Factor
## Data: data2
## AIC BIC logLik
## 345.1387 355.8438 -166.5693
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.8064014
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 103 4.076796 25.264938 0.0000
## Month_Factor36 4 2.536793 1.576794 0.1220
## Month_Factor42 7 3.409512 2.053080 0.0460
## Month_Factor48 9 3.976123 2.263512 0.0286
##
## Correlation:
## (Intr) Mn_F36 Mn_F42
## Month_Factor36 -0.311
## Month_Factor42 -0.418 0.672
## Month_Factor48 -0.488 0.526 0.774
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.699424e+00 -6.549862e-01 -2.012522e-15 7.080931e-01 1.841042e+00
##
## Residual standard error: 14.12244
## Degrees of freedom: 48 total; 44 residualanova(Model.GLS.AR)## Denom. DF: 44
## numDF F-value p-value
## (Intercept) 1 921.6425 <.0001
## Month_Factor 3 1.8041 0.1604getVarCov(Model.GLS.AR)## Marginal variance covariance matrix
## [,1] [,2] [,3] [,4]
## [1,] 199.44 160.83 129.69 104.59
## [2,] 160.83 199.44 160.83 129.69
## [3,] 129.69 160.83 199.44 160.83
## [4,] 104.59 129.69 160.83 199.44
## Standard Deviations: 14.122 14.122 14.122 14.122