CRAN Package Check Results for Package gamair

Last updated on 2022-01-27 19:50:59 CET.

Flavor Version Tinstall Tcheck Ttotal Status Flags
r-devel-linux-x86_64-debian-clang 1.0-2 2.03 92.15 94.18 OK
r-devel-linux-x86_64-debian-gcc 1.0-2 1.66 69.92 71.58 OK
r-devel-linux-x86_64-fedora-clang 1.0-2 113.19 OK
r-devel-linux-x86_64-fedora-gcc 1.0-2 103.85 OK
r-devel-windows-x86_64-new-UL 1.0-2 35.00 133.00 168.00 OK
r-devel-windows-x86_64-new-TK 1.0-2 ERROR
r-patched-linux-x86_64 1.0-2 2.75 87.43 90.18 OK
r-release-linux-x86_64 1.0-2 2.22 89.17 91.39 OK
r-release-macos-arm64 1.0-2 OK
r-release-macos-x86_64 1.0-2 OK
r-release-windows-ix86+x86_64 1.0-2 6.00 123.00 129.00 OK
r-oldrel-macos-x86_64 1.0-2 OK
r-oldrel-windows-ix86+x86_64 1.0-2 5.00 121.00 126.00 OK

Check Details

Version: 1.0-2
Check: package dependencies
Result: NOTE
    Package suggested but not available for checking: 'lme4'
Flavor: r-devel-windows-x86_64-new-TK

Version: 1.0-2
Check: examples
Result: ERROR
    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-38. 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)
    Error in library(lme4) : there is no package called 'lme4'
    Execution halted
Flavor: r-devel-windows-x86_64-new-TK