Last updated on 2024-06-14 01:52:34 CEST.
| Package | ERROR | NOTE | OK |
|---|---|---|---|
| gamair | 1 | 12 | |
| gamm4 | 1 | 12 | |
| mgcv | 2 | 11 |
Current CRAN status: ERROR: 1, OK: 12
Version: 1.0-2
Check: package dependencies
Result: NOTE
Package suggested but not available for checking: ‘lme4’
Flavor: r-devel-linux-x86_64-debian-gcc
Version: 1.0-2
Check: examples
Result: ERROR
Running examples in ‘gamair-Ex.R’ failed
The error most likely occurred in:
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### 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.9-1. 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 -1.1102e-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 249.2082 291.6289 -106.6041 1 vs 2 3.253456 0.0713
> 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 249.2082 291.6289 -106.6041 1 vs 2 12.54979 0.0137
> 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 249.2082 291.6289 -106.6041
m5 2 12 293.7081 321.9886 -134.8540 1 vs 2 56.49989 <.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)
Error in library(lme4) : there is no package called ‘lme4’
Execution halted
Flavor: r-devel-linux-x86_64-debian-gcc
Current CRAN status: ERROR: 1, OK: 12
Version: 0.2-6
Check: package dependencies
Result: ERROR
Package required but not available: ‘lme4’
See section ‘The DESCRIPTION file’ in the ‘Writing R Extensions’
manual.
Flavor: r-devel-linux-x86_64-debian-gcc
Current CRAN status: NOTE: 2, OK: 11
Version: 1.9-1
Check: installed package size
Result: NOTE
installed size is 5.3Mb
sub-directories of 1Mb or more:
R 2.0Mb
libs 2.1Mb
Flavors: r-release-macos-x86_64, r-oldrel-macos-x86_64