Running examples in 'gamair-Ex.R' failed
The error most likely occurred in:

> ### Name: ch2
> ### Title: Code for Chapter 2: Linear Mixed Models
> ### Aliases: ch2
> 
> ### ** Examples
> 
> library(gamair); library(mgcv)
Loading required package: nlme
This is mgcv 1.8-28. For overview type 'help("mgcv-package")'.
> 
> ## 2.1.1
> data(stomata)
> m1 <- lm(area ~ CO2 + tree,stomata)
> m0 <- lm(area ~ CO2,stomata)
> anova(m0,m1)
Analysis of Variance Table

Model 1: area ~ CO2
Model 2: area ~ CO2 + tree
  Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
1     22 2.1348                                
2     18 0.8604  4    1.2744 6.6654 0.001788 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> m2 <- lm(area ~ tree,stomata)
> anova(m2,m1)
Analysis of Variance Table

Model 1: area ~ tree
Model 2: area ~ CO2 + tree
  Res.Df    RSS Df  Sum of Sq F Pr(>F)
1     18 0.8604                       
2     18 0.8604  0 2.2204e-16         
> st <- aggregate(data.matrix(stomata),
+                 by=list(tree=stomata$tree),mean)
> st$CO2 <- as.factor(st$CO2);st
  tree     area CO2 tree
1    1 1.623374   1    1
2    2 1.598643   1    2
3    3 1.162961   1    3
4    4 2.789238   2    4
5    5 2.903544   2    5
6    6 2.329761   2    6
> m3 <- lm(area~CO2,st)
> anova(m3)
Analysis of Variance Table

Response: area
          Df  Sum Sq Mean Sq F value   Pr(>F)   
CO2        1 2.20531 2.20531  27.687 0.006247 **
Residuals  4 0.31861 0.07965                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> summary(m3)$sigma^2 - summary(m1)$sigma^2/4
[1] 0.06770177
> 
> ## 2.1.3
> library(nlme)  # load nlme `library', which contains data
> data(Rail)     # load data
> Rail
Grouped Data: travel ~ 1 | Rail
   Rail travel
1     1     55
2     1     53
3     1     54
4     2     26
5     2     37
6     2     32
7     3     78
8     3     91
9     3     85
10    4     92
11    4    100
12    4     96
13    5     49
14    5     51
15    5     50
16    6     80
17    6     85
18    6     83
> m1 <- lm(travel ~ Rail,Rail)
> anova(m1)
Analysis of Variance Table

Response: travel
          Df Sum Sq Mean Sq F value    Pr(>F)    
Rail       5 9310.5 1862.10  115.18 1.033e-09 ***
Residuals 12  194.0   16.17                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> rt <- aggregate(data.matrix(Rail),by=list(Rail$Rail),mean)
> rt
  Group.1 Rail   travel
1       2    1 31.66667
2       5    2 50.00000
3       1    3 54.00000
4       6    4 82.66667
5       3    5 84.66667
6       4    6 96.00000
> m0 <- lm(travel ~ 1,rt)   # fit model to aggregated data
> sigb <- (summary(m0)$sigma^2-summary(m1)$sigma^2/3)^0.5
> # sigb^2 is variance component for rail
> sig <- summary(m1)$sigma # sig^2 is resid. var. component
> sigb
[1] 24.80547
> sig
[1] 4.020779
> summary(m0)

Call:
lm(formula = travel ~ 1, data = rt)

Residuals:
     1      2      3      4      5      6 
-34.83 -16.50 -12.50  16.17  18.17  29.50 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)    66.50      10.17   6.538  0.00125 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 24.91 on 5 degrees of freedom

> 
> ## 2.1.4
> library(nlme)
> data(Machines)
> names(Machines)
[1] "Worker"  "Machine" "score"  
> attach(Machines)  # make data available without `Machines$'
> interaction.plot(Machine,Worker,score)
> m1 <- lm(score ~ Worker*Machine,Machines)
> m0 <- lm(score ~ Worker + Machine,Machines)
> anova(m0,m1)
Analysis of Variance Table

Model 1: score ~ Worker + Machine
Model 2: score ~ Worker * Machine
  Res.Df    RSS Df Sum of Sq     F    Pr(>F)    
1     46 459.82                                 
2     36  33.29 10    426.53 46.13 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> summary(m1)$sigma^2
[1] 0.9246296
> Mach <- aggregate(data.matrix(Machines),by=
+         list(Machines$Worker,Machines$Machine),mean)
> Mach$Worker <- as.factor(Mach$Worker)
> Mach$Machine <- as.factor(Mach$Machine)
> m0 <- lm(score ~ Worker + Machine,Mach)
> anova(m0)
Analysis of Variance Table

Response: score
          Df Sum Sq Mean Sq F value    Pr(>F)    
Worker     5 413.96  82.793  5.8232 0.0089495 ** 
Machine    2 585.09 292.544 20.5761 0.0002855 ***
Residuals 10 142.18  14.218                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> summary(m0)$sigma^2 - summary(m1)$sigma^2/3
[1] 13.90946
> M <- aggregate(data.matrix(Mach),by=list(Mach$Worker),mean)
> m00 <- lm(score ~ 1,M)
> summary(m00)$sigma^2 - (summary(m0)$sigma^2)/3
[1] 22.85844
> 
> ## 2.4.4
> llm <- function(theta,X,Z,y) {
+   ## untransform parameters...
+   sigma.b <- exp(theta[1])
+   sigma <- exp(theta[2])
+   ## extract dimensions...
+   n <- length(y); pr <- ncol(Z); pf <- ncol(X)
+   ## obtain \hat \beta, \hat b...
+   X1 <- cbind(X,Z)
+   ipsi <- c(rep(0,pf),rep(1/sigma.b^2,pr))
+   b1 <- solve(crossprod(X1)/sigma^2+diag(ipsi),
+               t(X1)%*%y/sigma^2)
+   ## compute log|Z'Z/sigma^2 + I/sigma.b^2|...
+   ldet <- sum(log(diag(chol(crossprod(Z)/sigma^2 + 
+               diag(ipsi[-(1:pf)])))))
+   ## compute log profile likelihood...
+   l <- (-sum((y-X1%*%b1)^2)/sigma^2 - sum(b1^2*ipsi) - 
+   n*log(sigma^2) - pr*log(sigma.b^2) - 2*ldet - n*log(2*pi))/2
+   attr(l,"b") <- as.numeric(b1) ## return \hat beta and \hat b
+   -l 
+ }
> library(nlme) ## for Rail data
> options(contrasts=c("contr.treatment","contr.treatment"))
> Z <- model.matrix(~Rail$Rail-1) ## r.e. model matrix
> X <- matrix(1,18,1)             ## fixed model matrix
> ## fit the model...
> rail.mod <- optim(c(0,0),llm,hessian=TRUE,
+                            X=X,Z=Z,y=Rail$travel)
> exp(rail.mod$par) ## variance components
[1] 22.629166  4.024072
> solve(rail.mod$hessian) ## approx cov matrix for theta 
              [,1]          [,2]
[1,]  0.0851408546 -0.0004397245
[2,] -0.0004397245  0.0417347933
> attr(llm(rail.mod$par,X,Z,Rail$travel),"b")
[1]  66.50000 -34.46999 -16.32789 -12.36961  15.99803  17.97717  29.19229
> 
> ## 2.5.1
> library(nlme)
> lme(travel~1,Rail,list(Rail=~1))
Linear mixed-effects model fit by REML
  Data: Rail 
  Log-restricted-likelihood: -61.0885
  Fixed: travel ~ 1 
(Intercept) 
       66.5 

Random effects:
 Formula: ~1 | Rail
        (Intercept) Residual
StdDev:    24.80547 4.020779

Number of Observations: 18
Number of Groups: 6 
> 
> ## 2.5.2
> 
> Loblolly$age <- Loblolly$age - mean(Loblolly$age)
> lmc <- lmeControl(niterEM=500,msMaxIter=100)
> m0 <- lme(height ~ age + I(age^2) + I(age^3),Loblolly,
+           random=list(Seed=~age+I(age^2)+I(age^3)),
+           correlation=corAR1(form=~age|Seed),control=lmc)
> plot(m0)
> m1 <- lme(height ~ age+I(age^2)+I(age^3)+I(age^4),Loblolly,
+           list(Seed=~age+I(age^2)+I(age^3)),
+           cor=corAR1(form=~age|Seed),control=lmc)
> plot(m1)
> m2 <- lme(height~age+I(age^2)+I(age^3)+I(age^4)+I(age^5),
+           Loblolly,list(Seed=~age+I(age^2)+I(age^3)),
+           cor=corAR1(form=~age|Seed),control=lmc)
> plot(m2)
> plot(m2,Seed~resid(.))
> qqnorm(m2,~resid(.))
> qqnorm(m2,~ranef(.))
> 
> m3 <- lme(height~age+I(age^2)+I(age^3)+I(age^4)+I(age^5),
+           Loblolly,list(Seed=~age+I(age^2)+I(age^3)),control=lmc)
> anova(m3,m2)
   Model df      AIC      BIC    logLik   Test  L.Ratio p-value
m3     1 17 250.4616 290.5257 -108.2308                        
m2     2 18 239.3576 281.7783 -101.6788 1 vs 2 13.10407   3e-04
> m4 <- lme(height~age+I(age^2)+I(age^3)+I(age^4)+I(age^5),
+           Loblolly,list(Seed=~age+I(age^2)),
+           correlation=corAR1(form=~age|Seed),control=lmc)
> anova(m4,m2)
   Model df      AIC      BIC    logLik   Test L.Ratio p-value
m4     1 14 253.7579 286.7519 -112.8790                       
m2     2 18 239.3576 281.7783 -101.6788 1 vs 2 22.4004   2e-04
> m5 <- lme(height~age+I(age^2)+I(age^3)+I(age^4)+I(age^5),
+           Loblolly,list(Seed=pdDiag(~age+I(age^2)+I(age^3))),
+           correlation=corAR1(form=~age|Seed),control=lmc)
> anova(m2,m5)
   Model df      AIC      BIC    logLik   Test L.Ratio p-value
m2     1 18 239.3576 281.7783 -101.6788                       
m5     2 12 293.7081 321.9886 -134.8540 1 vs 2 66.3505  <.0001
> plot(augPred(m2))
> 
> ## 2.5.3
> lme(score~Machine,Machines,list(Worker=~1,Machine=~1))
Linear mixed-effects model fit by REML
  Data: Machines 
  Log-restricted-likelihood: -107.8438
  Fixed: score ~ Machine 
(Intercept)    MachineB    MachineC 
  52.355556    7.966667   13.916667 

Random effects:
 Formula: ~1 | Worker
        (Intercept)
StdDev:     4.78105

 Formula: ~1 | Machine %in% Worker
        (Intercept)  Residual
StdDev:    3.729532 0.9615771

Number of Observations: 54
Number of Groups: 
             Worker Machine %in% Worker 
                  6                  18 
> 
> ## 2.5.4
> library(lme4)
Loading required package: Matrix
Error: package or namespace load failed for 'lme4' in library.dynam(lib, package, package.lib):
 DLL 'Rcpp' not found: maybe not installed for this architecture?
Execution halted