R Under development (unstable) (2025-02-20 r87772) -- "Unsuffered Consequences"
Copyright (C) 2025 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu

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 an HTML browser interface to help.
Type 'q()' to quit R.

> pkgname <- "datasets"
> source(file.path(R.home("share"), "R", "examples-header.R"))
> options(warn = 1)
> library('datasets')
> 
> base::assign(".oldSearch", base::search(), pos = 'CheckExEnv')
> base::assign(".old_wd", base::getwd(), pos = 'CheckExEnv')
> cleanEx()
> nameEx("AirPassengers")
> ### * AirPassengers
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: AirPassengers
> ### Title: Monthly Airline Passenger Numbers 1949-1960
> ### Aliases: AirPassengers
> ### Keywords: datasets
> 
> ### ** Examples
> 
> ## The classic 'airline model', by full ML
> (fit <- arima(log10(AirPassengers), c(0, 1, 1),
+               seasonal = list(order = c(0, 1, 1), period = 12)))

Call:
arima(x = log10(AirPassengers), order = c(0, 1, 1), seasonal = list(order = c(0, 
    1, 1), period = 12))

Coefficients:
          ma1     sma1
      -0.4018  -0.5569
s.e.   0.0896   0.0731

sigma^2 estimated as 0.0002543:  log likelihood = 353.96,  aic = -701.92
> update(fit, method = "CSS")

Call:
arima(x = log10(AirPassengers), order = c(0, 1, 1), seasonal = list(order = c(0, 
    1, 1), period = 12), method = "CSS")

Coefficients:
          ma1     sma1
      -0.3772  -0.5724
s.e.   0.0883   0.0704

sigma^2 estimated as 0.0002619:  part log likelihood = 354.32
> update(fit, x = window(log10(AirPassengers), start = 1954))

Call:
arima(x = window(log10(AirPassengers), start = 1954), order = c(0, 1, 1), seasonal = list(order = c(0, 
    1, 1), period = 12))

Coefficients:
          ma1     sma1
      -0.4797  -0.4460
s.e.   0.1000   0.1514

sigma^2 estimated as 0.0001603:  log likelihood = 208.02,  aic = -410.04
> pred <- predict(fit, n.ahead = 24)
> tl <- pred$pred - 1.96 * pred$se
> tu <- pred$pred + 1.96 * pred$se
> ts.plot(AirPassengers, 10^tl, 10^tu, log = "y", lty = c(1, 2, 2))
> 
> ## full ML fit is the same if the series is reversed, CSS fit is not
> ap0 <- rev(log10(AirPassengers))
> attributes(ap0) <- attributes(AirPassengers)
> arima(ap0, c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))

Call:
arima(x = ap0, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))

Coefficients:
          ma1     sma1
      -0.4018  -0.5569
s.e.   0.0896   0.0731

sigma^2 estimated as 0.0002543:  log likelihood = 353.95,  aic = -701.91
> arima(ap0, c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12),
+       method = "CSS")

Call:
arima(x = ap0, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12), 
    method = "CSS")

Coefficients:
          ma1     sma1
      -0.4188  -0.5224
s.e.   0.0875   0.0695

sigma^2 estimated as 0.0002673:  part log likelihood = 353
> 
> ## Structural Time Series
> ap <- log10(AirPassengers) - 2
> (fit <- StructTS(ap, type = "BSM"))

Call:
StructTS(x = ap, type = "BSM")

Variances:
    level      slope       seas    epsilon  
0.0001456  0.0000000  0.0002635  0.0000000  
> par(mfrow = c(1, 2))
> plot(cbind(ap, fitted(fit)), plot.type = "single")
> plot(cbind(ap, tsSmooth(fit)), plot.type = "single")
> 
> 
> 
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("BOD")
> ### * BOD
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: BOD
> ### Title: Biochemical Oxygen Demand
> ### Aliases: BOD
> ### Keywords: datasets
> 
> ### ** Examples
> 
> ## Don't show: 
> options(show.nls.convergence=FALSE)
> old <- options(digits = 5)
> ## End(Don't show)
> require(stats)
> # simplest form of fitting a first-order model to these data
> fm1 <- nls(demand ~ A*(1-exp(-exp(lrc)*Time)), data = BOD,
+            start = c(A = 20, lrc = log(.35)))
> coef(fm1)
       A      lrc 
19.14258 -0.63282 
> fm1
Nonlinear regression model
  model: demand ~ A * (1 - exp(-exp(lrc) * Time))
   data: BOD
     A    lrc 
19.143 -0.633 
 residual sum-of-squares: 26
> ## IGNORE_RDIFF_BEGIN
> # using the plinear algorithm  (trace o/p differs by platform)
> fm2 <- nls(demand ~ (1-exp(-exp(lrc)*Time)), data = BOD,
+            start = c(lrc = log(.35)), algorithm = "plinear", trace = TRUE)
32.946    (7.09e-01): par = (-1.0498 22.126)
25.992    (1.29e-02): par = (-0.62572 19.103)
25.990    (2.13e-04): par = (-0.6327 19.142)
25.990    (3.55e-06): par = (-0.63282 19.143)
> ## IGNORE_RDIFF_END
> # using a self-starting model
> fm3 <- nls(demand ~ SSasympOrig(Time, A, lrc), data = BOD)
> summary(fm3)

Formula: demand ~ SSasympOrig(Time, A, lrc)

Parameters:
    Estimate Std. Error t value Pr(>|t|)   
A     19.143      2.496    7.67   0.0016 **
lrc   -0.633      0.382   -1.65   0.1733   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.55 on 4 degrees of freedom

> ## Don't show: 
> options(old)
> ## End(Don't show)
> 
> 
> 
> cleanEx()
> nameEx("ChickWeight")
> ### * ChickWeight
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ChickWeight
> ### Title: Weight versus age of chicks on different diets
> ### Aliases: ChickWeight
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> coplot(weight ~ Time | Chick, data = ChickWeight,
+        type = "b", show.given = FALSE)
> 
> 
> 
> cleanEx()
> nameEx("DNase")
> ### * DNase
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: DNase
> ### Title: Elisa assay of DNase
> ### Aliases: DNase
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> ## Don't show: 
> options(show.nls.convergence=FALSE)
> ## End(Don't show)
> coplot(density ~ conc | Run, data = DNase,
+        show.given = FALSE, type = "b")
> coplot(density ~ log(conc) | Run, data = DNase,
+        show.given = FALSE, type = "b")
> ## fit a representative run
> fm1 <- nls(density ~ SSlogis( log(conc), Asym, xmid, scal ),
+     data = DNase, subset = Run == 1)
> ## compare with a four-parameter logistic
> fm2 <- nls(density ~ SSfpl( log(conc), A, B, xmid, scal ),
+     data = DNase, subset = Run == 1)
> summary(fm2)

Formula: density ~ SSfpl(log(conc), A, B, xmid, scal)

Parameters:
      Estimate Std. Error t value Pr(>|t|)    
A    -0.007897   0.017200  -0.459    0.654    
B     2.377239   0.109516  21.707 5.35e-11 ***
xmid  1.507403   0.102080  14.767 4.65e-09 ***
scal  1.062579   0.056996  18.643 3.16e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01981 on 12 degrees of freedom

> anova(fm1, fm2)
Analysis of Variance Table

Model 1: density ~ SSlogis(log(conc), Asym, xmid, scal)
Model 2: density ~ SSfpl(log(conc), A, B, xmid, scal)
  Res.Df Res.Sum Sq Df     Sum Sq F value Pr(>F)
1     13  0.0047896                             
2     12  0.0047073  1 8.2314e-05  0.2098 0.6551
> 
> 
> 
> cleanEx()
> nameEx("Formaldehyde")
> ### * Formaldehyde
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Formaldehyde
> ### Title: Determination of Formaldehyde
> ### Aliases: Formaldehyde
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> plot(optden ~ carb, data = Formaldehyde,
+      xlab = "Carbohydrate (ml)", ylab = "Optical Density",
+      main = "Formaldehyde data", col = 4, las = 1)
> abline(fm1 <- lm(optden ~ carb, data = Formaldehyde))
> summary(fm1)

Call:
lm(formula = optden ~ carb, data = Formaldehyde)

Residuals:
        1         2         3         4         5         6 
-0.006714  0.001029  0.002771  0.007143  0.007514 -0.011743 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.005086   0.007834   0.649    0.552    
carb        0.876286   0.013535  64.744 3.41e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.008649 on 4 degrees of freedom
Multiple R-squared:  0.999,	Adjusted R-squared:  0.9988 
F-statistic:  4192 on 1 and 4 DF,  p-value: 3.409e-07

> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0))
> plot(fm1)
> par(opar)
> 
> 
> 
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("HairEyeColor")
> ### * HairEyeColor
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: HairEyeColor
> ### Title: Hair and Eye Color of Statistics Students
> ### Aliases: HairEyeColor
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> ## Full mosaic
> mosaicplot(HairEyeColor)
> ## Aggregate over sex (as in Snee's original data)
> x <- apply(HairEyeColor, c(1, 2), sum)
> x
       Eye
Hair    Brown Blue Hazel Green
  Black    68   20    15     5
  Brown   119   84    54    29
  Red      26   17    14    14
  Blond     7   94    10    16
> mosaicplot(x, main = "Relation between hair and eye color")
> 
> 
> 
> cleanEx()
> nameEx("Harman23.cor")
> ### * Harman23.cor
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Harman23.cor
> ### Title: Harman Example 2.3
> ### Aliases: Harman23.cor
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats)
> (Harman23.FA <- factanal(factors = 1, covmat = Harman23.cor))

Call:
factanal(factors = 1, covmat = Harman23.cor)

Uniquenesses:
        height       arm.span        forearm      lower.leg         weight 
         0.158          0.135          0.190          0.187          0.760 
bitro.diameter    chest.girth    chest.width 
         0.829          0.877          0.801 

Loadings:
               Factor1
height         0.918  
arm.span       0.930  
forearm        0.900  
lower.leg      0.902  
weight         0.490  
bitro.diameter 0.413  
chest.girth    0.351  
chest.width    0.446  

               Factor1
SS loadings      4.064
Proportion Var   0.508

Test of the hypothesis that 1 factor is sufficient.
The chi square statistic is 611.44 on 20 degrees of freedom.
The p-value is 1.12e-116 
> for(factors in 2:4) print(update(Harman23.FA, factors = factors))

Call:
factanal(factors = factors, covmat = Harman23.cor)

Uniquenesses:
        height       arm.span        forearm      lower.leg         weight 
         0.170          0.107          0.166          0.199          0.089 
bitro.diameter    chest.girth    chest.width 
         0.364          0.416          0.537 

Loadings:
               Factor1 Factor2
height         0.865   0.287  
arm.span       0.927   0.181  
forearm        0.895   0.179  
lower.leg      0.859   0.252  
weight         0.233   0.925  
bitro.diameter 0.194   0.774  
chest.girth    0.134   0.752  
chest.width    0.278   0.621  

               Factor1 Factor2
SS loadings      3.335   2.617
Proportion Var   0.417   0.327
Cumulative Var   0.417   0.744

Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 75.74 on 13 degrees of freedom.
The p-value is 6.94e-11 

Call:
factanal(factors = factors, covmat = Harman23.cor)

Uniquenesses:
        height       arm.span        forearm      lower.leg         weight 
         0.127          0.005          0.193          0.157          0.090 
bitro.diameter    chest.girth    chest.width 
         0.359          0.411          0.490 

Loadings:
               Factor1 Factor2 Factor3
height          0.886   0.267  -0.130 
arm.span        0.937   0.195   0.280 
forearm         0.874   0.188         
lower.leg       0.877   0.230  -0.145 
weight          0.242   0.916  -0.106 
bitro.diameter  0.193   0.777         
chest.girth     0.137   0.755         
chest.width     0.261   0.646   0.159 

               Factor1 Factor2 Factor3
SS loadings      3.379   2.628   0.162
Proportion Var   0.422   0.329   0.020
Cumulative Var   0.422   0.751   0.771

Test of the hypothesis that 3 factors are sufficient.
The chi square statistic is 22.81 on 7 degrees of freedom.
The p-value is 0.00184 

Call:
factanal(factors = factors, covmat = Harman23.cor)

Uniquenesses:
        height       arm.span        forearm      lower.leg         weight 
         0.137          0.005          0.191          0.116          0.138 
bitro.diameter    chest.girth    chest.width 
         0.283          0.178          0.488 

Loadings:
               Factor1 Factor2 Factor3 Factor4
height          0.879   0.277          -0.115 
arm.span        0.937   0.194           0.277 
forearm         0.875   0.191                 
lower.leg       0.887   0.209   0.135  -0.188 
weight          0.246   0.882   0.111  -0.109 
bitro.diameter  0.187   0.822                 
chest.girth     0.117   0.729   0.526         
chest.width     0.263   0.644           0.141 

               Factor1 Factor2 Factor3 Factor4
SS loadings      3.382   2.595   0.323   0.165
Proportion Var   0.423   0.324   0.040   0.021
Cumulative Var   0.423   0.747   0.787   0.808

Test of the hypothesis that 4 factors are sufficient.
The chi square statistic is 4.63 on 2 degrees of freedom.
The p-value is 0.0988 
> 
> 
> 
> cleanEx()
> nameEx("Harman74.cor")
> ### * Harman74.cor
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Harman74.cor
> ### Title: Harman Example 7.4
> ### Aliases: Harman74.cor
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats)
> (Harman74.FA <- factanal(factors = 1, covmat = Harman74.cor))

Call:
factanal(factors = 1, covmat = Harman74.cor)

Uniquenesses:
      VisualPerception                  Cubes         PaperFormBoard 
                 0.677                  0.866                  0.830 
                 Flags     GeneralInformation  PargraphComprehension 
                 0.768                  0.487                  0.491 
    SentenceCompletion     WordClassification            WordMeaning 
                 0.500                  0.514                  0.474 
              Addition                   Code           CountingDots 
                 0.818                  0.731                  0.824 
StraightCurvedCapitals        WordRecognition      NumberRecognition 
                 0.681                  0.833                  0.863 
     FigureRecognition           ObjectNumber           NumberFigure 
                 0.775                  0.812                  0.778 
            FigureWord              Deduction       NumericalPuzzles 
                 0.816                  0.612                  0.676 
      ProblemReasoning       SeriesCompletion     ArithmeticProblems 
                 0.619                  0.524                  0.593 

Loadings:
                       Factor1
VisualPerception       0.569  
Cubes                  0.366  
PaperFormBoard         0.412  
Flags                  0.482  
GeneralInformation     0.716  
PargraphComprehension  0.713  
SentenceCompletion     0.707  
WordClassification     0.697  
WordMeaning            0.725  
Addition               0.426  
Code                   0.519  
CountingDots           0.419  
StraightCurvedCapitals 0.565  
WordRecognition        0.408  
NumberRecognition      0.370  
FigureRecognition      0.474  
ObjectNumber           0.434  
NumberFigure           0.471  
FigureWord             0.429  
Deduction              0.623  
NumericalPuzzles       0.569  
ProblemReasoning       0.617  
SeriesCompletion       0.690  
ArithmeticProblems     0.638  

               Factor1
SS loadings      7.438
Proportion Var   0.310

Test of the hypothesis that 1 factor is sufficient.
The chi square statistic is 622.91 on 252 degrees of freedom.
The p-value is 2.28e-33 
> for(factors in 2:5) print(update(Harman74.FA, factors = factors))

Call:
factanal(factors = factors, covmat = Harman74.cor)

Uniquenesses:
      VisualPerception                  Cubes         PaperFormBoard 
                 0.650                  0.864                  0.844 
                 Flags     GeneralInformation  PargraphComprehension 
                 0.778                  0.375                  0.316 
    SentenceCompletion     WordClassification            WordMeaning 
                 0.319                  0.503                  0.258 
              Addition                   Code           CountingDots 
                 0.670                  0.608                  0.581 
StraightCurvedCapitals        WordRecognition      NumberRecognition 
                 0.567                  0.832                  0.850 
     FigureRecognition           ObjectNumber           NumberFigure 
                 0.743                  0.770                  0.625 
            FigureWord              Deduction       NumericalPuzzles 
                 0.792                  0.629                  0.579 
      ProblemReasoning       SeriesCompletion     ArithmeticProblems 
                 0.634                  0.539                  0.553 

Loadings:
                       Factor1 Factor2
VisualPerception       0.506   0.306  
Cubes                  0.304   0.209  
PaperFormBoard         0.297   0.260  
Flags                  0.327   0.339  
GeneralInformation     0.240   0.753  
PargraphComprehension  0.171   0.809  
SentenceCompletion     0.163   0.809  
WordClassification     0.344   0.615  
WordMeaning            0.148   0.849  
Addition               0.563   0.115  
Code                   0.591   0.207  
CountingDots           0.647          
StraightCurvedCapitals 0.612   0.241  
WordRecognition        0.315   0.263  
NumberRecognition      0.328   0.205  
FigureRecognition      0.457   0.218  
ObjectNumber           0.431   0.209  
NumberFigure           0.601   0.116  
FigureWord             0.399   0.222  
Deduction              0.379   0.477  
NumericalPuzzles       0.604   0.237  
ProblemReasoning       0.390   0.462  
SeriesCompletion       0.486   0.474  
ArithmeticProblems     0.544   0.389  

               Factor1 Factor2
SS loadings      4.573   4.548
Proportion Var   0.191   0.190
Cumulative Var   0.191   0.380

Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 420.24 on 229 degrees of freedom.
The p-value is 2.01e-13 

Call:
factanal(factors = factors, covmat = Harman74.cor)

Uniquenesses:
      VisualPerception                  Cubes         PaperFormBoard 
                 0.500                  0.793                  0.662 
                 Flags     GeneralInformation  PargraphComprehension 
                 0.694                  0.352                  0.316 
    SentenceCompletion     WordClassification            WordMeaning 
                 0.300                  0.502                  0.256 
              Addition                   Code           CountingDots 
                 0.200                  0.586                  0.494 
StraightCurvedCapitals        WordRecognition      NumberRecognition 
                 0.569                  0.838                  0.848 
     FigureRecognition           ObjectNumber           NumberFigure 
                 0.643                  0.780                  0.635 
            FigureWord              Deduction       NumericalPuzzles 
                 0.788                  0.590                  0.580 
      ProblemReasoning       SeriesCompletion     ArithmeticProblems 
                 0.597                  0.498                  0.500 

Loadings:
                       Factor1 Factor2 Factor3
VisualPerception        0.176   0.656   0.198 
Cubes                   0.122   0.428         
PaperFormBoard          0.145   0.563         
Flags                   0.239   0.487   0.107 
GeneralInformation      0.745   0.191   0.237 
PargraphComprehension   0.780   0.249   0.118 
SentenceCompletion      0.802   0.175   0.160 
WordClassification      0.571   0.327   0.256 
WordMeaning             0.821   0.248         
Addition                0.162  -0.118   0.871 
Code                    0.198   0.219   0.572 
CountingDots                    0.179   0.688 
StraightCurvedCapitals  0.190   0.381   0.499 
WordRecognition         0.231   0.253   0.210 
NumberRecognition       0.158   0.299   0.195 
FigureRecognition       0.108   0.557   0.186 
ObjectNumber            0.178   0.267   0.342 
NumberFigure                    0.427   0.424 
FigureWord              0.167   0.355   0.240 
Deduction               0.392   0.472   0.181 
NumericalPuzzles        0.178   0.406   0.473 
ProblemReasoning        0.382   0.473   0.182 
SeriesCompletion        0.379   0.528   0.283 
ArithmeticProblems      0.377   0.226   0.554 

               Factor1 Factor2 Factor3
SS loadings      3.802   3.488   3.186
Proportion Var   0.158   0.145   0.133
Cumulative Var   0.158   0.304   0.436

Test of the hypothesis that 3 factors are sufficient.
The chi square statistic is 295.59 on 207 degrees of freedom.
The p-value is 5.12e-05 

Call:
factanal(factors = factors, covmat = Harman74.cor)

Uniquenesses:
      VisualPerception                  Cubes         PaperFormBoard 
                 0.438                  0.780                  0.644 
                 Flags     GeneralInformation  PargraphComprehension 
                 0.651                  0.352                  0.312 
    SentenceCompletion     WordClassification            WordMeaning 
                 0.283                  0.485                  0.257 
              Addition                   Code           CountingDots 
                 0.240                  0.551                  0.435 
StraightCurvedCapitals        WordRecognition      NumberRecognition 
                 0.491                  0.646                  0.696 
     FigureRecognition           ObjectNumber           NumberFigure 
                 0.549                  0.598                  0.593 
            FigureWord              Deduction       NumericalPuzzles 
                 0.762                  0.592                  0.583 
      ProblemReasoning       SeriesCompletion     ArithmeticProblems 
                 0.601                  0.497                  0.500 

Loadings:
                       Factor1 Factor2 Factor3 Factor4
VisualPerception        0.160   0.689   0.187   0.160 
Cubes                   0.117   0.436                 
PaperFormBoard          0.137   0.570           0.110 
Flags                   0.233   0.527                 
GeneralInformation      0.739   0.185   0.213   0.150 
PargraphComprehension   0.767   0.205           0.233 
SentenceCompletion      0.806   0.197   0.153         
WordClassification      0.569   0.339   0.242   0.132 
WordMeaning             0.806   0.201           0.227 
Addition                0.167  -0.118   0.831   0.166 
Code                    0.180   0.120   0.512   0.374 
CountingDots                    0.210   0.716         
StraightCurvedCapitals  0.188   0.438   0.525         
WordRecognition         0.197                   0.553 
NumberRecognition       0.122   0.116           0.520 
FigureRecognition               0.408           0.525 
ObjectNumber            0.142           0.219   0.574 
NumberFigure                    0.293   0.336   0.456 
FigureWord              0.148   0.239   0.161   0.365 
Deduction               0.378   0.402   0.118   0.301 
NumericalPuzzles        0.175   0.381   0.438   0.223 
ProblemReasoning        0.366   0.399   0.123   0.301 
SeriesCompletion        0.369   0.500   0.244   0.239 
ArithmeticProblems      0.370   0.158   0.496   0.304 

               Factor1 Factor2 Factor3 Factor4
SS loadings      3.647   2.872   2.657   2.290
Proportion Var   0.152   0.120   0.111   0.095
Cumulative Var   0.152   0.272   0.382   0.478

Test of the hypothesis that 4 factors are sufficient.
The chi square statistic is 226.68 on 186 degrees of freedom.
The p-value is 0.0224 

Call:
factanal(factors = factors, covmat = Harman74.cor)

Uniquenesses:
      VisualPerception                  Cubes         PaperFormBoard 
                 0.450                  0.781                  0.639 
                 Flags     GeneralInformation  PargraphComprehension 
                 0.649                  0.357                  0.288 
    SentenceCompletion     WordClassification            WordMeaning 
                 0.277                  0.485                  0.262 
              Addition                   Code           CountingDots 
                 0.215                  0.386                  0.444 
StraightCurvedCapitals        WordRecognition      NumberRecognition 
                 0.256                  0.639                  0.706 
     FigureRecognition           ObjectNumber           NumberFigure 
                 0.550                  0.614                  0.596 
            FigureWord              Deduction       NumericalPuzzles 
                 0.764                  0.521                  0.564 
      ProblemReasoning       SeriesCompletion     ArithmeticProblems 
                 0.580                  0.442                  0.478 

Loadings:
                       Factor1 Factor2 Factor3 Factor4 Factor5
VisualPerception        0.161   0.658   0.136   0.182   0.199 
Cubes                   0.113   0.435           0.107         
PaperFormBoard          0.135   0.562           0.107   0.116 
Flags                   0.231   0.533                         
GeneralInformation      0.736   0.188   0.192   0.162         
PargraphComprehension   0.775   0.187           0.251   0.113 
SentenceCompletion      0.809   0.208   0.136                 
WordClassification      0.568   0.348   0.223   0.131         
WordMeaning             0.800   0.215           0.224         
Addition                0.175  -0.100   0.844   0.176         
Code                    0.185           0.438   0.451   0.426 
CountingDots                    0.222   0.690   0.101   0.140 
StraightCurvedCapitals  0.186   0.425   0.458           0.559 
WordRecognition         0.197                   0.557         
NumberRecognition       0.121   0.130           0.508         
FigureRecognition               0.400           0.529         
ObjectNumber            0.145           0.208   0.562         
NumberFigure                    0.306   0.325   0.452         
FigureWord              0.147   0.242   0.145   0.364         
Deduction               0.370   0.452   0.139   0.287  -0.190 
NumericalPuzzles        0.170   0.402   0.439   0.230         
ProblemReasoning        0.358   0.423   0.126   0.302         
SeriesCompletion        0.360   0.549   0.256   0.223  -0.107 
ArithmeticProblems      0.371   0.185   0.502   0.307         

               Factor1 Factor2 Factor3 Factor4 Factor5
SS loadings      3.632   2.964   2.456   2.345   0.663
Proportion Var   0.151   0.124   0.102   0.098   0.028
Cumulative Var   0.151   0.275   0.377   0.475   0.503

Test of the hypothesis that 5 factors are sufficient.
The chi square statistic is 186.82 on 166 degrees of freedom.
The p-value is 0.128 
> Harman74.FA <- factanal(factors = 5, covmat = Harman74.cor,
+                         rotation = "promax")
> print(Harman74.FA$loadings, sort = TRUE)

Loadings:
                       Factor1 Factor2 Factor3 Factor4 Factor5
VisualPerception        0.831          -0.127           0.230 
Cubes                   0.534                                 
PaperFormBoard          0.736          -0.290           0.136 
Flags                   0.647                  -0.104         
SeriesCompletion        0.555   0.126   0.127                 
GeneralInformation              0.764                         
PargraphComprehension           0.845  -0.140   0.140         
SentenceCompletion              0.872          -0.140         
WordClassification      0.277   0.505   0.104                 
WordMeaning                     0.846  -0.108                 
Addition               -0.334           1.012                 
CountingDots            0.206  -0.200   0.722           0.185 
ArithmeticProblems              0.197   0.500   0.139         
WordRecognition        -0.126   0.127  -0.103   0.657         
NumberRecognition                               0.568         
FigureRecognition       0.399  -0.142  -0.207   0.562         
ObjectNumber           -0.108           0.107   0.613         
StraightCurvedCapitals  0.542           0.247           0.618 
Code                            0.112   0.288   0.486   0.424 
NumberFigure            0.255  -0.230   0.211   0.413         
FigureWord              0.187                   0.347         
Deduction               0.404   0.169           0.117  -0.203 
NumericalPuzzles        0.393           0.368                 
ProblemReasoning        0.381   0.188           0.169         

               Factor1 Factor2 Factor3 Factor4 Factor5
SS loadings      3.529   3.311   2.367   2.109   0.762
Proportion Var   0.147   0.138   0.099   0.088   0.032
Cumulative Var   0.147   0.285   0.384   0.471   0.503
> 
> 
> 
> cleanEx()
> nameEx("InsectSprays")
> ### * InsectSprays
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: InsectSprays
> ### Title: Effectiveness of Insect Sprays
> ### Aliases: InsectSprays
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> boxplot(count ~ spray, data = InsectSprays,
+         xlab = "Type of spray", ylab = "Insect count",
+         main = "InsectSprays data", varwidth = TRUE, col = "lightgray")
> fm1 <- aov(count ~ spray, data = InsectSprays)
> summary(fm1)
            Df Sum Sq Mean Sq F value Pr(>F)    
spray        5   2669   533.8    34.7 <2e-16 ***
Residuals   66   1015    15.4                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0))
> plot(fm1)
> fm2 <- aov(sqrt(count) ~ spray, data = InsectSprays)
> summary(fm2)
            Df Sum Sq Mean Sq F value Pr(>F)    
spray        5  88.44  17.688    44.8 <2e-16 ***
Residuals   66  26.06   0.395                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> plot(fm2)
> par(opar)
> 
> 
> 
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("JohnsonJohnson")
> ### * JohnsonJohnson
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: JohnsonJohnson
> ### Title: Quarterly Earnings per Johnson & Johnson Share
> ### Aliases: JohnsonJohnson
> ### Keywords: datasets
> 
> ### ** Examples
> 
> 
> cleanEx()
> nameEx("LifeCycleSavings")
> ### * LifeCycleSavings
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: LifeCycleSavings
> ### Title: Intercountry Life-Cycle Savings Data
> ### Aliases: LifeCycleSavings
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> pairs(LifeCycleSavings, panel = panel.smooth,
+       main = "LifeCycleSavings data")
> fm1 <- lm(sr ~ pop15 + pop75 + dpi + ddpi, data = LifeCycleSavings)
> summary(fm1)

Call:
lm(formula = sr ~ pop15 + pop75 + dpi + ddpi, data = LifeCycleSavings)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.2422 -2.6857 -0.2488  2.4280  9.7509 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 28.5660865  7.3545161   3.884 0.000334 ***
pop15       -0.4611931  0.1446422  -3.189 0.002603 ** 
pop75       -1.6914977  1.0835989  -1.561 0.125530    
dpi         -0.0003369  0.0009311  -0.362 0.719173    
ddpi         0.4096949  0.1961971   2.088 0.042471 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.803 on 45 degrees of freedom
Multiple R-squared:  0.3385,	Adjusted R-squared:  0.2797 
F-statistic: 5.756 on 4 and 45 DF,  p-value: 0.0007904

> 
> 
> 
> cleanEx()
> nameEx("Loblolly")
> ### * Loblolly
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Loblolly
> ### Title: Growth of Loblolly Pine Trees
> ### Aliases: Loblolly
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> plot(height ~ age, data = Loblolly, subset = Seed == 329,
+      xlab = "Tree age (yr)", las = 1,
+      ylab = "Tree height (ft)",
+      main = "Loblolly data and fitted curve (Seed 329 only)")
> fm1 <- nls(height ~ SSasymp(age, Asym, R0, lrc),
+            data = Loblolly, subset = Seed == 329)
> age <- seq(0, 30, length.out = 101)
> lines(age, predict(fm1, list(age = age)))
> 
> 
> 
> cleanEx()
> nameEx("Nile")
> ### * Nile
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Nile
> ### Title: Flow of the River Nile
> ### Aliases: Nile
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> par(mfrow = c(2, 2))
> plot(Nile)
> acf(Nile)
> pacf(Nile)
> ar(Nile) # selects order 2

Call:
ar(x = Nile)

Coefficients:
     1       2  
0.4081  0.1812  

Order selected 2  sigma^2 estimated as  21247
> cpgram(ar(Nile)$resid)
> par(mfrow = c(1, 1))
> arima(Nile, c(2, 0, 0))

Call:
arima(x = Nile, order = c(2, 0, 0))

Coefficients:
         ar1     ar2  intercept
      0.4096  0.1987   919.8397
s.e.  0.0974  0.0990    35.6410

sigma^2 estimated as 20291:  log likelihood = -637.98,  aic = 1283.96
> 
> ## Now consider missing values, following Durbin & Koopman
> NileNA <- Nile
> NileNA[c(21:40, 61:80)] <- NA
> arima(NileNA, c(2, 0, 0))

Call:
arima(x = NileNA, order = c(2, 0, 0))

Coefficients:
         ar1     ar2  intercept
      0.3622  0.1678   918.3103
s.e.  0.1273  0.1323    39.5037

sigma^2 estimated as 23676:  log likelihood = -387.7,  aic = 783.41
> plot(NileNA)
> pred <-
+    predict(arima(window(NileNA, 1871, 1890), c(2, 0, 0)), n.ahead = 20)
> lines(pred$pred, lty = 3, col = "red")
> lines(pred$pred + 2*pred$se, lty = 2, col = "blue")
> lines(pred$pred - 2*pred$se, lty = 2, col = "blue")
> pred <-
+    predict(arima(window(NileNA, 1871, 1930), c(2, 0, 0)), n.ahead = 20)
> lines(pred$pred, lty = 3, col = "red")
> lines(pred$pred + 2*pred$se, lty = 2, col = "blue")
> lines(pred$pred - 2*pred$se, lty = 2, col = "blue")
> 
> ## Structural time series models
> par(mfrow = c(3, 1))
> plot(Nile)
> ## local level model
> (fit <- StructTS(Nile, type = "level"))

Call:
StructTS(x = Nile, type = "level")

Variances:
  level  epsilon  
   1469    15099  
> lines(fitted(fit), lty = 2)              # contemporaneous smoothing
> lines(tsSmooth(fit), lty = 2, col = 4)   # fixed-interval smoothing
> plot(residuals(fit)); abline(h = 0, lty = 3)
> ## local trend model
> (fit2 <- StructTS(Nile, type = "trend")) ## constant trend fitted

Call:
StructTS(x = Nile, type = "trend")

Variances:
  level    slope  epsilon  
   1427        0    15047  
> pred <- predict(fit, n.ahead = 30)
> ## with 50% confidence interval
> ts.plot(Nile, pred$pred,
+         pred$pred + 0.67*pred$se, pred$pred -0.67*pred$se)
> 
> ## Now consider missing values
> plot(NileNA)
> (fit3 <- StructTS(NileNA, type = "level"))

Call:
StructTS(x = NileNA, type = "level")

Variances:
  level  epsilon  
  685.8  17899.8  
> lines(fitted(fit3), lty = 2)
> lines(tsSmooth(fit3), lty = 3)
> plot(residuals(fit3)); abline(h = 0, lty = 3)
> 
> 
> 
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("Orange")
> ### * Orange
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Orange
> ### Title: Growth of Orange Trees
> ### Aliases: Orange
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> coplot(circumference ~ age | Tree, data = Orange, show.given = FALSE)
> fm1 <- nls(circumference ~ SSlogis(age, Asym, xmid, scal),
+            data = Orange, subset = Tree == 3)
> plot(circumference ~ age, data = Orange, subset = Tree == 3,
+      xlab = "Tree age (days since 1968/12/31)",
+      ylab = "Tree circumference (mm)", las = 1,
+      main = "Orange tree data and fitted model (Tree 3 only)")
> age <- seq(0, 1600, length.out = 101)
> lines(age, predict(fm1, list(age = age)))
> 
> 
> 
> cleanEx()
> nameEx("OrchardSprays")
> ### * OrchardSprays
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: OrchardSprays
> ### Title: Potency of Orchard Sprays
> ### Aliases: OrchardSprays
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> pairs(OrchardSprays, main = "OrchardSprays data")
> 
> 
> 
> cleanEx()
> nameEx("PlantGrowth")
> ### * PlantGrowth
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: PlantGrowth
> ### Title: Results from an Experiment on Plant Growth
> ### Aliases: PlantGrowth
> ### Keywords: datasets
> 
> ### ** Examples
> 
> ## One factor ANOVA example from Dobson's book, cf. Table 7.4:
> require(stats); require(graphics)
> boxplot(weight ~ group, data = PlantGrowth, main = "PlantGrowth data",
+         ylab = "Dried weight of plants", col = "lightgray",
+         notch = TRUE, varwidth = TRUE)
Warning in (function (z, notch = FALSE, width = NULL, varwidth = FALSE,  :
  some notches went outside hinges ('box'): maybe set notch=FALSE
> anova(lm(weight ~ group, data = PlantGrowth))
Analysis of Variance Table

Response: weight
          Df  Sum Sq Mean Sq F value  Pr(>F)  
group      2  3.7663  1.8832  4.8461 0.01591 *
Residuals 27 10.4921  0.3886                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> 
> 
> cleanEx()
> nameEx("Puromycin")
> ### * Puromycin
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Puromycin
> ### Title: Reaction Velocity of an Enzymatic Reaction
> ### Aliases: Puromycin
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> ## Don't show: 
> options(show.nls.convergence=FALSE)
> ## End(Don't show)
> plot(rate ~ conc, data = Puromycin, las = 1,
+      xlab = "Substrate concentration (ppm)",
+      ylab = "Reaction velocity (counts/min/min)",
+      pch = as.integer(Puromycin$state),
+      col = as.integer(Puromycin$state),
+      main = "Puromycin data and fitted Michaelis-Menten curves")
> ## simplest form of fitting the Michaelis-Menten model to these data
> fm1 <- nls(rate ~ Vm * conc/(K + conc), data = Puromycin,
+            subset = state == "treated",
+            start = c(Vm = 200, K = 0.05))
> fm2 <- nls(rate ~ Vm * conc/(K + conc), data = Puromycin,
+            subset = state == "untreated",
+            start = c(Vm = 160, K = 0.05))
> summary(fm1)

Formula: rate ~ Vm * conc/(K + conc)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
Vm 2.127e+02  6.947e+00  30.615 3.24e-11 ***
K  6.412e-02  8.281e-03   7.743 1.57e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.93 on 10 degrees of freedom

> summary(fm2)

Formula: rate ~ Vm * conc/(K + conc)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
Vm 1.603e+02  6.480e+00  24.734 1.38e-09 ***
K  4.771e-02  7.782e-03   6.131 0.000173 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9.773 on 9 degrees of freedom

> ## add fitted lines to the plot
> conc <- seq(0, 1.2, length.out = 101)
> lines(conc, predict(fm1, list(conc = conc)), lty = 1, col = 1)
> lines(conc, predict(fm2, list(conc = conc)), lty = 2, col = 2)
> legend(0.8, 120, levels(Puromycin$state),
+        col = 1:2, lty = 1:2, pch = 1:2)
> 
> ## using partial linearity
> fm3 <- nls(rate ~ conc/(K + conc), data = Puromycin,
+            subset = state == "treated", start = c(K = 0.05),
+            algorithm = "plinear")
> 
> 
> 
> cleanEx()
> nameEx("Theoph")
> ### * Theoph
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Theoph
> ### Title: Pharmacokinetics of Theophylline
> ### Aliases: Theoph
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> ## Don't show: 
> options(show.nls.convergence=FALSE)
> ## End(Don't show)
> coplot(conc ~ Time | Subject, data = Theoph, show.given = FALSE)
> Theoph.4 <- subset(Theoph, Subject == 4)
> fm1 <- nls(conc ~ SSfol(Dose, Time, lKe, lKa, lCl),
+            data = Theoph.4)
> summary(fm1)

Formula: conc ~ SSfol(Dose, Time, lKe, lKa, lCl)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
lKe  -2.4365     0.2257 -10.797 4.77e-06 ***
lKa   0.1583     0.2297   0.689     0.51    
lCl  -3.2861     0.1448 -22.695 1.51e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8465 on 8 degrees of freedom

> plot(conc ~ Time, data = Theoph.4,
+      xlab = "Time since drug administration (hr)",
+      ylab = "Theophylline concentration (mg/L)",
+      main = "Observed concentrations and fitted model",
+      sub  = "Theophylline data - Subject 4 only",
+      las = 1, col = 4)
> xvals <- seq(0, par("usr")[2], length.out = 55)
> lines(xvals, predict(fm1, newdata = list(Time = xvals)),
+       col = 4)
> 
> 
> 
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("Titanic")
> ### * Titanic
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Titanic
> ### Title: Survival of passengers on the Titanic
> ### Aliases: Titanic
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> mosaicplot(Titanic, main = "Survival on the Titanic")
> ## Higher survival rates in children?
> apply(Titanic, c(3, 4), sum)
       Survived
Age       No Yes
  Child   52  57
  Adult 1438 654
> ## Higher survival rates in females?
> apply(Titanic, c(2, 4), sum)
        Survived
Sex        No Yes
  Male   1364 367
  Female  126 344
> ## Use loglm() in package 'MASS' for further analysis ...
> 
> 
> 
> cleanEx()
> nameEx("ToothGrowth")
> ### * ToothGrowth
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ToothGrowth
> ### Title: The Effect of Vitamin C on Tooth Growth in Guinea Pigs
> ### Aliases: ToothGrowth
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> coplot(len ~ dose | supp, data = ToothGrowth, panel = panel.smooth,
+        xlab = "ToothGrowth data: length vs dose, given type of supplement")
> 
> 
> 
> cleanEx()
> nameEx("UCBAdmissions")
> ### * UCBAdmissions
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: UCBAdmissions
> ### Title: Student Admissions at UC Berkeley
> ### Aliases: UCBAdmissions
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> ## Data aggregated over departments
> apply(UCBAdmissions, c(1, 2), sum)
          Gender
Admit      Male Female
  Admitted 1198    557
  Rejected 1493   1278
> mosaicplot(apply(UCBAdmissions, c(1, 2), sum),
+            main = "Student admissions at UC Berkeley")
> ## Data for individual departments
> opar <- par(mfrow = c(2, 3), oma = c(0, 0, 2, 0))
> for(i in 1:6)
+   mosaicplot(UCBAdmissions[,,i],
+     xlab = "Admit", ylab = "Sex",
+     main = paste("Department", LETTERS[i]))
> mtext(expression(bold("Student admissions at UC Berkeley")),
+       outer = TRUE, cex = 1.5)
> par(opar)
> 
> 
> 
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("UKDriverDeaths")
> ### * UKDriverDeaths
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: UKDriverDeaths
> ### Title: Road Casualties in Great Britain 1969-84
> ### Aliases: UKDriverDeaths Seatbelts
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> ## work with pre-seatbelt period to identify a model, use logs
> work <- window(log10(UKDriverDeaths), end = 1982+11/12)
> par(mfrow = c(3, 1))
> plot(work); acf(work); pacf(work)
> par(mfrow = c(1, 1))
> (fit <- arima(work, c(1, 0, 0), seasonal = list(order = c(1, 0, 0))))

Call:
arima(x = work, order = c(1, 0, 0), seasonal = list(order = c(1, 0, 0)))

Coefficients:
         ar1    sar1  intercept
      0.4378  0.6281     3.2274
s.e.  0.0764  0.0637     0.0131

sigma^2 estimated as 0.00157:  log likelihood = 300.85,  aic = -593.7
> z <- predict(fit, n.ahead = 24)
> ts.plot(log10(UKDriverDeaths), z$pred, z$pred+2*z$se, z$pred-2*z$se,
+         lty = c(1, 3, 2, 2), col = c("black", "red", "blue", "blue"))
> 
> ## now see the effect of the explanatory variables
> X <- Seatbelts[, c("kms", "PetrolPrice", "law")]
> X[, 1] <- log10(X[, 1]) - 4
> arima(log10(Seatbelts[, "drivers"]), c(1, 0, 0),
+       seasonal = list(order = c(1, 0, 0)), xreg = X)

Call:
arima(x = log10(Seatbelts[, "drivers"]), order = c(1, 0, 0), seasonal = list(order = c(1, 
    0, 0)), xreg = X)

Coefficients:
         ar1    sar1  intercept     kms  PetrolPrice      law
      0.3348  0.6672     3.3539  0.0082      -1.2224  -0.0963
s.e.  0.0775  0.0612     0.0441  0.0902       0.3839   0.0166

sigma^2 estimated as 0.001476:  log likelihood = 349.73,  aic = -685.46
> 
> 
> 
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("UKLungDeaths")
> ### * UKLungDeaths
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: UKLungDeaths
> ### Title: Monthly Deaths from Lung Diseases in the UK
> ### Aliases: UKLungDeaths ldeaths fdeaths mdeaths
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics) # for time
> plot(ldeaths)
> plot(mdeaths, fdeaths)
> ## Better labels:
> yr <- floor(tt <- time(mdeaths))
> plot(mdeaths, fdeaths,
+      xy.labels = paste(month.abb[12*(tt - yr)], yr-1900, sep = "'"))
> 
> 
> 
> cleanEx()
> nameEx("UKgas")
> ### * UKgas
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: UKgas
> ### Title: UK Quarterly Gas Consumption
> ### Aliases: UKgas
> ### Keywords: datasets
> 
> ### ** Examples
> 
> ## maybe str(UKgas) ; plot(UKgas) ...
> 
> 
> 
> cleanEx()
> nameEx("USArrests")
> ### * USArrests
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: USArrests
> ### Title: Violent Crime Rates by US State
> ### Aliases: USArrests
> ### Keywords: datasets
> 
> ### ** Examples
> 
> summary(USArrests)
     Murder          Assault         UrbanPop          Rape      
 Min.   : 0.800   Min.   : 45.0   Min.   :32.00   Min.   : 7.30  
 1st Qu.: 4.075   1st Qu.:109.0   1st Qu.:54.50   1st Qu.:15.07  
 Median : 7.250   Median :159.0   Median :66.00   Median :20.10  
 Mean   : 7.788   Mean   :170.8   Mean   :65.54   Mean   :21.23  
 3rd Qu.:11.250   3rd Qu.:249.0   3rd Qu.:77.75   3rd Qu.:26.18  
 Max.   :17.400   Max.   :337.0   Max.   :91.00   Max.   :46.00  
> 
> require(graphics)
> pairs(USArrests, panel = panel.smooth, main = "USArrests data")
> 
> ## Difference between 'USArrests' and its correction
> USArrests["Maryland", "UrbanPop"] # 67 -- the transcription error
[1] 67
> UA.C <- USArrests
> UA.C["Maryland", "UrbanPop"] <- 76.6
> 
> ## also +/- 0.5 to restore the original  <n>.5  percentages
> s5u <- c("Colorado", "Florida", "Mississippi", "Wyoming")
> s5d <- c("Nebraska", "Pennsylvania")
> UA.C[s5u, "UrbanPop"] <- UA.C[s5u, "UrbanPop"] + 0.5
> UA.C[s5d, "UrbanPop"] <- UA.C[s5d, "UrbanPop"] - 0.5
> 
> ## ==> UA.C  is now a *C*orrected version of  USArrests
> 
> 
> 
> cleanEx()
> nameEx("USJudgeRatings")
> ### * USJudgeRatings
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: USJudgeRatings
> ### Title: Lawyers' Ratings of State Judges in the US Superior Court
> ### Aliases: USJudgeRatings
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> pairs(USJudgeRatings, main = "USJudgeRatings data")
> 
> 
> 
> cleanEx()
> nameEx("USPersonalExpenditure")
> ### * USPersonalExpenditure
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: USPersonalExpenditure
> ### Title: Personal Expenditure Data
> ### Aliases: USPersonalExpenditure
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats) # for medpolish
> USPersonalExpenditure
                      1940   1945  1950 1955  1960
Food and Tobacco    22.200 44.500 59.60 73.2 86.80
Household Operation 10.500 15.500 29.00 36.5 46.20
Medical and Health   3.530  5.760  9.71 14.0 21.10
Personal Care        1.040  1.980  2.45  3.4  5.40
Private Education    0.341  0.974  1.80  2.6  3.64
> medpolish(log10(USPersonalExpenditure))
1: 1.126317
2: 1.032421
Final: 1.032421

Median Polish Results (Dataset: "log10(USPersonalExpenditure)")

Overall: 0.9872192

Row Effects:
   Food and Tobacco Household Operation  Medical and Health       Personal Care 
          0.7880270           0.4327608           0.0000000          -0.5606543 
  Private Education 
         -0.7319467 

Column Effects:
      1940       1945       1950       1955       1960 
-0.4288933 -0.2267967  0.0000000  0.1423128  0.3058289 

Residuals:
                         1940       1945      1950      1955      1960
Food and Tobacco     0.000000  0.0999105  0.000000 -0.053048 -0.142555
Household Operation  0.030103 -0.0028516  0.042418  0.000000 -0.061167
Medical and Health  -0.010551  0.0000000  0.000000  0.016596  0.031234
Personal Care        0.019362  0.0968971 -0.037399 -0.037399  0.000000
Private Education   -0.293625 -0.0399168  0.000000  0.017388  0.000000

> 
> 
> 
> cleanEx()
> nameEx("VADeaths")
> ### * VADeaths
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: VADeaths
> ### Title: Death Rates in Virginia (1940)
> ### Aliases: VADeaths
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> n <- length(dr <- c(VADeaths))
> nam <- names(VADeaths)
> d.VAD <- data.frame(
+  Drate = dr,
+  age = rep(ordered(rownames(VADeaths)), length.out = n),
+  gender = gl(2, 5, n, labels = c("M", "F")),
+  site =  gl(2, 10, labels = c("rural", "urban")))
> coplot(Drate ~ as.numeric(age) | gender * site, data = d.VAD,
+        panel = panel.smooth, xlab = "VADeaths data - Given: gender")
> summary(aov.VAD <- aov(Drate ~ .^2, data = d.VAD))
            Df Sum Sq Mean Sq F value   Pr(>F)    
age          4   6288  1572.1 590.858 8.55e-06 ***
gender       1    648   647.5 243.361 9.86e-05 ***
site         1     77    76.8  28.876  0.00579 ** 
age:gender   4     86    21.6   8.100  0.03358 *  
age:site     4     43    10.6   3.996  0.10414    
gender:site  1     73    73.0  27.422  0.00636 ** 
Residuals    4     11     2.7                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0))
> plot(aov.VAD)
> par(opar)
> 
> 
> 
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("WWWusage")
> ### * WWWusage
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: WWWusage
> ### Title: Internet Usage per Minute
> ### Aliases: WWWusage
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> work <- diff(WWWusage)
> par(mfrow = c(2, 1)); plot(WWWusage); plot(work)
> ## Not run: 
> ##D require(stats)
> ##D aics <- matrix(, 6, 6, dimnames = list(p = 0:5, q = 0:5))
> ##D for(q in 1:5) aics[1, 1+q] <- arima(WWWusage, c(0, 1, q),
> ##D     optim.control = list(maxit = 500))$aic
> ##D for(p in 1:5)
> ##D    for(q in 0:5) aics[1+p, 1+q] <- arima(WWWusage, c(p, 1, q),
> ##D        optim.control = list(maxit = 500))$aic
> ##D round(aics - min(aics, na.rm = TRUE), 2)
> ## End(Not run)
> 
> 
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("WorldPhones")
> ### * WorldPhones
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: WorldPhones
> ### Title: The World's Telephones
> ### Aliases: WorldPhones
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> matplot(rownames(WorldPhones), WorldPhones, type = "b", log = "y",
+         xlab = "Year", ylab = "Number of telephones (1000's)")
> legend(1951.5, 80000, colnames(WorldPhones), col = 1:6, lty = 1:5,
+        pch = rep(21, 7))
> title(main = "World phones data: log scale for response")
> 
> 
> 
> cleanEx()
> nameEx("ability.cov")
> ### * ability.cov
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ability.cov
> ### Title: Ability and Intelligence Tests
> ### Aliases: ability.cov
> ### Keywords: datasets
> 
> ### ** Examples
> 
> 
> cleanEx()
> nameEx("airmiles")
> ### * airmiles
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: airmiles
> ### Title: Passenger Miles on Commercial US Airlines, 1937-1960
> ### Aliases: airmiles
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> plot(airmiles, main = "airmiles data",
+      xlab = "Passenger-miles flown by U.S. commercial airlines", col = 4)
> 
> 
> 
> cleanEx()
> nameEx("airquality")
> ### * airquality
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: airquality
> ### Title: New York Air Quality Measurements
> ### Aliases: airquality
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> pairs(airquality, panel = panel.smooth, main = "airquality data")
> 
> 
> 
> cleanEx()
> nameEx("anscombe")
> ### * anscombe
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: anscombe
> ### Title: Anscombe's Quartet of 'Identical' Simple Linear Regressions
> ### Aliases: anscombe
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> summary(anscombe)
       x1             x2             x3             x4           y1        
 Min.   : 4.0   Min.   : 4.0   Min.   : 4.0   Min.   : 8   Min.   : 4.260  
 1st Qu.: 6.5   1st Qu.: 6.5   1st Qu.: 6.5   1st Qu.: 8   1st Qu.: 6.315  
 Median : 9.0   Median : 9.0   Median : 9.0   Median : 8   Median : 7.580  
 Mean   : 9.0   Mean   : 9.0   Mean   : 9.0   Mean   : 9   Mean   : 7.501  
 3rd Qu.:11.5   3rd Qu.:11.5   3rd Qu.:11.5   3rd Qu.: 8   3rd Qu.: 8.570  
 Max.   :14.0   Max.   :14.0   Max.   :14.0   Max.   :19   Max.   :10.840  
       y2              y3              y4        
 Min.   :3.100   Min.   : 5.39   Min.   : 5.250  
 1st Qu.:6.695   1st Qu.: 6.25   1st Qu.: 6.170  
 Median :8.140   Median : 7.11   Median : 7.040  
 Mean   :7.501   Mean   : 7.50   Mean   : 7.501  
 3rd Qu.:8.950   3rd Qu.: 7.98   3rd Qu.: 8.190  
 Max.   :9.260   Max.   :12.74   Max.   :12.500  
> 
> ##-- now some "magic" to do the 4 regressions in a loop:
> ff <- y ~ x
> mods <- setNames(as.list(1:4), paste0("lm", 1:4))
> for(i in 1:4) {
+   ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
+   ## or   ff[[2]] <- as.name(paste0("y", i))
+   ##      ff[[3]] <- as.name(paste0("x", i))
+   mods[[i]] <- lmi <- lm(ff, data = anscombe)
+   print(anova(lmi))
+ }
Analysis of Variance Table

Response: y1
          Df Sum Sq Mean Sq F value  Pr(>F)   
x1         1 27.510 27.5100   17.99 0.00217 **
Residuals  9 13.763  1.5292                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Analysis of Variance Table

Response: y2
          Df Sum Sq Mean Sq F value   Pr(>F)   
x2         1 27.500 27.5000  17.966 0.002179 **
Residuals  9 13.776  1.5307                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Analysis of Variance Table

Response: y3
          Df Sum Sq Mean Sq F value   Pr(>F)   
x3         1 27.470 27.4700  17.972 0.002176 **
Residuals  9 13.756  1.5285                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Analysis of Variance Table

Response: y4
          Df Sum Sq Mean Sq F value   Pr(>F)   
x4         1 27.490 27.4900  18.003 0.002165 **
Residuals  9 13.742  1.5269                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> ## See how close they are (numerically!)
> sapply(mods, coef)
                  lm1      lm2       lm3       lm4
(Intercept) 3.0000909 3.000909 3.0024545 3.0017273
x1          0.5000909 0.500000 0.4997273 0.4999091
> lapply(mods, function(fm) coef(summary(fm)))
$lm1
             Estimate Std. Error  t value    Pr(>|t|)
(Intercept) 3.0000909  1.1247468 2.667348 0.025734051
x1          0.5000909  0.1179055 4.241455 0.002169629

$lm2
            Estimate Std. Error  t value    Pr(>|t|)
(Intercept) 3.000909  1.1253024 2.666758 0.025758941
x2          0.500000  0.1179637 4.238590 0.002178816

$lm3
             Estimate Std. Error  t value    Pr(>|t|)
(Intercept) 3.0024545  1.1244812 2.670080 0.025619109
x3          0.4997273  0.1178777 4.239372 0.002176305

$lm4
             Estimate Std. Error  t value    Pr(>|t|)
(Intercept) 3.0017273  1.1239211 2.670763 0.025590425
x4          0.4999091  0.1178189 4.243028 0.002164602

> 
> ## Now, do what you should have done in the first place: PLOTS
> op <- par(mfrow = c(2, 2), mar = 0.1+c(4,4,1,1), oma =  c(0, 0, 2, 0))
> for(i in 1:4) {
+   ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
+   plot(ff, data = anscombe, col = "red", pch = 21, bg = "orange", cex = 1.2,
+        xlim = c(3, 19), ylim = c(3, 13))
+   abline(mods[[i]], col = "blue")
+ }
> mtext("Anscombe's 4 Regression data sets", outer = TRUE, cex = 1.5)
> par(op)
> 
> 
> 
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("attenu")
> ### * attenu
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: attenu
> ### Title: The Joyner-Boore Attenuation Data
> ### Aliases: attenu
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> ## check the data class of the variables
> sapply(attenu, data.class)
    event       mag   station      dist     accel 
"numeric" "numeric"  "factor" "numeric" "numeric" 
> summary(attenu)
     event            mag           station         dist       
 Min.   : 1.00   Min.   :5.000   117    :  5   Min.   :  0.50  
 1st Qu.: 9.00   1st Qu.:5.300   1028   :  4   1st Qu.: 11.32  
 Median :18.00   Median :6.100   113    :  4   Median : 23.40  
 Mean   :14.74   Mean   :6.084   112    :  3   Mean   : 45.60  
 3rd Qu.:20.00   3rd Qu.:6.600   135    :  3   3rd Qu.: 47.55  
 Max.   :23.00   Max.   :7.700   (Other):147   Max.   :370.00  
                                 NA's   : 16                   
     accel        
 Min.   :0.00300  
 1st Qu.:0.04425  
 Median :0.11300  
 Mean   :0.15422  
 3rd Qu.:0.21925  
 Max.   :0.81000  
                  
> pairs(attenu, main = "attenu data")
> coplot(accel ~ dist | as.factor(event), data = attenu, show.given = FALSE)
> coplot(log(accel) ~ log(dist) | as.factor(event),
+        data = attenu, panel = panel.smooth, show.given = FALSE)
> 
> 
> 
> cleanEx()
> nameEx("attitude")
> ### * attitude
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: attitude
> ### Title: The Chatterjee-Price Attitude Data
> ### Aliases: attitude
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> pairs(attitude, main = "attitude data")
> summary(attitude)
     rating        complaints     privileges       learning         raises     
 Min.   :40.00   Min.   :37.0   Min.   :30.00   Min.   :34.00   Min.   :43.00  
 1st Qu.:58.75   1st Qu.:58.5   1st Qu.:45.00   1st Qu.:47.00   1st Qu.:58.25  
 Median :65.50   Median :65.0   Median :51.50   Median :56.50   Median :63.50  
 Mean   :64.63   Mean   :66.6   Mean   :53.13   Mean   :56.37   Mean   :64.63  
 3rd Qu.:71.75   3rd Qu.:77.0   3rd Qu.:62.50   3rd Qu.:66.75   3rd Qu.:71.00  
 Max.   :85.00   Max.   :90.0   Max.   :83.00   Max.   :75.00   Max.   :88.00  
    critical        advance     
 Min.   :49.00   Min.   :25.00  
 1st Qu.:69.25   1st Qu.:35.00  
 Median :77.50   Median :41.00  
 Mean   :74.77   Mean   :42.93  
 3rd Qu.:80.00   3rd Qu.:47.75  
 Max.   :92.00   Max.   :72.00  
> summary(fm1 <- lm(rating ~ ., data = attitude))

Call:
lm(formula = rating ~ ., data = attitude)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.9418  -4.3555   0.3158   5.5425  11.5990 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 10.78708   11.58926   0.931 0.361634    
complaints   0.61319    0.16098   3.809 0.000903 ***
privileges  -0.07305    0.13572  -0.538 0.595594    
learning     0.32033    0.16852   1.901 0.069925 .  
raises       0.08173    0.22148   0.369 0.715480    
critical     0.03838    0.14700   0.261 0.796334    
advance     -0.21706    0.17821  -1.218 0.235577    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.068 on 23 degrees of freedom
Multiple R-squared:  0.7326,	Adjusted R-squared:  0.6628 
F-statistic:  10.5 on 6 and 23 DF,  p-value: 1.24e-05

> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0),
+             mar = c(4.1, 4.1, 2.1, 1.1))
> plot(fm1)
> summary(fm2 <- lm(rating ~ complaints, data = attitude))

Call:
lm(formula = rating ~ complaints, data = attitude)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.8799  -5.9905   0.1783   6.2978   9.6294 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 14.37632    6.61999   2.172   0.0385 *  
complaints   0.75461    0.09753   7.737 1.99e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.993 on 28 degrees of freedom
Multiple R-squared:  0.6813,	Adjusted R-squared:  0.6699 
F-statistic: 59.86 on 1 and 28 DF,  p-value: 1.988e-08

> plot(fm2)
> par(opar)
> 
> 
> 
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("beavers")
> ### * beavers
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: beavers
> ### Title: Body Temperature Series of Two Beavers
> ### Aliases: beavers beaver1 beaver2
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> (yl <- range(beaver1$temp, beaver2$temp))
[1] 36.33 38.35
> 
> beaver.plot <- function(bdat, ...) {
+   nam <- deparse(substitute(bdat))
+   with(bdat, {
+     # Hours since start of day:
+     hours <- time %/% 100 + 24*(day - day[1]) + (time %% 100)/60
+     plot (hours, temp, type = "l", ...,
+           main = paste(nam, "body temperature"))
+     abline(h = 37.5, col = "gray", lty = 2)
+     is.act <- activ == 1
+     points(hours[is.act], temp[is.act], col = 2, cex = .8)
+   })
+ }
> op <- par(mfrow = c(2, 1), mar = c(3, 3, 4, 2), mgp = 0.9 * 2:0)
>  beaver.plot(beaver1, ylim = yl)
>  beaver.plot(beaver2, ylim = yl)
> par(op)
> 
> 
> 
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("cars")
> ### * cars
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: cars
> ### Title: Speed and Stopping Distances of Cars
> ### Aliases: cars
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)",
+      las = 1)
> lines(lowess(cars$speed, cars$dist, f = 2/3, iter = 3), col = "red")
> title(main = "cars data")
> plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)",
+      las = 1, log = "xy")
> title(main = "cars data (logarithmic scales)")
> lines(lowess(cars$speed, cars$dist, f = 2/3, iter = 3), col = "red")
> summary(fm1 <- lm(log(dist) ~ log(speed), data = cars))

Call:
lm(formula = log(dist) ~ log(speed), data = cars)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.00215 -0.24578 -0.02898  0.20717  0.88289 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.7297     0.3758  -1.941   0.0581 .  
log(speed)    1.6024     0.1395  11.484 2.26e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4053 on 48 degrees of freedom
Multiple R-squared:  0.7331,	Adjusted R-squared:  0.7276 
F-statistic: 131.9 on 1 and 48 DF,  p-value: 2.259e-15

> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0),
+             mar = c(4.1, 4.1, 2.1, 1.1))
> plot(fm1)
> par(opar)
> 
> ## An example of polynomial regression
> plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)",
+     las = 1, xlim = c(0, 25))
> d <- seq(0, 25, length.out = 200)
> for(degree in 1:4) {
+   fm <- lm(dist ~ poly(speed, degree), data = cars)
+   assign(paste("cars", degree, sep = "."), fm)
+   lines(d, predict(fm, data.frame(speed = d)), col = degree)
+ }
> anova(cars.1, cars.2, cars.3, cars.4)
Analysis of Variance Table

Model 1: dist ~ poly(speed, degree)
Model 2: dist ~ poly(speed, degree)
Model 3: dist ~ poly(speed, degree)
Model 4: dist ~ poly(speed, degree)
  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1     48 11354                           
2     47 10825  1    528.81 2.3108 0.1355
3     46 10634  1    190.35 0.8318 0.3666
4     45 10298  1    336.55 1.4707 0.2316
> 
> 
> 
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("chickwts")
> ### * chickwts
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: chickwts
> ### Title: Chicken Weights by Feed Type
> ### Aliases: chickwts
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> boxplot(weight ~ feed, data = chickwts, col = "lightgray",
+     varwidth = TRUE, notch = TRUE, main = "chickwt data",
+     ylab = "Weight at six weeks (gm)")
Warning in (function (z, notch = FALSE, width = NULL, varwidth = FALSE,  :
  some notches went outside hinges ('box'): maybe set notch=FALSE
> anova(fm1 <- lm(weight ~ feed, data = chickwts))
Analysis of Variance Table

Response: weight
          Df Sum Sq Mean Sq F value    Pr(>F)    
feed       5 231129   46226  15.365 5.936e-10 ***
Residuals 65 195556    3009                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0),
+             mar = c(4.1, 4.1, 2.1, 1.1))
> plot(fm1)
> par(opar)
> 
> 
> 
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("co2")
> ### * co2
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: co2
> ### Title: Mauna Loa Atmospheric CO2 Concentration
> ### Aliases: co2
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> plot(co2, ylab = expression("Atmospheric concentration of CO"[2]),
+      las = 1)
> title(main = "co2 data set")
> 
> 
> 
> cleanEx()
> nameEx("crimtab")
> ### * crimtab
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: crimtab
> ### Title: Student's 3000 Criminals Data
> ### Aliases: crimtab
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats)
> dim(crimtab)
[1] 42 22
> utils::str(crimtab)
 'table' int [1:42, 1:22] 0 0 0 0 0 0 1 0 0 0 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:42] "9.4" "9.5" "9.6" "9.7" ...
  ..$ : chr [1:22] "142.24" "144.78" "147.32" "149.86" ...
> ## for nicer printing:
> local({cT <- crimtab
+        colnames(cT) <- substring(colnames(cT), 2, 3)
+        print(cT, zero.print = " ")
+ })
     42 44 47 49 52 54 57 60 62 65 67 70 72 75 77 80 82 85 87 90 93 95
9.4                                                                   
9.5                  1                                                
9.6                                                                   
9.7                                                                   
9.8                     1                                             
9.9         1     1     1                                             
10    1        1  2     2        1                                    
10.1           1  3  1     1  1                                       
10.2        2  2  2  1     2     1                                    
10.3     1  1  3  2  2  3  5                                          
10.4        1  1  2  3  3  4  3  3                                    
10.5           1  3  7  6  4  3  1  3  1     1                        
10.6           1  4  5  9 14  6  3  1        1                        
10.7        1  2  4  9 14 16 15  7  3  1  2                           
10.8           2  5  6 14 27 10  7  1  2  1                           
10.9              2  6 14 24 27 14 10  4  1                           
11             2  6 12 15 31 37 27 17 10  6                           
11.1           3  3 12 22 26 24 26 24  7  4  1                        
11.2           3  2  7 21 30 38 29 27 20  4  1                       1
11.3           1     5 10 24 26 39 26 24  7  2                        
11.4              3  4  9 29 56 58 26 22 10 11                        
11.5                 5 11 17 33 57 38 34 25 11  2                     
11.6              2  1  4 13 37 39 48 38 27 12  2  2     1            
11.7                 2  9 17 30 37 48 45 24  9  9  2                  
11.8              1     2 11 15 35 41 34 29 10  5  1                  
11.9              1  1  2 12 10 27 32 35 19 10  9  3  1               
12                      1  4  8 19 42 39 22 16  8  2  2               
12.1                       2  4 13 22 28 15 27 10  4  1               
12.2                    1  2  5  6 23 17 16 11  8  1  1               
12.3                          4  8 10 13 20 23  6  5                  
12.4                    1  1  1  2  7 12  4  7  7  1        1         
12.5                       1     1  3 12 11  8  6  8     2            
12.6                             1     3  5  7  8  6  3  1  1         
12.7                             1  1  7  5  5  8  2  2               
12.8                                1  2  3  1  8  5  3  1  1         
12.9                                   1  2  2     1  1               
13                                  3     1     1     2  1            
13.1                                   1  1                           
13.2                                1  1     1     3                  
13.3                                                  1     1         
13.4                                                                  
13.5                                                     1            
> 
> ## Repeat Student's experiment:
> 
> # 1) Reconstitute 3000 raw data for heights in inches and rounded to
> #    nearest integer as in Student's paper:
> 
> (heIn <- round(as.numeric(colnames(crimtab)) / 2.54))
 [1] 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
> d.hei <- data.frame(height = rep(heIn, colSums(crimtab)))
> 
> # 2) shuffle the data:
> 
> set.seed(1)
> d.hei <- d.hei[sample(1:3000), , drop = FALSE]
> 
> # 3) Make 750 samples each of size 4:
> 
> d.hei$sample <- as.factor(rep(1:750, each = 4))
> 
> # 4) Compute the means and standard deviations (n) for the 750 samples:
> 
> h.mean <- with(d.hei, tapply(height, sample, FUN = mean))
> h.sd   <- with(d.hei, tapply(height, sample, FUN = sd)) * sqrt(3/4)
> 
> # 5) Compute the difference between the mean of each sample and
> #    the mean of the population and then divide by the
> #    standard deviation of the sample:
> 
> zobs <- (h.mean - mean(d.hei[,"height"]))/h.sd
> 
> # 6) Replace infinite values by +/- 6 as in Student's paper:
> 
> zobs[infZ <- is.infinite(zobs)] # none of them 
named numeric(0)
> zobs[infZ] <- 6 * sign(zobs[infZ])
> 
> # 7) Plot the distribution:
> 
> require(grDevices); require(graphics)
> hist(x = zobs, probability = TRUE, xlab = "Student's z",
+      col = grey(0.8), border = grey(0.5),
+      main = "Distribution of Student's z score  for 'crimtab' data")
> 
> 
> 
> cleanEx()
> nameEx("discoveries")
> ### * discoveries
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: discoveries
> ### Title: Yearly Numbers of Important Discoveries
> ### Aliases: discoveries
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> plot(discoveries, ylab = "Number of important discoveries",
+      las = 1)
> title(main = "discoveries data set")
> 
> 
> 
> cleanEx()
> nameEx("esoph")
> ### * esoph
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: esoph
> ### Title: Smoking, Alcohol and (O)esophageal Cancer
> ### Aliases: esoph
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats)
> require(graphics) # for mosaicplot
> summary(esoph)
   agegp          alcgp         tobgp        ncases         ncontrols     
 25-34:15   0-39g/day:23   0-9g/day:24   Min.   : 0.000   Min.   : 0.000  
 35-44:15   40-79    :23   10-19   :24   1st Qu.: 0.000   1st Qu.: 1.000  
 45-54:16   80-119   :21   20-29   :20   Median : 1.000   Median : 4.000  
 55-64:16   120+     :21   30+     :20   Mean   : 2.273   Mean   : 8.807  
 65-74:15                                3rd Qu.: 4.000   3rd Qu.:10.000  
 75+  :11                                Max.   :17.000   Max.   :60.000  
> ## effects of alcohol, tobacco and interaction, age-adjusted
> model1 <- glm(cbind(ncases, ncontrols) ~ agegp + tobgp * alcgp,
+               data = esoph, family = binomial())
> anova(model1)
Analysis of Deviance Table

Model: binomial, link: logit

Response: cbind(ncases, ncontrols)

Terms added sequentially (first to last)


            Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                           87     367.95              
agegp        5  121.045        82     246.91 < 2.2e-16 ***
tobgp        3   36.639        79     210.27 5.485e-08 ***
alcgp        3  127.933        76      82.34 < 2.2e-16 ***
tobgp:alcgp  9    5.451        67      76.89    0.7934    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> ## Try a linear effect of alcohol and tobacco
> model2 <- glm(cbind(ncases, ncontrols) ~ agegp + unclass(tobgp)
+                                          + unclass(alcgp),
+               data = esoph, family = binomial())
> summary(model2)

Call:
glm(formula = cbind(ncases, ncontrols) ~ agegp + unclass(tobgp) + 
    unclass(alcgp), family = binomial(), data = esoph)

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -4.92067    0.37344 -13.177  < 2e-16 ***
agegp.L         3.81892    0.67862   5.627 1.83e-08 ***
agegp.Q        -1.49473    0.60671  -2.464   0.0138 *  
agegp.C         0.07923    0.46318   0.171   0.8642    
agegp^4         0.12136    0.32203   0.377   0.7063    
agegp^5        -0.24856    0.21153  -1.175   0.2400    
unclass(tobgp)  0.43955    0.09623   4.568 4.93e-06 ***
unclass(alcgp)  1.06766    0.10493  10.175  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 367.953  on 87  degrees of freedom
Residual deviance:  91.121  on 80  degrees of freedom
AIC: 222.18

Number of Fisher Scoring iterations: 6

> ## Re-arrange data for a mosaic plot
> ttt <- table(esoph$agegp, esoph$alcgp, esoph$tobgp)
> o <- with(esoph, order(tobgp, alcgp, agegp))
> ttt[ttt == 1] <- esoph$ncases[o]
> tt1 <- table(esoph$agegp, esoph$alcgp, esoph$tobgp)
> tt1[tt1 == 1] <- esoph$ncontrols[o]
> tt <- array(c(ttt, tt1), c(dim(ttt),2),
+             c(dimnames(ttt), list(c("Cancer", "control"))))
> mosaicplot(tt, main = "esoph data set", color = TRUE)
> 
> 
> 
> cleanEx()
> nameEx("euro")
> ### * euro
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: euro
> ### Title: Conversion Rates of Euro Currencies
> ### Aliases: euro euro.cross
> ### Keywords: datasets
> 
> ### ** Examples
> 
> cbind(euro)
           euro
ATS   13.760300
BEF   40.339900
DEM    1.955830
ESP  166.386000
FIM    5.945730
FRF    6.559570
IEP    0.787564
ITL 1936.270000
LUF   40.339900
NLG    2.203710
PTE  200.482000
> 
> ## These relations hold:
> euro == signif(euro, 6) # [6 digit precision in Euro's definition]
 ATS  BEF  DEM  ESP  FIM  FRF  IEP  ITL  LUF  NLG  PTE 
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE 
> all(euro.cross == outer(1/euro, euro))
[1] TRUE
> 
> ## Convert 20 Euro to Belgian Franc
> 20 * euro["BEF"]
    BEF 
806.798 
> ## Convert 20 Austrian Schilling to Euro
> 20 / euro["ATS"]
     ATS 
1.453457 
> ## Convert 20 Spanish Pesetas to Italian Lira
> 20 * euro.cross["ESP", "ITL"]
[1] 232.7443
> 
> require(graphics)
> dotchart(euro,
+          main = "euro data: 1 Euro in currency unit")
> dotchart(1/euro,
+          main = "euro data: 1 currency unit in Euros")
> dotchart(log(euro, 10),
+          main = "euro data: log10(1 Euro in currency unit)")
> 
> 
> 
> cleanEx()
> nameEx("faithful")
> ### * faithful
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: faithful
> ### Title: Old Faithful Geyser Data
> ### Aliases: faithful
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> f.tit <-  "faithful data: Eruptions of Old Faithful"
> 
> ne60 <- round(e60 <- 60 * faithful$eruptions)
> all.equal(e60, ne60)             # relative diff. ~ 1/10000
[1] "Mean relative difference: 9.515332e-05"
> table(zapsmall(abs(e60 - ne60))) # 0, 0.02 or 0.04

   0 0.02 0.04 
 106  163    3 
> faithful$better.eruptions <- ne60 / 60
> te <- table(ne60)
> te[te >= 4]                      # (too) many multiples of 5 !
ne60
105 108 110 112 113 120 216 230 240 245 249 250 255 260 261 262 265 270 272 275 
  6   4   7   8   4   4   4   5   6   5   4   4   4   5   4   4   4   8   5   4 
276 282 288 
  4   6   6 
> plot(names(te), te, type = "h", main = f.tit, xlab = "Eruption time (sec)")
> 
> plot(faithful[, -3], main = f.tit,
+      xlab = "Eruption time (min)",
+      ylab = "Waiting time to next eruption (min)")
> lines(lowess(faithful$eruptions, faithful$waiting, f = 2/3, iter = 3),
+       col = "red")
> 
> 
> 
> cleanEx()
> nameEx("freeny")
> ### * freeny
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: freeny
> ### Title: Freeny's Revenue Data
> ### Aliases: freeny freeny.x freeny.y
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> summary(freeny)
       y         lag.quarterly.revenue  price.index     income.level  
 Min.   :8.791   Min.   :8.791         Min.   :4.278   Min.   :5.821  
 1st Qu.:9.045   1st Qu.:9.020         1st Qu.:4.392   1st Qu.:5.948  
 Median :9.314   Median :9.284         Median :4.510   Median :6.061  
 Mean   :9.306   Mean   :9.281         Mean   :4.496   Mean   :6.039  
 3rd Qu.:9.591   3rd Qu.:9.561         3rd Qu.:4.605   3rd Qu.:6.139  
 Max.   :9.794   Max.   :9.775         Max.   :4.710   Max.   :6.200  
 market.potential
 Min.   :12.97   
 1st Qu.:13.01   
 Median :13.07   
 Mean   :13.07   
 3rd Qu.:13.12   
 Max.   :13.17   
> pairs(freeny, main = "freeny data")
> # gives warning: freeny$y has class "ts"
> 
> summary(fm1 <- lm(y ~ ., data = freeny))

Call:
lm(formula = y ~ ., data = freeny)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0259426 -0.0101033  0.0003824  0.0103236  0.0267124 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)           -10.4726     6.0217  -1.739   0.0911 .  
lag.quarterly.revenue   0.1239     0.1424   0.870   0.3904    
price.index            -0.7542     0.1607  -4.693 4.28e-05 ***
income.level            0.7675     0.1339   5.730 1.93e-06 ***
market.potential        1.3306     0.5093   2.613   0.0133 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01473 on 34 degrees of freedom
Multiple R-squared:  0.9981,	Adjusted R-squared:  0.9978 
F-statistic:  4354 on 4 and 34 DF,  p-value: < 2.2e-16

> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0),
+             mar = c(4.1, 4.1, 2.1, 1.1))
> plot(fm1)
> par(opar)
> 
> 
> 
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("gait")
> ### * gait
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: gait
> ### Title: Hip and Knee Angle while Walking
> ### Aliases: gait
> ### Keywords: datasets
> 
> ### ** Examples
> 
> plot(gait[, 1, ], type = "b",
+      xlim = range(gait[,,1]), ylim = range(gait[,,2]), 
+      xlab = "Hip Angle", ylab = "Knee Angle", main = "'gait' data : Boy 1")
> mtext("all other boys", col = "thistle"); grid()
> matlines(gait[, -1, 1], gait[, -1, 2], type = "l", lty = 1, col = adjustcolor("thistle", 1/3))
> 
> ## The data array, two matrices :
> op <- options(width = 128) # on a wide console
> aperm(gait, c(2:1, 3))
, , Variable = Hip Angle

       Time
Subject 0.025 0.075 0.125 0.175 0.225 0.275 0.325 0.375 0.425 0.475 0.525 0.575 0.625 0.675 0.725 0.775 0.825 0.875 0.925 0.975
  boy1     37    36    33    29    23    18    15    12     9     6     4     6    13    22    31    38    43    44    40    35
  boy2     47    46    42    34    27    21    16    12     4    -1    -2     2    12    24    34    42    47    48    45    43
  boy3     46    44    39    34    33    27    21    15     7     1    -3    -1    10    22    34    43    51    53    50    49
  boy4     37    36    27    20    15    15    10     6     3     0     3     4    10    14    23    31    33    35    34    32
  boy5     20    18    11     8     7     5     1     0    -2    -4    -4     1    10    18    27    36    34    33    32    28
  boy6     57    48    44    35    31    27    21    18    14    10     8     8    11    21    32    40    49    55    56    55
  boy7     46    38    33    25    18    15    14    10     8     4     0     3     9    19    31    37    46    45    46    41
  boy8     46    46    43    40    36    30    24    17    12     7     1     3    12    24    35    44    49    49    47    43
  boy9     46    42    37    34    31    25    19    17    13    10     5     7    13    21    29    39    47    51    52    47
  boy10    35    34    29    28    19    15    13     6     0    -1     2    11    20    30    39    45    48    48    43    39
  boy11    38    37    33    29    26    20    14     9     3    -2    -1     5    15    24    36    44    47    45    42    40
  boy12    35    31    29    26    22    19    15    11     6     3     0     4    16    27    35    41    44    45    40    39
  boy13    34    31    27    23    19    15    13    10     6     1     0     3     8    16    27    34    40    41    39    36
  boy14    43    41    36    31    26    20    16    13    10     8     4    12    19    27    37    41    44    44    41    37
  boy15    43    37    35    28    26    21    20    15    13     7     3     8    16    25    35    43    49    50    45    46
  boy16    40    41    36    32    27    20    17    10     5    -2    -4    -2     7    24    38    49    57    59    54    46
  boy17    51    49    45    39    31    23    16     9     4     0    -5     1    10    21    33    44    51    55    56    51
  boy18    52    46    41    35    31    24    18    12     6     3     2     5    12    22    32    42    46    48    49    46
  boy19    36    33    28    22    18    13    11     7     1    -2    -3    -3     2    13    22    32    39    41    38    34
  boy20    35    37    33    27    22    14     9     6     4    -1     1     2    11    17    34    35    40    43    43    42
  boy21    46    38    30    23    17    13     7     3     0    -3    -3     1    11    20    36    48    55    57    56    50
  boy22    43    41    37    30    24    16    12     8     5     2     1     8    16    24    33    42    48    48    48    46
  boy23    55    51    47    41    35    30    26    23    20    11     8    10    19    35    41    52    57    61    63    58
  boy24    39    38    31    27    21    14     9     9     6     4     2     3    13    30    37    47    53    53    49    44
  boy25    37    34    30    27    26    19    15    10     4     0     0     4    12    22    31    40    43    43    38    38
  boy26    36    33    28    22    18    13    11     7     1    -2    -3    -3     2    13    22    32    39    41    38    34
  boy27    36    33    30    28    21    15     8     1    -5   -11   -12    -7     4    16    26    37    44    47    44    37
  boy28    42    40    40    34    23    15    11     7     5     6     8    12    22    33    43    51    57    58    54    46
  boy29    38    34    30    23    17    12     8     4     1     0    -4    -4     1    10    22    32    38    41    41    40
  boy30    46    47    44    37    29    23    19    14     8     3    -2     1    12    26    39    48    52    48    43    42
  boy31    54    48    44    37    30    27    21    18    15    11    12    20    30    41    50    56    61    59    57    58
  boy32    52    44    44    33    28    27    23    24    19    15    15    16    24    32    43    52    58    59    57    52
  boy33    32    28    26    22    19    13     8     5    -1    -6    -5     0    12    22    30    36    39    36    30    29
  boy34    46    41    38    31    25    20    13     7     1    -4    -3    -2     8    22    34    45    53    57    55    43
  boy35    46    44    40    35    31    25    19    15    10     5     3     6    14    27    40    48    53    53    50    47
  boy36    48    42    42    35    30    23    19    14     9     4    -1     3     9    25    37    45    52    53    52    46
  boy37    44    41    38    32    24    18    10     6     3     0    -2    -1     6    20    36    44    49    46    38    35
  boy38    55    56    51    46    41    36    30    25    21    15     8     5     9    19    31    43    52    56    59    59
  boy39    48    50    47    42    37    29    22    14     8     5     8    15    24    36    51    59    63    64    61    55

, , Variable = Knee Angle

       Time
Subject 0.025 0.075 0.125 0.175 0.225 0.275 0.325 0.375 0.425 0.475 0.525 0.575 0.625 0.675 0.725 0.775 0.825 0.875 0.925 0.975
  boy1     10    15    18    18    15    14    13    13    15    21    26    39    52    64    70    66    57    40    22    11
  boy2     16    25    28    25    18    12    10     9     7     9    17    32    51    67    70    66    55    39    23    16
  boy3     18    27    32    32    28    23    16    11     8    11    18    33    56    72    79    77    67    46    22    14
  boy4      5    14    16    17    10     8     6     5     5     9    20    40    57    66    67    61    43    18     3     0
  boy5      2     6     6     6     5     6     6     7     8    15    24    41    59    69    71    66    44    23     5     4
  boy6     15    17    23    23    20    19    15    13     9    10    14    25    40    58    65    66    61    45    21     9
  boy7     13    16    22    17    12     9     8    10    14    16    22    36    51    64    71    68    65    40    17     3
  boy8     14    17    18    19    19    15    11     8     7     8    12    26    47    68    79    76    61    35    13     5
  boy9     15    20    23    26    25    21    16    13     8     8    10    22    40    56    65    67    61    47    26    10
  boy10     9    22    25    21    10     9     9     9    11    19    33    55    69    75    77    70    55    36    17    16
  boy11    13    24    27    23    18    13    12    15    16    19    27    43    60    70    75    69    50    25    13     6
  boy12     7     8    11    12     8     6     4     5     5    10    19    36    58    70    72    66    51    30     8     5
  boy13     9    14    16    15    15    12    13    14    14    15    20    31    45    60    71    70    62    45    26    12
  boy14    15    20    22    22    21    19    19    19    19    20    25    41    58    70    73    68    53    32    16    10
  boy15     6    11    20    18    13     9     9     8    11    12    20    38    56    67    71    65    48    25     5     8
  boy16    11    19    30    28    25    17    15    10     7     7    13    25    46    68    79    77    67    48    23    13
  boy17    24    32    35    33    29    24    18    14    12    14    20    36    55    69    75    72    60    43    26    16
  boy18    16    20    21    20    18    14    13    11     9    11    19    34    53    66    73    72    61    43    22    11
  boy19    16    20    22    21    20    20    20    19    16    16    21    31    49    69    74    71    60    41    23    15
  boy20     7    13    14    17    14     8     7     7     7     8    17    30    49    60    71    65    58    39    20    10
  boy21    21    24    25    21    16     9     4     4     5     9    19    36    56    70    80    80    68    48    27    16
  boy22    11    14    14    11     8     5     6     7     8    10    18    36    53    65    72    73    64    46    28    13
  boy23    12    17    20    20    18    12    11     9    11    10    16    30    49    68    69    69    63    44    28     8
  boy24     8    12    14    13    12     9     7    10     9    13    21    39    59    76    76    71    55    26     8    12
  boy25    11    20    22    21    21    17    16    17    15    19    28    42    59    70    73    66    49    24     6     9
  boy26    16    20    22    21    20    20    20    19    16    16    21    31    49    69    74    71    60    41    23    15
  boy27    19    26    28    28    24    18    12     8     7     8    16    33    52    68    78    80    74    57    35    18
  boy28    13    23    30    28    19    10     6     5    10    19    33    50    68    79    81    75    61    42    25    16
  boy29    11    15    19    20    18    17    15    11     9     9    12    23    40    58    71    70    56    36    19    11
  boy30    17    25    30    30    27    22    17    11     7     7    13    31    52    69    75    68    46    21     6     8
  boy31    20    20    22    16    10    10    10    12    13    17    28    47    64    72    72    63    47    23     9    19
  boy32    18    18    25    23    18    19    21    25    24    25    27    39    59    72    77    69    60    38    17    14
  boy33     9    12    16    15    14    11     8     6     5     9    18    39    56    68    71    67    54    30     9     6
  boy34     8    10    17    16    12    10     8     4     2     4    12    27    48    66    76    77    71    54    26     3
  boy35     9    18    19    19    19    15    12    11    10    12    20    36    56    72    80    75    59    34    10     9
  boy36    13    18    27    26    25    18    12    10     9    13    16    33    53    73    77    75    64    44    23    16
  boy37    19    23    26    25    21    18    14    13    14    17    24    37    56    73    82    76    59    31     7     7
  boy38    16    23    28    28    25    21    19    17    16    13    13    19    37    57    69    71    62    45    28    20
  boy39    14    25    32    34    30    20    12     7     6    11    24    45    67    79    82    76    65    46    25    15

> options(op)
> 
> 
> 
> cleanEx()
> nameEx("infert")
> ### * infert
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: infert
> ### Title: Infertility after Spontaneous and Induced Abortion
> ### Aliases: infert
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats)
> model1 <- glm(case ~ spontaneous+induced, data = infert, family = binomial())
> summary(model1)

Call:
glm(formula = case ~ spontaneous + induced, family = binomial(), 
    data = infert)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.7079     0.2677  -6.380 1.78e-10 ***
spontaneous   1.1972     0.2116   5.657 1.54e-08 ***
induced       0.4181     0.2056   2.033    0.042 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 316.17  on 247  degrees of freedom
Residual deviance: 279.61  on 245  degrees of freedom
AIC: 285.61

Number of Fisher Scoring iterations: 4

> ## adjusted for other potential confounders:
> summary(model2 <- glm(case ~ age+parity+education+spontaneous+induced,
+                      data = infert, family = binomial()))

Call:
glm(formula = case ~ age + parity + education + spontaneous + 
    induced, family = binomial(), data = infert)

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -1.14924    1.41220  -0.814   0.4158    
age               0.03958    0.03120   1.269   0.2046    
parity           -0.82828    0.19649  -4.215 2.49e-05 ***
education6-11yrs -1.04424    0.79255  -1.318   0.1876    
education12+ yrs -1.40321    0.83416  -1.682   0.0925 .  
spontaneous       2.04591    0.31016   6.596 4.21e-11 ***
induced           1.28876    0.30146   4.275 1.91e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 316.17  on 247  degrees of freedom
Residual deviance: 257.80  on 241  degrees of freedom
AIC: 271.8

Number of Fisher Scoring iterations: 4

> ## Really should be analysed by conditional logistic regression
> ## which is in the survival package
> 
> 
> 
> cleanEx()
> nameEx("iris")
> ### * iris
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: iris
> ### Title: Edgar Anderson's Iris Data
> ### Aliases: iris iris3
> ### Keywords: datasets
> 
> ### ** Examples
> 
> summary(iris)
  Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
 Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
 1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
 Median :5.800   Median :3.000   Median :4.350   Median :1.300  
 Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
 3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
 Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
       Species  
 setosa    :50  
 versicolor:50  
 virginica :50  
                
                
                
> 
> ## Fisher's (1936) research question: whether (compound measurements of)
> ## Iris versicolor "differs twice as much from I. setosa as from I. virginica"
> pairs(iris[1:4], col = iris$Species)
> legend(0.5, 1, levels(iris$Species), fill = 1:3, bty = "n",
+        horiz = TRUE, xjust = 0.5, yjust = 0, xpd = TRUE)
> 
> ## equivalence of legacy array (iris3) and data.frame (iris) representation
> dni3 <- dimnames(iris3)
> ii <- data.frame(matrix(aperm(iris3, c(1,3,2)), ncol = 4,
+                         dimnames = list(NULL, sub(" L.",".Length",
+                                         sub(" W.",".Width", dni3[[2]])))),
+     Species = gl(3, 50, labels = sub("S", "s", sub("V", "v", dni3[[3]]))))
> stopifnot(all.equal(ii, iris))
> 
> 
> 
> cleanEx()
> nameEx("islands")
> ### * islands
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: islands
> ### Title: Areas of the World's Major Landmasses
> ### Aliases: islands
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> dotchart(log(islands, 10),
+    main = "islands data: log10(area) (log10(sq. miles))")
> dotchart(log(islands[order(islands)], 10),
+    main = "islands data: log10(area) (log10(sq. miles))")
> 
> 
> 
> cleanEx()
> nameEx("longley")
> ### * longley
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: longley
> ### Title: Longley's Economic Regression Data
> ### Aliases: longley
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> ## give the data set in the form it was used in S-PLUS:
> longley.x <- data.matrix(longley[, 1:6])
> longley.y <- longley[, "Employed"]
> pairs(longley, main = "longley data")
> summary(fm1 <- lm(Employed ~ ., data = longley))

Call:
lm(formula = Employed ~ ., data = longley)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.41011 -0.15767 -0.02816  0.10155  0.45539 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -3.482e+03  8.904e+02  -3.911 0.003560 ** 
GNP.deflator  1.506e-02  8.492e-02   0.177 0.863141    
GNP          -3.582e-02  3.349e-02  -1.070 0.312681    
Unemployed   -2.020e-02  4.884e-03  -4.136 0.002535 ** 
Armed.Forces -1.033e-02  2.143e-03  -4.822 0.000944 ***
Population   -5.110e-02  2.261e-01  -0.226 0.826212    
Year          1.829e+00  4.555e-01   4.016 0.003037 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3049 on 9 degrees of freedom
Multiple R-squared:  0.9955,	Adjusted R-squared:  0.9925 
F-statistic: 330.3 on 6 and 9 DF,  p-value: 4.984e-10

> opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1.1, 0),
+             mar = c(4.1, 4.1, 2.1, 1.1))
> plot(fm1)
> par(opar)
> 
> 
> 
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("morley")
> ### * morley
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: morley
> ### Title: Michelson Speed of Light Data
> ### Aliases: morley
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> michelson <- transform(morley,
+                        Expt = factor(Expt), Run = factor(Run))
> xtabs(~ Expt + Run, data = michelson)  # 5 x 20 balanced (two-way)
    Run
Expt 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
   1 1 1 1 1 1 1 1 1 1  1  1  1  1  1  1  1  1  1  1  1
   2 1 1 1 1 1 1 1 1 1  1  1  1  1  1  1  1  1  1  1  1
   3 1 1 1 1 1 1 1 1 1  1  1  1  1  1  1  1  1  1  1  1
   4 1 1 1 1 1 1 1 1 1  1  1  1  1  1  1  1  1  1  1  1
   5 1 1 1 1 1 1 1 1 1  1  1  1  1  1  1  1  1  1  1  1
> plot(Speed ~ Expt, data = michelson,
+      main = "Speed of Light Data", xlab = "Experiment No.")
> fm <- aov(Speed ~ Run + Expt, data = michelson)
> summary(fm)
            Df Sum Sq Mean Sq F value  Pr(>F)   
Run         19 113344    5965   1.105 0.36321   
Expt         4  94514   23629   4.378 0.00307 **
Residuals   76 410166    5397                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> fm0 <- update(fm, . ~ . - Run)
> anova(fm0, fm)
Analysis of Variance Table

Model 1: Speed ~ Expt
Model 2: Speed ~ Run + Expt
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     95 523510                           
2     76 410166 19    113344 1.1053 0.3632
> 
> 
> 
> cleanEx()
> nameEx("mtcars")
> ### * mtcars
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: mtcars
> ### Title: Motor Trend Car Road Tests
> ### Aliases: mtcars
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> pairs(mtcars, main = "mtcars data", gap = 1/4)
> coplot(mpg ~ disp | as.factor(cyl), data = mtcars,
+        panel = panel.smooth, rows = 1)
> ## possibly more meaningful, e.g., for summary() or bivariate plots:
> mtcars2 <- within(mtcars, {
+    vs <- factor(vs, labels = c("V", "S"))
+    am <- factor(am, labels = c("automatic", "manual"))
+    cyl  <- ordered(cyl)
+    gear <- ordered(gear)
+    carb <- ordered(carb)
+ })
> summary(mtcars2)
      mpg        cyl         disp             hp             drat      
 Min.   :10.40   4:11   Min.   : 71.1   Min.   : 52.0   Min.   :2.760  
 1st Qu.:15.43   6: 7   1st Qu.:120.8   1st Qu.: 96.5   1st Qu.:3.080  
 Median :19.20   8:14   Median :196.3   Median :123.0   Median :3.695  
 Mean   :20.09          Mean   :230.7   Mean   :146.7   Mean   :3.597  
 3rd Qu.:22.80          3rd Qu.:326.0   3rd Qu.:180.0   3rd Qu.:3.920  
 Max.   :33.90          Max.   :472.0   Max.   :335.0   Max.   :4.930  
       wt             qsec       vs             am     gear   carb  
 Min.   :1.513   Min.   :14.50   V:18   automatic:19   3:15   1: 7  
 1st Qu.:2.581   1st Qu.:16.89   S:14   manual   :13   4:12   2:10  
 Median :3.325   Median :17.71                         5: 5   3: 3  
 Mean   :3.217   Mean   :17.85                                4:10  
 3rd Qu.:3.610   3rd Qu.:18.90                                6: 1  
 Max.   :5.424   Max.   :22.90                                8: 1  
> 
> 
> 
> cleanEx()
> nameEx("nhtemp")
> ### * nhtemp
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: nhtemp
> ### Title: Average Yearly Temperatures in New Haven
> ### Aliases: nhtemp
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> plot(nhtemp, main = "nhtemp data",
+   ylab = "Mean annual temperature in New Haven, CT (deg. F)")
> 
> 
> 
> cleanEx()
> nameEx("nottem")
> ### * nottem
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: nottem
> ### Title: Average Monthly Temperatures at Nottingham, 1920-1939
> ### Aliases: nottem
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> nott <- window(nottem, end = c(1936,12))
> fit <- arima(nott, order = c(1,0,0), list(order = c(2,1,0), period = 12))
> nott.fore <- predict(fit, n.ahead = 36)
> ts.plot(nott, nott.fore$pred, nott.fore$pred+2*nott.fore$se,
+         nott.fore$pred-2*nott.fore$se, gpars = list(col = c(1,1,4,4)))
> 
> 
> 
> cleanEx()
> nameEx("npk")
> ### * npk
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: npk
> ### Title: Classical N, P, K Factorial Experiment
> ### Aliases: npk
> ### Keywords: datasets
> 
> ### ** Examples
> 
> 
> cleanEx()
> nameEx("occupationalStatus")
> ### * occupationalStatus
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: occupationalStatus
> ### Title: Occupational Status of Fathers and their Sons
> ### Aliases: occupationalStatus
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> 
> plot(occupationalStatus)
> 
> ##  Fit a uniform association model separating diagonal effects
> Diag <- as.factor(diag(1:8))
> Rscore <- scale(as.numeric(row(occupationalStatus)), scale = FALSE)
> Cscore <- scale(as.numeric(col(occupationalStatus)), scale = FALSE)
> modUnif <- glm(Freq ~ origin + destination + Diag + Rscore:Cscore,
+                family = poisson, data = occupationalStatus)
> 
> summary(modUnif)

Call:
glm(formula = Freq ~ origin + destination + Diag + Rscore:Cscore, 
    family = poisson, data = occupationalStatus)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.568592   0.183358   3.101 0.001929 ** 
origin2       0.431314   0.149415   2.887 0.003893 ** 
origin3       1.461862   0.131141  11.147  < 2e-16 ***
origin4       1.788731   0.126588  14.130  < 2e-16 ***
origin5       0.441011   0.144844   3.045 0.002329 ** 
origin6       2.491735   0.121219  20.556  < 2e-16 ***
origin7       1.127536   0.129032   8.738  < 2e-16 ***
origin8       0.796445   0.131863   6.040 1.54e-09 ***
destination2  0.873202   0.166902   5.232 1.68e-07 ***
destination3  1.813992   0.153861  11.790  < 2e-16 ***
destination4  2.082515   0.150887  13.802  < 2e-16 ***
destination5  1.366383   0.155590   8.782  < 2e-16 ***
destination6  2.816369   0.146054  19.283  < 2e-16 ***
destination7  1.903918   0.147810  12.881  < 2e-16 ***
destination8  1.398585   0.151658   9.222  < 2e-16 ***
Diag1         1.665495   0.237383   7.016 2.28e-12 ***
Diag2         0.959681   0.212122   4.524 6.06e-06 ***
Diag3         0.021750   0.156530   0.139 0.889490    
Diag4         0.226399   0.124364   1.820 0.068689 .  
Diag5         0.808646   0.229754   3.520 0.000432 ***
Diag6         0.132277   0.077191   1.714 0.086597 .  
Diag7         0.506709   0.115936   4.371 1.24e-05 ***
Diag8         0.221880   0.134803   1.646 0.099771 .  
Rscore:Cscore 0.136974   0.007489  18.289  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 4679.004  on 63  degrees of freedom
Residual deviance:   58.436  on 40  degrees of freedom
AIC: 428.78

Number of Fisher Scoring iterations: 4

> plot(modUnif) # 4 plots, with warning about  h_ii ~= 1
Warning: not plotting observations with leverage one:
  1, 10, 19, 28, 37, 46, 55, 64
> 
> 
> 
> cleanEx()
> nameEx("penguins")
> ### * penguins
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: penguins
> ### Title: Measurements of Penguins near Palmer Station, Antarctica
> ### Aliases: penguins penguins_raw
> ### Keywords: datasets
> 
> ### ** Examples
> 
> ## view summaries
> summary(penguins)
      species          island       bill_len        bill_dep    
 Adelie   :152   Biscoe   :168   Min.   :32.10   Min.   :13.10  
 Chinstrap: 68   Dream    :124   1st Qu.:39.23   1st Qu.:15.60  
 Gentoo   :124   Torgersen: 52   Median :44.45   Median :17.30  
                                 Mean   :43.92   Mean   :17.15  
                                 3rd Qu.:48.50   3rd Qu.:18.70  
                                 Max.   :59.60   Max.   :21.50  
                                 NA's   :2       NA's   :2      
  flipper_len      body_mass        sex           year     
 Min.   :172.0   Min.   :2700   female:165   Min.   :2007  
 1st Qu.:190.0   1st Qu.:3550   male  :168   1st Qu.:2007  
 Median :197.0   Median :4050   NA's  : 11   Median :2008  
 Mean   :200.9   Mean   :4202                Mean   :2008  
 3rd Qu.:213.0   3rd Qu.:4750                3rd Qu.:2009  
 Max.   :231.0   Max.   :6300                Max.   :2009  
 NA's   :2       NA's   :2                                 
> summary(penguins_raw) # not useful for character vectors
  studyName         Sample Number      Species             Region         
 Length:344         Min.   :  1.00   Length:344         Length:344        
 Class :character   1st Qu.: 29.00   Class :character   Class :character  
 Mode  :character   Median : 58.00   Mode  :character   Mode  :character  
                    Mean   : 63.15                                        
                    3rd Qu.: 95.25                                        
                    Max.   :152.00                                        
                                                                          
    Island             Stage           Individual ID      Clutch Completion 
 Length:344         Length:344         Length:344         Length:344        
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
                                                                            
    Date Egg          Culmen Length (mm) Culmen Depth (mm) Flipper Length (mm)
 Min.   :2007-11-09   Min.   :32.10      Min.   :13.10     Min.   :172.0      
 1st Qu.:2007-11-28   1st Qu.:39.23      1st Qu.:15.60     1st Qu.:190.0      
 Median :2008-11-09   Median :44.45      Median :17.30     Median :197.0      
 Mean   :2008-11-27   Mean   :43.92      Mean   :17.15     Mean   :200.9      
 3rd Qu.:2009-11-16   3rd Qu.:48.50      3rd Qu.:18.70     3rd Qu.:213.0      
 Max.   :2009-12-01   Max.   :59.60      Max.   :21.50     Max.   :231.0      
                      NA's   :2          NA's   :2         NA's   :2          
 Body Mass (g)      Sex            Delta 15 N (o/oo) Delta 13 C (o/oo)
 Min.   :2700   Length:344         Min.   : 7.632    Min.   :-27.02   
 1st Qu.:3550   Class :character   1st Qu.: 8.300    1st Qu.:-26.32   
 Median :4050   Mode  :character   Median : 8.652    Median :-25.83   
 Mean   :4202                      Mean   : 8.733    Mean   :-25.69   
 3rd Qu.:4750                      3rd Qu.: 9.172    3rd Qu.:-25.06   
 Max.   :6300                      Max.   :10.025    Max.   :-23.79   
 NA's   :2                         NA's   :14        NA's   :13       
   Comments        
 Length:344        
 Class :character  
 Mode  :character  
                   
                   
                   
                   
> ## convert character vectors to factors first
> dFactor <- function(dat) {
+   dat[] <- lapply(dat, \(.) if (is.character(.)) as.factor(.) else .)
+   dat
+ }
> summary(dFactor(penguins_raw))
   studyName   Sample Number                                         Species   
 PAL0708:110   Min.   :  1.00   Adelie Penguin (Pygoscelis adeliae)      :152  
 PAL0809:114   1st Qu.: 29.00   Chinstrap penguin (Pygoscelis antarctica): 68  
 PAL0910:120   Median : 58.00   Gentoo penguin (Pygoscelis papua)        :124  
               Mean   : 63.15                                                  
               3rd Qu.: 95.25                                                  
               Max.   :152.00                                                  
                                                                               
    Region          Island                   Stage     Individual ID
 Anvers:344   Biscoe   :168   Adult, 1 Egg Stage:344   N13A1  :  3  
              Dream    :124                            N13A2  :  3  
              Torgersen: 52                            N18A1  :  3  
                                                       N18A2  :  3  
                                                       N21A1  :  3  
                                                       N21A2  :  3  
                                                       (Other):326  
 Clutch Completion    Date Egg          Culmen Length (mm) Culmen Depth (mm)
 No : 36           Min.   :2007-11-09   Min.   :32.10      Min.   :13.10    
 Yes:308           1st Qu.:2007-11-28   1st Qu.:39.23      1st Qu.:15.60    
                   Median :2008-11-09   Median :44.45      Median :17.30    
                   Mean   :2008-11-27   Mean   :43.92      Mean   :17.15    
                   3rd Qu.:2009-11-16   3rd Qu.:48.50      3rd Qu.:18.70    
                   Max.   :2009-12-01   Max.   :59.60      Max.   :21.50    
                                        NA's   :2          NA's   :2        
 Flipper Length (mm) Body Mass (g)      Sex      Delta 15 N (o/oo)
 Min.   :172.0       Min.   :2700   FEMALE:165   Min.   : 7.632   
 1st Qu.:190.0       1st Qu.:3550   MALE  :168   1st Qu.: 8.300   
 Median :197.0       Median :4050   NA's  : 11   Median : 8.652   
 Mean   :200.9       Mean   :4202                Mean   : 8.733   
 3rd Qu.:213.0       3rd Qu.:4750                3rd Qu.: 9.172   
 Max.   :231.0       Max.   :6300                Max.   :10.025   
 NA's   :2           NA's   :2                   NA's   :14       
 Delta 13 C (o/oo)                                  Comments  
 Min.   :-27.02    Nest never observed with full clutch.: 34  
 1st Qu.:-26.32    Not enough blood for isotopes.       :  7  
 Median :-25.83    Sexing primers did not amplify.      :  4  
 Mean   :-25.69    No blood sample obtained for sexing. :  2  
 3rd Qu.:-25.06    No blood sample obtained.            :  2  
 Max.   :-23.79    (Other)                              :  5  
 NA's   :13        NA's                                 :290  
> 
> ## visualise distribution across factors
> plot(island ~ species, data = penguins)
> plot(sex ~ interaction(island, species, sep = "\n"), data = penguins)
> 
> ## bill depth vs. length by species (color) and sex (symbol):
> ## positive correlations for all species, males tend to have bigger bills
> sym <- c(1, 16)
> pal <- c("darkorange","purple","cyan4")
> plot(bill_dep ~ bill_len, data = penguins, pch = sym[sex], col = pal[species])
> 
> ## simplified sex dimorphism analysis for Adelie species:
> ## proportion of males increases with several size measurements
> adelie <- subset(penguins, species == "Adelie")
> plot(sex ~ bill_len, data = adelie)
> plot(sex ~ bill_dep, data = adelie)
> plot(sex ~ body_mass, data = adelie)
> m <- glm(sex ~ bill_len + bill_dep + body_mass, data = adelie, family = binomial)
> summary(m)

Call:
glm(formula = sex ~ bill_len + bill_dep + body_mass, family = binomial, 
    data = adelie)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -85.088438  17.811711  -4.777 1.78e-06 ***
bill_len      0.840401   0.237656   3.536 0.000406 ***
bill_dep      1.305989   0.423911   3.081 0.002064 ** 
body_mass     0.007790   0.001969   3.957 7.59e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 202.399  on 145  degrees of freedom
Residual deviance:  54.024  on 142  degrees of freedom
  (6 observations deleted due to missingness)
AIC: 62.024

Number of Fisher Scoring iterations: 8

> 
> ## Produce the long variable names as from {palmerpenguins} pkg:
> long_nms <- sub("len", "length_mm",
+                 sub("dep","depth_mm",
+                     sub("mass", "mass_g", colnames(penguins))))
> ## compare long and short names:
> noquote(rbind(long_nms, nms = colnames(penguins)))
         [,1]    [,2]   [,3]           [,4]          [,5]             
long_nms species island bill_length_mm bill_depth_mm flipper_length_mm
nms      species island bill_len       bill_dep      flipper_len      
         [,6]        [,7] [,8]
long_nms body_mass_g sex  year
nms      body_mass   sex  year
> 
> ## Not run: 
> ##D  # << keeping shorter 'penguins' names in this example:
> ##D     colnames(penguins) <- long_nms
> ## End(Not run)
> 
> 
> 
> cleanEx()
> nameEx("precip")
> ### * precip
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: precip
> ### Title: Annual Precipitation in US Cities
> ### Aliases: precip
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> dotchart(precip[order(precip)], main = "precip data")
> title(sub = "Average annual precipitation (in.)")
> 
> ## Old ("wrong") version of dataset (just name change):
> precip.O <- local({
+    p <- precip; names(p)[names(p) == "Cincinnati"] <- "Cincinati" ; p })
> stopifnot(all(precip == precip.O),
+ 	  match("Cincinnati", names(precip)) == 46,
+ 	  identical(names(precip)[-46], names(precip.O)[-46]))
> 
> 
> 
> cleanEx()
> nameEx("presidents")
> ### * presidents
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: presidents
> ### Title: Quarterly Approval Ratings of US Presidents
> ### Aliases: presidents
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> plot(presidents, las = 1, ylab = "Approval rating (%)",
+      main = "presidents data")
> 
> 
> 
> cleanEx()
> nameEx("pressure")
> ### * pressure
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: pressure
> ### Title: Vapor Pressure of Mercury as a Function of Temperature
> ### Aliases: pressure
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> plot(pressure, xlab = "Temperature (deg C)",
+      ylab = "Pressure (mm of Hg)",
+      main = "pressure data: Vapor Pressure of Mercury")
> plot(pressure, xlab = "Temperature (deg C)",  log = "y",
+      ylab = "Pressure (mm of Hg)",
+      main = "pressure data: Vapor Pressure of Mercury")
> 
> 
> 
> cleanEx()
> nameEx("quakes")
> ### * quakes
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: quakes
> ### Title: Locations of Earthquakes off Fiji
> ### Aliases: quakes
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> pairs(quakes, main = "Fiji Earthquakes, N = 1000", cex.main = 1.2, pch = ".")
> 
> 
> 
> cleanEx()
> nameEx("randu")
> ### * randu
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: randu
> ### Title: Random Numbers from Congruential Generator RANDU
> ### Aliases: randu
> ### Keywords: datasets
> 
> ### ** Examples
> 
> ## We could re-generate the dataset by the following R code
> seed <- as.double(1)
> RANDU <- function() {
+     seed <<- ((2^16 + 3) * seed) %% (2^31)
+     seed/(2^31)
+ }
> myrandu <- matrix(NA_real_, 400, 3, dimnames = list(NULL, c("x","y","z")))
> for(i in 1:400) {
+     U <- c(RANDU(), RANDU(), RANDU(), RANDU(), RANDU())
+     myrandu[i,] <- round(U[1:3], 6)
+ }
> stopifnot(all.equal(randu, as.data.frame(myrandu), tolerance = 1e-5))
> 
> 
> 
> cleanEx()
> nameEx("sleep")
> ### * sleep
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sleep
> ### Title: Student's Sleep Data
> ### Aliases: sleep
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats)
> ## Student's paired t-test
> with(sleep,
+      t.test(extra[group == 1],
+             extra[group == 2], paired = TRUE))

	Paired t-test

data:  extra[group == 1] and extra[group == 2]
t = -4.0621, df = 9, p-value = 0.002833
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 -2.4598858 -0.7001142
sample estimates:
mean difference 
          -1.58 

> 
> ## The sleep *prolongations*
> sleep1 <- with(sleep, extra[group == 2] - extra[group == 1])
> summary(sleep1)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00    1.05    1.30    1.58    1.70    4.60 
> stripchart(sleep1, method = "stack", xlab = "hours",
+            main = "Sleep prolongation (n = 10)")
> boxplot(sleep1, horizontal = TRUE, add = TRUE,
+         at = .6, pars = list(boxwex = 0.5, staplewex = 0.25))
> 
> 
> 
> cleanEx()
> nameEx("stackloss")
> ### * stackloss
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: stackloss
> ### Title: Brownlee's Stack Loss Plant Data
> ### Aliases: stackloss stack.loss stack.x
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats)
> summary(lm.stack <- lm(stack.loss ~ stack.x))

Call:
lm(formula = stack.loss ~ stack.x)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.2377 -1.7117 -0.4551  2.3614  5.6978 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -39.9197    11.8960  -3.356  0.00375 ** 
stack.xAir.Flow     0.7156     0.1349   5.307  5.8e-05 ***
stack.xWater.Temp   1.2953     0.3680   3.520  0.00263 ** 
stack.xAcid.Conc.  -0.1521     0.1563  -0.973  0.34405    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.243 on 17 degrees of freedom
Multiple R-squared:  0.9136,	Adjusted R-squared:  0.8983 
F-statistic:  59.9 on 3 and 17 DF,  p-value: 3.016e-09

> 
> 
> 
> cleanEx()
> nameEx("state")
> ### * state
> 
> flush(stderr()); flush(stdout())
> 
> ### Encoding: UTF-8
> ### Name: state
> ### Title: US State Facts and Figures
> ### Aliases: state state.abb state.area state.center state.division
> ###   state.name state.region state.x77
> ### Keywords: datasets
> 
> ### ** Examples
> 
> (dst <- dxy <- data.frame(state.center, row.names=state.abb))
           x       y
AL  -86.7509 32.5901
AK -127.2500 49.2500
AZ -111.6250 34.2192
AR  -92.2992 34.7336
CA -119.7730 36.5341
CO -105.5130 38.6777
CT  -72.3573 41.5928
DE  -74.9841 38.6777
FL  -81.6850 27.8744
GA  -83.3736 32.3329
HI -126.2500 31.7500
ID -113.9300 43.5648
IL  -89.3776 40.0495
IN  -86.0808 40.0495
IA  -93.3714 41.9358
KS  -98.1156 38.4204
KY  -84.7674 37.3915
LA  -92.2724 30.6181
ME  -68.9801 45.6226
MD  -76.6459 39.2778
MA  -71.5800 42.3645
MI  -84.6870 43.1361
MN  -94.6043 46.3943
MS  -89.8065 32.6758
MO  -92.5137 38.3347
MT -109.3200 46.8230
NE  -99.5898 41.3356
NV -116.8510 39.1063
NH  -71.3924 43.3934
NJ  -74.2336 39.9637
NM -105.9420 34.4764
NY  -75.1449 43.1361
NC  -78.4686 35.4195
ND -100.0990 47.2517
OH  -82.5963 40.2210
OK  -97.1239 35.5053
OR -120.0680 43.9078
PA  -77.4500 40.9069
RI  -71.1244 41.5928
SC  -80.5056 33.6190
SD  -99.7238 44.3365
TN  -86.4560 35.6767
TX  -98.7857 31.3897
UT -111.3300 39.1063
VT  -72.5450 44.2508
VA  -78.2005 37.5630
WA -119.7460 47.4231
WV  -80.6665 38.4204
WI  -89.9941 44.5937
WY -107.2560 43.0504
> ## Alaska and Hawaii are placed just off the West Coast (for compact map drawing):
> dst[c("AK", "HI"),]
         x     y
AK -127.25 49.25
HI -126.25 31.75
> ## state.center2 := version of state.center with "correct" coordinates for AK & HI:
> ## From https://pubs.usgs.gov/gip/Elevations-Distances/elvadist.html#Geographic%20Centers
> ##   Alaska   63°50' N., 152°00' W., 60 miles northwest of Mount McKinley
> ##   Hawaii   20°15' N., 156°20' W., off Maui Island
> dxy["AK",] <- c(-152.  , 63.83) # or  c(-152.11, 65.17)
> dxy["HI",] <- c(-156.33, 20.25) # or  c(-156.69, 20.89)
> state.center2 <- as.list(dxy)
> 
> plot(dxy, asp=1.2, pch=3, col=2)
> text(state.center2, state.abb, cex=1/2, pos=4, offset=1/4)
> i <- c("AK","HI")
> do.call(arrows, c(setNames(c(dst[i,], dxy[i,]), c("x0","y0", "x1","y1")),
+                   col=adjustcolor(4, .7), length=1/8))
> points(dst[i,], col=2)
> if(FALSE) { # if(require("maps")) {
+    map("state", interior = FALSE,          add = TRUE)
+    map("state", boundary = FALSE, lty = 2, add = TRUE)
+ }
> 
> 
> 
> cleanEx()
> nameEx("sunspot.month")
> ### * sunspot.month
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sunspot.month
> ### Title: Monthly Sunspot Data, from 1749 to "Present"
> ### Aliases: sunspot.m2014 sunspot.month
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> ## Compare the monthly series
> plot (sunspot.month,
+       main="sunspot.month & sunspots [package 'datasets']", col=2)
> lines(sunspots) # -> clear recalibration (to *larger* values)
> at. <- seq(1750, 2030, by=10)
> atyr <- at.[at. %% 50 != 0] #%
> Axis(time(sunspot.month), at = atyr, side = 1,
+      tck= -1/100, padj = -1.5, cex.axis = 3/4)
> 
> ## Now look at the difference :
> all(tsp(sunspots)     [c(1,3)] ==
+     tsp(sunspot.month)[c(1,3)]) ## Start & Periodicity are the same
[1] TRUE
> n1 <- length(sunspots)
> table(eq <- sunspots == sunspot.m2014[1:n1]) #>  143  are different !

FALSE  TRUE 
  143  2677 
> i <- which(!eq)
> rug(time(eq)[i])
> s1 <- sunspots[i] ; s2 <- sunspot.m2014[i]
> cbind(i = i, time = time(sunspots)[i], sunspots = s1, ss.month = s2,
+       perc.diff = round(100*2*abs(s1-s2)/(s1+s2), 1))
          i     time sunspots ss.month perc.diff
  [1,]   55 1753.500     22.2     22.0       0.9
  [2,]  838 1818.750     31.7     31.6       0.3
  [3,]  841 1819.000     32.5     32.8       0.9
  [4,]  862 1820.750      9.0      8.9       1.1
  [5,]  864 1820.917      9.7      9.1       6.4
  [6,]  866 1821.083      4.3      4.2       2.4
  [7,]  876 1821.917      0.0      0.2     200.0
  [8,]  901 1824.000     21.6     21.7       0.5
  [9,]  917 1825.333     15.4     15.5       0.6
 [10,]  920 1825.583     25.4     25.7       1.2
 [11,]  943 1827.500     42.9     42.3       1.4
 [12,]  946 1827.750     57.2     56.1       1.9
 [13,]  955 1828.500     54.3     54.2       0.2
 [14,]  960 1828.917     46.6     46.9       0.6
 [15,]  965 1829.333     67.5     67.4       0.1
 [16,]  968 1829.583     78.3     77.6       0.9
 [17,]  976 1830.250    107.1    106.3       0.7
 [18,]  988 1831.250     54.6     54.5       0.2
 [19,]  992 1831.583     54.9     55.0       0.2
 [20,]  994 1831.750     46.2     46.3       0.2
 [21,]  998 1832.083     55.5     55.6       0.2
 [22,] 1003 1832.500     13.9     14.0       0.7
 [23,] 1047 1836.167     98.1     98.2       0.1
 [24,] 1061 1837.333    111.3    111.7       0.4
 [25,] 1081 1839.000    107.6    105.6       1.9
 [26,] 1087 1839.500     84.7     84.8       0.1
 [27,] 1090 1839.750     90.8     90.9       0.1
 [28,] 1092 1839.917     63.6     63.7       0.2
 [29,] 1095 1840.167     55.5     67.8      20.0
 [30,] 1102 1840.750     49.8     55.0       9.9
 [31,] 1105 1841.000     24.0     24.1       0.4
 [32,] 1108 1841.250     42.6     40.2       5.8
 [33,] 1109 1841.333     67.4     67.5       0.1
 [34,] 1113 1841.667     35.1     36.5       3.9
 [35,] 1124 1842.583     26.5     26.6       0.4
 [36,] 1125 1842.667     18.5     18.4       0.5
 [37,] 1132 1843.250      8.8      9.5       7.7
 [38,] 1145 1844.333     12.0     11.6       3.4
 [39,] 1149 1844.667      6.9      7.0       1.4
 [40,] 1156 1845.250     56.9     57.0       0.2
 [41,] 1168 1846.250     69.2     69.3       0.1
 [42,] 1185 1847.667    161.2    160.9       0.2
 [43,] 1191 1848.167    108.9    108.6       0.3
 [44,] 1194 1848.417    123.8    129.0       4.1
 [45,] 1196 1848.583    132.5    132.6       0.1
 [46,] 1200 1848.917    159.9    159.5       0.3
 [47,] 1201 1849.000    156.7    157.0       0.2
 [48,] 1202 1849.083    131.7    131.8       0.1
 [49,] 1203 1849.167     96.5     96.2       0.3
 [50,] 1206 1849.417     81.2     81.1       0.1
 [51,] 1208 1849.583     61.3     67.7       9.9
 [52,] 1211 1849.833     99.7     99.0       0.7
 [53,] 1224 1850.917     60.0     61.0       1.7
 [54,] 1235 1851.833     50.9     51.0       0.2
 [55,] 1238 1852.083     67.5     66.4       1.6
 [56,] 1243 1852.500     42.0     42.1       0.2
 [57,] 1256 1853.583     50.4     50.5       0.2
 [58,] 1258 1853.750     42.3     42.4       0.2
 [59,] 1264 1854.250     26.4     26.5       0.4
 [60,] 1270 1854.750     12.7     12.6       0.8
 [61,] 1272 1854.917     21.4     21.6       0.9
 [62,] 1282 1855.750      9.7      9.6       1.0
 [63,] 1283 1855.833      4.3      4.2       2.4
 [64,] 1290 1856.417      5.0      5.2       3.9
 [65,] 1301 1857.333     29.2     28.5       2.4
 [66,] 1333 1860.000     81.5     82.4       1.1
 [67,] 1334 1860.083     88.0     88.3       0.3
 [68,] 1346 1861.083     77.8     77.7       0.1
 [69,] 1350 1861.417     87.8     88.1       0.3
 [70,] 1366 1862.750     42.0     41.9       0.2
 [71,] 1407 1866.167     24.6     24.5       0.4
 [72,] 1424 1867.583      4.9      4.8       2.1
 [73,] 1427 1867.833      9.3      9.6       3.2
 [74,] 1429 1868.000     15.6     15.5       0.6
 [75,] 1430 1868.083     15.8     15.7       0.6
 [76,] 1435 1868.500     28.6     29.0       1.4
 [77,] 1437 1868.667     43.8     47.2       7.5
 [78,] 1438 1868.750     61.7     61.6       0.2
 [79,] 1442 1869.083     59.3     59.9       1.0
 [80,] 1445 1869.333    104.0    103.9       0.1
 [81,] 1450 1869.750     59.4     59.3       0.2
 [82,] 1451 1869.833     77.4     78.1       0.9
 [83,] 1452 1869.917    104.3    104.4       0.1
 [84,] 1455 1870.167    159.4    157.5       1.2
 [85,] 1472 1871.583    110.0    110.1       0.1
 [86,] 1476 1871.917     90.3     90.4       0.1
 [87,] 1486 1872.750    103.5    102.6       0.9
 [88,] 1497 1873.667     47.5     47.1       0.8
 [89,] 1498 1873.750     47.4     47.1       0.6
 [90,] 1514 1875.083     22.2     21.5       3.2
 [91,] 1527 1876.167     31.2     30.6       1.9
 [92,] 1539 1877.167     11.7     11.9       1.7
 [93,] 1541 1877.333     21.2     21.6       1.9
 [94,] 1542 1877.417     13.4     14.2       5.8
 [95,] 1543 1877.500      5.9      6.0       1.7
 [96,] 1545 1877.667     16.4     16.9       3.0
 [97,] 1547 1877.833     14.5     14.2       2.1
 [98,] 1548 1877.917      2.3      2.2       4.4
 [99,] 1550 1878.083      6.0      6.6       9.5
[100,] 1553 1878.333      5.8      5.9       1.7
[101,] 1561 1879.000      0.8      1.0      22.2
[102,] 1571 1879.833     12.9     13.1       1.5
[103,] 1572 1879.917      7.2      7.3       1.4
[104,] 1574 1880.083     27.5     27.2       1.1
[105,] 1575 1880.167     19.5     19.3       1.0
[106,] 1576 1880.250     19.3     19.5       1.0
[107,] 1588 1881.250     51.7     51.6       0.2
[108,] 1592 1881.583     58.0     58.4       0.7
[109,] 1594 1881.750     64.0     64.4       0.6
[110,] 1598 1882.083     69.3     69.5       0.3
[111,] 1599 1882.167     67.5     66.8       1.0
[112,] 1613 1883.333     32.1     31.5       1.9
[113,] 1614 1883.417     76.5     76.3       0.3
[114,] 1623 1884.167     86.8     87.5       0.8
[115,] 1643 1885.833     33.3     30.9       7.5
[116,] 1656 1886.917     12.4     13.0       4.7
[117,] 1663 1887.500     23.3     23.4       0.4
[118,] 1683 1889.167      7.0      6.7       4.4
[119,] 1687 1889.500      9.7      9.4       3.1
[120,] 1712 1891.583     33.2     33.0       0.6
[121,] 1716 1891.917     32.3     32.5       0.6
[122,] 1723 1892.500     76.8     76.5       0.4
[123,] 1734 1893.417     88.2     89.9       1.9
[124,] 1735 1893.500     88.8     88.6       0.2
[125,] 1738 1893.750     79.7     80.0       0.4
[126,] 1774 1896.750     28.4     28.7       1.1
[127,] 1837 1902.000      5.2      5.5       5.6
[128,] 2126 1926.083     70.0     69.9       0.1
[129,] 2151 1928.167     85.4     85.5       0.1
[130,] 2153 1928.333     76.9     77.0       0.1
[131,] 2162 1929.083     64.1     62.8       2.0
[132,] 2174 1930.083     49.2     49.9       1.4
[133,] 2233 1935.000     18.9     18.6       1.6
[134,] 2315 1941.833     38.3     38.4       0.3
[135,] 2329 1943.000     12.4     12.5       0.8
[136,] 2378 1947.083    113.4    133.4      16.2
[137,] 2427 1951.167     59.9     55.9       6.9
[138,] 2498 1957.083    130.2    130.3       0.1
[139,] 2552 1961.583     55.9     55.8       0.2
[140,] 2556 1961.917     40.0     39.9       0.3
[141,] 2594 1965.083     14.2     14.3       0.7
[142,] 2790 1981.417     90.0     90.9       1.0
[143,] 2819 1983.833     33.3     33.4       0.3
> 
> ## How to recreate the "old" sunspot.month (R <= 3.0.3) =: sunspot.month.0
> .sunspot.diff <- cbind(
+     i = c(1202L, 1256L, 1258L, 1301L, 1407L, 1429L, 1452L, 1455L,
+           1663L, 2151L, 2329L, 2498L, 2594L, 2694L, 2819L),
+     res10 = c(1L, 1L, 1L, -1L, -1L, -1L, 1L, -1L,
+           1L, 1L, 1L, 1L, 1L, 20L, 1L))
> ssm0 <- sunspot.m2014[1:2988]
> with(as.data.frame(.sunspot.diff), ssm0[i] <<- ssm0[i] - res10/10)
> sunspot.month.0 <- ts(ssm0, start = 1749, frequency = 12)
> stopifnot(length(sunspot.month.0) == 2988)
> 
> 
> 
> cleanEx()
> nameEx("sunspot.year")
> ### * sunspot.year
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sunspot.year
> ### Title: Yearly Sunspot Data, 1700-1988
> ### Aliases: sunspot.year
> ### Keywords: datasets
> 
> ### ** Examples
> 
> utils::str(sm <- sunspots)# the monthly version we keep unchanged
 Time-Series [1:2820] from 1749 to 1984: 58 62.6 70 55.7 85 83.5 94.8 66.3 75.9 75.5 ...
> utils::str(sy <- sunspot.year)
 Time-Series [1:289] from 1700 to 1988: 5 11 16 23 36 58 29 20 10 8 ...
> ## The common time interval
> (t1 <- c(max(start(sm), start(sy)),     1)) # Jan 1749
[1] 1749    1
> (t2 <- c(min(  end(sm)[1],end(sy)[1]), 12)) # Dec 1983 (will not be updated!)
[1] 1983   12
> s.m <- window(sm, start=t1, end=t2)
> s.y <- window(sy, start=t1, end=t2[1])
> stopifnot(length(s.y) * 12 == length(s.m),
+           ## The yearly series *is* close to the averages of the monthly one:
+           all.equal(s.y, aggregate(s.m, FUN = mean), tolerance = 0.0020))
> ## NOTE: Strangely, correctly weighting the number of days per month
> ##       (using 28.25 for February) is *not* closer than the simple mean:
> ndays <- c(31, 28.25, rep(c(31,30, 31,30, 31), 2))
> all.equal(s.y, aggregate(s.m, FUN = mean))                     # 0.0013
[1] "Mean relative difference: 0.001312539"
> all.equal(s.y, aggregate(s.m, FUN = weighted.mean, w = ndays)) # 0.0017
[1] "Mean relative difference: 0.001692215"
> 
> 
> 
> cleanEx()
> nameEx("sunspots")
> ### * sunspots
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sunspots
> ### Title: Monthly Sunspot Numbers, 1749-1983
> ### Aliases: sunspots
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> plot(sunspots, main = "sunspots data", xlab = "Year",
+      ylab = "Monthly sunspot numbers")
> 
> 
> 
> cleanEx()
> nameEx("swiss")
> ### * swiss
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: swiss
> ### Title: Swiss Fertility and Socioeconomic Indicators (1888) Data
> ### Aliases: swiss
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> pairs(swiss, panel = panel.smooth, main = "swiss data",
+       col = 3 + (swiss$Catholic > 50))
> summary(lm(Fertility ~ . , data = swiss))

Call:
lm(formula = Fertility ~ ., data = swiss)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.2743  -5.2617   0.5032   4.1198  15.3213 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      66.91518   10.70604   6.250 1.91e-07 ***
Agriculture      -0.17211    0.07030  -2.448  0.01873 *  
Examination      -0.25801    0.25388  -1.016  0.31546    
Education        -0.87094    0.18303  -4.758 2.43e-05 ***
Catholic          0.10412    0.03526   2.953  0.00519 ** 
Infant.Mortality  1.07705    0.38172   2.822  0.00734 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.165 on 41 degrees of freedom
Multiple R-squared:  0.7067,	Adjusted R-squared:  0.671 
F-statistic: 19.76 on 5 and 41 DF,  p-value: 5.594e-10

> 
> 
> 
> cleanEx()
> nameEx("trees")
> ### * trees
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: trees
> ### Title: Diameter, Height and Volume for Black Cherry Trees
> ### Aliases: trees
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> pairs(trees, panel = panel.smooth, main = "trees data")
> plot(Volume ~ Girth, data = trees, log = "xy")
> coplot(log(Volume) ~ log(Girth) | Height, data = trees,
+        panel = panel.smooth)
> summary(fm1 <- lm(log(Volume) ~ log(Girth), data = trees))

Call:
lm(formula = log(Volume) ~ log(Girth), data = trees)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.205999 -0.068702  0.001011  0.072585  0.247963 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -2.35332    0.23066  -10.20 4.18e-11 ***
log(Girth)   2.19997    0.08983   24.49  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.115 on 29 degrees of freedom
Multiple R-squared:  0.9539,	Adjusted R-squared:  0.9523 
F-statistic: 599.7 on 1 and 29 DF,  p-value: < 2.2e-16

> summary(fm2 <- update(fm1, ~ . + log(Height), data = trees))

Call:
lm(formula = log(Volume) ~ log(Girth) + log(Height), data = trees)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.168561 -0.048488  0.002431  0.063637  0.129223 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -6.63162    0.79979  -8.292 5.06e-09 ***
log(Girth)   1.98265    0.07501  26.432  < 2e-16 ***
log(Height)  1.11712    0.20444   5.464 7.81e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.08139 on 28 degrees of freedom
Multiple R-squared:  0.9777,	Adjusted R-squared:  0.9761 
F-statistic: 613.2 on 2 and 28 DF,  p-value: < 2.2e-16

> step(fm2)
Start:  AIC=-152.69
log(Volume) ~ log(Girth) + log(Height)

              Df Sum of Sq    RSS      AIC
<none>                     0.1855 -152.685
- log(Height)  1    0.1978 0.3832 -132.185
- log(Girth)   1    4.6275 4.8130  -53.743

Call:
lm(formula = log(Volume) ~ log(Girth) + log(Height), data = trees)

Coefficients:
(Intercept)   log(Girth)  log(Height)  
     -6.632        1.983        1.117  

> ## i.e., Volume ~= c * Height * Girth^2  seems reasonable
> 
> 
> 
> cleanEx()
> nameEx("uspop")
> ### * uspop
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: uspop
> ### Title: Populations Recorded by the US Census
> ### Aliases: uspop
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> plot(uspop, log = "y", main = "uspop data", xlab = "Year",
+      ylab = "U.S. Population (millions)")
> 
> 
> 
> cleanEx()
> nameEx("volcano")
> ### * volcano
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: volcano
> ### Title: Topographic Information on Auckland's Maunga Whau Volcano
> ### Aliases: volcano
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(grDevices); require(graphics)
> filled.contour(volcano, color.palette = terrain.colors, asp = 1)
> title(main = "volcano data: filled contour map")
> 
> 
> 
> cleanEx()
> nameEx("warpbreaks")
> ### * warpbreaks
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: warpbreaks
> ### Title: The Number of Breaks in Yarn during Weaving
> ### Aliases: warpbreaks
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> summary(warpbreaks)
     breaks      wool   tension
 Min.   :10.00   A:27   L:18   
 1st Qu.:18.25   B:27   M:18   
 Median :26.00          H:18   
 Mean   :28.15                 
 3rd Qu.:34.00                 
 Max.   :70.00                 
> opar <- par(mfrow = c(1, 2), oma = c(0, 0, 1.1, 0))
> plot(breaks ~ tension, data = warpbreaks, col = "lightgray",
+      varwidth = TRUE, subset = wool == "A", main = "Wool A")
> plot(breaks ~ tension, data = warpbreaks, col = "lightgray",
+      varwidth = TRUE, subset = wool == "B", main = "Wool B")
> mtext("warpbreaks data", side = 3, outer = TRUE)
> par(opar)
> summary(fm1 <- lm(breaks ~ wool*tension, data = warpbreaks))

Call:
lm(formula = breaks ~ wool * tension, data = warpbreaks)

Residuals:
     Min       1Q   Median       3Q      Max 
-19.5556  -6.8889  -0.6667   7.1944  25.4444 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)      44.556      3.647  12.218 2.43e-16 ***
woolB           -16.333      5.157  -3.167 0.002677 ** 
tensionM        -20.556      5.157  -3.986 0.000228 ***
tensionH        -20.000      5.157  -3.878 0.000320 ***
woolB:tensionM   21.111      7.294   2.895 0.005698 ** 
woolB:tensionH   10.556      7.294   1.447 0.154327    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.94 on 48 degrees of freedom
Multiple R-squared:  0.3778,	Adjusted R-squared:  0.3129 
F-statistic: 5.828 on 5 and 48 DF,  p-value: 0.0002772

> anova(fm1)
Analysis of Variance Table

Response: breaks
             Df Sum Sq Mean Sq F value    Pr(>F)    
wool          1  450.7  450.67  3.7653 0.0582130 .  
tension       2 2034.3 1017.13  8.4980 0.0006926 ***
wool:tension  2 1002.8  501.39  4.1891 0.0210442 *  
Residuals    48 5745.1  119.69                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> 
> 
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("women")
> ### * women
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: women
> ### Title: Average Heights and Weights for American Women
> ### Aliases: women
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(graphics)
> plot(women, xlab = "Height (in)", ylab = "Weight (lb)",
+      main = "women data: American women aged 30-39")
> 
> 
> 
> cleanEx()
> nameEx("zCO2")
> ### * zCO2
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: CO2
> ### Title: Carbon Dioxide Uptake in Grass Plants
> ### Aliases: CO2
> ### Keywords: datasets
> 
> ### ** Examples
> 
> require(stats); require(graphics)
> ## Don't show: 
> options(show.nls.convergence=FALSE)
> ## End(Don't show)
> coplot(uptake ~ conc | Plant, data = CO2, show.given = FALSE, type = "b")
> ## fit the data for the first plant
> fm1 <- nls(uptake ~ SSasymp(conc, Asym, lrc, c0),
+    data = CO2, subset = Plant == "Qn1")
> summary(fm1)

Formula: uptake ~ SSasymp(conc, Asym, lrc, c0)

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
Asym  38.1398     0.9164  41.620 1.99e-06 ***
lrc  -34.2766    18.9661  -1.807    0.145    
c0    -4.3806     0.2042 -21.457 2.79e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.663 on 4 degrees of freedom

> ## fit each plant separately
> fmlist <- list()
> for (pp in levels(CO2$Plant)) {
+   fmlist[[pp]] <- nls(uptake ~ SSasymp(conc, Asym, lrc, c0),
+       data = CO2, subset = Plant == pp)
+ }
> ## check the coefficients by plant
> print(sapply(fmlist, coef), digits = 3)
        Qn1    Qn2    Qn3   Qc1    Qc3    Qc2    Mn3    Mn2   Mn1   Mc2     Mc3
Asym  38.14  42.87  44.23 36.43  40.68  39.82  28.48  32.13 34.08 13.56   18.54
lrc  -34.28 -29.66 -37.63 -9.90 -11.54 -51.53 -17.37 -29.04 -8.81 -1.98 -136.11
c0    -4.38  -4.67  -4.49 -4.86  -4.95  -4.46  -4.59  -4.47 -5.06 -4.56   -3.47
       Mc1
Asym 21.79
lrc   2.45
c0   -5.14
> 
> 
> 
> ### * <FOOTER>
> ###
> cleanEx()
> options(digits = 7L)
> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
Time elapsed:  2.05 0.135 2.209 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')