5 Chapter 5

#-*- R -*-

# initialization

library(nlme)
options(width = 65, digits = 5)
options(contrasts = c(unordered = "contr.helmert", ordered = "contr.poly"))
#pdf(file = "ch05.pdf")

# Chapter 5    Extending the Basic Linear Mixed-Effects Models

# 5.1 General Formulation of the Extended Model

vf1Fixed <- varFixed(~ age)
vf1Fixed <- Initialize(vf1Fixed, data = Orthodont)
varWeights(vf1Fixed)

##   [1] 0.35355 0.31623 0.28868 0.26726 0.35355 0.31623 0.28868
##   [8] 0.26726 0.35355 0.31623 0.28868 0.26726 0.35355 0.31623
##  [15] 0.28868 0.26726 0.35355 0.31623 0.28868 0.26726 0.35355
##  [22] 0.31623 0.28868 0.26726 0.35355 0.31623 0.28868 0.26726
##  [29] 0.35355 0.31623 0.28868 0.26726 0.35355 0.31623 0.28868
##  [36] 0.26726 0.35355 0.31623 0.28868 0.26726 0.35355 0.31623
##  [43] 0.28868 0.26726 0.35355 0.31623 0.28868 0.26726 0.35355
##  [50] 0.31623 0.28868 0.26726 0.35355 0.31623 0.28868 0.26726
##  [57] 0.35355 0.31623 0.28868 0.26726 0.35355 0.31623 0.28868
##  [64] 0.26726 0.35355 0.31623 0.28868 0.26726 0.35355 0.31623
##  [71] 0.28868 0.26726 0.35355 0.31623 0.28868 0.26726 0.35355
##  [78] 0.31623 0.28868 0.26726 0.35355 0.31623 0.28868 0.26726
##  [85] 0.35355 0.31623 0.28868 0.26726 0.35355 0.31623 0.28868
##  [92] 0.26726 0.35355 0.31623 0.28868 0.26726 0.35355 0.31623
##  [99] 0.28868 0.26726 0.35355 0.31623 0.28868 0.26726 0.35355
## [106] 0.31623 0.28868 0.26726

vf1Ident <- varIdent(c(Female = 0.5), ~ 1 | Sex)
vf1Ident <- Initialize(vf1Ident, Orthodont)
varWeights(vf1Ident)

##   Male   Male   Male   Male   Male   Male   Male   Male   Male 
##      1      1      1      1      1      1      1      1      1 
##   Male   Male   Male   Male   Male   Male   Male   Male   Male 
##      1      1      1      1      1      1      1      1      1 
##   Male   Male   Male   Male   Male   Male   Male   Male   Male 
##      1      1      1      1      1      1      1      1      1 
##   Male   Male   Male   Male   Male   Male   Male   Male   Male 
##      1      1      1      1      1      1      1      1      1 
##   Male   Male   Male   Male   Male   Male   Male   Male   Male 
##      1      1      1      1      1      1      1      1      1 
##   Male   Male   Male   Male   Male   Male   Male   Male   Male 
##      1      1      1      1      1      1      1      1      1 
##   Male   Male   Male   Male   Male   Male   Male   Male   Male 
##      1      1      1      1      1      1      1      1      1 
##   Male Female Female Female Female Female Female Female Female 
##      1      2      2      2      2      2      2      2      2 
## Female Female Female Female Female Female Female Female Female 
##      2      2      2      2      2      2      2      2      2 
## Female Female Female Female Female Female Female Female Female 
##      2      2      2      2      2      2      2      2      2 
## Female Female Female Female Female Female Female Female Female 
##      2      2      2      2      2      2      2      2      2 
## Female Female Female Female Female Female Female Female Female 
##      2      2      2      2      2      2      2      2      2

vf2Ident <- varIdent(form =  ~ 1 | Sex, fixed = c(Female = 0.5))
vf2Ident <- Initialize(vf2Ident, Orthodont)
varWeights(vf2Ident)

##   Male   Male   Male   Male   Male   Male   Male   Male   Male 
##      1      1      1      1      1      1      1      1      1 
##   Male   Male   Male   Male   Male   Male   Male   Male   Male 
##      1      1      1      1      1      1      1      1      1 
##   Male   Male   Male   Male   Male   Male   Male   Male   Male 
##      1      1      1      1      1      1      1      1      1 
##   Male   Male   Male   Male   Male   Male   Male   Male   Male 
##      1      1      1      1      1      1      1      1      1 
##   Male   Male   Male   Male   Male   Male   Male   Male   Male 
##      1      1      1      1      1      1      1      1      1 
##   Male   Male   Male   Male   Male   Male   Male   Male   Male 
##      1      1      1      1      1      1      1      1      1 
##   Male   Male   Male   Male   Male   Male   Male   Male   Male 
##      1      1      1      1      1      1      1      1      1 
##   Male Female Female Female Female Female Female Female Female 
##      1      2      2      2      2      2      2      2      2 
## Female Female Female Female Female Female Female Female Female 
##      2      2      2      2      2      2      2      2      2 
## Female Female Female Female Female Female Female Female Female 
##      2      2      2      2      2      2      2      2      2 
## Female Female Female Female Female Female Female Female Female 
##      2      2      2      2      2      2      2      2      2 
## Female Female Female Female Female Female Female Female Female 
##      2      2      2      2      2      2      2      2      2

vf3Ident <- varIdent(form = ~ 1 | Sex * age)
vf3Ident <- Initialize(vf3Ident, Orthodont)
varWeights(vf3Ident)

##    Male*8   Male*10   Male*12   Male*14    Male*8   Male*10 
##         1         1         1         1         1         1 
##   Male*12   Male*14    Male*8   Male*10   Male*12   Male*14 
##         1         1         1         1         1         1 
##    Male*8   Male*10   Male*12   Male*14    Male*8   Male*10 
##         1         1         1         1         1         1 
##   Male*12   Male*14    Male*8   Male*10   Male*12   Male*14 
##         1         1         1         1         1         1 
##    Male*8   Male*10   Male*12   Male*14    Male*8   Male*10 
##         1         1         1         1         1         1 
##   Male*12   Male*14    Male*8   Male*10   Male*12   Male*14 
##         1         1         1         1         1         1 
##    Male*8   Male*10   Male*12   Male*14    Male*8   Male*10 
##         1         1         1         1         1         1 
##   Male*12   Male*14    Male*8   Male*10   Male*12   Male*14 
##         1         1         1         1         1         1 
##    Male*8   Male*10   Male*12   Male*14    Male*8   Male*10 
##         1         1         1         1         1         1 
##   Male*12   Male*14    Male*8   Male*10   Male*12   Male*14 
##         1         1         1         1         1         1 
##    Male*8   Male*10   Male*12   Male*14  Female*8 Female*10 
##         1         1         1         1         1         1 
## Female*12 Female*14  Female*8 Female*10 Female*12 Female*14 
##         1         1         1         1         1         1 
##  Female*8 Female*10 Female*12 Female*14  Female*8 Female*10 
##         1         1         1         1         1         1 
## Female*12 Female*14  Female*8 Female*10 Female*12 Female*14 
##         1         1         1         1         1         1 
##  Female*8 Female*10 Female*12 Female*14  Female*8 Female*10 
##         1         1         1         1         1         1 
## Female*12 Female*14  Female*8 Female*10 Female*12 Female*14 
##         1         1         1         1         1         1 
##  Female*8 Female*10 Female*12 Female*14  Female*8 Female*10 
##         1         1         1         1         1         1 
## Female*12 Female*14  Female*8 Female*10 Female*12 Female*14 
##         1         1         1         1         1         1

vf1Power <- varPower(1)
formula(vf1Power)

## ~fitted(.)
## <environment: 0x00000000145abd40>

vf2Power <- varPower(fixed = 0.5)
vf3Power <- varPower(form = ~ fitted(.) | Sex,
  fixed = list(Male = 0.5, Female = 0))
vf1Exp <- varExp(form = ~ age | Sex, fixed = c(Female = 0))
vf1ConstPower <- varConstPower(power = 0.5,
      fixed = list(const = 1))
vf1Comb <- varComb(varIdent(c(Female = 0.5), ~ 1 | Sex),
                     varExp(1, ~ age))
vf1Comb <- Initialize(vf1Comb, Orthodont)
varWeights(vf1Comb)

##   [1] 3.3546e-04 4.5400e-05 6.1442e-06 8.3153e-07 3.3546e-04
##   [6] 4.5400e-05 6.1442e-06 8.3153e-07 3.3546e-04 4.5400e-05
##  [11] 6.1442e-06 8.3153e-07 3.3546e-04 4.5400e-05 6.1442e-06
##  [16] 8.3153e-07 3.3546e-04 4.5400e-05 6.1442e-06 8.3153e-07
##  [21] 3.3546e-04 4.5400e-05 6.1442e-06 8.3153e-07 3.3546e-04
##  [26] 4.5400e-05 6.1442e-06 8.3153e-07 3.3546e-04 4.5400e-05
##  [31] 6.1442e-06 8.3153e-07 3.3546e-04 4.5400e-05 6.1442e-06
##  [36] 8.3153e-07 3.3546e-04 4.5400e-05 6.1442e-06 8.3153e-07
##  [41] 3.3546e-04 4.5400e-05 6.1442e-06 8.3153e-07 3.3546e-04
##  [46] 4.5400e-05 6.1442e-06 8.3153e-07 3.3546e-04 4.5400e-05
##  [51] 6.1442e-06 8.3153e-07 3.3546e-04 4.5400e-05 6.1442e-06
##  [56] 8.3153e-07 3.3546e-04 4.5400e-05 6.1442e-06 8.3153e-07
##  [61] 3.3546e-04 4.5400e-05 6.1442e-06 8.3153e-07 6.7093e-04
##  [66] 9.0800e-05 1.2288e-05 1.6631e-06 6.7093e-04 9.0800e-05
##  [71] 1.2288e-05 1.6631e-06 6.7093e-04 9.0800e-05 1.2288e-05
##  [76] 1.6631e-06 6.7093e-04 9.0800e-05 1.2288e-05 1.6631e-06
##  [81] 6.7093e-04 9.0800e-05 1.2288e-05 1.6631e-06 6.7093e-04
##  [86] 9.0800e-05 1.2288e-05 1.6631e-06 6.7093e-04 9.0800e-05
##  [91] 1.2288e-05 1.6631e-06 6.7093e-04 9.0800e-05 1.2288e-05
##  [96] 1.6631e-06 6.7093e-04 9.0800e-05 1.2288e-05 1.6631e-06
## [101] 6.7093e-04 9.0800e-05 1.2288e-05 1.6631e-06 6.7093e-04
## [106] 9.0800e-05 1.2288e-05 1.6631e-06

fm1Dial.lme <-
    lme(rate ~(pressure + I(pressure^2) + I(pressure^3) + I(pressure^4))*QB,
        Dialyzer, ~ pressure + I(pressure^2))
fm1Dial.lme

## Linear mixed-effects model fit by REML
##   Data: Dialyzer 
##   Log-restricted-likelihood: -326.39
##   Fixed: rate ~ (pressure + I(pressure^2) + I(pressure^3) + I(pressure^4)) *      QB 
##       (Intercept)          pressure     I(pressure^2) 
##        -16.597990         88.673253        -42.731978 
##     I(pressure^3)     I(pressure^4)               QB1 
##          9.216490         -0.775615         -0.631728 
##      pressure:QB1 I(pressure^2):QB1 I(pressure^3):QB1 
##          0.310398          1.574171          0.050908 
## I(pressure^4):QB1 
##         -0.085964 
## 
## Random effects:
##  Formula: ~pressure + I(pressure^2) | Subject
##  Structure: General positive-definite, Log-Cholesky parametrization
##               StdDev Corr         
## (Intercept)   1.4987 (Intr) pressr
## pressure      4.9072 -0.507       
## I(pressure^2) 1.4739  0.311 -0.944
## Residual      1.8214              
## 
## Number of Observations: 140
## Number of Groups: 20

plot(fm1Dial.lme, resid(.) ~ pressure, abline = 0)

fm2Dial.lme <- update(fm1Dial.lme,
                        weights = varPower(form = ~ pressure))
fm2Dial.lme

## Linear mixed-effects model fit by REML
##   Data: Dialyzer 
##   Log-restricted-likelihood: -309.51
##   Fixed: rate ~ (pressure + I(pressure^2) + I(pressure^3) + I(pressure^4)) *      QB 
##       (Intercept)          pressure     I(pressure^2) 
##         -17.67989          93.71144         -49.18679 
##     I(pressure^3)     I(pressure^4)               QB1 
##          12.24531          -1.24264          -0.92080 
##      pressure:QB1 I(pressure^2):QB1 I(pressure^3):QB1 
##           1.35310           0.48033           0.49135 
## I(pressure^4):QB1 
##          -0.14627 
## 
## Random effects:
##  Formula: ~pressure + I(pressure^2) | Subject
##  Structure: General positive-definite, Log-Cholesky parametrization
##               StdDev Corr         
## (Intercept)   1.8571 (Intr) pressr
## pressure      5.3281 -0.522       
## I(pressure^2) 1.6483  0.362 -0.954
## Residual      1.2626              
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~pressure 
##  Parameter estimates:
##   power 
## 0.74933 
## Number of Observations: 140
## Number of Groups: 20

anova(fm1Dial.lme, fm2Dial.lme)

##             Model df    AIC    BIC  logLik   Test L.Ratio
## fm1Dial.lme     1 17 686.78 735.53 -326.39               
## fm2Dial.lme     2 18 655.01 706.63 -309.51 1 vs 2   33.77
##             p-value
## fm1Dial.lme        
## fm2Dial.lme  <.0001

plot(fm2Dial.lme, resid(., type = "p") ~ pressure,
     abline = 0)

intervals(fm2Dial.lme)

## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                       lower      est.     upper
## (Intercept)       -19.14806 -17.67989 -16.21172
## pressure           87.23117  93.71144 100.19171
## I(pressure^2)     -57.61614 -49.18679 -40.75744
## I(pressure^3)       7.96723  12.24531  16.52338
## I(pressure^4)      -1.95269  -1.24264  -0.53260
## QB1                -2.47755  -0.92080   0.63596
## pressure:QB1       -5.12717   1.35310   7.83337
## I(pressure^2):QB1  -7.94902   0.48033   8.90968
## I(pressure^3):QB1  -3.78673   0.49135   4.76942
## I(pressure^4):QB1  -0.85631  -0.14627   0.56378
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: Subject 
##                                   lower     est.     upper
## sd((Intercept))                 1.25583  1.85706  2.746115
## sd(pressure)                    3.62470  5.32810  7.832007
## sd(I(pressure^2))               1.09141  1.64827  2.489272
## cor((Intercept),pressure)      -0.80266 -0.52242 -0.053226
## cor((Intercept),I(pressure^2)) -0.16526  0.36227  0.728629
## cor(pressure,I(pressure^2))    -0.98500 -0.95390 -0.862778
## 
##  Variance function:
##         lower    est.   upper
## power 0.50801 0.74933 0.99065
## attr(,"label")
## [1] "Variance function:"
## 
##  Within-group standard error:
##  lower   est.  upper 
## 1.0607 1.2626 1.5029

plot(fm2Dial.lme, resid(.) ~ pressure|QB, abline = 0)

fm3Dial.lme <- update(fm2Dial.lme,
                      weights=varPower(form = ~ pressure | QB))
fm3Dial.lme

## Linear mixed-effects model fit by REML
##   Data: Dialyzer 
##   Log-restricted-likelihood: -309.15
##   Fixed: rate ~ (pressure + I(pressure^2) + I(pressure^3) + I(pressure^4)) *      QB 
##       (Intercept)          pressure     I(pressure^2) 
##         -17.69504          93.75943         -49.23064 
##     I(pressure^3)     I(pressure^4)               QB1 
##          12.26045          -1.24439          -1.01692 
##      pressure:QB1 I(pressure^2):QB1 I(pressure^3):QB1 
##           1.83955          -0.19437           0.82722 
## I(pressure^4):QB1 
##          -0.20035 
## 
## Random effects:
##  Formula: ~pressure + I(pressure^2) | Subject
##  Structure: General positive-definite, Log-Cholesky parametrization
##               StdDev Corr         
## (Intercept)   1.8240 (Intr) pressr
## pressure      5.2423 -0.502       
## I(pressure^2) 1.6372  0.338 -0.951
## Residual      1.2590              
## 
## Variance function:
##  Structure: Power of variance covariate, different strata
##  Formula: ~pressure | QB 
##  Parameter estimates:
##     200     300 
## 0.64775 0.83780 
## Number of Observations: 140
## Number of Groups: 20

anova(fm2Dial.lme, fm3Dial.lme)

##             Model df    AIC    BIC  logLik   Test L.Ratio
## fm2Dial.lme     1 18 655.01 706.63 -309.51               
## fm3Dial.lme     2 19 656.30 710.78 -309.15 1 vs 2 0.71091
##             p-value
## fm2Dial.lme        
## fm3Dial.lme  0.3991

fm4Dial.lme <- update(fm2Dial.lme,
                      weights = varConstPower(form = ~ pressure))
anova(fm2Dial.lme, fm4Dial.lme)

##             Model df    AIC    BIC  logLik   Test L.Ratio
## fm2Dial.lme     1 18 655.01 706.63 -309.51               
## fm4Dial.lme     2 19 656.85 711.34 -309.43 1 vs 2 0.15915
##             p-value
## fm2Dial.lme        
## fm4Dial.lme  0.6899

plot(augPred(fm2Dial.lme), grid = TRUE)

anova(fm2Dial.lme)

##                  numDF denDF F-value p-value
## (Intercept)          1   112  552.98  <.0001
## pressure             1   112 2328.82  <.0001
## I(pressure^2)        1   112 1174.68  <.0001
## I(pressure^3)        1   112  359.92  <.0001
## I(pressure^4)        1   112   12.47  0.0006
## QB                   1    18    4.83  0.0414
## pressure:QB          1   112   80.11  <.0001
## I(pressure^2):QB     1   112    1.35  0.2476
## I(pressure^3):QB     1   112    2.24  0.1370
## I(pressure^4):QB     1   112    0.17  0.6839

anova(fm2Dial.lme, Terms = 8:10)

## F-test for: I(pressure^2):QB, I(pressure^3):QB, I(pressure^4):QB 
##   numDF denDF F-value p-value
## 1     3   112  1.2535  0.2939

options(contrasts = c("contr.treatment", "contr.poly"))
fm1BW.lme <- lme(weight ~ Time * Diet, BodyWeight,
                   random = ~ Time)
fm1BW.lme

## Linear mixed-effects model fit by REML
##   Data: BodyWeight 
##   Log-restricted-likelihood: -575.86
##   Fixed: weight ~ Time * Diet 
## (Intercept)        Time       Diet2       Diet3  Time:Diet2 
##   251.65165     0.35964   200.66549   252.07168     0.60584 
##  Time:Diet3 
##     0.29834 
## 
## Random effects:
##  Formula: ~Time | Rat
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr  
## (Intercept) 36.93907 (Intr)
## Time         0.24841 -0.149
## Residual     4.44361       
## 
## Number of Observations: 176
## Number of Groups: 16

fm2BW.lme <- update(fm1BW.lme, weights = varPower())
fm2BW.lme

## Linear mixed-effects model fit by REML
##   Data: BodyWeight 
##   Log-restricted-likelihood: -570.96
##   Fixed: weight ~ Time * Diet 
## (Intercept)        Time       Diet2       Diet3  Time:Diet2 
##   251.60215     0.36109   200.77697   252.17017     0.60183 
##  Time:Diet3 
##     0.29524 
## 
## Random effects:
##  Formula: ~Time | Rat
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr  
## (Intercept) 36.89802 (Intr)
## Time         0.24365 -0.145
## Residual     0.17517       
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~fitted(.) 
##  Parameter estimates:
##   power 
## 0.54284 
## Number of Observations: 176
## Number of Groups: 16

anova(fm1BW.lme, fm2BW.lme)

##           Model df    AIC    BIC  logLik   Test L.Ratio p-value
## fm1BW.lme     1 10 1171.7 1203.1 -575.86                       
## fm2BW.lme     2 11 1163.9 1198.4 -570.96 1 vs 2  9.7983  0.0017

summary(fm2BW.lme)

## Linear mixed-effects model fit by REML
##  Data: BodyWeight 
##      AIC    BIC  logLik
##   1163.9 1198.4 -570.96
## 
## Random effects:
##  Formula: ~Time | Rat
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr  
## (Intercept) 36.89802 (Intr)
## Time         0.24365 -0.145
## Residual     0.17517       
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~fitted(.) 
##  Parameter estimates:
##   power 
## 0.54284 
## Fixed effects: weight ~ Time * Diet 
##               Value Std.Error  DF t-value p-value
## (Intercept) 251.602   13.0673 157 19.2543  0.0000
## Time          0.361    0.0884 157  4.0857  0.0001
## Diet2       200.777   22.6562  13  8.8619  0.0000
## Diet3       252.170   22.6616  13 11.1276  0.0000
## Time:Diet2    0.602    0.1554 157  3.8722  0.0002
## Time:Diet3    0.295    0.1559 157  1.8938  0.0601
##  Correlation: 
##            (Intr) Time   Diet2  Diet3  Tm:Dt2
## Time       -0.152                            
## Diet2      -0.577  0.088                     
## Diet3      -0.577  0.088  0.333              
## Time:Diet2  0.087 -0.569 -0.157 -0.050       
## Time:Diet3  0.086 -0.567 -0.050 -0.158  0.322
## 
## Standardized Within-Group Residuals:
##       Min        Q1       Med        Q3       Max 
## -2.937368 -0.443866  0.079864  0.580848  2.264898 
## 
## Number of Observations: 176
## Number of Groups: 16

anova(fm2BW.lme, L = c("Time:Diet2" = 1, "Time:Diet3" = -1))

## F-test for linear combination(s)
## Time:Diet2 Time:Diet3 
##          1         -1 
##   numDF denDF F-value p-value
## 1     1   157  2.8626  0.0926

cs1CompSymm <- corCompSymm(value = 0.3, form = ~ 1 | Subject)
cs2CompSymm <- corCompSymm(value = 0.3, form = ~ age | Subject)
cs1CompSymm <- Initialize(cs1CompSymm, data = Orthodont)
corMatrix(cs1CompSymm)

## $M01
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $M02
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $M03
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $M04
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $M05
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $M06
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $M07
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $M08
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $M09
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $M10
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $M11
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $M12
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $M13
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $M14
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $M15
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $M16
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $F01
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $F02
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $F03
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $F04
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $F05
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $F06
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $F07
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $F08
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $F09
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $F10
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0
## 
## $F11
##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.3  0.3  0.3
## [2,]  0.3  1.0  0.3  0.3
## [3,]  0.3  0.3  1.0  0.3
## [4,]  0.3  0.3  0.3  1.0

cs1Symm <- corSymm(value = c(0.2, 0.1, -0.1, 0, 0.2, 0),
                   form = ~ 1 | Subject)
cs1Symm <- Initialize(cs1Symm, data = Orthodont)
corMatrix(cs1Symm)

## $M01
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $M02
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $M03
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $M04
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $M05
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $M06
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $M07
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $M08
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $M09
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $M10
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $M11
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $M12
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $M13
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $M14
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $M15
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $M16
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $F01
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $F02
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $F03
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $F04
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $F05
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $F06
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $F07
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $F08
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $F09
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $F10
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00
## 
## $F11
##      [,1]       [,2]        [,3]        [,4]
## [1,]  1.0 2.0000e-01  1.0000e-01 -1.0000e-01
## [2,]  0.2 1.0000e+00  9.0206e-17  2.0000e-01
## [3,]  0.1 9.0206e-17  1.0000e+00 -1.0408e-16
## [4,] -0.1 2.0000e-01 -1.0408e-16  1.0000e+00

cs1AR1 <- corAR1(0.8, form = ~ 1 | Subject)
cs1AR1 <- Initialize(cs1AR1, data = Orthodont)
corMatrix(cs1AR1)

## $M01
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $M02
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $M03
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $M04
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $M05
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $M06
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $M07
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $M08
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $M09
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $M10
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $M11
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $M12
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $M13
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $M14
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $M15
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $M16
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $F01
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $F02
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $F03
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $F04
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $F05
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $F06
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $F07
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $F08
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $F09
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $F10
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000
## 
## $F11
##       [,1] [,2] [,3]  [,4]
## [1,] 1.000 0.80 0.64 0.512
## [2,] 0.800 1.00 0.80 0.640
## [3,] 0.640 0.80 1.00 0.800
## [4,] 0.512 0.64 0.80 1.000

cs1ARMA <- corARMA(0.4, form = ~ 1 | Subject, q = 1)
cs1ARMA <- Initialize(cs1ARMA, data = Orthodont)
corMatrix(cs1ARMA)

## $M01
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $M02
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $M03
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $M04
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $M05
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $M06
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $M07
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $M08
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $M09
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $M10
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $M11
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $M12
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $M13
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $M14
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $M15
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $M16
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $F01
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $F02
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $F03
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $F04
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $F05
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $F06
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $F07
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $F08
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $F09
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $F10
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000
## 
## $F11
##         [,1]    [,2]    [,3]    [,4]
## [1,] 1.00000 0.34483 0.00000 0.00000
## [2,] 0.34483 1.00000 0.34483 0.00000
## [3,] 0.00000 0.34483 1.00000 0.34483
## [4,] 0.00000 0.00000 0.34483 1.00000

cs2ARMA <- corARMA(c(0.8, 0.4), form = ~ 1 | Subject, p=1, q=1)
cs2ARMA <- Initialize(cs2ARMA, data = Orthodont)
corMatrix(cs2ARMA)

## $M01
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $M02
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $M03
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $M04
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $M05
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $M06
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $M07
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $M08
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $M09
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $M10
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $M11
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $M12
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $M13
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $M14
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $M15
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $M16
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $F01
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $F02
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $F03
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $F04
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $F05
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $F06
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $F07
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $F08
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $F09
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $F10
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000
## 
## $F11
##        [,1]  [,2]  [,3]   [,4]
## [1,] 1.0000 0.880 0.704 0.5632
## [2,] 0.8800 1.000 0.880 0.7040
## [3,] 0.7040 0.880 1.000 0.8800
## [4,] 0.5632 0.704 0.880 1.0000

spatDat <- data.frame(x = (0:4)/4, y = (0:4)/4)
cs1Exp <- corExp(1, form = ~ x + y)
cs1Exp <- Initialize(cs1Exp, spatDat)
corMatrix(cs1Exp)

##         [,1]    [,2]    [,3]    [,4]    [,5]
## [1,] 1.00000 0.70219 0.49307 0.34623 0.24312
## [2,] 0.70219 1.00000 0.70219 0.49307 0.34623
## [3,] 0.49307 0.70219 1.00000 0.70219 0.49307
## [4,] 0.34623 0.49307 0.70219 1.00000 0.70219
## [5,] 0.24312 0.34623 0.49307 0.70219 1.00000

cs2Exp <- corExp(1, form = ~ x + y, metric = "man")
cs2Exp <- Initialize(cs2Exp, spatDat)
corMatrix(cs2Exp)

##         [,1]    [,2]    [,3]    [,4]    [,5]
## [1,] 1.00000 0.60653 0.36788 0.22313 0.13534
## [2,] 0.60653 1.00000 0.60653 0.36788 0.22313
## [3,] 0.36788 0.60653 1.00000 0.60653 0.36788
## [4,] 0.22313 0.36788 0.60653 1.00000 0.60653
## [5,] 0.13534 0.22313 0.36788 0.60653 1.00000

cs3Exp <- corExp(c(1, 0.2), form = ~ x + y, nugget = TRUE)
cs3Exp <- Initialize(cs3Exp, spatDat)
corMatrix(cs3Exp)

##         [,1]    [,2]    [,3]    [,4]    [,5]
## [1,] 1.00000 0.56175 0.39445 0.27698 0.19449
## [2,] 0.56175 1.00000 0.56175 0.39445 0.27698
## [3,] 0.39445 0.56175 1.00000 0.56175 0.39445
## [4,] 0.27698 0.39445 0.56175 1.00000 0.56175
## [5,] 0.19449 0.27698 0.39445 0.56175 1.00000

fm1Ovar.lme <- lme(follicles ~ sin(2*pi*Time) + cos(2*pi*Time),
                   data = Ovary, random = pdDiag(~sin(2*pi*Time)))
fm1Ovar.lme

## Linear mixed-effects model fit by REML
##   Data: Ovary 
##   Log-restricted-likelihood: -813.04
##   Fixed: follicles ~ sin(2 * pi * Time) + cos(2 * pi * Time) 
##        (Intercept) sin(2 * pi * Time) cos(2 * pi * Time) 
##           12.18202           -3.29854           -0.86237 
## 
## Random effects:
##  Formula: ~sin(2 * pi * Time) | Mare
##  Structure: Diagonal
##         (Intercept) sin(2 * pi * Time) Residual
## StdDev:      3.0521             2.0793   3.1129
## 
## Number of Observations: 308
## Number of Groups: 11

ACF(fm1Ovar.lme)

##    lag       ACF
## 1    0  1.000000
## 2    1  0.379480
## 3    2  0.179722
## 4    3  0.035693
## 5    4  0.059779
## 6    5  0.002097
## 7    6  0.064327
## 8    7  0.071635
## 9    8  0.048578
## 10   9  0.027782
## 11  10 -0.034276
## 12  11 -0.077204
## 13  12 -0.161132
## 14  13 -0.196030
## 15  14 -0.289337

plot(ACF(fm1Ovar.lme,  maxLag = 10), alpha = 0.01)

fm2Ovar.lme <- update(fm1Ovar.lme, correlation = corAR1())
anova(fm1Ovar.lme, fm2Ovar.lme)

##             Model df    AIC    BIC  logLik   Test L.Ratio
## fm1Ovar.lme     1  6 1638.1 1660.4 -813.04               
## fm2Ovar.lme     2  7 1563.5 1589.5 -774.72 1 vs 2  76.634
##             p-value
## fm1Ovar.lme        
## fm2Ovar.lme  <.0001

if (interactive()) intervals(fm2Ovar.lme)
fm3Ovar.lme <- update(fm1Ovar.lme, correlation = corARMA(q = 2))
fm3Ovar.lme

## Linear mixed-effects model fit by REML
##   Data: Ovary 
##   Log-restricted-likelihood: -777.62
##   Fixed: follicles ~ sin(2 * pi * Time) + cos(2 * pi * Time) 
##        (Intercept) sin(2 * pi * Time) cos(2 * pi * Time) 
##           12.19419           -3.11453           -0.86911 
## 
## Random effects:
##  Formula: ~sin(2 * pi * Time) | Mare
##  Structure: Diagonal
##         (Intercept) sin(2 * pi * Time) Residual
## StdDev:      2.9741             1.6732   3.2442
## 
## Correlation Structure: ARMA(0,2)
##  Formula: ~1 | Mare 
##  Parameter estimate(s):
##  Theta1  Theta2 
## 0.47524 0.25701 
## Number of Observations: 308
## Number of Groups: 11

anova(fm2Ovar.lme, fm3Ovar.lme, test = F)

##             Model df    AIC    BIC  logLik
## fm2Ovar.lme     1  7 1563.5 1589.5 -774.72
## fm3Ovar.lme     2  8 1571.2 1601.0 -777.62

fm4Ovar.lme <- update(fm1Ovar.lme,
                       correlation = corCAR1(form = ~Time))
anova(fm2Ovar.lme, fm4Ovar.lme, test = F)

##             Model df    AIC    BIC  logLik
## fm2Ovar.lme     1  7 1563.5 1589.5 -774.72
## fm4Ovar.lme     2  7 1565.5 1591.6 -775.77

(fm5Ovar.lme <- update(fm1Ovar.lme,
                       corr = corARMA(p = 1, q = 1)))

## Linear mixed-effects model fit by REML
##   Data: Ovary 
##   Log-restricted-likelihood: -771.95
##   Fixed: follicles ~ sin(2 * pi * Time) + cos(2 * pi * Time) 
##        (Intercept) sin(2 * pi * Time) cos(2 * pi * Time) 
##           12.12487           -2.91983           -0.84871 
## 
## Random effects:
##  Formula: ~sin(2 * pi * Time) | Mare
##  Structure: Diagonal
##         (Intercept) sin(2 * pi * Time) Residual
## StdDev:      2.6144             1.0049   3.7334
## 
## Correlation Structure: ARMA(1,1)
##  Formula: ~1 | Mare 
##  Parameter estimate(s):
##     Phi1   Theta1 
##  0.78689 -0.27936 
## Number of Observations: 308
## Number of Groups: 11

anova(fm2Ovar.lme, fm5Ovar.lme)

##             Model df    AIC    BIC  logLik   Test L.Ratio
## fm2Ovar.lme     1  7 1563.5 1589.5 -774.72               
## fm5Ovar.lme     2  8 1559.9 1589.7 -771.95 1 vs 2  5.5538
##             p-value
## fm2Ovar.lme        
## fm5Ovar.lme  0.0184

plot(ACF(fm5Ovar.lme,  maxLag = 10, resType = "n"), alpha = 0.01)

Variogram(fm2BW.lme, form = ~ Time)

##     variog dist n.pairs
## 1  0.34506    1      16
## 2  0.99326    6      16
## 3  0.76201    7     144
## 4  0.68493    8      16
## 5  0.68184   13      16
## 6  0.95119   14     128
## 7  0.89961   15      16
## 8  1.69441   20      16
## 9  1.12511   21     112
## 10 1.08828   22      16
## 11 0.89696   28      96
## 12 0.93233   29      16
## 13 0.85146   35      80
## 14 0.75460   36      16
## 15 1.08217   42      64
## 16 1.56669   43      16
## 17 0.64376   49      48
## 18 0.67355   56      32
## 19 0.58680   63      16

plot(Variogram(fm2BW.lme, form = ~ Time, maxDist = 42))

fm3BW.lme <- update(fm2BW.lme,
                    correlation = corExp(form = ~ Time))
intervals(fm3BW.lme)

## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                  lower      est.     upper
## (Intercept) 225.637859 251.48677 277.33568
## Time          0.192914   0.36259   0.53226
## Diet2       151.729596 200.78555 249.84150
## Diet3       203.512665 252.58990 301.66713
## Time:Diet2    0.321775   0.62402   0.92627
## Time:Diet3    0.002626   0.30656   0.61050
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: Rat 
##                          lower     est.    upper
## sd((Intercept))       25.02845 36.91913 54.45891
## sd(Time)               0.14747  0.23281  0.36756
## cor((Intercept),Time) -0.62699 -0.14661  0.41455
## 
##  Correlation structure:
##        lower   est.  upper
## range 2.4616 4.8864 9.6996
## attr(,"label")
## [1] "Correlation structure:"
## 
##  Variance function:
##         lower    est.   upper
## power 0.24463 0.59434 0.94405
## attr(,"label")
## [1] "Variance function:"
## 
##  Within-group standard error:
##    lower     est.    upper 
## 0.018105 0.138387 1.057789

anova(fm2BW.lme, fm3BW.lme)

##           Model df    AIC    BIC  logLik   Test L.Ratio p-value
## fm2BW.lme     1 11 1163.9 1198.4 -570.96                       
## fm3BW.lme     2 12 1145.1 1182.8 -560.57 1 vs 2  20.781  <.0001

fm4BW.lme <-
      update(fm3BW.lme, correlation = corExp(form =  ~ Time,
                        nugget = TRUE))
anova(fm3BW.lme, fm4BW.lme)

##           Model df    AIC    BIC  logLik   Test L.Ratio p-value
## fm3BW.lme     1 12 1145.1 1182.8 -560.57                       
## fm4BW.lme     2 13 1137.6 1178.4 -555.82 1 vs 2  9.5025  0.0021

plot(Variogram(fm3BW.lme, form = ~ Time, maxDist = 42))

plot(Variogram(fm3BW.lme, form = ~ Time, maxDist = 42,
               resType = "n", robust = TRUE))

fm5BW.lme <- update(fm3BW.lme, correlation = corRatio(form = ~ Time))
fm6BW.lme <- update(fm3BW.lme, correlation = corSpher(form = ~ Time))
fm7BW.lme <- update(fm3BW.lme, correlation = corLin(form = ~ Time))
fm8BW.lme <- update(fm3BW.lme, correlation = corGaus(form = ~ Time))
anova(fm3BW.lme, fm5BW.lme, fm6BW.lme, fm7BW.lme, fm8BW.lme)

##           Model df    AIC    BIC  logLik
## fm3BW.lme     1 12 1145.1 1182.8 -560.57
## fm5BW.lme     2 12 1148.8 1186.4 -562.38
## fm6BW.lme     3 12 1150.8 1188.4 -563.39
## fm7BW.lme     4 12 1150.8 1188.4 -563.39
## fm8BW.lme     5 12 1150.8 1188.4 -563.39

fm1Orth.gls <- gls(distance ~ Sex * I(age - 11), Orthodont,
                   correlation = corSymm(form = ~ 1 | Subject),
                   weights = varIdent(form = ~ 1 | age))
fm1Orth.gls

## Generalized least squares fit by REML
##   Model: distance ~ Sex * I(age - 11) 
##   Data: Orthodont 
##   Log-restricted-likelihood: -212.27
## 
## Coefficients:
##           (Intercept)             SexFemale 
##              24.93712              -2.27174 
##           I(age - 11) SexFemale:I(age - 11) 
##               0.82680              -0.35044 
## 
## Correlation Structure: General
##  Formula: ~1 | Subject 
##  Parameter estimate(s):
##  Correlation: 
##   1     2     3    
## 2 0.568            
## 3 0.659 0.581      
## 4 0.522 0.725 0.740
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | age 
##  Parameter estimates:
##       8      10      12      14 
## 1.00000 0.87888 1.07446 0.95869 
## Degrees of freedom: 108 total; 104 residual
## Residual standard error: 2.3292

intervals(fm1Orth.gls)

## Approximate 95% confidence intervals
## 
##  Coefficients:
##                          lower     est.     upper
## (Intercept)           23.99941 24.93712 25.874835
## SexFemale             -3.74086 -2.27174 -0.802630
## I(age - 11)            0.66376  0.82680  0.989845
## SexFemale:I(age - 11) -0.60587 -0.35044 -0.095003
## attr(,"label")
## [1] "Coefficients:"
## 
##  Correlation structure:
##            lower    est.   upper
## cor(1,2) 0.25295 0.56820 0.77437
## cor(1,3) 0.38394 0.65895 0.82658
## cor(1,4) 0.18478 0.52204 0.74930
## cor(2,3) 0.27267 0.58061 0.78064
## cor(2,4) 0.48205 0.72492 0.86433
## cor(3,4) 0.51333 0.73962 0.86975
## attr(,"label")
## [1] "Correlation structure:"
## 
##  Variance function:
##      lower    est.  upper
## 10 0.63223 0.87888 1.2217
## 12 0.79935 1.07446 1.4443
## 14 0.68137 0.95869 1.3489
## attr(,"label")
## [1] "Variance function:"
## 
##  Residual standard error:
##  lower   est.  upper 
## 1.7597 2.3292 3.0831

fm2Orth.gls <-
   update(fm1Orth.gls, corr = corCompSymm(form = ~ 1 | Subject))
anova(fm1Orth.gls, fm2Orth.gls)

##             Model df    AIC    BIC  logLik   Test L.Ratio
## fm1Orth.gls     1 14 452.55 489.57 -212.27               
## fm2Orth.gls     2  9 449.97 473.77 -215.99 1 vs 2  7.4256
##             p-value
## fm1Orth.gls        
## fm2Orth.gls  0.1909

intervals(fm2Orth.gls)

## Approximate 95% confidence intervals
## 
##  Coefficients:
##                          lower     est.     upper
## (Intercept)           23.92956 24.86832 25.807077
## SexFemale             -3.66815 -2.19740 -0.726641
## I(age - 11)            0.64160  0.79431  0.947029
## SexFemale:I(age - 11) -0.55482 -0.31556 -0.076301
## attr(,"label")
## [1] "Coefficients:"
## 
##  Correlation structure:
##       lower    est.   upper
## Rho 0.44608 0.63529 0.77825
## attr(,"label")
## [1] "Correlation structure:"
## 
##  Variance function:
##      lower    est.  upper
## 10 0.63545 0.86236 1.1703
## 12 0.76905 1.03395 1.3901
## 14 0.67823 0.92037 1.2490
## attr(,"label")
## [1] "Variance function:"
## 
##  Residual standard error:
##  lower   est.  upper 
## 1.8023 2.3868 3.1609

fm3Orth.gls <- update(fm2Orth.gls, weights = NULL)
anova(fm2Orth.gls, fm3Orth.gls)

##             Model df    AIC    BIC  logLik   Test L.Ratio
## fm2Orth.gls     1  9 449.97 473.77 -215.99               
## fm3Orth.gls     2  6 445.76 461.62 -216.88 1 vs 2  1.7849
##             p-value
## fm2Orth.gls        
## fm3Orth.gls  0.6182

plot(fm3Orth.gls, resid(., type = "n") ~ age | Sex)

fm4Orth.gls <- update(fm3Orth.gls,
                      weights = varIdent(form = ~ 1 | Sex))
anova(fm3Orth.gls, fm4Orth.gls)

##             Model df    AIC    BIC  logLik   Test L.Ratio
## fm3Orth.gls     1  6 445.76 461.62 -216.88               
## fm4Orth.gls     2  7 436.19 454.70 -211.09 1 vs 2  11.569
##             p-value
## fm3Orth.gls        
## fm4Orth.gls   7e-04

qqnorm(fm4Orth.gls, ~resid(., type = "n"))

# not in book but needed for the following command
fm3Orth.lme <-
    lme(distance~Sex*I(age-11), data = Orthodont,
        random = ~ I(age-11) | Subject,
        weights = varIdent(form = ~ 1 | Sex))
anova(fm3Orth.lme, fm4Orth.gls, test = FALSE)

##             Model df    AIC    BIC  logLik
## fm3Orth.lme     1  9 429.52 453.32 -205.76
## fm4Orth.gls     2  7 436.19 454.70 -211.09

fm1Dial.gls <-
  gls(rate ~(pressure + I(pressure^2) + I(pressure^3) + I(pressure^4))*QB,
      Dialyzer)
plot(fm1Dial.gls, resid(.) ~ pressure, abline = 0)

fm2Dial.gls <- update(fm1Dial.gls,
                      weights = varPower(form = ~ pressure))
anova(fm1Dial.gls, fm2Dial.gls)

##             Model df    AIC    BIC  logLik   Test L.Ratio
## fm1Dial.gls     1 11 761.17 792.71 -369.59               
## fm2Dial.gls     2 12 738.22 772.63 -357.11 1 vs 2   24.95
##             p-value
## fm1Dial.gls        
## fm2Dial.gls  <.0001

ACF(fm2Dial.gls, form = ~ 1 | Subject)

##   lag      ACF
## 1   0 1.000000
## 2   1 0.770851
## 3   2 0.632301
## 4   3 0.408306
## 5   4 0.200737
## 6   5 0.073116
## 7   6 0.077802

plot(ACF(fm2Dial.gls, form = ~ 1 | Subject), alpha = 0.01)

(fm3Dial.gls <- update(fm2Dial.gls,
                      corr = corAR1(0.771, form = ~ 1 | Subject)))

## Generalized least squares fit by REML
##   Model: rate ~ (pressure + I(pressure^2) + I(pressure^3) + I(pressure^4)) *      QB 
##   Data: Dialyzer 
##   Log-restricted-likelihood: -308.34
## 
## Coefficients:
##         (Intercept)            pressure       I(pressure^2) 
##           -16.81845            92.33424           -49.26516 
##       I(pressure^3)       I(pressure^4)               QB300 
##            11.39968            -1.01964            -1.59419 
##      pressure:QB300 I(pressure^2):QB300 I(pressure^3):QB300 
##             1.70543             2.12680             0.47971 
## I(pressure^4):QB300 
##            -0.22064 
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | Subject 
##  Parameter estimate(s):
##    Phi 
## 0.7526 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~pressure 
##  Parameter estimates:
##   power 
## 0.51824 
## Degrees of freedom: 140 total; 130 residual
## Residual standard error: 3.0463

intervals(fm3Dial.gls)

## Approximate 95% confidence intervals
## 
##  Coefficients:
##                        lower      est.      upper
## (Intercept)         -18.8968 -16.81845 -14.740093
## pressure             81.9144  92.33424 102.754095
## I(pressure^2)       -63.1040 -49.26516 -35.426272
## I(pressure^3)         4.5648  11.39968  18.234542
## I(pressure^4)        -2.1248  -1.01964   0.085558
## QB300                -4.7565  -1.59419   1.568148
## pressure:QB300      -13.6410   1.70543  17.051845
## I(pressure^2):QB300 -17.9484   2.12680  22.201987
## I(pressure^3):QB300  -9.3503   0.47971  10.309706
## I(pressure^4):QB300  -1.8021  -0.22064   1.360833
## attr(,"label")
## [1] "Coefficients:"
## 
##  Correlation structure:
##       lower   est.   upper
## Phi 0.62844 0.7526 0.83937
## attr(,"label")
## [1] "Correlation structure:"
## 
##  Variance function:
##        lower    est.   upper
## power 0.3806 0.51824 0.65587
## attr(,"label")
## [1] "Variance function:"
## 
##  Residual standard error:
##  lower   est.  upper 
## 2.5013 3.0463 3.7101

anova(fm2Dial.gls, fm3Dial.gls)

##             Model df    AIC    BIC  logLik   Test L.Ratio
## fm2Dial.gls     1 12 738.22 772.63 -357.11               
## fm3Dial.gls     2 13 642.67 679.95 -308.34 1 vs 2  97.546
##             p-value
## fm2Dial.gls        
## fm3Dial.gls  <.0001

anova(fm3Dial.gls, fm2Dial.lme, test = FALSE)

##             Model df    AIC    BIC  logLik
## fm3Dial.gls     1 13 642.67 679.95 -308.34
## fm2Dial.lme     2 18 655.01 706.63 -309.51

fm1Wheat2 <- gls(yield ~ variety - 1, Wheat2)
Variogram(fm1Wheat2, form = ~ latitude + longitude)

##     variog    dist n.pairs
## 1  0.37026  4.3000    1143
## 2  0.39572  5.6080    1259
## 3  0.47019  8.3863    1263
## 4  0.50774  9.3231    1241
## 5  0.54502 10.5190    1242
## 6  0.63965 12.7472    1241
## 7  0.61183 13.3929    1283
## 8  0.65711 14.7635    1252
## 9  0.73844 16.1818    1221
## 10 0.72832 17.3666    1261
## 11 0.75117 18.4567    1288
## 12 0.87492 20.2428    1254
## 13 0.80516 21.6335    1256
## 14 0.87118 22.6736    1182
## 15 0.86812 24.6221    1257
## 16 0.85929 26.2427    1264
## 17 0.97088 28.5601    1235
## 18 0.99293 30.7877    1226
## 19 1.09617 34.5879    1263
## 20 1.34146 39.3641    1234

plot(Variogram(fm1Wheat2, form = ~ latitude + longitude,
      maxDist = 32), xlim = c(0,32))

fm2Wheat2 <- update(fm1Wheat2, corr = corSpher(c(28, 0.2),
                               form = ~ latitude + longitude,
                               nugget = TRUE))
fm2Wheat2

## Generalized least squares fit by REML
##   Model: yield ~ variety - 1 
##   Data: Wheat2 
##   Log-restricted-likelihood: -533.93
## 
## Coefficients:
##   varietyARAPAHOE      varietyBRULE   varietyBUCKSKIN 
##            26.659            25.850            34.848 
##    varietyCENTURA  varietyCENTURK78   varietyCHEYENNE 
##            25.095            26.334            24.673 
##       varietyCODY       varietyCOLT       varietyGAGE 
##            22.489            25.184            24.321 
##  varietyHOMESTEAD   varietyKS831374     varietyLANCER 
##            21.669            26.856            23.333 
##    varietyLANCOTA    varietyNE83404    varietyNE83406 
##            21.291            24.049            25.314 
##    varietyNE83407    varietyNE83432    varietyNE83498 
##            25.192            21.818            28.705 
##    varietyNE83T12    varietyNE84557    varietyNE85556 
##            22.145            21.781            27.973 
##    varietyNE85623    varietyNE86482    varietyNE86501 
##            23.924            24.974            25.018 
##    varietyNE86503    varietyNE86507    varietyNE86509 
##            27.224            27.453            22.408 
##    varietyNE86527    varietyNE86582    varietyNE86606 
##            25.882            22.587            26.796 
##    varietyNE86607   varietyNE86T666    varietyNE87403 
##            25.910            16.787            21.490 
##    varietyNE87408    varietyNE87409    varietyNE87446 
##            24.319            26.320            22.219 
##    varietyNE87451    varietyNE87457    varietyNE87463 
##            24.191            23.536            23.226 
##    varietyNE87499    varietyNE87512    varietyNE87513 
##            22.173            22.575            21.825 
##    varietyNE87522    varietyNE87612    varietyNE87613 
##            19.483            27.391            27.562 
##    varietyNE87615    varietyNE87619    varietyNE87627 
##            23.819            28.485            18.480 
##     varietyNORKAN    varietyREDLAND varietyROUGHRIDER 
##            22.145            27.983            25.748 
##    varietySCOUT66  varietySIOUXLAND     varietyTAM107 
##            26.878            25.656            22.770 
##     varietyTAM200       varietyVONA 
##            18.764            24.782 
## 
## Correlation Structure: Spherical spatial correlation
##  Formula: ~latitude + longitude 
##  Parameter estimate(s):
##    range   nugget 
## 27.45748  0.20931 
## Degrees of freedom: 224 total; 168 residual
## Residual standard error: 7.4106

fm3Wheat2 <- update(fm1Wheat2,
                    corr = corRatio(c(12.5, 0.2),
                    form = ~ latitude + longitude, nugget = TRUE))
fm3Wheat2

## Generalized least squares fit by REML
##   Model: yield ~ variety - 1 
##   Data: Wheat2 
##   Log-restricted-likelihood: -532.64
## 
## Coefficients:
##   varietyARAPAHOE      varietyBRULE   varietyBUCKSKIN 
##            26.546            26.284            35.037 
##    varietyCENTURA  varietyCENTURK78   varietyCHEYENNE 
##            24.867            26.695            24.430 
##       varietyCODY       varietyCOLT       varietyGAGE 
##            23.366            25.176            24.459 
##  varietyHOMESTEAD   varietyKS831374     varietyLANCER 
##            21.517            26.549            23.046 
##    varietyLANCOTA    varietyNE83404    varietyNE83406 
##            21.177            24.592            25.734 
##    varietyNE83407    varietyNE83432    varietyNE83498 
##            25.527            21.817            29.084 
##    varietyNE83T12    varietyNE84557    varietyNE85556 
##            21.576            21.321            27.917 
##    varietyNE85623    varietyNE86482    varietyNE86501 
##            23.739            24.373            24.879 
##    varietyNE86503    varietyNE86507    varietyNE86509 
##            27.300            27.366            22.174 
##    varietyNE86527    varietyNE86582    varietyNE86606 
##            24.989            23.274            27.317 
##    varietyNE86607   varietyNE86T666    varietyNE87403 
##            25.710            17.305            21.837 
##    varietyNE87408    varietyNE87409    varietyNE87446 
##            24.657            26.314            22.113 
##    varietyNE87451    varietyNE87457    varietyNE87463 
##            24.383            23.624            23.435 
##    varietyNE87499    varietyNE87512    varietyNE87513 
##            21.900            22.743            21.644 
##    varietyNE87522    varietyNE87612    varietyNE87613 
##            19.632            28.330            27.732 
##    varietyNE87615    varietyNE87619    varietyNE87627 
##            24.013            28.661            19.099 
##     varietyNORKAN    varietyREDLAND varietyROUGHRIDER 
##            22.684            27.670            25.619 
##    varietySCOUT66  varietySIOUXLAND     varietyTAM107 
##            26.289            25.740            22.476 
##     varietyTAM200       varietyVONA 
##            18.693            25.046 
## 
## Correlation Structure: Rational quadratic spatial correlation
##  Formula: ~latitude + longitude 
##  Parameter estimate(s):
##    range   nugget 
## 13.46134  0.19358 
## Degrees of freedom: 224 total; 168 residual
## Residual standard error: 8.8468

anova(fm2Wheat2, fm3Wheat2)

##           Model df    AIC    BIC  logLik
## fm2Wheat2     1 59 1185.9 1370.2 -533.93
## fm3Wheat2     2 59 1183.3 1367.6 -532.64

anova(fm1Wheat2, fm3Wheat2)

##           Model df    AIC    BIC  logLik   Test L.Ratio p-value
## fm1Wheat2     1 57 1354.7 1532.8 -620.37                       
## fm3Wheat2     2 59 1183.3 1367.6 -532.64 1 vs 2  175.46  <.0001

plot(Variogram(fm3Wheat2, resType = "n"))

plot(fm3Wheat2, resid(., type = "n") ~ fitted(.), abline = 0)

qqnorm(fm3Wheat2, ~ resid(., type = "n"))

fm4Wheat2 <- update(fm3Wheat2, model = yield ~ variety)
anova(fm4Wheat2)

## Denom. DF: 168 
##             numDF F-value p-value
## (Intercept)     1 30.3994  <.0001
## variety        55  1.8509  0.0015

anova(fm3Wheat2, L = c(-1, 0, 1))

## Denom. DF: 168 
##  F-test for linear combination(s)
## varietyARAPAHOE varietyBUCKSKIN 
##              -1               1 
##   numDF F-value p-value
## 1     1  7.6967  0.0062

# cleanup

proc.time()

##    user  system elapsed 
##   55.47    2.37   62.01