R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("SASmixed-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('SASmixed')
Loading required package: lme4
Loading required package: Matrix
Loading required package: lattice
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "Animal"
> 
> ### * Animal
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Animal
> ### Title: Animal breeding experiment
> ### Aliases: Animal
> ### Keywords: datasets
> 
> ### ** Examples
> 
> str(Animal)
`data.frame':	20 obs. of  3 variables:
 $ Sire        : Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ...
 $ Dam         : Factor w/ 2 levels "1","2": 1 1 2 2 1 1 2 2 1 1 ...
 $ AvgDailyGain: num  2.24 1.85 2.05 2.41 1.99 1.93 2.72 2.32 2.33 2.68 ...
 - attr(*, "ginfo")=List of 7
  ..$ formula     :Class 'formula' length 3 AvgDailyGain ~ 1 | Sire/Dam
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ order.groups:List of 2
  .. ..$ Sire: logi TRUE
  .. ..$ Dam : logi TRUE
  ..$ FUN         :function (x)  
  ..$ outer       : NULL
  ..$ inner       : NULL
  ..$ labels      :List of 1
  .. ..$ AvgDailyGain: chr "Average Daily Weight Gain"
  ..$ units       : list()
> 
> 
> 
> cleanEx(); ..nameEx <- "AvgDailyGain"
> 
> ### * AvgDailyGain
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: AvgDailyGain
> ### Title: Average daily weight gain of steers on different diets
> ### Aliases: AvgDailyGain
> ### Keywords: datasets
> 
> ### ** Examples
> 
> options(contrasts = c(unordered = "contr.SAS", ordered = "contr.poly"))
> ## plot of adg versus Treatment by Block
> if (require("latticeExtra", quietly = TRUE)) gplot(AvgDailyGain)
> fm1Adg <- lmer(adg ~ InitWt * Treatment - 1 + (1 | Block), AvgDailyGain)
> summary(fm1Adg) # compare with output 5.1, p. 178
Linear mixed-effects model fit by REML
Formula: adg ~ InitWt * Treatment - 1 + (1 | Block) 
   Data: AvgDailyGain 
      AIC     BIC    logLik MLdeviance REMLdeviance
 85.32685 99.9842 -32.66342   10.09817     65.32685
Random effects:
 Groups   Name        Variance Std.Dev.
 Block    (Intercept) 0.25930  0.50922 
 Residual             0.04943  0.22233 
# of obs: 32, groups: Block, 8

Fixed effects:
                     Estimate Std. Error DF t value Pr(>|t|)  
InitWt              0.0044480  0.0020816 24  2.1368  0.04301 *
Treatment0          0.4391279  0.7110925 24  0.6175  0.54269  
Treatment10         1.4261132  0.6375493 24  2.2369  0.03485 *
Treatment20         0.4796212  0.5488892 24  0.8738  0.39089  
Treatment30         0.2001150  0.7752039 24  0.2581  0.79850  
InitWt:Treatment0  -0.0021543  0.0027863 24 -0.7732  0.44697  
InitWt:Treatment10 -0.0033651  0.0025148 24 -1.3381  0.19341  
InitWt:Treatment20 -0.0010823  0.0024876 24 -0.4351  0.66739  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Correlation of Fixed Effects:
            InitWt Trtmn0 Trtm10 Trtm20 Trtm30 InW:T0 IW:T10
Treatment0   0.050                                          
Treatment10 -0.032  0.039                                   
Treatment20  0.035  0.080  0.334                            
Treatment30 -0.967  0.011  0.097  0.043                     
IntWt:Trtm0 -0.780 -0.640  0.046 -0.024  0.754              
IntWt:Trt10 -0.808 -0.021 -0.535 -0.178  0.781  0.617       
IntWt:Trt20 -0.856 -0.040 -0.106 -0.512  0.828  0.666  0.775
> anova(fm1Adg)   # checking significance of terms
Analysis of Variance Table
                 Df  Sum Sq Mean Sq   Denom F value    Pr(>F)    
InitWt            1  4.5320  4.5320 24.0000 91.6859 1.139e-09 ***
Treatment         4  1.7425  0.4356 24.0000  8.8130 0.0001587 ***
InitWt:Treatment  3  0.1381  0.0460 24.0000  0.9312 0.4408869    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> fm2Adg <- lmer(adg ~ InitWt + Treatment + (1 | Block), AvgDailyGain)
> summary(fm2Adg)
Linear mixed-effects model fit by REML
Formula: adg ~ InitWt + Treatment + (1 | Block) 
   Data: AvgDailyGain 
      AIC      BIC    logLik MLdeviance REMLdeviance
 50.33733 60.59748 -18.16866   13.62304     36.33733
Random effects:
 Groups   Name        Variance Std.Dev.
 Block    (Intercept) 0.24084  0.49076 
 Residual             0.05008  0.22379 
# of obs: 32, groups: Block, 8

Fixed effects:
               Estimate  Std. Error DF t value  Pr(>|t|)    
(Intercept)  0.80110753  0.35566101 27  2.2524  0.032628 *  
InitWt       0.00277972  0.00083335 27  3.3356  0.002486 ** 
Treatment0  -0.55207371  0.11481324 27 -4.8084 5.097e-05 ***
Treatment10 -0.06856620  0.11896910 27 -0.5763  0.569162    
Treatment20 -0.08812918  0.11628794 27 -0.7579  0.455103    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Correlation of Fixed Effects:
            (Intr) InitWt Trtmn0 Trtm10
InitWt      -0.844                     
Treatment0   0.036 -0.224              
Treatment10  0.139 -0.340  0.534       
Treatment20  0.079 -0.272  0.530  0.545
> anova(fm2Adg)
Analysis of Variance Table
          Df  Sum Sq Mean Sq   Denom F value    Pr(>F)    
InitWt     1  0.5146  0.5146 27.0000  10.275 0.0034525 ** 
Treatment  3  1.5267  0.5089 27.0000  10.162 0.0001185 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> summary(lmer(adg ~ InitWt + Treatment - 1 + (1 | Block), AvgDailyGain))
Linear mixed-effects model fit by REML
Formula: adg ~ InitWt + Treatment - 1 + (1 | Block) 
   Data: AvgDailyGain 
      AIC      BIC    logLik MLdeviance REMLdeviance
 50.33733 60.59748 -18.16866   13.62304     36.33733
Random effects:
 Groups   Name        Variance Std.Dev.
 Block    (Intercept) 0.24084  0.49076 
 Residual             0.05008  0.22379 
# of obs: 32, groups: Block, 8

Fixed effects:
              Estimate Std. Error DF t value Pr(>|t|)   
InitWt      2.7797e-03 8.3335e-04 27  3.3356 0.002486 **
Treatment0  2.4903e-01 3.7763e-01 27  0.6595 0.515185   
Treatment10 7.3254e-01 3.9038e-01 27  1.8765 0.071437 . 
Treatment20 7.1298e-01 3.8277e-01 27  1.8627 0.073421 . 
Treatment30 8.0111e-01 3.5566e-01 27  2.2524 0.032628 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Correlation of Fixed Effects:
            InitWt Trtmn0 Trtm10 Trtm20
Treatment0  -0.863                     
Treatment10 -0.873  0.957              
Treatment20 -0.867  0.957  0.958       
Treatment30 -0.844  0.953  0.953  0.953
> 
> 
> 
> options(contrasts = c(unordered = "contr.treatment",ordered = "contr.poly"))
> cleanEx(); ..nameEx <- "BIB"
> 
> ### * BIB
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: BIB
> ### Title: Data from a balanced incomplete block design
> ### Aliases: BIB
> ### Keywords: datasets
> 
> ### ** Examples
> 
> options(contrasts = c(unordered = "contr.SAS", ordered = "contr.poly"))
> if (require("latticeExtra", quietly = TRUE)) gplot(BIB)
> fm1BIB <- lmer(y ~ Treatment * x + (1 | Block), BIB)
> summary(fm1BIB)   # compare with Output 5.7, p. 188
Linear mixed-effects model fit by REML
Formula: y ~ Treatment * x + (1 | Block) 
   Data: BIB 
      AIC     BIC    logLik MLdeviance REMLdeviance
 124.8945 136.675 -52.44723    93.4988     104.8945
Random effects:
 Groups   Name        Variance Std.Dev.
 Block    (Intercept) 18.238   4.2706  
 Residual              1.201   1.0959  
# of obs: 24, groups: Block, 8

Fixed effects:
              Estimate Std. Error DF t value  Pr(>|t|)    
(Intercept)  22.368136   3.102127 16  7.2106 2.077e-06 ***
Treatment1    4.429357   3.365812 16  1.3160  0.206723    
Treatment2   -0.437474   2.933883 16 -0.1491  0.883329    
Treatment3    6.278356   3.282780 16  1.9125  0.073883 .  
x             0.442537   0.087082 16  5.0818  0.000111 ***
Treatment1:x -0.223761   0.106107 16 -2.1088  0.051072 .  
Treatment2:x  0.053385   0.097165 16  0.5494  0.590298    
Treatment3:x -0.179166   0.115736 16 -1.5481  0.141158    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Correlation of Fixed Effects:
            (Intr) Trtmn1 Trtmn2 Trtmn3 x      Trtm1: Trtm2:
Treatment1  -0.729                                          
Treatment2  -0.778  0.797                                   
Treatment3  -0.797  0.827  0.826                            
x           -0.859  0.797  0.865  0.886                     
Treatmnt1:x  0.709 -0.979 -0.774 -0.797 -0.799              
Treatmnt2:x  0.722 -0.731 -0.965 -0.763 -0.829  0.729       
Treatmnt3:x  0.769 -0.789 -0.790 -0.976 -0.879  0.777  0.748
> anova(fm1BIB)     # strong evidence of different slopes
Analysis of Variance Table
            Df  Sum Sq Mean Sq   Denom  F value    Pr(>F)    
Treatment    3  23.448   7.816  16.000   6.5078  0.004376 ** 
x            1 136.808 136.808  16.000 113.9129 1.102e-08 ***
Treatment:x  3  18.426   6.142  16.000   5.1142  0.011368 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> fm2BIB <- lmer(y ~ Treatment + x : Grp + (1 | Block), BIB)
> summary(fm2BIB)   # compare with Output 5.9, p. 193
Linear mixed-effects model fit by REML
Formula: y ~ Treatment + x:Grp + (1 | Block) 
   Data: BIB 
      AIC      BIC    logLik MLdeviance REMLdeviance
 115.1770 124.6015 -49.58851   94.08929     99.17702
Random effects:
 Groups   Name        Variance Std.Dev.
 Block    (Intercept) 18.5214  4.3036  
 Residual              1.0380  1.0188  
# of obs: 24, groups: Block, 8

Fixed effects:
             Estimate Std. Error DF t value  Pr(>|t|)    
(Intercept) 20.945232   2.062233 18 10.1566 7.028e-09 ***
Treatment1   5.341392   1.975836 18  2.7034  0.014548 *  
Treatment2   1.135550   0.714037 18  1.5903  0.129171    
Treatment3   8.180984   1.770218 18  4.6215  0.000212 ***
x:Grp13      0.239519   0.042966 18  5.5746 2.724e-05 ***
x:Grp24      0.489228   0.044125 18 11.0874 1.783e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Correlation of Fixed Effects:
           (Intr) Trtmn1 Trtmn2 Trtmn3 x:Gr13
Treatment1 -0.501                            
Treatment2 -0.431  0.559                     
Treatment3 -0.527  0.942  0.581              
x:Grp13     0.027 -0.663 -0.165 -0.605       
x:Grp24    -0.639  0.651  0.452  0.688  0.042
> anova(fm2BIB)
Analysis of Variance Table
          Df  Sum Sq Mean Sq   Denom F value    Pr(>F)    
Treatment  3  23.424   7.808  18.000  7.5225  0.001820 ** 
x:Grp      2 154.733  77.366  18.000 74.5363 1.956e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> 
> 
> 
> options(contrasts = c(unordered = "contr.treatment",ordered = "contr.poly"))
> cleanEx(); ..nameEx <- "Bond"
> 
> ### * Bond
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Bond
> ### Title: Strengths of metal bonds
> ### Aliases: Bond
> ### Keywords: datasets
> 
> ### ** Examples
> 
> options(contrasts = c(unordered = "contr.SAS", ordered = "contr.poly"))
> fm1Bond <- lmer(pressure ~ Metal + (1|Ingot), Bond)
> summary(fm1Bond)               # compare with output 1.1 on p. 6
Linear mixed-effects model fit by REML
Formula: pressure ~ Metal + (1 | Ingot) 
   Data: Bond 
      AIC      BIC   logLik MLdeviance REMLdeviance
 117.7902 123.0128 -53.8951   115.7074     107.7902
Random effects:
 Groups   Name        Variance Std.Dev.
 Ingot    (Intercept) 11.448   3.3835  
 Residual             10.372   3.2205  
# of obs: 21, groups: Ingot, 7

Fixed effects:
            Estimate Std. Error DF t value Pr(>|t|)    
(Intercept) 71.10000    1.76552 18 40.2715  < 2e-16 ***
Metalc      -0.91429    1.72143 18 -0.5311  0.60183    
Metali       4.80000    1.72143 18  2.7884  0.01213 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Correlation of Fixed Effects:
       (Intr) Metalc
Metalc -0.488       
Metali -0.488  0.500
> anova(fm1Bond)
Analysis of Variance Table
      Df Sum Sq Mean Sq  Denom F value   Pr(>F)   
Metal  2 131.90   65.95  18.00  6.3588 0.008147 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> 
> 
> 
> options(contrasts = c(unordered = "contr.treatment",ordered = "contr.poly"))
> cleanEx(); ..nameEx <- "Cultivation"
> 
> ### * Cultivation
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Cultivation
> ### Title: Bacterial innoculation applied to grass cultivars
> ### Aliases: Cultivation
> ### Keywords: datasets
> 
> ### ** Examples
> 
> options(contrasts = c(unordered = "contr.SAS", ordered = "contr.poly"))
> str(Cultivation)
`data.frame':	24 obs. of  4 variables:
 $ Block: Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 2 2 2 2 ...
 $ Cult : Factor w/ 2 levels "a","b": 1 1 1 2 2 2 1 1 1 2 ...
 $ Inoc : Factor w/ 3 levels "con","dea","liv": 1 2 3 1 2 3 1 2 3 1 ...
 $ drywt: num  27.4 29.7 34.5 29.4 32.5 34.4 28.9 28.7 33.4 28.7 ...
 - attr(*, "ginfo")=List of 7
  ..$ formula     :Class 'formula' length 3 drywt ~ 1 | Block/Cult
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ order.groups:List of 2
  .. ..$ Block: logi TRUE
  .. ..$ Cult : logi TRUE
  ..$ FUN         :function (x)  
  ..$ outer       : NULL
  ..$ inner       :List of 1
  .. ..$ Cult:Class 'formula' length 2 ~Inoc
  .. .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ labels      :List of 1
  .. ..$ drywt: chr "Yield"
  ..$ units       : list()
> xtabs(~Block+Cult, Cultivation)
     Cult
Block a b
    1 3 3
    2 3 3
    3 3 3
    4 3 3
> fm1Cult <- lmer(drywt ~ Inoc * Cult + (1|Block) + (1|Cult), Cultivation)
> summary(fm1Cult)   # compare with Output 2.10, page 58
Linear mixed-effects model fit by REML
Formula: drywt ~ Inoc * Cult + (1 | Block) + (1 | Cult) 
   Data: Cultivation 
      AIC     BIC    logLik MLdeviance REMLdeviance
 86.48742 97.0899 -34.24371   74.94172     68.48742
Random effects:
 Groups   Name        Variance Std.Dev.
 Block    (Intercept) 1.20737  1.0988  
 Cult     (Intercept) 0.26585  0.5156  
 Residual             1.19632  1.0938  
# of obs: 24, groups: Block, 4; Cult, 2

Fixed effects:
              Estimate Std. Error DF t value  Pr(>|t|)    
(Intercept)   33.52500    0.93100 18 36.0095 < 2.2e-16 ***
Inoccon       -5.50000    0.77341 18 -7.1114 1.256e-06 ***
Inocdea       -2.87500    0.77341 18 -3.7173  0.001577 ** 
Culta         -0.37500    1.06295 18 -0.3528  0.728341    
Inoccon:Culta  0.25000    1.09376 18  0.2286  0.821781    
Inocdea:Culta -1.02500    1.09376 18 -0.9371  0.361096    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Correlation of Fixed Effects:
            (Intr) Inoccn Inocde Culta  Incc:C
Inoccon     -0.415                            
Inocdea     -0.415  0.500                     
Culta       -0.571  0.364  0.364              
Inoccon:Clt  0.294 -0.707 -0.354 -0.514       
Inocdea:Clt  0.294 -0.354 -0.707 -0.514  0.500
> anova(fm1Cult)
Analysis of Variance Table
          Df  Sum Sq Mean Sq   Denom F value   Pr(>F)    
Inoc       2 118.176  59.088  18.000 49.3915 4.91e-08 ***
Cult       1   0.656   0.656  18.000  0.5487   0.4684    
Inoc:Cult  2   1.826   0.913  18.000  0.7631   0.4807    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> fm2Cult <- lmer(drywt ~ Inoc + Cult + (1|Block) + (1|Cult), Cultivation)
> anova(fm2Cult)
Analysis of Variance Table
     Df  Sum Sq Mean Sq   Denom F value    Pr(>F)    
Inoc  2 118.176  59.088  20.000 50.8073 1.447e-08 ***
Cult  1   0.656   0.656  20.000  0.5644    0.4612    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> fm3Cult <- lmer(drywt ~ Inoc + (1|Block) + (1|Cult), Cultivation)
> anova(fm3Cult)
Analysis of Variance Table
     Df  Sum Sq Mean Sq   Denom F value    Pr(>F)    
Inoc  2 118.176  59.088  21.000  50.807 8.988e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> summary(fm3Cult)
Linear mixed-effects model fit by REML
Formula: drywt ~ Inoc + (1 | Block) + (1 | Cult) 
   Data: Cultivation 
      AIC      BIC    logLik MLdeviance REMLdeviance
 87.67784 94.74616 -37.83892   77.32082     75.67784
Random effects:
 Groups   Name        Variance Std.Dev.
 Block    (Intercept) 1.21283  1.10129 
 Cult     (Intercept) 0.10364  0.32193 
 Residual             1.16299  1.07842 
# of obs: 24, groups: Block, 4; Cult, 2

Fixed effects:
            Estimate Std. Error DF t value  Pr(>|t|)    
(Intercept) 33.33750    0.70739 21 47.1274 < 2.2e-16 ***
Inoccon     -5.37500    0.53921 21 -9.9683 2.048e-09 ***
Inocdea     -3.38750    0.53921 21 -6.2823 3.134e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Correlation of Fixed Effects:
        (Intr) Inoccn
Inoccon -0.381       
Inocdea -0.381  0.500
> 
> 
> 
> options(contrasts = c(unordered = "contr.treatment",ordered = "contr.poly"))
> cleanEx(); ..nameEx <- "Demand"
> 
> ### * Demand
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Demand
> ### Title: Per-capita demand deposits by state and year
> ### Aliases: Demand
> ### Keywords: datasets
> 
> ### ** Examples
> 
> str(Demand)
`data.frame':	77 obs. of  7 variables:
 $ State: Factor w/ 7 levels "CA","DC","FL",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Year : num  1949 1950 1951 1952 1953 ...
 $ d    : num  533 603 669 651 609 634 665 676 642 678 ...
 $ y    : num  1347 1464 1608 1636 1669 ...
 $ rd   : num  0.343 0.364 0.367 0.369 0.41 0.499 0.496 0.533 0.63 0.667 ...
 $ rt   : num  1.11 1.16 1.49 1.57 1.59 ...
 $ rs   : num  2.90 2.94 3.09 3.07 3.36 ...
 - attr(*, "ginfo")=List of 7
  ..$ formula     :Class 'formula' length 3 d ~ Year | State
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ order.groups: logi TRUE
  ..$ FUN         :function (x)  
  ..$ outer       : NULL
  ..$ inner       : NULL
  ..$ labels      :List of 1
  .. ..$ d: chr "per capita demand deposits"
  ..$ units       : list()
> fm1Demand <-
+   lmer(log(d) ~ log(y) + log(rd) + log(rt) + log(rs) + (1|State) + (1|Year),
+        Demand)
> summary(fm1Demand)        # compare to output 3.13, p. 132
Linear mixed-effects model fit by REML
Formula: log(d) ~ log(y) + log(rd) + log(rt) + log(rs) + (1 | State) +      (1 | Year) 
   Data: Demand 
       AIC       BIC   logLik MLdeviance REMLdeviance
 -224.1652 -205.4148 120.0826  -260.5234    -240.1652
Random effects:
 Groups   Name        Variance   Std.Dev.
 Year     (Intercept) 0.00026461 0.016267
 State    (Intercept) 0.02944038 0.171582
 Residual             0.00111723 0.033425
# of obs: 77, groups: Year, 11; State, 7

Fixed effects:
             Estimate Std. Error DF t value  Pr(>|t|)    
(Intercept) -1.284710   0.723390 72 -1.7760  0.079965 .  
log(y)       1.069890   0.103922 72 10.2951 8.515e-16 ***
log(rd)     -0.295393   0.052460 72 -5.6308 3.248e-07 ***
log(rt)      0.039888   0.027891 72  1.4301  0.157004    
log(rs)     -0.326755   0.114390 72 -2.8565  0.005595 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Correlation of Fixed Effects:
        (Intr) log(y) lg(rd) lg(rt)
log(y)  -0.976                     
log(rd)  0.383 -0.227              
log(rt)  0.077 -0.062 -0.337       
log(rs)  0.444 -0.600 -0.270 -0.323
> 
> 
> 
> cleanEx(); ..nameEx <- "Genetics"
> 
> ### * Genetics
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Genetics
> ### Title: Heritability data
> ### Aliases: Genetics
> ### Keywords: datasets
> 
> ### ** Examples
> 
> options(contrasts = c(unordered = "contr.SAS", ordered = "contr.poly"))
> str(Genetics)
`data.frame':	60 obs. of  4 variables:
 $ Location: Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
 $ Block   : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...
 $ Family  : Factor w/ 5 levels "1","2","3","4",..: 1 1 1 2 2 2 3 3 3 4 ...
 $ Yield   : num  268 279 261 242 261 258 242 245 234 225 ...
 - attr(*, "ginfo")=List of 7
  ..$ formula     :Class 'formula' length 3 Yield ~ 1 | Location/Block
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ order.groups:List of 2
  .. ..$ Location: logi TRUE
  .. ..$ Block   : logi TRUE
  ..$ FUN         :function (x)  
  ..$ outer       : NULL
  ..$ inner       : NULL
  ..$ labels      :List of 1
  .. ..$ Yield: chr "Wheat yield"
  ..$ units       : list()
> Genetics$LocFam <- with(Genetics, Location:Family)
> Genetics$LocBloc <- with(Genetics, Location:Block)
> ## Not run: 
> ##D (fm1Gen <- lmer(Yield ~ 1 + (1|LocBloc) + (1|LocFam) + (1|Location), Genetics))
> ##D (fm2Gen <- lmer(Yield ~ Family + (1|LocBloc) + (1|LocFam) + (1|Location), Genetics))
> ## End(Not run)
> 
> 
> 
> options(contrasts = c(unordered = "contr.treatment",ordered = "contr.poly"))
> cleanEx(); ..nameEx <- "HR"
> 
> ### * HR
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: HR
> ### Title: Heart rates of patients on different drug treatments
> ### Aliases: HR
> ### Keywords: datasets
> 
> ### ** Examples
> 
> options(contrasts = c(unordered = "contr.SAS", ordered = "contr.poly"))
> if (require("latticeExtra", quietly = TRUE)) gplot(HR)
> (fm1HR <- lmer(HR ~ Time * Drug + baseHR + (Time|Patient), HR))  # linear trend in time
Linear mixed-effects model fit by REML
Formula: HR ~ Time * Drug + baseHR + (Time | Patient) 
   Data: HR 
     AIC      BIC    logLik MLdeviance REMLdeviance
 789.607 820.2694 -383.8035   788.1222      767.607
Random effects:
 Groups   Name        Variance Std.Dev. Corr   
 Patient  (Intercept) 60.633   7.7867          
          Time        37.784   6.1469   -0.563 
 Residual             24.361   4.9357          
# of obs: 120, groups: Patient, 24

Fixed effects:
             Estimate Std. Error  DF t value  Pr(>|t|)    
(Intercept)  33.97767   10.28299 113  3.3043  0.001276 ** 
Time         -3.19704    3.08492 113 -1.0363  0.302254    
Druga         3.59918    4.23139 113  0.8506  0.396795    
Drugb         7.09122    4.20941 113  1.6846  0.094824 .  
baseHR        0.54342    0.11615 113  4.6787 8.065e-06 ***
Time:Druga   -7.50131    4.36274 113 -1.7194  0.088279 .  
Time:Drugb   -3.98942    4.36274 113 -0.9144  0.362438    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> anova(fm1HR)
Analysis of Variance Table
          Df Sum Sq Mean Sq  Denom F value    Pr(>F)    
Time       1 379.23  379.23 113.00 15.5671 0.0001386 ***
Drug       2  92.88   46.44 113.00  1.9064 0.1533651    
baseHR     1 533.27  533.27 113.00 21.8904 8.065e-06 ***
Time:Drug  2  72.12   36.06 113.00  1.4802 0.2319775    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> ## Not run: 
> ##D fm2HR <- update(fm1HR, weights = varPower(0.5)) # use power-of-mean variance
> ##D summary(fm2HR)
> ##D intervals(fm2HR)             # variance function does not seem significant
> ##D anova(fm1HR, fm2HR)         # confirm with likelihood ratio
> ## End(Not run)
> (fm3HR <- lmer(HR ~ Time + Drug + baseHR + (Time|Patient), HR))
Linear mixed-effects model fit by REML
Formula: HR ~ Time + Drug + baseHR + (Time | Patient) 
   Data: HR 
      AIC      BIC    logLik MLdeviance REMLdeviance
 797.8283 822.9158 -389.9142   791.2093     779.8283
Random effects:
 Groups   Name        Variance Std.Dev. Corr   
 Patient  (Intercept) 61.560   7.8460          
          Time        40.963   6.4002   -0.571 
 Residual             24.361   4.9357          
# of obs: 120, groups: Patient, 24

Fixed effects:
             Estimate Std. Error  DF t value  Pr(>|t|)    
(Intercept)  36.04641   10.19449 115  3.5359 0.0005868 ***
Time         -7.02729    1.81789 115 -3.8656 0.0001839 ***
Druga        -0.45237    3.51456 115 -0.1287 0.8978084    
Drugb         4.93648    3.48807 115  1.4152 0.1596980    
baseHR        0.54342    0.11615 115  4.6787 7.937e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> anova(fm3HR)
Analysis of Variance Table
       Df Sum Sq Mean Sq  Denom F value    Pr(>F)    
Time    1 364.03  364.03 115.00 14.9431 0.0001839 ***
Drug    2  92.88   46.44 115.00  1.9064 0.1532828    
baseHR  1 533.27  533.27 115.00 21.8905 7.937e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> (fm4HR <- lmer(HR ~ Time + baseHR + (Time|Patient), HR))  # remove Drug term
Linear mixed-effects model fit by REML
Formula: HR ~ Time + baseHR + (Time | Patient) 
   Data: HR 
      AIC      BIC    logLik MLdeviance REMLdeviance
 805.1481 824.6605 -395.5740   794.2833     791.1481
Random effects:
 Groups   Name        Variance Std.Dev. Corr   
 Patient  (Intercept) 63.024   7.9388          
          Time        40.965   6.4004   -0.553 
 Residual             24.361   4.9357          
# of obs: 120, groups: Patient, 24

Fixed effects:
             Estimate Std. Error  DF t value  Pr(>|t|)    
(Intercept)  36.93167    9.90111 117  3.7301 0.0002968 ***
Time         -7.02729    1.81791 117 -3.8656 0.0001825 ***
baseHR        0.55077    0.11754 117  4.6859  7.59e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> anova(fm4HR)
Analysis of Variance Table
       Df Sum Sq Mean Sq  Denom F value    Pr(>F)    
Time    1 364.02  364.02 117.00  14.943 0.0001825 ***
baseHR  1 534.90  534.90 117.00  21.957  7.59e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> 
> 
> 
> options(contrasts = c(unordered = "contr.treatment",ordered = "contr.poly"))
> cleanEx(); ..nameEx <- "IncBlk"
> 
> ### * IncBlk
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: IncBlk
> ### Title: An unbalanced incomplete block experiment
> ### Aliases: IncBlk
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(IncBlk)
> str(IncBlk)
`data.frame':	24 obs. of  4 variables:
 $ Block    : Factor w/ 12 levels "1","10","11",..: 1 1 5 5 6 6 7 7 8 8 ...
 $ Treatment: Factor w/ 4 levels "1","2","3","4": 1 2 1 2 1 2 1 2 1 2 ...
 $ y        : num  0.62 0.91 0.41 0.48 0.41 0.49 0.26 0.28 0.29 0.37 ...
 $ x        : num  0.078 0.01 0.032 0.05 0 0.015 0.01 0.016 0.053 0.069 ...
 - attr(*, "ginfo")=List of 7
  ..$ formula     :Class 'formula' length 3 y ~ x | Block
  .. .. ..- attr(*, ".Environment")=length 3 <environment> 
  ..$ order.groups: logi TRUE
  ..$ FUN         :function (x)  
  ..$ outer       : NULL
  ..$ inner       :Class 'formula' length 2 ~Treatment
  .. .. ..- attr(*, ".Environment")=length 3 <environment> 
  ..$ labels      : list()
  ..$ units       : list()
> 
> 
> 
> cleanEx(); ..nameEx <- "Mississippi"
> 
> ### * Mississippi
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Mississippi
> ### Title: Nitrogen concentrations in the Mississippi River
> ### Aliases: Mississippi
> ### Keywords: datasets
> 
> ### ** Examples
> 
> options(contrasts = c(unordered = "contr.SAS", ordered = "contr.poly"))
> str(Mississippi)
`data.frame':	37 obs. of  3 variables:
 $ influent: Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 2 ...
 $ y       : num  21 27 29 17 19 12 29 20 20 21 ...
 $ Type    : Factor w/ 3 levels "1","2","3": 2 2 2 2 2 2 2 2 2 2 ...
 - attr(*, "ginfo")=List of 7
  ..$ formula     :Class 'formula' length 3 y ~ 1 | influent
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ order.groups: logi TRUE
  ..$ FUN         :function (x)  
  ..$ outer       :Class 'formula' length 2 ~Type
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ inner       : NULL
  ..$ labels      :List of 1
  .. ..$ y: chr "Nitrogen concentration in Mississippi River"
  ..$ units       :List of 1
  .. ..$ y: chr "(ppm)"
> #plot(Mississippi)
> fm1Miss <- lmer(y ~ 1 + (1|influent), Mississippi)
> summary(fm1Miss)        # compare with output 4.1, p. 142
Linear mixed-effects model fit by REML
Formula: y ~ 1 + (1 | influent) 
   Data: Mississippi 
      AIC      BIC    logLik MLdeviance REMLdeviance
 258.3511 263.1839 -126.1756   256.6398     252.3511
Random effects:
 Groups   Name        Variance Std.Dev.
 influent (Intercept) 63.324   7.9576  
 Residual             42.658   6.5313  
# of obs: 37, groups: influent, 6

Fixed effects:
            Estimate Std. Error DF t value  Pr(>|t|)    
(Intercept)   21.223      3.429 36  6.1892 3.885e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> fm1MLMiss <- update(fm1Miss, method = "ML")
> summary(fm1MLMiss)        # compare with output 4.2, p. 143
Linear mixed-effects model fit by maximum likelihood
Formula: y ~ 1 + (1 | influent) 
   Data: Mississippi 
     AIC      BIC    logLik MLdeviance REMLdeviance
 262.557 267.3898 -128.2785    256.557     252.4286
Random effects:
 Groups   Name        Variance Std.Dev.
 influent (Intercept) 52.679   7.2580  
 Residual             43.883   6.6245  
# of obs: 37, groups: influent, 6

Fixed effects:
            Estimate Std. Error DF t value  Pr(>|t|)    
(Intercept)   21.217      3.122 36   6.796 6.089e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> ranef(fm1MLMiss)          # BLUP's of random effects on p. 144
$influent
  (Intercept)
1   0.3097833
2  -6.5772278
3  -3.7862748
4   2.8826711
5  -5.8435210
6  13.0145691

attr(,"varFac")
attr(,"varFac")$influent
, , 1

          [,1]
[1,] 0.1016979

, , 2

          [,1]
[1,] 0.1276643

, , 3

          [,1]
[1,] 0.1714372

, , 4

          [,1]
[1,] 0.1463477

, , 5

          [,1]
[1,] 0.1714372

, , 6

          [,1]
[1,] 0.1714372


attr(,"stdErr")
[1] 6.534319
attr(,"class")
[1] "lmer.ranef"
attr(,"class")attr(,"package")
[1] "Matrix"
> ranef(fm1Miss)            # BLUP's of random effects on p. 142
$influent
  (Intercept)
1    0.309286
2   -6.719335
3   -3.897948
4    2.946106
5   -6.012988
6   13.374879

attr(,"varFac")
attr(,"varFac")$influent
, , 1

          [,1]
[1,] 0.1033736

, , 2

          [,1]
[1,] 0.1303161

, , 3

          [,1]
[1,] 0.1762533

, , 4

         [,1]
[1,] 0.149843

, , 5

          [,1]
[1,] 0.1762533

, , 6

          [,1]
[1,] 0.1762533


attr(,"stdErr")
[1] 6.531315
attr(,"class")
[1] "lmer.ranef"
attr(,"class")attr(,"package")
[1] "Matrix"
> #intervals(fm1Miss)       # interval estimates of variance components
> fm2Miss <- lmer(y ~ Type+(1|influent), Mississippi, method = "REML")
> summary(fm2Miss)         # compare to output 4.8 and 4.9, pp. 150-152
Linear mixed-effects model fit by REML
Formula: y ~ Type + (1 | influent) 
   Data: Mississippi 
      AIC      BIC    logLik MLdeviance REMLdeviance
 244.5246 252.5792 -117.2623   247.4686     234.5246
Random effects:
 Groups   Name        Variance Std.Dev.
 influent (Intercept) 14.970   3.8691  
 Residual             42.514   6.5202  
# of obs: 37, groups: influent, 6

Fixed effects:
            Estimate Std. Error DF t value  Pr(>|t|)    
(Intercept)  36.4000     4.8449 34  7.5131 1.011e-08 ***
Type1       -20.8000     5.9338 34 -3.5054  0.001302 ** 
Type2       -16.4619     5.5168 34 -2.9840  0.005238 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Correlation of Fixed Effects:
      (Intr) Type1 
Type1 -0.816       
Type2 -0.878  0.717
> anova(fm2Miss)
Analysis of Variance Table
     Df Sum Sq Mean Sq  Denom F value   Pr(>F)   
Type  2 541.76  270.88  34.00  6.3716 0.004466 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> 
> 
> 
> options(contrasts = c(unordered = "contr.treatment",ordered = "contr.poly"))
> cleanEx(); ..nameEx <- "Multilocation"
> 
> ### * Multilocation
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Multilocation
> ### Title: A multilocation trial
> ### Aliases: Multilocation
> ### Keywords: datasets
> 
> ### ** Examples
> 
> options(contrasts = c(unordered = "contr.SAS", ordered = "contr.poly"))
> str(Multilocation)
`data.frame':	108 obs. of  7 variables:
 $ obs     : num  3 4 6 7 9 10 12 16 19 20 ...
 $ Location: Factor w/ 9 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Block   : Factor w/ 3 levels "1","2","3": 1 1 1 1 2 2 2 2 3 3 ...
 $ Trt     : Factor w/ 4 levels "1","2","3","4": 3 4 2 1 2 1 3 4 1 2 ...
 $ Adj     : num  3.16 3.12 3.16 3.25 2.71 ...
 $ Fe      : num  7.10 6.68 6.83 6.53 8.25 ...
 $ Grp     : Factor w/ 27 levels "A/1","A/2","A/3",..: 1 1 1 1 2 2 2 2 3 3 ...
 - attr(*, "ginfo")=List of 7
  ..$ formula     :Class 'formula' length 3 Adj ~ 1 | Location/Block
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ order.groups:List of 2
  .. ..$ Location: logi TRUE
  .. ..$ Block   : logi TRUE
  ..$ FUN         :function (x)  
  ..$ outer       : NULL
  ..$ inner       :List of 1
  .. ..$ Block:Class 'formula' length 2 ~Trt
  .. .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ labels      :List of 1
  .. ..$ Adj: chr "Adjusted yield"
  ..$ units       : list()
> ### Create a Block 
> Multilocation$Grp <- with(Multilocation, Block:Location)
> (fm1Mult <- lmer(Adj ~ Location * Trt + (1|Grp), Multilocation))
Linear mixed-effects model fit by REML
Formula: Adj ~ Location * Trt + (1 | Grp) 
   Data: Multilocation 
      AIC      BIC    logLik MLdeviance REMLdeviance
 86.64621 188.5672 -5.323106  -87.14598     10.64621
Random effects:
 Groups   Name        Variance  Std.Dev.
 Grp      (Intercept) 0.0056193 0.074962
 Residual             0.0345787 0.185953
# of obs: 108, groups: Grp, 27

Fixed effects:
                Estimate Std. Error DF t value  Pr(>|t|)    
(Intercept)     2.359233   0.115755 72 20.3812 < 2.2e-16 ***
LocationA       0.649300   0.163703 72  3.9663 0.0001705 ***
LocationB       0.066433   0.163703 72  0.4058 0.6860811    
LocationC       0.545333   0.163703 72  3.3312 0.0013667 ** 
LocationD       0.374133   0.163703 72  2.2854 0.0252337 *  
LocationE       0.550000   0.163703 72  3.3597 0.0012505 ** 
LocationF       0.998100   0.163703 72  6.0970 4.861e-08 ***
LocationG       0.360567   0.163703 72  2.2026 0.0308276 *  
LocationH       1.014033   0.163703 72  6.1943 3.252e-08 ***
Trt1            0.227200   0.151830 72  1.4964 0.1389186    
Trt2           -0.001400   0.151830 72 -0.0092 0.9926685    
Trt3            0.423233   0.151830 72  2.7875 0.0067874 ** 
LocationA:Trt1 -0.188533   0.214721 72 -0.8780 0.3828425    
LocationB:Trt1 -0.275233   0.214721 72 -1.2818 0.2040178    
LocationC:Trt1 -0.040000   0.214721 72 -0.1863 0.8527423    
LocationD:Trt1 -0.535133   0.214721 72 -2.4922 0.0149969 *  
LocationE:Trt1 -0.262967   0.214721 72 -1.2247 0.2246830    
LocationF:Trt1 -0.271533   0.214721 72 -1.2646 0.2100968    
LocationG:Trt1  0.203233   0.214721 72  0.9465 0.3470587    
LocationH:Trt1 -0.149533   0.214721 72 -0.6964 0.4884150    
LocationA:Trt2 -0.093467   0.214721 72 -0.4353 0.6646509    
LocationB:Trt2 -0.322733   0.214721 72 -1.5030 0.1372028    
LocationC:Trt2  0.089600   0.214721 72  0.4173 0.6777105    
LocationD:Trt2 -0.296933   0.214721 72 -1.3829 0.1709748    
LocationE:Trt2 -0.306933   0.214721 72 -1.4295 0.1571983    
LocationF:Trt2 -0.309933   0.214721 72 -1.4434 0.1532374    
LocationG:Trt2 -0.108600   0.214721 72 -0.5058 0.6145606    
LocationH:Trt2 -0.330600   0.214721 72 -1.5397 0.1280231    
LocationA:Trt3 -0.402467   0.214721 72 -1.8744 0.0649358 .  
LocationB:Trt3 -0.565500   0.214721 72 -2.6337 0.0103329 *  
LocationC:Trt3 -0.122467   0.214721 72 -0.5704 0.5702135    
LocationD:Trt3 -0.548400   0.214721 72 -2.5540 0.0127654 *  
LocationE:Trt3 -0.328633   0.214721 72 -1.5305 0.1302711    
LocationF:Trt3 -0.462567   0.214721 72 -2.1543 0.0345659 *  
LocationG:Trt3 -0.252967   0.214721 72 -1.1781 0.2426279    
LocationH:Trt3 -0.372033   0.214721 72 -1.7326 0.0874414 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> anova(fm1Mult)
Analysis of Variance Table
             Df Sum Sq Mean Sq  Denom F value    Pr(>F)    
Location      8  6.947   0.868 72.000 25.1147 < 2.2e-16 ***
Trt           3  1.222   0.407 72.000 11.7774 2.307e-06 ***
Location:Trt 24  0.997   0.042 72.000  1.2008    0.2710    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> (fm2Mult <- lmer(Adj ~ Location + Trt + (1|Grp), Multilocation))
Linear mixed-effects model fit by REML
Formula: Adj ~ Location + Trt + (1 | Grp) 
   Data: Multilocation 
      AIC      BIC   logLik MLdeviance REMLdeviance
 21.99894 59.54877 3.000531  -51.21968    -6.001063
Random effects:
 Groups   Name        Variance  Std.Dev.
 Grp      (Intercept) 0.0050851 0.07131 
 Residual             0.0367154 0.19161 
# of obs: 108, groups: Grp, 27

Fixed effects:
             Estimate Std. Error DF t value  Pr(>|t|)    
(Intercept)  2.532965   0.075990 96 33.3327 < 2.2e-16 ***
LocationA    0.478183   0.097516 96  4.9037 3.828e-06 ***
LocationB   -0.224433   0.097516 96 -2.3015 0.0235251 *  
LocationC    0.527117   0.097516 96  5.4055 4.710e-07 ***
LocationD    0.029017   0.097516 96  0.2976 0.7666828    
LocationE    0.325367   0.097516 96  3.3366 0.0012075 ** 
LocationF    0.737092   0.097516 96  7.5587 2.411e-11 ***
LocationG    0.320983   0.097516 96  3.2916 0.0013947 ** 
LocationH    0.800992   0.097516 96  8.2140 9.996e-13 ***
Trt1         0.058344   0.052150 96  1.1188 0.2660283    
Trt2        -0.188022   0.052150 96 -3.6054 0.0004966 ***
Trt3         0.083785   0.052150 96  1.6066 0.1114247    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> fm3Mult <- lmer(Adj ~ Location + (1|Grp), Multilocation)
> fm4Mult <- lmer(Adj ~ Trt + (1|Grp), Multilocation)
> fm5Mult <- lmer(Adj ~ 1 + (1|Grp), Multilocation)
> anova(fm2Mult)
Analysis of Variance Table
         Df Sum Sq Mean Sq  Denom F value    Pr(>F)    
Location  8  7.377   0.922 96.000  25.115 < 2.2e-16 ***
Trt       3  1.222   0.407 96.000  11.092 2.571e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> anova(fm1Mult, fm2Mult, fm3Mult, fm4Mult, fm5Mult)
Data: Multilocation
Models:
fm5Mult: Adj ~ 1 + (1 | Grp)
fm4Mult: Adj ~ Trt + (1 | Grp)
fm3Mult: Adj ~ Location + (1 | Grp)
fm2Mult: Adj ~ Location + Trt + (1 | Grp)
fm1Mult: Adj ~ Location * Trt + (1 | Grp)
        Df     AIC     BIC  logLik  Chisq Chi Df Pr(>Chisq)    
fm5Mult  3  53.327  61.374 -23.664                             
fm4Mult  6  43.506  59.598 -15.753 15.822      3  0.0012336 ** 
fm3Mult 11  31.820  61.324  -4.910 21.685      5  0.0006009 ***
fm2Mult 14  21.999  59.549   3.001 15.822      3  0.0012336 ** 
fm1Mult 38  86.646 188.567  -5.323  0.000     24  1.0000000    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> Multilocation$LocTrt <- with(Multilocation, Location:Trt)
> ### Treating the location as a random effect
> ## Not run: 
> ##D (fm1MultR <- lmer(Adj ~ Trt + (1|LocTrt) + (1|Location) + (1|Grp), Multilocation))
> ##D #intervals(fm1MultR)
> ##D anova(fm1MultR)
> ##D fm2MultR <- lmer(Adj ~ Trt + (Trt - 1|Location) + (1|Block), Multilocation)
> ##D anova(fm1MultR, fm2MultR)
> ## End(Not run)
> 
> 
> 
> options(contrasts = c(unordered = "contr.treatment",ordered = "contr.poly"))
> cleanEx(); ..nameEx <- "PBIB"
> 
> ### * PBIB
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: PBIB
> ### Title: A partially balanced incomplete block experiment
> ### Aliases: PBIB
> ### Keywords: datasets
> 
> ### ** Examples
> 
> options(contrasts = c(unordered = "contr.SAS", ordered = "contr.poly"))
> str(PBIB)
`data.frame':	60 obs. of  3 variables:
 $ response : num  2.4 2.5 2.6 2 2.7 2.8 2.4 2.7 2.6 2.8 ...
 $ Treatment: Factor w/ 15 levels "1","10","11",..: 7 15 1 5 11 13 14 1 2 1 ...
 $ Block    : Factor w/ 15 levels "1","10","11",..: 1 1 1 1 8 8 8 8 9 9 ...
 - attr(*, "ginfo")=List of 7
  ..$ formula     :Class 'formula' length 3 response ~ Treatment | Block
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ order.groups: logi TRUE
  ..$ FUN         :function (x)  
  ..$ outer       : NULL
  ..$ inner       : NULL
  ..$ labels      : list()
  ..$ units       : list()
> fm1PBIB <- lmer(response ~ Treatment + (1|Block), PBIB)
> summary(fm1PBIB)     # compare with output 1.7  pp. 24-25
Linear mixed-effects model fit by REML
Formula: response ~ Treatment + (1 | Block) 
   Data: PBIB 
     AIC      BIC    logLik MLdeviance REMLdeviance
 85.9849 121.5888 -25.99245   22.82831     51.98489
Random effects:
 Groups   Name        Variance Std.Dev.
 Block    (Intercept) 0.046522 0.21569 
 Residual             0.085559 0.29250 
# of obs: 60, groups: Block, 15

Fixed effects:
              Estimate Std. Error DF t value Pr(>|t|)    
(Intercept)  2.8913111  0.1664127 45 17.3743  < 2e-16 ***
Treatment1  -0.0737886  0.2220608 45 -0.3323  0.74121    
Treatment10 -0.4002495  0.2220608 45 -1.8024  0.07818 .  
Treatment11  0.0073879  0.2220608 45  0.0333  0.97361    
Treatment12  0.1615103  0.2220608 45  0.7273  0.47079    
Treatment13 -0.2735419  0.2220608 45 -1.2318  0.22441    
Treatment14 -0.4000000  0.2272003 45 -1.7606  0.08511 .  
Treatment15 -0.0320781  0.2220608 45 -0.1445  0.88579    
Treatment2  -0.4859962  0.2220608 45 -2.1886  0.03386 *  
Treatment3  -0.4363680  0.2220608 45 -1.9651  0.05560 .  
Treatment4  -0.1074807  0.2272003 45 -0.4731  0.63845    
Treatment5  -0.0864131  0.2220608 45 -0.3891  0.69901    
Treatment6   0.0193828  0.2220608 45  0.0873  0.93083    
Treatment7  -0.1023261  0.2220608 45 -0.4608  0.64716    
Treatment8  -0.1097056  0.2220608 45 -0.4940  0.62369    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Correlation of Fixed Effects:
            (Intr) Trtmn1 Trtm10 Trtm11 Trtm12 Trtm13 Trtm14 Trtm15 Trtmn2
Treatment1  -0.667                                                        
Treatment10 -0.667  0.500                                                 
Treatment11 -0.667  0.477  0.500                                          
Treatment12 -0.667  0.500  0.500  0.500                                   
Treatment13 -0.667  0.500  0.500  0.500  0.500                            
Treatment14 -0.683  0.512  0.512  0.512  0.512  0.512                     
Treatment15 -0.667  0.500  0.477  0.500  0.500  0.500  0.512              
Treatment2  -0.667  0.500  0.500  0.500  0.477  0.500  0.512  0.500       
Treatment3  -0.667  0.500  0.500  0.500  0.500  0.477  0.512  0.500  0.500
Treatment4  -0.683  0.512  0.512  0.512  0.512  0.512  0.500  0.512  0.512
Treatment5  -0.667  0.500  0.477  0.500  0.500  0.500  0.512  0.477  0.500
Treatment6  -0.667  0.477  0.500  0.477  0.500  0.500  0.512  0.500  0.500
Treatment7  -0.667  0.500  0.500  0.500  0.477  0.500  0.512  0.500  0.477
Treatment8  -0.667  0.500  0.500  0.500  0.500  0.477  0.512  0.500  0.500
            Trtmn3 Trtmn4 Trtmn5 Trtmn6 Trtmn7
Treatment1                                    
Treatment10                                   
Treatment11                                   
Treatment12                                   
Treatment13                                   
Treatment14                                   
Treatment15                                   
Treatment2                                    
Treatment3                                    
Treatment4   0.512                            
Treatment5   0.500  0.512                     
Treatment6   0.500  0.512  0.500              
Treatment7   0.500  0.512  0.500  0.500       
Treatment8   0.477  0.512  0.500  0.500  0.500
> anova(fm1PBIB)
Analysis of Variance Table
          Df Sum Sq Mean Sq  Denom F value Pr(>F)
Treatment 14  1.834   0.131 45.000  1.5312 0.1388
> 
> 
> 
> options(contrasts = c(unordered = "contr.treatment",ordered = "contr.poly"))
> cleanEx(); ..nameEx <- "SIMS"
> 
> ### * SIMS
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: SIMS
> ### Title: Second International Mathematics Study data
> ### Aliases: SIMS
> ### Keywords: datasets
> 
> ### ** Examples
> 
> str(SIMS)
`data.frame':	3691 obs. of  3 variables:
 $ Pretot: num  29 38 31 31 29 23 23 33 30 32 ...
 $ Gain  : num  2 0 6 6 5 9 7 2 1 3 ...
 $ Class : Factor w/ 190 levels "1","10","100",..: 1 1 1 1 1 1 1 1 1 1 ...
 - attr(*, "ginfo")=List of 7
  ..$ formula     :Class 'formula' length 3 Gain ~ Pretot | Class
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ order.groups: logi TRUE
  ..$ FUN         :function (x)  
  ..$ outer       : NULL
  ..$ inner       : NULL
  ..$ labels      :List of 2
  .. ..$ Pretot: chr "Sum of pre-test core item scores"
  .. ..$ Gain  : chr "Gain in mathematics achievement score"
  ..$ units       : list()
> fm1SIMS <- lmer(Gain ~ Pretot + (Pretot | Class), data = SIMS,
+                 control = list(msVerbose = TRUE))
iter    0 value 22383.012955
final  value 22383.012947 
converged
> summary(fm1SIMS)              # compare to output 7.4, p. 262
Linear mixed-effects model fit by REML
Formula: Gain ~ Pretot + (Pretot | Class) 
   Data: SIMS 
      AIC      BIC    logLik MLdeviance REMLdeviance
 22395.01 22432.29 -11191.51   22375.65     22383.01
Random effects:
 Groups   Name        Variance   Std.Dev. Corr   
 Class    (Intercept) 11.7351668 3.425663        
          Pretot       0.0098813 0.099405 -0.417 
 Residual             22.1942169 4.711074        
# of obs: 3691, groups: Class, 190

Fixed effects:
               Estimate  Std. Error   DF t value  Pr(>|t|)    
(Intercept)    7.062263    0.345898 3689  20.417 < 2.2e-16 ***
Pretot        -0.191333    0.016492 3689 -11.601 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Correlation of Fixed Effects:
       (Intr)
Pretot -0.713
> anova(fm1SIMS)
Analysis of Variance Table
       Df Sum Sq Mean Sq  Denom F value    Pr(>F)    
Pretot  1 2987.2  2987.2 3689.0  134.59 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> 
> 
> 
> cleanEx(); ..nameEx <- "Semi2"
> 
> ### * Semi2
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Semi2
> ### Title: Oxide layer thicknesses on semiconductors
> ### Aliases: Semi2
> ### Keywords: datasets
> 
> ### ** Examples
> 
> options(contrasts = c(unordered = "contr.SAS", ordered = "contr.poly"))
> str(Semi2)
`data.frame':	72 obs. of  5 variables:
 $ Source   : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
 $ Lot      : Factor w/ 8 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 2 ...
 $ Wafer    : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 3 3 3 1 ...
 $ Site     : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...
 $ Thickness: num  2006 1999 2007 1980 1988 ...
 - attr(*, "ginfo")=List of 7
  ..$ formula     :Class 'formula' length 3 Thickness ~ 1 | Lot/Wafer
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ order.groups:List of 2
  .. ..$ Lot  : logi TRUE
  .. ..$ Wafer: logi TRUE
  ..$ FUN         :function (x)  
  ..$ outer       :Class 'formula' length 2 ~Source
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ inner       : NULL
  ..$ labels      :List of 1
  .. ..$ Thickness: chr "Thickness of oxide layer"
  ..$ units       : list()
> xtabs(~Lot + Wafer, Semi2)
   Wafer
Lot 1 2 3
  1 3 3 3
  2 3 3 3
  3 3 3 3
  4 3 3 3
  5 3 3 3
  6 3 3 3
  7 3 3 3
  8 3 3 3
> Semi2$LotWafer <- with(Semi2, Lot:Wafer)
> fm1Semi2 <- lmer(Thickness ~ 1 + (1|LotWafer) + (1|Lot), Semi2)
> summary(fm1Semi2)       # compare with output 4.13, p. 156
Linear mixed-effects model fit by REML
Formula: Thickness ~ 1 + (1 | LotWafer) + (1 | Lot) 
   Data: Semi2 
      AIC      BIC    logLik MLdeviance REMLdeviance
 462.0221 471.1287 -227.0110   458.7378     454.0221
Random effects:
 Groups   Name        Variance Std.Dev.
 LotWafer (Intercept)  35.866   5.9888 
 Lot      (Intercept) 129.857  11.3955 
 Residual              12.570   3.5454 
# of obs: 72, groups: LotWafer, 24; Lot, 8

Fixed effects:
            Estimate Std. Error DF t value  Pr(>|t|)    
(Intercept) 2000.153      4.231 71  472.74 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> fm2Semi2 <- lmer(Thickness ~ Source + (1|LotWafer) + (1|Lot), Semi2)
> summary(fm2Semi2)       # compare with output 4.15, p. 159
Linear mixed-effects model fit by REML
Formula: Thickness ~ Source + (1 | LotWafer) + (1 | Lot) 
   Data: Semi2 
      AIC      BIC    logLik MLdeviance REMLdeviance
 456.4779 467.8612 -223.2389   457.1361     446.4779
Random effects:
 Groups   Name        Variance Std.Dev.
 LotWafer (Intercept)  35.866   5.9889 
 Lot      (Intercept) 119.825  10.9465 
 Residual              12.570   3.5454 
# of obs: 72, groups: LotWafer, 24; Lot, 8

Fixed effects:
             Estimate Std. Error DF  t value Pr(>|t|)    
(Intercept) 2005.1944     5.7701 70 347.5128   <2e-16 ***
Source1      -10.0833     8.1602 70  -1.2357   0.2207    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Correlation of Fixed Effects:
        (Intr)
Source1 -0.707
> anova(fm2Semi2)
Analysis of Variance Table
       Df Sum Sq Mean Sq  Denom F value Pr(>F)
Source  1 19.193  19.193 70.000  1.5269 0.2207
> Semi2$LotSrc <- with(Semi2, Lot:Source)
> fm3Semi2 <- lmer(Thickness ~ Source + (1|LotWafer) + (1|LotSrc), Semi2)
> summary(fm3Semi2)       # compare with output 4.17, p. 163
Linear mixed-effects model fit by REML
Formula: Thickness ~ Source + (1 | LotWafer) + (1 | LotSrc) 
   Data: Semi2 
      AIC      BIC    logLik MLdeviance REMLdeviance
 456.4779 467.8612 -223.2389   457.1361     446.4779
Random effects:
 Groups   Name        Variance Std.Dev.
 LotWafer (Intercept)  35.866   5.9889 
 LotSrc   (Intercept) 119.825  10.9465 
 Residual              12.570   3.5454 
# of obs: 72, groups: LotWafer, 24; LotSrc, 8

Fixed effects:
             Estimate Std. Error DF  t value Pr(>|t|)    
(Intercept) 2005.1944     5.7701 70 347.5128   <2e-16 ***
Source1      -10.0833     8.1602 70  -1.2357   0.2207    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Correlation of Fixed Effects:
        (Intr)
Source1 -0.707
> ## This is not the same as the SAS model.
> 
> 
> 
> options(contrasts = c(unordered = "contr.treatment",ordered = "contr.poly"))
> cleanEx(); ..nameEx <- "Semiconductor"
> 
> ### * Semiconductor
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Semiconductor
> ### Title: Semiconductor split-plot experiment
> ### Aliases: Semiconductor
> ### Keywords: datasets
> 
> ### ** Examples
> 
> options(contrasts = c(unordered = "contr.SAS", ordered = "contr.poly"))
> str(Semiconductor)
`data.frame':	48 obs. of  5 variables:
 $ resistance: num  5.22 5.61 6.11 6.33 6.13 6.14 5.6 5.91 5.49 4.6 ...
 $ ET        : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
 $ Wafer     : Factor w/ 3 levels "1","2","3": 1 1 1 1 2 2 2 2 3 3 ...
 $ position  : Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
 $ Grp       : Factor w/ 12 levels "1/1","1/2","1/3",..: 1 1 1 1 2 2 2 2 3 3 ...
 - attr(*, "ginfo")=List of 7
  ..$ formula     :Class 'formula' length 3 resistance ~ 1 | Grp
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ order.groups: logi TRUE
  ..$ FUN         :function (x)  
  ..$ outer       :Class 'formula' length 2 ~ET
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ inner       :Class 'formula' length 2 ~position
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ labels      : list()
  ..$ units       : list()
> (fm1Semi <- lmer(resistance ~ ET * position + (1|Grp), Semiconductor))
Linear mixed-effects model fit by REML
Formula: resistance ~ ET * position + (1 | Grp) 
   Data: Semiconductor 
     AIC      BIC    logLik MLdeviance REMLdeviance
 86.6505 120.3321 -25.32525   30.14672      50.6505
Random effects:
 Groups   Name        Variance Std.Dev.
 Grp      (Intercept) 0.10579  0.32525 
 Residual             0.11115  0.33339 
# of obs: 48, groups: Grp, 12

Fixed effects:
               Estimate Std. Error DF t value Pr(>|t|)    
(Intercept)    6.540000   0.268911 32 24.3203  < 2e-16 ***
ET1           -0.653333   0.380298 32 -1.7180  0.09547 .  
ET2           -0.623333   0.380298 32 -1.6391  0.11100    
ET3           -0.446667   0.380298 32 -1.1745  0.24885    
position1     -0.200000   0.272212 32 -0.7347  0.46786    
position2      0.013333   0.272212 32  0.0490  0.96124    
position3     -0.643333   0.272212 32 -2.3634  0.02435 *  
ET1:position1 -0.073333   0.384966 32 -0.1905  0.85013    
ET2:position1  0.276667   0.384966 32  0.7187  0.47755    
ET3:position1  0.243333   0.384966 32  0.6321  0.53182    
ET1:position2 -0.450000   0.384966 32 -1.1689  0.25106    
ET2:position2  0.256667   0.384966 32  0.6667  0.50973    
ET3:position2  0.240000   0.384966 32  0.6234  0.53742    
ET1:position3  0.310000   0.384966 32  0.8053  0.42661    
ET2:position3  0.493333   0.384966 32  1.2815  0.20923    
ET3:position3  0.323333   0.384966 32  0.8399  0.40720    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> anova(fm1Semi)
Analysis of Variance Table
            Df Sum Sq Mean Sq  Denom F value  Pr(>F)  
ET           3  0.647   0.216 32.000  1.9415 0.14273  
position     3  1.129   0.376 32.000  3.3855 0.02991 *
ET:position  9  0.809   0.090 32.000  0.8092 0.61127  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> (fm2Semi <- lmer(resistance ~ ET + position + (1|Grp), Semiconductor))
Linear mixed-effects model fit by REML
Formula: resistance ~ ET + position + (1 | Grp) 
   Data: Semiconductor 
     AIC    BIC   logLik MLdeviance REMLdeviance
 71.0862 87.927 -26.5431   40.11902      53.0862
Random effects:
 Groups   Name        Variance Std.Dev.
 Grp      (Intercept) 0.10724  0.32747 
 Residual             0.10537  0.32460 
# of obs: 48, groups: Grp, 12

Fixed effects:
             Estimate Std. Error DF t value Pr(>|t|)    
(Intercept)  6.438750   0.226078 41 28.4802  < 2e-16 ***
ET1         -0.706667   0.298415 41 -2.3681  0.02268 *  
ET2         -0.366667   0.298415 41 -1.2287  0.22619    
ET3         -0.245000   0.298415 41 -0.8210  0.41639    
position1   -0.088333   0.132518 41 -0.6666  0.50878    
position2    0.025000   0.132518 41  0.1887  0.85129    
position3   -0.361667   0.132518 41 -2.7292  0.00931 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> anova(fm2Semi)
Analysis of Variance Table
         Df Sum Sq Mean Sq  Denom F value  Pr(>F)  
ET        3  0.614   0.205 41.000  1.9415 0.13796  
position  3  1.129   0.376 41.000  3.5714 0.02198 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> 
> 
> 
> options(contrasts = c(unordered = "contr.treatment",ordered = "contr.poly"))
> cleanEx(); ..nameEx <- "TeachingI"
> 
> ### * TeachingI
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: TeachingI
> ### Title: Teaching Methods I
> ### Aliases: TeachingI
> ### Keywords: datasets
> 
> ### ** Examples
> 
> str(TeachingI)
`data.frame':	96 obs. of  7 variables:
 $ Method    : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
 $ Teacher   : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 2 2 ...
 $ Gender    : Factor w/ 2 levels "f","m": 1 1 1 1 2 2 2 2 1 1 ...
 $ Student   : Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
 $ score     : num  15 17 16 16 17 16 17 17 18 17 ...
 $ Experience: num  11 11 11 11 11 11 11 11 8 8 ...
 $ uTeacher  : Factor w/ 12 levels "1/1","1/2","1/3",..: 1 1 1 1 1 1 1 1 2 2 ...
 - attr(*, "ginfo")=List of 7
  ..$ formula     :Class 'formula' length 3 score ~ 1 | uTeacher
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ order.groups: logi TRUE
  ..$ FUN         :function (x)  
  ..$ outer       : NULL
  ..$ inner       :Class 'formula' length 2 ~Gender
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ labels      : list()
  ..$ units       : list()
> 
> 
> 
> cleanEx(); ..nameEx <- "TeachingII"
> 
> ### * TeachingII
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: TeachingII
> ### Title: Teaching Methods II
> ### Aliases: TeachingII
> ### Keywords: datasets
> 
> ### ** Examples
> 
> str(TeachingII)
`data.frame':	96 obs. of  6 variables:
 $ Method  : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
 $ Teacher : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 2 2 ...
 $ Gender  : Factor w/ 2 levels "f","m": 1 1 1 1 2 2 2 2 1 1 ...
 $ IQ      : num  89 105 108 116 95 103 91 82 83 103 ...
 $ score   : num  54 55 54 64 59 58 42 48 48 56 ...
 $ uTeacher: Factor w/ 12 levels "1/1","1/2","1/3",..: 1 1 1 1 1 1 1 1 2 2 ...
 - attr(*, "ginfo")=List of 7
  ..$ formula     :Class 'formula' length 3 score ~ IQ | uTeacher
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ order.groups: logi TRUE
  ..$ FUN         :function (x)  
  ..$ outer       : NULL
  ..$ inner       :Class 'formula' length 2 ~Gender
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ labels      :List of 1
  .. ..$ IQ: chr "Intelligence Quotient"
  ..$ units       : list()
> 
> 
> 
> cleanEx(); ..nameEx <- "WWheat"
> 
> ### * WWheat
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: WWheat
> ### Title: Winter wheat
> ### Aliases: WWheat
> ### Keywords: datasets
> 
> ### ** Examples
> 
> str(WWheat)
`data.frame':	60 obs. of  3 variables:
 $ Variety : Factor w/ 10 levels "1","10","2","3",..: 1 1 1 1 1 1 3 3 3 3 ...
 $ Yield   : num  41 69 53 66 64 64 49 44 44 46 ...
 $ Moisture: num  10 57 32 52 47 48 30 21 20 26 ...
 - attr(*, "ginfo")=List of 7
  ..$ formula     :Class 'formula' length 3 Yield ~ Moisture | Variety
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ order.groups: logi TRUE
  ..$ FUN         :function (x)  
  ..$ outer       : NULL
  ..$ inner       : NULL
  ..$ labels      :List of 2
  .. ..$ Moisture: chr "Preplant moisture content of top 36 inches of soil"
  .. ..$ Yield   : chr "Winter Wheat Yield"
  ..$ units       :List of 1
  .. ..$ Yield: chr "Bu/Acre"
> 
> 
> 
> cleanEx(); ..nameEx <- "WaferTypes"
> 
> ### * WaferTypes
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: WaferTypes
> ### Title: Data on different types of silicon wafers
> ### Aliases: WaferTypes
> ### Keywords: datasets
> 
> ### ** Examples
> 
> str(WaferTypes)
`data.frame':	144 obs. of  8 variables:
 $ Group      : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
 $ Temperature: Factor w/ 3 levels "1000","1100",..: 3 3 3 3 3 3 3 3 3 3 ...
 $ Type       : Factor w/ 2 levels "A","B": 1 1 1 1 1 1 2 2 2 2 ...
 $ Wafer      : num  1 1 1 2 2 2 1 1 1 2 ...
 $ Site       : num  1 2 3 1 2 3 1 2 3 1 ...
 $ delta      : num  291 295 294 318 315 315 349 348 345 332 ...
 $ Thick      : num  1919 1919 1919 2113 2113 ...
 $ uWafer     : Factor w/ 48 levels "1/1000/A/1","1/1000/A/2",..: 9 9 9 10 10 10 11 11 11 12 ...
 - attr(*, "ginfo")=List of 7
  ..$ formula     :Class 'formula' length 3 delta ~ 1 | uWafer
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ order.groups: logi TRUE
  ..$ FUN         :function (x)  
  ..$ outer       : NULL
  ..$ inner       : NULL
  ..$ labels      :List of 1
  .. ..$ delta: chr "Rate of deposition of polysilicon"
  ..$ units       : list()
> 
> 
> 
> cleanEx(); ..nameEx <- "Weights"
> 
> ### * Weights
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Weights
> ### Title: Data from a weight-lifting program
> ### Aliases: Weights
> ### Keywords: datasets
> 
> ### ** Examples
> 
> options(contrasts = c(unordered = "contr.SAS", ordered = "contr.poly"))
> str(Weights)
`data.frame':	399 obs. of  5 variables:
 $ strength: num  85 85 86 85 87 86 87 80 79 79 ...
 $ Subject : Factor w/ 21 levels "1","10","11",..: 1 1 1 1 1 1 1 12 12 12 ...
 $ Program : Factor w/ 3 levels "CONT","RI","WI": 1 1 1 1 1 1 1 1 1 1 ...
 $ Subj    : Factor w/ 57 levels "CONT/1","CONT/10",..: 1 1 1 1 1 1 1 12 12 12 ...
 $ Time    : num  1 3 5 7 9 11 13 1 3 5 ...
 - attr(*, "ginfo")=List of 7
  ..$ formula     :Class 'formula' length 3 strength ~ Time | Subj
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ order.groups: logi TRUE
  ..$ FUN         :function (x)  
  ..$ outer       :Class 'formula' length 2 ~Program
  .. .. ..- attr(*, ".Environment")=length 2 <environment> 
  ..$ inner       : NULL
  ..$ labels      :List of 2
  .. ..$ Time    : chr "Time since beginning of study"
  .. ..$ strength: chr "Strength"
  ..$ units       :List of 1
  .. ..$ Time: chr "(days)"
> fm1Weight <- lmer(strength ~ Program * Time + (1|Subj), Weights)
Warning in "LMEoptimize<-"(`*tmp*`, value = list(maxIter = 200, msMaxIter = 200,  : 
	 optim or nlminb returned message ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 

> summary(fm1Weight)               # compare with output 3.1, p. 91
Linear mixed-effects model fit by REML
Formula: strength ~ Program * Time + (1 | Subj) 
   Data: Weights 
      AIC      BIC    logLik MLdeviance REMLdeviance
 1455.363 1487.275 -719.6815   1425.918     1439.363
Random effects:
 Groups   Name        Variance Std.Dev.
 Subj     (Intercept) 9.6047   3.0991  
 Residual             1.1876   1.0897  
# of obs: 399, groups: Subj, 57

Fixed effects:
                   Estimate Std. Error  DF  t value  Pr(>|t|)    
(Intercept)       81.105442   0.700132 393 115.8432 < 2.2e-16 ***
ProgramCONT       -1.115264   1.002436 393  -1.1126    0.2666    
ProgramRI         -1.045174   1.064683 393  -0.9817    0.3269    
Time               0.159864   0.022470 393   7.1145 5.357e-12 ***
ProgramCONT:Time  -0.183971   0.032172 393  -5.7183 2.133e-08 ***
ProgramRI:Time    -0.054953   0.034170 393  -1.6082    0.1086    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Correlation of Fixed Effects:
            (Intr) PrCONT PrgrRI Time   PCONT:
ProgramCONT -0.698                            
ProgramRI   -0.658  0.459                     
Time        -0.225  0.157  0.148              
PrgrmCONT:T  0.157 -0.225 -0.103 -0.698       
ProgrmRI:Tm  0.148 -0.103 -0.225 -0.658  0.459
> anova(fm1Weight)
Analysis of Variance Table
             Df Sum Sq Mean Sq  Denom F value    Pr(>F)    
Program       2   7.28    3.64 393.00  3.0651   0.04776 *  
Time          1  40.74   40.74 393.00 34.3079 9.951e-09 ***
Program:Time  2  40.39   20.20 393.00 17.0061 8.253e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> fm2Weight <- lmer(strength ~ Program * Time + (Time|Subj), Weights)
> anova(fm1Weight, fm2Weight)
Data: Weights
Models:
fm1Weight: strength ~ Program * Time + (1 | Subj)
fm2Weight: strength ~ Program * Time + (Time | Subj)
          Df     AIC     BIC  logLik  Chisq Chi Df Pr(>Chisq)    
fm1Weight  8 1455.36 1487.27 -719.68                             
fm2Weight 10 1343.40 1383.29 -661.70 115.97      2  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> summary(fm2Weight)
Linear mixed-effects model fit by REML
Formula: strength ~ Program * Time + (Time | Subj) 
   Data: Weights 
      AIC      BIC    logLik MLdeviance REMLdeviance
 1343.396 1383.285 -661.6979   1313.470     1323.396
Random effects:
 Groups   Name        Variance Std.Dev. Corr   
 Subj     (Intercept) 9.033464 3.00557         
          Time        0.031084 0.17631  -0.118 
 Residual             0.633018 0.79562         
# of obs: 399, groups: Subj, 57

Fixed effects:
                   Estimate Std. Error  DF  t value  Pr(>|t|)    
(Intercept)       81.105442   0.669073 393 121.2206 < 2.2e-16 ***
ProgramCONT       -1.115264   0.957967 393  -1.1642 0.2450500    
ProgramRI         -1.045174   1.017454 393  -1.0272 0.3049368    
Time               0.159864   0.041825 393   3.8222 0.0001537 ***
ProgramCONT:Time  -0.183971   0.059884 393  -3.0721 0.0022735 ** 
ProgramRI:Time    -0.054953   0.063603 393  -0.8640 0.3881129    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Correlation of Fixed Effects:
            (Intr) PrCONT PrgrRI Time   PCONT:
ProgramCONT -0.698                            
ProgramRI   -0.658  0.459                     
Time        -0.174  0.121  0.114              
PrgrmCONT:T  0.121 -0.174 -0.080 -0.698       
ProgrmRI:Tm  0.114 -0.080 -0.174 -0.658  0.459
> ## Not run: 
> ##D intervals(fm2Weight)
> ##D fm3Weight <- update(fm2Weight, correlation = corAR1())
> ##D anova(fm2Weight, fm3Weight)
> ##D fm4Weight <- update(fm3Weight, strength ~ Program * (Time + I(Time^2)),
> ##D                     random = ~Time|Subj)
> ##D summary(fm4Weight)
> ##D anova(fm4Weight)
> ##D intervals(fm4Weight)
> ## End(Not run)
> 
> 
> 
> options(contrasts = c(unordered = "contr.treatment",ordered = "contr.poly"))
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  22.72 0.37 24 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')