R version 3.1.1 RC (2014-07-04 r66081) -- "Sock it to Me"
Copyright (C) 2014 The R Foundation for Statistical Computing
Platform: x86_64-unknown-linux-gnu (64-bit)

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.

  Natural language support but running in an English locale

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')
> cleanEx()
> nameEx("AirPassengers")
> ### * AirPassengers
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: AirPassengers
> ### Title: Monthly Airline Passenger Numbers 1949-1960
> ### Aliases: AirPassengers
> ### Keywords: datasets
> 
> ### ** Examples
> 
> ## Not run: 
> ##D ## These are quite slow and so not run by example(AirPassengers)
> ##D 
> ##D ## The classic 'airline model', by full ML
> ##D (fit <- arima(log10(AirPassengers), c(0, 1, 1),
> ##D               seasonal = list(order = c(0, 1, 1), period = 12)))
> ##D update(fit, method = "CSS")
> ##D update(fit, x = window(log10(AirPassengers), start = 1954))
> ##D pred <- predict(fit, n.ahead = 24)
> ##D tl <- pred$pred - 1.96 * pred$se
> ##D tu <- pred$pred + 1.96 * pred$se
> ##D ts.plot(AirPassengers, 10^tl, 10^tu, log = "y", lty = c(1, 2, 2))
> ##D 
> ##D ## full ML fit is the same if the series is reversed, CSS fit is not
> ##D ap0 <- rev(log10(AirPassengers))
> ##D attributes(ap0) <- attributes(AirPassengers)
> ##D arima(ap0, c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))
> ##D arima(ap0, c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12),
> ##D       method = "CSS")
> ##D 
> ##D ## Structural Time Series
> ##D ap <- log10(AirPassengers) - 2
> ##D (fit <- StructTS(ap, type = "BSM"))
> ##D par(mfrow = c(1, 2))
> ##D plot(cbind(ap, fitted(fit)), plot.type = "single")
> ##D plot(cbind(ap, tsSmooth(fit)), plot.type = "single")
> ## End(Not run)
> 
> 
> 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
> # using the plinear algorithm
> fm2 <- nls(demand ~ (1-exp(-exp(lrc)*Time)), data = BOD,
+    start = c(lrc = log(.35)), algorithm = "plinear", trace = TRUE)
32.946 :  -1.0498 22.1260
25.992 :  -0.62572 19.10319
25.99 :  -0.6327 19.1419
25.99 :  -0.63282 19.14256
> # 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
> 
> 
> 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 bxp(list(stats = c(4.17, 4.53, 5.155, 5.33, 6.11, 3.59, 4.17,  :
  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
> 
> require(graphics)
> pairs(USArrests, panel = panel.smooth, main = "USArrests data")
> 
> 
> 
> 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 bxp(list(stats = c(216, 271.5, 342, 373.5, 404, 108, 136, 151.5,  :
  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)] # 3 of them
  73  312  674 
-Inf  Inf -Inf 
> 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.   : 1.00  
 35-44:15   40-79    :23   10-19   :24   1st Qu.: 0.000   1st Qu.: 3.00  
 45-54:16   80-119   :21   20-29   :20   Median : 1.000   Median : 6.00  
 55-64:16   120+     :21   30+     :20   Mean   : 2.273   Mean   :11.08  
 65-74:15                                3rd Qu.: 4.000   3rd Qu.:14.00  
 75+  :11                                Max.   :17.000   Max.   :60.00  
> ## 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
NULL                           87    227.241
agegp        5   88.128        82    139.112
tobgp        3   19.085        79    120.028
alcgp        3   66.054        76     53.973
tobgp:alcgp  9    6.489        67     47.484
> ## 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)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.7628  -0.6426  -0.2709   0.3043   2.0421  

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -4.01097    0.31224 -12.846  < 2e-16 ***
agegp.L         2.96113    0.65092   4.549 5.39e-06 ***
agegp.Q        -1.33735    0.58918  -2.270  0.02322 *  
agegp.C         0.15292    0.44792   0.341  0.73281    
agegp^4         0.06668    0.30776   0.217  0.82848    
agegp^5        -0.20288    0.19523  -1.039  0.29872    
unclass(tobgp)  0.26162    0.08198   3.191  0.00142 ** 
unclass(alcgp)  0.65308    0.08452   7.727 1.10e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 227.241  on 87  degrees of freedom
Residual deviance:  59.277  on 80  degrees of freedom
AIC: 222.76

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("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)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6678  -0.8360  -0.5772   0.9030   1.9362  

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)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.7603  -0.8162  -0.4956   0.8349   2.6536  

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
> 
> 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]]))))
> all.equal(ii, iris) # TRUE
[1] TRUE
> 
> 
> 
> 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 is 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")
> coplot(mpg ~ disp | as.factor(cyl), data = mtcars,
+        panel = panel.smooth, rows = 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
> 
> 
> 
> cleanEx()
> nameEx("npk")
> ### * npk
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: npk
> ### Title: Classical N, P, K Factorial Experiment
> ### Aliases: npk
> ### Keywords: datasets
> 
> ### ** Examples
> 
> 
> base::options(contrasts = c(unordered = "contr.treatment",ordered = "contr.poly"))
> 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)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.6521  -0.6267   0.0000   0.5913   2.0964  

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
Warning: not plotting observations with leverage one:
  1, 10, 19, 28, 37, 46, 55, 64
> 
> 
> 
> 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.)")
> 
> 
> 
> 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
> 
> 
> 
> 
> 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 difference in means is not equal to 0
95 percent confidence interval:
 -2.4598858 -0.7001142
sample estimates:
mean of the differences 
                  -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("sunspot.month")
> ### * sunspot.month
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sunspot.month
> ### Title: Monthly Sunspot Data, from 1749 to "Present"
> ### Aliases: 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) # -> faint differences where they overlap
> 
> ## 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.month[1:n1]) #>  132  are different !

FALSE  TRUE 
  143  2677 
> i <- which(!eq)
> rug(time(eq)[i])
> s1 <- sunspots[i] ; s2 <- sunspot.month[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.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.month[1:2988]
> with(as.data.frame(.sunspot.diff), ssm0[i] <<- ssm0[i] - res10/10)
> sunspot.month.0 <- ts(ssm0, start = 1749, frequency = 12)
> 
> 
> 
> 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
[1] 1983   12
> s.m <- window(sm, start=t1, end=t2)
> s.y <- window(sy, start=t1, end=t2[1]) # {irrelevant warning}
> 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), tol = 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: Girth, 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>
> ###
> options(digits = 7L)
> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
Time elapsed:  1.783 0.04 1.831 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')