                             M A C A N O V A   4.07

       An Interactive Program for Statistical Analysis and Matrix Algebra
            For information on major features, type 'help(macanova)'
          For information on linear models and GLM's, type 'help(glm)'
              For latest information on changes, type 'help(news)'
               For information on Unix version, type 'help(unix)'
                        Version of 03/25/99 (HP/UX [gcc])
             Type 'help(copyright)' for copyright and warranty info
        Copyright (C) 1994 - 1999 Gary W. Oehlert and Christopher Bingham

       MacAnova will automatically attempt to read undefined macros from a
           file.  See macros, addmacrofile, getmacros, and macroread()

Cmd> # File to test commands in Chapter 10 of MacAnova Users' Guide

Cmd> # Version of 980813

Cmd> # Set default configuration

Cmd> startTime <- gettime(quiet:T,keep:T)

Cmd> setoptions(default:T)

Cmd> DEGPERRAD <- 180/PI; E <- exp(1)

Cmd> PLATFORM <- vecread(string:VERSION,char:T)[5]

Cmd> MACINTOSH <- PLATFORM == "(Mac" || PLATFORM == "(Macintosh" ||\
	PLATFORM == "(Power"

Cmd> DOS <- PLATFORM == "(DOS" || PLATFORM == "(Win32s"

Cmd> if (MACINTOSH || DOS){DATAPATH <- DATAPATHS <- HOME;;}

Cmd> PATHSEP <- if (MACINTOSH){":"}else{"/"}

Cmd> DATAFILE <- "macanova.dat"; MACROFILE <-  "macanova.mac"

Cmd> MACROFILES <- vector(MACROFILE,"tser.mac","design.mac")

Cmd> setoptions(errors:20,dumbplot:T,height:25,width:80)

Cmd> #10.  Examples

Cmd> #10.1 Introduction

Cmd> #10.2  Simple descriptive statistics

Cmd> rem <- vector(.5,1,2.4,5,15,19,48,83)

Cmd> rem <- vector(rem,5,5.5,6.7,13.5,31,40,47,47,59,62,68,72,\
78,84,89,105,180)

Cmd> secretion <- factor(vector(rep(2,8),rep(1,17)))

Cmd> boxplot(split(rem,secretion),\
xlab:"Variable rem split by secretion",\
title:"Data from Example 3.19, page 100 of Devore & Peck, Ed. 2")
                Data from Example 3.19, page 100 of Devore & Peck, Ed. 2
           +-------+------+------+------+-------+------+------+------+-------++
           |                                                                  |
           |                                                                  |
           |                                                                  |
           |                                                                  |
           |                                                                  |
           |                                                                  |
 B         |+--+-------+                                                      |
 o         ||  |       |                                                      |
 x        2||  |       |-----+           *                                    |
           |+--+-------+                                                      |
 N         |                                                                  |
 u         |           +---------+------+                                     |
 m        1| +---------|         |      |--------+                           *|
 b         |           |         |      |                                     |
 e         |           +---------+------+                                     |
 r         |                                                                  |
           |                                                                  |
           |                                                                  |
           |                                                                  |
           |                                                                  |
           +-------+------+------+------+-------+------+------+------+-------++
           0      20     40     60     80      100    120    140    160    180
                             Variable rem split by secretion

Cmd> stemleaf(rem[secretion==1], title:\
"Data from Example 3.19 Devore & Peck 2nd Ed., normal secretion")
Data from Example 3.19 Devore & Peck 2nd Ed., normal secretion
    3    +0|556
    4     1|3
    4     2|
    5     3|1
    8     4|077
  ( 1)    5|9
    8     6|28
    6     7|28
    4     8|49
    2     9|
    2    10|5
 High  180
          1|1 represents 11  Leaf digit unit = 1

Cmd> stemleaf(rem[secretion==2],title:\
"Data from Example 3.19 Devore & Peck 2nd Ed., hyper secretion")
Data from Example 3.19 Devore & Peck 2nd Ed., hyper secretion
    4    +0|0125
    4     1|59
    2     2|
    2     3|
    2     4|8
 High  83
          1|1 represents 11  Leaf digit unit = 1

Cmd> describe(rem[secretion==1]) # stats for hyper secretion
component: n
(1)           17
component: min
(1)            5
component: q1
(1)           31
component: median
(1)           59
component: q3
(1)           78
component: max
(1)          180
component: mean
(1)       58.394
component: var
(1)         1932

Cmd> tt <- t2val(rem[secretion==1],rem[secretion==2])#2 sample t

Cmd> vector(tt,1 - cumstu(tt,8+17-2),2*(1 - cumstu(abs(tt),8+17-2)))
(1)       2.1333     0.021897     0.043793

Cmd> #10.3 Simple regression

Cmd> x <- vector(.3,1,5,10,20,25,25,30,40,50,75,100) # distance

Cmd> y <- vector(62.75,37.51,29.70,20.71,17.65,15.41,14.15,\
13.50,12.11,11.40, 10.85,10.85) # lead content

Cmd> plot(Distance:x,Lead:y,title:\
"Data from Devore & Peck Ex. 18, p. 153") # It's pretty curved
                         Data from Devore & Peck Ex. 18, p. 153
           +------------+------------+------------+------------+-------------++
           |*                                                                 |
           |                                                                  |
         60+                                                                  +
           |                                                                  |
           |                                                                  |
         50+                                                                  +
           |                                                                  |
           |                                                                  |
           |                                                                  |
 L       40+                                                                  +
 e         |*                                                                 |
 a         |                                                                  |
 d       30+                                                                  +
           |   *                                                              |
           |                                                                  |
           |                                                                  |
         20+      *                                                           +
           |            *                                                     |
           |                *                                                 |
           |                   *     *      *               *                *|
         10+------------+------------+------------+------------+-------------++
           0           20           40           60           80           100
                                        Distance

Cmd> plot(log10(x),log10(y),\
xlab:"Log distance",ylab:"Log lead",\
title:"Log-log plot of Data from Devore & Peck Ex. 4.19, p. 177")
                Log-log plot of Data from Devore & Peck Ex. 4.19, p. 177
           ++------------+------------+------------+------------+------------++
        1.8+*            :                                                    +
           |             :                                                    |
           |             :                                                    |
        1.7+             :                                                    +
           |             :                                                    |
        1.6+             :                                                    +
           |             *                                                    |
 L         |             :                                                    |
 o      1.5+             :                                                    +
 g         |             :                 *                                  |
        1.4+             :                                                    +
 l         |             :                                                    |
 e         |             :                                                    |
 a      1.3+             :                         *                          +
 d         |             :                                *                   |
        1.2+             :                                                    +
           |             :                                   *                |
           |             :                                   * *              |
        1.1+             :                                        *           +
           |             :                                           *   *   *|
          1++------------+------------+------------+------------+------------++
          -0.5           0           0.5           1           1.5           2
                                      Log distance

Cmd> # log vs log plot is much straighter; regress log(y) on log(x)

Cmd> logx <- log10(x); logy <- log10(y)

Cmd> regress("logy=logx") # Do simple linear regression of logs
Model used is logy=logx
                 Coef       StdErr            t
CONSTANT       1.6225     0.020599       78.767
logx         -0.31472     0.015039      -20.928

N: 12, MSE: 0.0014366, DF: 10, R^2: 0.97768
Regression F(1,10): 437.97, Durbin-Watson: 1.8778
To see the ANOVA table type 'anova()'

Cmd> anova() # gives ANOVA table for preceding regression
Model used is logy=logx
WARNING: summaries are sequential
                DF           SS           MS
CONSTANT         1       18.969       18.969
logx             1      0.62921      0.62921
ERROR1          10     0.014366    0.0014366

Cmd> #Now use macros resvshyat & resvsrankits to look at residuals

Cmd> resvsyhat(title:\
"Residuals vs predicted, log10(lead) on log10(distance)")
                 Residuals vs predicted, log10(lead) on log10(distance)
           ++-------+-------+--------+-------+-------+-------+-------+--------+
          2+                                 *                                +
           |                                                                  |
 S         |                                                                  |
 t      1.5+                                                                  +
 u         |*                                                                 |
 d         |                                                                  |
 e        1+                  *                                               +
 n         |                                                                  |
 t         |                                                                  |
 i      0.5+                                                                 *+
 z         |               *         *                                        |
 e        0+...*..............................................................+
 d         |                                                                  |
           |                                                                  |
 R     -0.5+                                                                  +
 e         |             *                                                    |
 s         |       *       *                                                  |
 i       -1+          *                                                       +
 d         |                                                                  |
 s     -1.5+                                                   *              +
           ++-------+-------+--------+-------+-------+-------+-------+--------+
            1      1.1     1.2      1.3     1.4     1.5     1.6     1.7
                                          Yhat

Cmd> resvsrankits(title:\
"Residuals vs rankits, log10(lead) on log10(distance)")
                  Residuals vs rankits, log10(lead) on log10(distance)
           +--+---------+---------+---------+---------+---------+---------+---+
          2+                                :                                *+
           |                                :                                 |
 S         |                                :                                 |
 t      1.5+                                :                                 +
 u         |                                :                     *           |
 d         |                                :                                 |
 e        1+                                :               *                 +
 n         |                                :                                 |
 t         |                                :                                 |
 i      0.5+                                :          *                      +
 z         |                                : *   *                           |
 e        0+..............................*.:.................................+
 d         |                                :                                 |
           |                                :                                 |
 R     -0.5+                                :                                 +
 e         |                          *     :                                 |
 s         |                *    *          :                                 |
 i       -1+          *                     :                                 +
 d         |                                :                                 |
 s     -1.5+*                               :                                 +
           +--+---------+---------+---------+---------+---------+---------+---+
            -1.5       -1       -0.5        0        0.5        1        1.5
                                         Rankits

Cmd> tstats <- vector(coefs()/secoefs(coefs:F))

Cmd> twotailt(tstats,DF[3]) # Sec. 2.
(1)   2.6645e-15    1.377e-09

Cmd> fstat <- (SS[2]/DF[2])/(SS[3]/DF[3])

Cmd> vector(fstat,1 - cumF(fstat,DF[2],DF[3]))
(1)       437.97    1.377e-09

Cmd> regress("{log10(y)}={log10(x)}",pvals:T)
Model used is {log10(y)}={log10(x)}
                 Coef       StdErr            t      P-Value
CONSTANT       1.6225     0.020599       78.767   2.6645e-15
log10(x)     -0.31472     0.015039      -20.928    1.377e-09

N: 12, MSE: 0.0014366, DF: 10, R^2: 0.97768
Regression F(1,10): 437.97, P-value: 1.377e-09, Durbin-Watson: 1.8778
To see the ANOVA table type 'anova()'

Cmd> anova(,fstat:T) # to suppress P values, use pvals:F
Model used is {log10(y)}={log10(x)}
WARNING: summaries are sequential
                DF           SS           MS            F      P-value
CONSTANT         1       18.969       18.969  13203.46834            0
{log10(x)}       1      0.62921      0.62921    437.96893    1.377e-09
ERROR1          10     0.014366    0.0014366

Cmd> #setoptions(pvals:T,fstats:T)

Cmd> #10.4 Multiple linear regression

Cmd> iron <- vector(61,175,111,124,130,173,169,\
169,160,244,257,333,199)

Cmd> aluminum <- vector(13,21,24,23,64,38,33,61,39,71,112,88,54)

Cmd> adsorption <- vector(4,18,14,18,26,26,21,30,28,36,65,62,40)

Cmd> regress("adsorption = iron + aluminum")
Model used is adsorption = iron + aluminum
                 Coef       StdErr            t
CONSTANT      -7.3507       3.4847      -2.1094
iron          0.11273     0.029691       3.7969
aluminum        0.349     0.071306       4.8944

N: 13, MSE: 19.179, DF: 10, R^2: 0.94847
Regression F(2,10): 92.026, Durbin-Watson: 2.6339
To see the ANOVA table type 'anova()'

Cmd> anova() # display ANOVA table, too
Model used is adsorption = iron + aluminum
WARNING: summaries are sequential
                DF           SS           MS
CONSTANT         1        11580        11580
iron             1       3070.5       3070.5
aluminum         1       459.43       459.43
ERROR1          10       191.79       19.179

Cmd> print(DF,SS,ms:SS/DF)#print side effect variables; see Sec. 3.6
DF:
(1)            1            1            1           10
SS:
     CONSTANT         iron     aluminum       ERROR1
        11580       3070.5       459.43       191.79
ms:
     CONSTANT         iron     aluminum       ERROR1
        11580       3070.5       459.43       19.179

Cmd> COEF # side effect variable containing coefficients
     CONSTANT         iron     aluminum
      -7.3507      0.11273        0.349

Cmd> regcoefs() # coefficients in readable format; see Sec. 3.13.1
                 Coef       StdErr            t
CONSTANT      -7.3507       3.4847      -2.1094
iron          0.11273     0.029691       3.7969
aluminum        0.349     0.071306       4.8944

Cmd> dfe <- DF[4] ; mse <- SS[4]/dfe;mse # compute error mean square
       ERROR1
       19.179

Cmd> # Predict absorption when Iron==200,Al==60; see Sec. 3..18

Cmd> regpred(vector(200,60))
component: estimate
(1)       36.136
component: SEest
(1)       1.3017
component: SEpred
(1)       4.5687

Cmd> HII # Look at leverages; see Sec. 3.6
 (1)      0.30754      0.27815      0.15449      0.14506      0.35893
 (6)      0.10313      0.12684      0.13338      0.08757      0.15238
(11)      0.53456      0.53144     0.086537

Cmd> studres <- RESIDUALS/sqrt(mse*(1 - HII))

Cmd> studres # internally studentized residuals
 (1)    -0.017302     -0.45867      0.11455      0.82601      -1.0383
 (6)      0.14126     -0.54206     -0.73346      0.88505      -2.2161
(11)       1.4359      0.36646       1.4504

Cmd> # Now compute externally studentized residuals = t-statistics

Cmd> studres*(sqrt((dfe-1)/(dfe-studres*studres)))
 (1)    -0.016414     -0.43978      0.10875      0.81181      -1.0428
 (6)      0.13414     -0.52196     -0.71533      0.87459      -2.9471
(11)        1.529      0.35002       1.5484

Cmd> studres*studres*HII/(1-HII)/3 # Cook's distance;3 = DF in model
 (1)   4.4318e-05     0.027022   0.00079928     0.038588      0.20119
 (6)   0.00076483     0.014227     0.027598      0.02506       0.2943
(11)      0.78938     0.050772      0.06643

Cmd> # Macro resid computes all these with one command; see Sec.3.17

Cmd> resid()
           Depvar    StdResids          HII     Cook's D      t-stats
(1)             4    -0.017302      0.30754   4.4318e-05    -0.016414
(2)            18     -0.45867      0.27815     0.027022     -0.43978
(3)            14      0.11455      0.15449   0.00079928      0.10875
(4)            18      0.82601      0.14506     0.038588      0.81181
(5)            26      -1.0383      0.35893      0.20119      -1.0428
(6)            26      0.14126      0.10313   0.00076483      0.13414
(7)            21     -0.54206      0.12684     0.014227     -0.52196
(8)            30     -0.73346      0.13338     0.027598     -0.71533
(9)            28      0.88505      0.08757      0.02506      0.87459
(10)           36      -2.2161      0.15238       0.2943      -2.9471
(11)           65       1.4359      0.53456      0.78938        1.529
(12)           62      0.36646      0.53144     0.050772      0.35002
(13)           40       1.4504     0.086537      0.06643       1.5484

Cmd> #10.5  One way ANOVA

Cmd> strength <- vector(7,7,15,11,9, 12,17,12,18,18,\
 14,18,18,19,19, 19,25,22,19,23, 7,10,11,15,11)

Cmd> percent <- rep(run(15,35,5),rep(5,5))

Cmd> grp <- makefactor(percent) # use macro to create factor

Cmd> print(format:"2.0f",percent,grp) # print using changed format
percent:
(1) 15 15 15 15 15 20 20 20 20 20 25 25 25 25 25 30 30 30 30 30 35 35 35 35 35
grp:
(1)  1  1  1  1  1  2  2  2  2  2  3  3  3  3  3  4  4  4  4  4  5  5  5  5  5

Cmd> boxplot(split(strength,grp),xlab:"Strength",title:\
"Fiber strength split by percent cotton (15, 20, 25, 30, 35)")
               Fiber strength split by percent cotton (15, 20, 25, 30, 35)
           +---+------+-------+------+------+------+------+-------+------+----+
           |                                                                  |
           |          +---+                                                   |
          5|*         |   |             o                                     |
           |          |   |                                                   |
           |          +---+                                                   |
           |                                           +----------+--+        |
 B        4|                                           |          |  |-------+|
 o         |                                           +----------+--+        |
 x         |                                                                  |
           |                                       +---+                      |
 N        3|                         o             |   |                      |
 u         |                                       +---+                      |
 m         |                                                                  |
 b         |                  +-----------------+--+                          |
 e        2|                  |                 |  |                          |
 r         |                  +-----------------+--+                          |
           |+------+------+                                                   |
           ||      |      |                                                   |
          1||      |      |-------------+                                     |
           |+------+------+                                                   |
           +---+------+-------+------+------+------+------+-------+------+----+
               8     10      12     14     16     18     20      22     24
                                        Strength

Cmd> # Clear differences in strength as percent cotton varies.

Cmd> # See Sec. 2.12.2 for boxplot()

Cmd> anova("strength=grp") # One-way ANOVA with 5 groups
Model used is strength=grp
                DF           SS           MS
CONSTANT         1         5655         5655
grp              4       475.76       118.94
ERROR1          20        161.2         8.06

Cmd> contrast("grp",vector(-2,-1,0,1,2)) # linear effect of percent
component: estimate
(1)          8.2
component: ss
(1)        33.62
component: se
(1)        4.015

Cmd> # See Sec. 3.16 for contrast()

Cmd> # Compare with the polynomial regression model in Sec. 10.6

Cmd> tabs(strength, grp) # Cell by cell statistics; see Sec. 3.12
component: mean
(1)          9.8         15.4         17.6         21.6         10.8
component: var
(1)         11.2          9.8          4.3          6.8          8.2
component: count
(1)            5            5            5            5            5

Cmd> #10.6  Polynomial regression

Cmd> percent2 <- percent*percent  # quadratic term

Cmd> percent3 <- percent2*percent # cubic term

Cmd> percent4 <- percent3*percent # quartic term

Cmd> anova("strength=percent+percent2+percent3+percent4",fstat:T)
Model used is strength=percent+percent2+percent3+percent4
WARNING: summaries are sequential
                DF           SS           MS            F      P-value
CONSTANT         1         5655         5655    701.61787            0
percent          1        33.62        33.62      4.17122     0.054524
percent2         1       343.21       343.21     42.58242   2.3255e-06
percent3         1        64.98        64.98      8.06203     0.010133
percent4         1       33.946       33.946      4.21163     0.053469
ERROR1          20        161.2         8.06

Cmd> regcoefs() # Coefficients, etc. see Sec. 3.13.1
                 Coef       StdErr            t
CONSTANT       -406.4       231.52      -1.7554
percent        73.777        40.63       1.8158
percent2      -4.8077       2.5851      -1.8598
percent3      0.13773     0.070868       1.9435
percent4   -0.0014533   0.00070817      -2.0522

Cmd> # Notice SS for 'percent' is linear contrast SS from above

Cmd> # Compute F-statistics and P-values "by hand"

Cmd> f <- (SS[run(2,5)]/DF[run(2,5)])/(SS[6]/DF[6]) # F-stats

Cmd> f ; 1-cumF(f,DF[run(2,5)],DF[6]); invF(1-.05,1,DF[6])
      percent     percent2     percent3     percent4
       4.1712       42.582        8.062       4.2116
(1)     0.054524   2.3255e-06     0.010133     0.053469
(1)       4.3512

Cmd> regress("strength=P3(percent)",pvals:T)
Model used is strength=P3(percent)
                    Coef       StdErr            t      P-Value
CONSTANT          62.611       39.757       1.5748      0.13024
percent          -9.0114       5.1966      -1.7341     0.097558
(percent)^2      0.48143      0.21605       2.2284     0.036915
(percent)^3      -0.0076     0.002874      -2.6444     0.015164

N: 25, MSE: 9.2927, DF: 21, R^2: 0.69363
Regression F(3,21): 15.848, P-value: 1.2953e-05, Durbin-Watson: 2.2203
To see the ANOVA table type 'anova()'

Cmd> anova(,fstat:T)
Model used is strength=P3(percent)
WARNING: summaries are sequential
                DF           SS           MS            F      P-value
CONSTANT         1         5655         5655    608.54957            0
{percent}        1        33.62        33.62      3.61791     0.070966
{(percent)^2}    1       343.21       343.21     36.93394   4.9628e-06
{(percent)^3}    1        64.98        64.98      6.99262     0.015164
ERROR1          21       195.15       9.2927

Cmd> #10.7 Variance stabilizing transformations

Cmd> discharge <- vector(.34,.12,1.23,.7,1.75,.12,.91,2.94,2.14,\
2.36,2.86,4.55,6.31,8.37,9.75,6.09,9.82,7.24,17.15,11.82,\
10.95,17.2,14.35,16.82)

Cmd> method <- factor(rep(run(4),rep(6,4)))#6 1's,6 2's,6 3's,6 4's

Cmd> anova("discharge=method")
Model used is discharge=method
                DF           SS           MS
CONSTANT         1       1012.6       1012.6
method           3       708.35       236.12
ERROR1          20       62.081       3.1041

Cmd> tabs(discharge,method) # Group statistics; see Sec. 3.12
component: mean
(1)         0.71       2.6267         7.93       14.715
component: var
(1)      0.43704       1.4213       2.7128        7.845
component: count
(1)            6            6            6            6

Cmd> # Let's look at residuals; see Sec.3.17

Cmd> resvsyhat(title:\
"Residuals vs predicted, Montgomery Example 4-2")
                     Residuals vs predicted, Montgomery Example 4-2
           +-----+---------+--------+--------+---------+--------+--------+----+
           |                                                                 *|
        1.5+                                                                 *+
 S         |        *                        *                               *|
 t        1+                                 *                                +
 u         |                                                                  |
 d         |*                                                                 |
 e      0.5+*                                                                 +
 n         |        *                        *                                |
 t        0+*.................................................................+
 i         |*       *                                                        *|
 z     -0.5+*       *                        *                                +
 e         |                                                                  |
 d         |                                                                  |
         -1+        *                        *                                +
 R         |                                 *                                |
 e     -1.5+                                                                  +
 s         |                                                                  |
 i         |                                                                 *|
 d       -2+                                                                  +
 s         |                                                                 *|
       -2.5+-----+---------+--------+--------+---------+--------+--------+----+
                 2         4        6        8        10       12       14
                                          Yhat

Cmd> @tmp <- tabs(discharge,method);means <- @tmp$mean;\
			vars <- @tmp$var

Cmd> plot(log(means),log(vars),xlab:"Log cell means",\
ylab:"Log cell vars",title:"Montgomery Example 4-2")
                                 Montgomery Example 4-2
           +-------+----------+---------+----------+----------+---------+-----+
           |       :                                                         *|
          2+       :                                                          +
           |       :                                                          |
           |       :                                                          |
 L      1.5+       :                                                          +
 o         |       :                                                          |
 g         |       :                                                          |
          1+       :                                           *              +
 c         |       :                                                          |
 e         |       :                                                          |
 l         |       :                                                          |
 l      0.5+       :                                                          +
           |       :                    *                                     |
 v         |       :                                                          |
 a        0+.......:..........................................................+
 r         |       :                                                          |
 s         |       :                                                          |
       -0.5+       :                                                          +
           |       :                                                          |
           |*      :                                                          |
           +-------+----------+---------+----------+----------+---------+-----+
                   0         0.5        1         1.5         2        2.5
                                     Log cell means

Cmd> regress("{log10(sqrt(vars))}={log10(means)}",pvals:T)
Model used is {log10(sqrt(vars))}={log10(means)}
                     Coef       StdErr            t      P-Value
CONSTANT         -0.12077     0.043562      -2.7724       0.1092
log10(means)      0.44647     0.056589       7.8897     0.015688

N: 4, MSE: 0.0032157, DF: 2, R^2: 0.96887
Regression F(1,2): 62.248, P-value: 0.015688, Durbin-Watson: 2.7556
To see the ANOVA table type 'anova()'

Cmd> anova("{discharge^.5}=method")
Model used is {discharge^.5}=method
                DF           SS           MS
CONSTANT         1       120.52       120.52
method           3       32.684       10.895
ERROR1          20       2.6884      0.13442

Cmd> tabs(discharge^.5,method) # Much better on this scale
component: mean
(1)      0.75742        1.582       2.8033       3.8208
component: var
(1)      0.16358      0.14879     0.085844      0.13947
component: count
(1)            6            6            6            6

Cmd> powers <- run(0,1,.05); ss <- 0*powers

Cmd> for(@i,run(length(powers))){
	@tmp <- boxcox(discharge,powers[@i])
	anova("@tmp=method",silent:T) # silent:T suppresses anova() output
	ss[@i] <- SS[3];;}

Cmd> lineplot(Power:powers,RSS:ss,\
title:"Residual SS vs Box-Cox power")
                              Residual SS vs Box-Cox power
           ++------------+------------+------------+------------+------------++
           |.                                                                 |
         90+:.                                                                +
           |:.                                                                |
           |: .                                                               |
           |: .                                                               |
         80+:  .                                                              +
           |:   .                                                             |
           |:    .                                                            |
         70+:    .                                                            +
 R         |:     .                                                           |
 S         |:      .                                                          |
 S       60+:       .                                                       ..+
           |:        ..                                                   ..  |
           |:          ..                                               ..    |
           |:            .                                            ..      |
         50+:             ..                                        ..        +
           |:               ..                                    ..          |
           |:                 ...                             ....            |
         40+:                    ....                   ......                +
           |:                        ...................                      |
           ++------------+------------+------------+------------+------------++
            0           0.2          0.4          0.6          0.8           1
                                          Power

Cmd> getmacros(boxcoxvec,quiet:T)

Cmd> help(boxcoxvec, file:"design.hlp") # see Sec. 8.6.1
boxcoxvec(rhs_model,y,powers:pow) computes the error SS for y
transformed to the boxcox powers given in pow and modeled using the
explanatory variables in CHARACTER scalar or quoted string rhs_model.
The returned value is a structure with named components 'power' and 'SS'
giving the boxcox powers and the error SS.

rhs_model is the right hand side of the anova model (the part following
'=') as a CHARACTER scalar, for example "x + a + b".  If powers:pow is
omitted, the default powers are run(-1,2,.25)

Examples:
   Cmd> y <- rnorm(20);ey <- exp(y)
   NOTE: random number seeds set to 59622139 and 172924584

   Cmd> a <- factor(rep(run(5),4))

   Cmd> boxcoxvec("a",ey)  # default powers
   component: power
    (1)           -1        -0.75         -0.5        -0.25            0
    (6)         0.25          0.5         0.75            1         1.25
   (11)          1.5         1.75            2
   component: SS
    (1)       179.93       72.889       34.971       21.487       17.979
    (6)       20.624       31.191       58.929       131.81       332.94
   (11)       915.71         2673       8140.8

   Cmd> boxcoxvec("a",ey,powers:run(-.25,.25,1/16))#explore powers ~= 0
   component: power
   (1)        -0.25      -0.1875       -0.125      -0.0625            0
   (6)       0.0625        0.125       0.1875         0.25
   component: SS
   (1)       21.487       19.938       18.875       18.236       17.979
   (6)       18.083       18.544        19.38       20.624


Cmd> boxcoxvec("method",discharge,powers:powers)
component: power
 (1)            0         0.05          0.1         0.15          0.2
 (6)         0.25          0.3         0.35          0.4         0.45
(11)          0.5         0.55          0.6         0.65          0.7
(16)         0.75          0.8         0.85          0.9         0.95
(21)            1
component: SS
 (1)       91.958       78.439       67.658       59.084       52.303
 (6)       46.989       42.891       39.817       37.617       36.179
(11)       35.424       35.294       35.756       36.792       38.404
(16)       40.607       43.432       46.923       51.142       56.163
(21)       62.081

Cmd> ss[run(5)] # some of values computed previously
(1)       91.958       78.439       67.658       59.084       52.303

Cmd> #10.8 Randomized complete blocks

Cmd> tiptype <-factor(rep(run(4),rep(4,4)))

Cmd> # tiptype is vector(1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4)

Cmd> specimen <- factor(rep(run(4),4))# vector(1,2,3,4,1,2,3,4,....)

Cmd> depth <- vector(9.3,9.4,9.6,10,9.4,9.3,9.8,9.9,9.2,9.4,\
9.5,9.7,9.7,9.6,10,10.2)

Cmd> anova("depth=specimen+tiptype")
Model used is depth=specimen+tiptype
                DF           SS           MS
CONSTANT         1       1482.2       1482.2
specimen         3        0.825        0.275
tiptype          3        0.385      0.12833
ERROR1           9         0.08    0.0088889

Cmd> anova("depth=tiptype+specimen") # different factor order
Model used is depth=tiptype+specimen
                DF           SS           MS
CONSTANT         1       1482.2       1482.2
tiptype          3        0.385      0.12833
specimen         3        0.825        0.275
ERROR1           9         0.08    0.0088889

Cmd> coefs("tiptype") # determine the treatment effects; Sec. 3.13
(1)        -0.05       -0.025       -0.175         0.25

Cmd> # These do not depend on order of factors in model

Cmd> # Compute a contrast in tiptype.  See Sec. 3.16

Cmd> contrast("tiptype",vector(-2,3,1,-2))
component: estimate
(1)        -0.65
component: ss
(1)     0.093889
component: se
(1)          0.2

Cmd> sum(vector(-2,3,1,-2)*coefs("tiptype")) # confirmation of value
(1)        -0.65

Cmd> contrast("tiptype",vector(-2,3,1,-2),"specimen")# Note by-var.
component: estimate
(1)         -0.6         -0.7         -0.3           -1
component: ss
(1)         0.02     0.027222        0.005     0.055556
component: se
(1)          0.4          0.4          0.4          0.4

Cmd> #10.8.1 Multiple comparisons.

Cmd> getmacros(pairedcomp)
pairedcomp   MACRO
) pairedcomp(termname,lev [,method:T] [,error:term] ) does paired comparisons
) of the means of the levels of the factor given in termname at level of
) signficance lev.  If only these two arguments are used, then the paired
) comparisons are done using the Bonferroni method.  Other techniques can be
) requested by using a method keyword.  Currently supported methods are: lsd
) (least significant difference), bsd (Bonferroni), snk (Student-Newman-
) Keuls), hsd (Tukey's honest significant difference, also called the
) studentized range procedure), regwb (regw with Bonferroni test), and regwr
) (regw with studentized range).
)
)
) By default, the error in these tests is taken from the last error term of
) the current model (the last line of the ANOVA table).  You may specify a
) different error term with a keyword error:term.  term may be a number,
) indicating a line in the ANOVA table, or the name of a term in the ANOVA
) table, for example, error:4 or error:"ERROR1".
)
) An alternative form is pairedcomp(termname,critval:val)  In this form, all
) paired t-test statistics are compared to the critical value given in val.
) The same output is printed as before.

Cmd> pairedcomp("tiptype",.05) # BSD at .05 level
  |    3   -0.175
  |    1    -0.05
  |    2   -0.025
       4     0.25

Cmd> pairedcomp("tiptype",.01,hsd:T) #HSD at .01 level
  |      3   -0.175
  |      1    -0.05
  | |    2   -0.025
    |    4     0.25

Cmd> # From Steel and Torrie, the 1% critical value for the

Cmd> # Studentized range with 4 groups and 9 error d.f. is 5.96

Cmd> pairedcomp("tiptype",critval:5.96) # HSD @ 1% level
  |      3   -0.175
  | |    1    -0.05
  | |    2   -0.025
    |    4     0.25

Cmd> usage(invstudrng,orig:T)
invstudrng(P, ngroup, errorDf [,epsilon:eps]), elements of P between 0
  and 1, elements of ngroup integers >= 2, elements of errorDf >= 1,
  eps > 0 small

Cmd> Q <- invstudrng(1-.01, 4,DF[4]); Q
(1)       5.9576

Cmd> pairedcomp("tiptype",critval:5.96/sqrt(2))
  |      3   -0.175
  |      1    -0.05
  | |    2   -0.025
    |    4     0.25

Cmd> effects <- coefs(tiptype);effects # effects of tiptype
(1)        -0.05       -0.025       -0.175         0.25

Cmd> stderr <- sqrt(2*(SS[4]/DF[4])/4)

Cmd> (effects[-1] - effects[1])/stderr # t-statistics
(1)        0.375       -1.875          4.5

Cmd> invstu(1 - .025/3, DF[4]) # 2 tail critical value
(1)       2.9333

Cmd> invdunnett(1-.05, 4, DF[4])# 4 groups, DF[4] error d.f.
(1)       2.8116

Cmd> #10.8.2 Randomized block with data missing

Cmd> depth[1] <- ?; anova("depth=specimen+tiptype")
Model used is depth=specimen+tiptype
WARNING: cases with missing values deleted
WARNING: summaries are sequential
                DF           SS           MS
CONSTANT         1       1395.9       1395.9
specimen         3      0.72567      0.24189
tiptype          3      0.37611      0.12537
ERROR1           8     0.075556    0.0094444

Cmd> contrast("tiptype",vector(-2,3,1,-2))
component: estimate
(1)     -0.69444
component: ss
(1)     0.097534
component: se
(1)       0.2161

Cmd> contrast("tiptype",vector(-2,3,1,-2),"specimen") #by-variable
component: estimate
(1)      MISSING         -0.7         -0.3           -1
component: ss
(1)      MISSING     0.027222        0.005     0.055556
component: se
(1)      MISSING      0.41231      0.41231      0.41231

Cmd> pairedcomp("tiptype",.05) # BSD at .05 level
  |    3   -0.181
  |    1  -0.0333
  |    2  -0.0306
       4    0.244

Cmd> anova("depth=tiptype+specimen")
Model used is depth=tiptype+specimen
WARNING: cases with missing values deleted
WARNING: summaries are sequential
                DF           SS           MS
CONSTANT         1       1395.9       1395.9
tiptype          3      0.37317      0.12439
specimen         3      0.72861      0.24287
ERROR1           8     0.075556    0.0094444

Cmd> #10.9  Latin squares

Cmd> force <- vector(24,20,19,24,24, 17,24,30,27,36, 18,38,26,27,21,\
	26,31,26,23,22, 22,30,20,29,31)

Cmd> operators <- factor(rep(run(5),5)) # 1,2,3,4,5,1,2,3,4,5,...

Cmd> batches <- factor(rep(run(5),rep(5,5))) # 1,1,1,1,1,2,2,2,...

Cmd> mix <- factor(vector(1,2,3,4,5, 2,3,4,5,1, 3,4,5,1,2,\
	4,5,1,2,3, 5,1,2,3,4))

Cmd> anova("force=batches+operators+mix",fstat:T)
Model used is force=batches+operators+mix
                DF           SS           MS            F      P-value
CONSTANT         1        16129        16129   1512.09375   5.3957e-14
batches          4           68           17      1.59375      0.23906
operators        4          150         37.5      3.51562     0.040373
mix              4          330         82.5      7.73438    0.0025365
ERROR1          12          128       10.667

Cmd> #10.10 Balanced incomplete blocks

Cmd> time <- vector(73,73,75,74,75,75,67,68,72,71,72,75)

Cmd> batches <- factor(rep(run(4),rep(3,4)))#vector(1,1,1,2,2,2,...)

Cmd> catalyst <- factor(vector(1,3,4,1,2,3,2,3,4,1,2,4))

Cmd> anova("time=batches+catalyst",fstat:T)
Model used is time=batches+catalyst
WARNING: summaries are sequential
                DF           SS           MS            F      P-value
CONSTANT         1        63075        63075  97038.46154   6.4698e-12
batches          3           55       18.333     28.20513    0.0014678
catalyst         3        22.75       7.5833     11.66667     0.010739
ERROR1           5         3.25         0.65

Cmd> secoefs("catalyst") # adjusted treatment effects
component: coefs
(1)       -1.125       -0.875         -0.5          2.5
component: se
(1)      0.42757      0.42757      0.42757      0.42757

Cmd> contrast("catalyst",vector(1,1,1,-3))#compare last with 1st 3
component: estimate
(1)          -10
component: ss
(1)       22.222
component: se
(1)       1.7103

Cmd> temp <- (catalyst==1) + (catalyst==2) + (catalyst==3) - \
3*(catalyst==4) # create vector equivalent to contrast

Cmd> print(format:"3.0f",temp)
temp:
(1)   1   1  -3   1   1   1   1   1  -3   1   1  -3

Cmd> anova("time=batches+temp+catalyst",fstat:T)
Model used is time=batches+temp+catalyst
WARNING: summaries are sequential
                DF           SS           MS            F      P-value
CONSTANT         1        63075        63075  97038.46154   6.4698e-12
batches          3           55       18.333     28.20513    0.0014678
temp             1       22.222       22.222     34.18803    0.0020716
catalyst         2      0.52778      0.26389      0.40598      0.68646
ERROR1           5         3.25         0.65

Cmd> 

Cmd> #10.11  Analysis of covariance

Cmd> strength <- vector(36,41,39,42,49,40,48,39,45,44,35,37,42,34,32)

Cmd> diameter <- vector(20,25,24,25,32,22,28,22,30,28,21,23,26,21,15)

Cmd> machine <-factor(rep(run(3),rep(5,3)))#vector(1,1,1,1,1,2,2,...,)

Cmd> anova("strength=machine",fstat:T)# analysis ignoring diameter
Model used is strength=machine
                DF           SS           MS            F      P-value
CONSTANT         1        24241        24241   1412.07379   8.1046e-14
machine          2        140.4         70.2      4.08932     0.044232
ERROR1          12          206       17.167

Cmd> predtable(seest:T) # unadjusted treatment means & std errors
component: estimate
(1)         41.4         43.2           36
component: SEest
(1)       1.8529       1.8529       1.8529

Cmd> secoefs("machine") # unadjusted machine effects
component: coefs
(1)          1.2            3         -4.2
component: se
(1)       1.5129       1.5129       1.5129

Cmd> contrast("machine",vector(1,-1,0)) # unadjusted contrast
component: estimate
(1)         -1.8
component: ss
(1)          8.1
component: se
(1)       2.6204

Cmd> contrast("machine",vector(1,1,-2)) # unadjusted contrast
component: estimate
(1)         12.6
component: ss
(1)        132.3
component: se
(1)       4.5387

Cmd> anova("strength=diameter+machine",fstat:T)
Model used is strength=diameter+machine
WARNING: summaries are sequential
                DF           SS           MS            F      P-value
CONSTANT         1        24241        24241   9527.89438            0
diameter         1       305.13       305.13    119.93304   2.9601e-07
machine          2       13.284       6.6419      2.61064      0.11808
ERROR1          11       27.986       2.5442

Cmd> predtable(seest:T) # adjusted treatment means and SE's
component: estimate
(1)       40.382       41.419       38.798
component: SEest
(1)      0.72363      0.74442      0.78788

Cmd> secoefs("machine"); secoefs("diameter")
component: coefs
(1)      0.18241       1.2192      -1.4016
component: se
(1)        0.595      0.62012      0.67167
component: coefs
(1)      0.95399
component: se
(1)      0.11405

Cmd> contrast("machine", vector(1,-1,0))
component: estimate
(1)      -1.0368
component: ss
(1)       2.6656
component: se
(1)       1.0129

Cmd> anova("strength=machine+diameter+diameter.machine",fstat:T)
Model used is strength=machine+diameter+diameter.machine
WARNING: summaries are sequential
                   DF           SS           MS            F      P-value
CONSTANT            1        24241        24241   8640.65457            0
machine             2        140.4         70.2     25.02306   0.00021073
diameter            1       178.01       178.01     63.45381   2.2905e-05
machine.diameter    2       2.7372       1.3686      0.48784      0.62929
ERROR1              9       25.249       2.8054

Cmd> coefs(diameter)# different from after strength=diameter+machine
(1)      0.94187

Cmd> coefs("machine.diameter") # deviations from .94187
(1,1)      0.16241
(2,1)     -0.08473
(3,1)    -0.077675

Cmd> anova("strength=machine+diameter.machine")
Model used is strength=machine+diameter.machine
WARNING: summaries are sequential
                   DF           SS           MS
CONSTANT            1        24241        24241
machine             2        140.4         70.2
machine.diameter    3       180.75        60.25
ERROR1              9       25.249       2.8054

Cmd> coefs("machine.diameter")
(1,1)       1.1043
(2,1)      0.85714
(3,1)       0.8642

Cmd> coefs("machine.diameter") - 0.94187 #deviations from ave slope
(1,1)      0.16241
(2,1)    -0.084727
(3,1)    -0.077672

Cmd> #10.12 Factorial models

Cmd> volume <- vector(-3,-1,0,1,5,4,-1,0,2,1,7,6,-1,0,2,3,7,9,1,1,\
6,5,10,11)

Cmd> percent <- factor(rep(rep(run(3),rep(2,3)),4))

Cmd> pressure <- factor(rep(run(2),rep(12,2)))

Cmd> speed <- factor(rep(rep(run(2),rep(6,2)),2))

Cmd> print(format:"2.0f",percent,pressure,speed)
percent:
(1)  1  1  2  2  3  3  1  1  2  2  3  3  1  1  2  2  3  3  1  1  2  2  3  3
pressure:
(1)  1  1  1  1  1  1  1  1  1  1  1  1  2  2  2  2  2  2  2  2  2  2  2  2
speed:
(1)  1  1  1  1  1  1  2  2  2  2  2  2  1  1  1  1  1  1  2  2  2  2  2  2

Cmd> anova("volume=percent*pressure*speed") # Full factorial model
Model used is volume=percent*pressure*speed
                         DF           SS           MS
CONSTANT                  1       234.38       234.38
percent                   2       252.75       126.38
pressure                  1       45.375       45.375
percent.pressure          2         5.25        2.625
speed                     1       22.042       22.042
percent.speed             2      0.58333      0.29167
pressure.speed            1       1.0417       1.0417
percent.pressure.speed    2       1.0833      0.54167
ERROR1                   12          8.5      0.70833

Cmd> anova("volume=percent*pressure*speed
	-percent.pressure.speed",fstat:T)
Model used is volume=percent*pressure*speed
	-percent.pressure.speed
                   DF           SS           MS            F      P-value
CONSTANT            1       234.38       234.38    342.39130    3.077e-11
percent             2       252.75       126.38    184.61739   8.6824e-11
pressure            1       45.375       45.375     66.28696   1.1157e-06
percent.pressure    2         5.25        2.625      3.83478     0.046983
speed               1       22.042       22.042     32.20000   5.7379e-05
percent.speed       2      0.58333      0.29167      0.42609      0.66125
pressure.speed      1       1.0417       1.0417      1.52174      0.23767
ERROR1             14       9.5833      0.68452

Cmd> contrast("pressure",vector(-1,1),"percent") # percent is byvar
component: estimate
(1)          1.5            3         3.75
component: ss
(1)          4.5           18       28.125
component: se
(1)      0.58503      0.58503      0.58503

Cmd> colplot(coefs("percent.pressure"),\
title:"Interaction plot of percent.pressure",\
xlab:"Level of factor percent",ylab:"Interaction effects")
                          Interaction plot of percent.pressure
           ++---------------+---------------+---------------+----------------++
           |1.                                                                |
 I      0.6+  ...                                                             +
 n         |     ..                                                        --2|
 t         |       ...                                               ------   |
 e      0.4+          ...                                       -----         +
 r         |             ..                               ------              |
 a         |               ...                       -----                    |
 c      0.2+                  ...              ------                         +
 t         |                     ..        -2--                               |
 i         |                       ...  ---                                   |
 o        0+.........................---......................................+
 n         |                      ---    ..                                   |
           |                   ---         .1..                               |
 e     -0.2+                ---                ......                         +
 f         |              --                         .....                    |
 f         |           ---                                ......              |
 e     -0.4+        ---                                         .....         +
 c         |     ---                                                 ......   |
 t         |  ---                                                          ..1|
 s     -0.6+2-                                                                +
           ++---------------+---------------+---------------+----------------++
            1              1.5              2              2.5               3
                                 Level of factor percent

Cmd> colplot(tabs(volume,percent,pressure)$mean,title:\
"Volume averaged over speed for each level of percent.pressure",\
xlab:"Level of factor percent",ylab:"Volume average")
              Volume averaged over speed for each level of percent.pressure
           ++---------------+---------------+---------------+----------------++
           |                                                                -2|
           |                                                             ---  |
           |                                                         ----     |
          8+                                                      ---         +
 V         |                                                   ---            |
 o         |                                               ----               |
 l        6+                                            ---                   +
 u         |                                         ---                   ..1|
 m         |                                     ----                  ....   |
 e         |                                  ---                  ....       |
          4+                              --2-                 ....           +
 a         |                          ----                 ....               |
 v         |                     -----                 ....                   |
 e         |                -----                  ....                       |
 r        2+            ----                   ....                           +
 a         |       -----                ....1..                               |
 g         |   ----             ........                                      |
 e        0+2--...............................................................+
           |    ........                                                      |
           |1...                                                              |
           ++---------------+---------------+---------------+----------------++
            1              1.5              2              2.5               3
                                 Level of factor percent

Cmd> coefs()[vector(2,3,4)] # get coefficients for 3 terms
percent
(1)       -3.625       -0.625         4.25
pressure
(1)       -1.375        1.375
percent.pressure
(1,1)        0.625       -0.625
(2,1)       -0.125        0.125
(3,1)         -0.5          0.5

Cmd> concoefs <- outer(vector(-1,0,1),vector(-1,1)); concoefs
(1,1)            1           -1
(2,1)            0            0
(3,1)           -1            1

Cmd> contrast("percent.pressure",concoefs)
component: estimate
(1)         2.25
component: ss
(1)       5.0625
component: se
(1)      0.82736

Cmd> #10.13  Factorial designs with confounding

Cmd> y <- vector(-3,2,2,1, 0,-1,-1,6, -1,0,3,5, 1,0,1,1)

Cmd> a <- factor(vector(1,2,2,1, 2,1,1,2, 1,1,2,2, 2,1,2,1))

Cmd> b <- factor(vector(1,2,1,2, 1,2,1,2, 1,1,2,2, 1,2,1,2))

Cmd> c <- factor(vector(1,1,2,2, 1,1,2,2, 1,2,1,2, 1,1,2,2))

Cmd> repl <- factor(rep(run(2),rep(8,2)))

Cmd> blk <- factor(rep(vector(1,2,1,2),rep(4,4)))

Cmd> print(format:"3.0f",hconcat(a,b,c,repl,blk,y)')
MATRIX:
(1,1)   1   2   2   1   2   1   1   2   1   1   2   2   2   1   2   1
(2,1)   1   2   1   2   1   2   1   2   1   1   2   2   1   2   1   2
(3,1)   1   1   2   2   1   1   2   2   1   2   1   2   1   1   2   2
(4,1)   1   1   1   1   1   1   1   1   2   2   2   2   2   2   2   2
(5,1)   1   1   1   1   2   2   2   2   1   1   1   1   2   2   2   2
(6,1)  -3   2   2   1   0  -1  -1   6  -1   0   3   5   1   0   1   1

Cmd> anova("y=repl/blk+a*b*c# blocks nested within replicates",\
	fstats:T)
Model used is y=repl/blk+a*b*c# blocks nested within replicates
WARNING: summaries are sequential
                DF           SS           MS            F      P-value
CONSTANT         1           16           16     21.33333    0.0057416
repl             1            1            1      1.33333       0.3004
repl.blk         2          2.5         1.25      1.66667      0.27885
a                1           36           36     48.00000   0.00096125
b                1        20.25        20.25     27.00000    0.0034782
a.b              1          0.5          0.5      0.66667      0.45135
c                1        12.25        12.25     16.33333    0.0099085
a.c              1         0.25         0.25      0.33333      0.58872
b.c              1            1            1      1.33333       0.3004
a.b.c            1          0.5          0.5      0.66667      0.45135
ERROR1           5         3.75         0.75

Cmd> secoefs(coefs:F)[vector(4,5,6)]^2 # variances of some effects
component: a
(1)     0.046875     0.046875
component: b
(1)     0.046875     0.046875
component: a.b
(1,1)      0.09375      0.09375
(2,1)      0.09375      0.09375

Cmd> #10.14  Fractional factorial designs

Cmd> a <- factor(rep(run(2),4))

Cmd> b <- factor(rep(rep(run(2),rep(2,2)),2))

Cmd> c <- factor(rep(run(2),rep(4,2)))

Cmd> d <- factor(vector(1,2,2,1,2,1,1,2))

Cmd> y <- vector(45,100,45,65,75,60,80,96)

Cmd> print(format:"3.0f",hconcat(a,b,c,d,y)')
MATRIX:
(1,1)   1   2   1   2   1   2   1   2
(2,1)   1   1   2   2   1   1   2   2
(3,1)   1   1   1   1   2   2   2   2
(4,1)   1   2   2   1   2   1   1   2
(5,1)  45 100  45  65  75  60  80  96

Cmd> anova("y=a+b+c+d # Main effects model",fstat:T)
Model used is y=a+b+c+d # Main effects model
                DF           SS           MS            F      P-value
CONSTANT         1        40044        40044     85.29180    0.0026858
a                1          722          722      1.53781       0.3031
b                1          4.5          4.5      0.00958      0.92819
c                1          392          392      0.83493      0.42823
d                1        544.5        544.5      1.15974       0.3604
ERROR1           3       1408.5        469.5

Cmd> anova("y=d+a*b*c") # all main effects + a.b + a.c +b.c + a.b.c
Model used is y=d+a*b*c
WARNING: summaries are sequential
                DF           SS           MS
CONSTANT         1        40044        40044
d                1        544.5        544.5
a                1          722          722
b                1          4.5          4.5
a.b              1            2            2
c                1          392          392
a.c              1        684.5        684.5
b.c              1          722          722
a.b.c            0            0    undefined
ERROR1           0            0    undefined

Cmd> anova("y=d+(a+b+c)^3") # same model n a different order
Model used is y=d+(a+b+c)^3
WARNING: summaries are sequential
                DF           SS           MS
CONSTANT         1        40044        40044
d                1        544.5        544.5
a                1          722          722
b                1          4.5          4.5
c                1          392          392
a.b              1            2            2
a.c              1        684.5        684.5
b.c              1          722          722
a.b.c            0            0    undefined
ERROR1           0            0    undefined

Cmd> #10.15 Split plot designs

Cmd> repl <- factor(rep(run(2),rep(12,2))) # replicate

Cmd> whole <- factor(rep(rep(run(3),rep(4,3)),2))# whole plot factor

Cmd> sub <- factor(rep(run(4),6)) # subplot factor

Cmd> y <- vector(30,35,37,36,34,41,38,42,29,26,33,36,\
28,32,40,41,31,36,42,40,31,30,32,40)

Cmd> anova("y=repl+whole+E(repl.whole)+sub+whole.sub",fstat:T)
Model used is y=repl+whole+E(repl.whole)+sub+whole.sub
                DF           SS           MS            F      P-value
CONSTANT         1        29400        29400   4126.31579   0.00024226
repl             1          1.5          1.5      0.21053      0.69139
whole            2       138.25       69.125      9.70175     0.093443
ERROR1           2        14.25        7.125      1.20423      0.34401
sub              3       266.33       88.778     15.00469   0.00075634
whole.sub        6       58.417       9.7361      1.64554      0.24071
ERROR2           9        53.25       5.9167

Cmd> secoefs("whole",errorterm:"ERROR1") # use whole plot error
component: coefs
(1)       -0.125            3       -2.875
component: se
(1)      0.77055      0.77055      0.77055

Cmd> secoefs("whole",coefs:F) # by default is uses subplot error
(1)      0.70218      0.70218      0.70218

Cmd> contrast("whole",vector(0,1,-1),errorterm:"ERROR1")
component: estimate
(1)        5.875
component: ss
(1)       138.06
component: se
(1)       1.3346

Cmd> secoefs("sub") # uses subplot error by default which is o.k.
component: coefs
(1)         -4.5      -1.6667            2       4.1667
component: se
(1)      0.85999      0.85999      0.85999      0.85999

Cmd> #10.16  Multivariate analysis of variance

Cmd> y <- matread("macanova.dat","irisdata")
irisdata       150     5 FORMAT
) Data from R. A. Fisher, The use of multiple measurements in taxonomic
) problems, Annals of Eugenics 7 (1936) 376-386
) Col. 1: Variety number (1 = I. Setosa, 2 = I. Versicolor, 3 = I.Virginica)
) Col. 2: X1 = Sepal length
) Col. 3: X2 = Sepal width
) Col. 4: X3 = Petal length
) Col. 5: X4 = Petal width
) Rows 1-50:    Group 1 = Iris Setosa
) Rows 51-100:  Group 2 = Iris Versicolor
) Rows 101-150: Group 3 = Iris Virginica

Cmd> varieties <- factor(y[,1]); y <- y[,-1]

Cmd> y <- matrix(y,labels:structure("@",vector("Sepal Len",\
"Sepal Wid","Petal Len", "Petal Wid"))) # add labels, Sec. 8.4.1

Cmd> manova("y=varieties")
Model used is y=varieties
WARNING: summaries are sequential
                          SS and SP Matrices
                DF
CONSTANT         1
             Sepal Len    Sepal Wid    Petal Len    Petal Wid
Sepal Len       5121.7       2679.8       3293.9       1051.2
Sepal Wid       2679.8       1402.1       1723.4       550.01
Petal Len       3293.9       1723.4       2118.4       676.06
Petal Wid       1051.2       550.01       676.06       215.76
varieties        2
             Sepal Len    Sepal Wid    Petal Len    Petal Wid
Sepal Len       63.212      -19.953       165.25       71.279
Sepal Wid      -19.953       11.345       -57.24      -22.933
Petal Len       165.25       -57.24        437.1       186.77
Petal Wid       71.279      -22.933       186.77       80.413
ERROR1         147
             Sepal Len    Sepal Wid    Petal Len    Petal Wid
Sepal Len       38.956        13.63       24.625        5.645
Sepal Wid        13.63       16.962       8.1208       4.8084
Petal Len       24.625       8.1208       27.223       6.2718
Petal Wid        5.645       4.8084       6.2718       6.1566

Cmd> h <- SS[2,,];e <- SS[3,,]#hypothesis (between) & error (within)

Cmd> f <- (diag(h)/DF[2])/(diag(e)/DF[3])

Cmd> f # f-statistics
(1)       119.26        49.16       1180.2       960.01

Cmd> 1 - cumF(f,DF[2],DF[3]) # P values
(1)            0            0            0            0

Cmd> manova(,sssp:F) # uses most recent model
Model used is y=varieties
WARNING: summaries are sequential
                          SS and SP Matrices
                DF
CONSTANT         1
                    Type 'SS[1,,]' to see SS/SP matrix
varieties        2
                    Type 'SS[2,,]' to see SS/SP matrix
ERROR1         147
                    Type 'SS[3,,]' to see SS/SP matrix

Cmd> SS[2,,] # varieties matrix
                       Sepal Len    Sepal Wid    Petal Len    Petal Wid
varieties Sepal Len       63.212      -19.953       165.25       71.279
          Sepal Wid      -19.953       11.345       -57.24      -22.933
          Petal Len       165.25       -57.24        437.1       186.77
          Petal Wid       71.279      -22.933       186.77       80.413

Cmd> manova(,fstat:T)# to suppress P values, use pvals:F
Model used is y=varieties
WARNING: summaries are sequential
                DF           SS           MS            F      P-value
CONSTANT         1
 Sepal Len               5121.7       5121.7  19326.50528            0
 Sepal Wid               1402.1       1402.1  12151.14260            0
 Petal Len               2118.4       2118.4  11439.11809            0
 Petal Wid               215.76       215.76   5151.66322            0
varieties        2
 Sepal Len               63.212       31.606    119.26450            0
 Sepal Wid               11.345       5.6725     49.16004            0
 Petal Len                437.1       218.55   1180.16118            0
 Petal Wid               80.413       40.207    960.00715            0
ERROR1         147
 Sepal Len               38.956      0.26501
 Sepal Wid               16.962      0.11539
 Petal Len               27.223      0.18519
 Petal Wid               6.1566     0.041882

Cmd> manova("y=varieties",byvar:T,fstat:T)
Model used is y=varieties
WARNING: summaries are sequential
                               Sepal Len
                DF           SS           MS            F      P-value
CONSTANT         1       5121.7       5121.7  19326.50528            0
varieties        2       63.212       31.606    119.26450            0
ERROR1         147       38.956      0.26501
                               Sepal Wid
                DF           SS           MS            F      P-value
CONSTANT         1       1402.1       1402.1  12151.14260            0
varieties        2       11.345       5.6725     49.16004            0
ERROR1         147       16.962      0.11539
                               Petal Len
                DF           SS           MS            F      P-value
CONSTANT         1       2118.4       2118.4  11439.11809            0
varieties        2        437.1       218.55   1180.16118            0
ERROR1         147       27.223      0.18519
                               Petal Wid
                DF           SS           MS            F      P-value
CONSTANT         1       215.76       215.76   5151.66322            0
varieties        2       80.413       40.207    960.00715            0
ERROR1         147       6.1566     0.041882

Cmd> #10.16.1 Multivariate test statistics

Cmd> lambda <- det(e)/det(e+h); lambda# Wilk's lambda
(1)     0.023439

Cmd> fe <- DF[3]; fh <- DF[2]; p <- ncols(y); vector(fe, fh, p)
(1)          147            2            4

Cmd> m1 <- fe - (p - fh + 1)/2

Cmd> w <- -m1 * log(lambda); vector(m1,w)
(1)        145.5       546.12

Cmd> df <- p*fh

Cmd> vector(df,1-cumchi(w,df))#-m1*log(lambda)~chisq(p*fh)
(1)            8            0

Cmd> m2 <- fe - p + 1; m3 <- fe + fh; vector(m2, m3)
(1)          144          149

Cmd> m2*trace(solve(e, h)) # T
(1)       4676.7

Cmd> m3*trace(solve(h+e,h)) # V
(1)       177.59

Cmd> eigs <- releigen(h,e); eigs # relative eigen things
component: values
(1)       32.192      0.28539   3.4751e-15  -1.0354e-15
component: vectors
                   (1)          (2)          (3)          (4)
Sepal Len    -0.068406    0.0019879      0.22042      0.14317
Sepal Wid     -0.12656      0.17853     -0.22026     0.020649
Petal Len      0.18155    -0.076864     -0.24547     0.038476
Petal Wid       0.2318      0.23417      0.31196     -0.21039

Cmd> m1 * sum(log(1 + eigs$values)) # log likelihood ratio
(1)       546.12

Cmd> m2*sum(eigs$values) # T
(1)       4676.7

Cmd> m3*sum(eigs$values/(1+eigs$values)) # V
(1)       177.59

Cmd> #10.16.2 MANOVA canonical variables

Cmd> z <- y %*% eigs$vectors # compute canonical variables

Cmd> chplot(CanVar_1:z[,1],CanVar_2:z[,2],varieties,\
title:"MANOVA Canonical Variable plot for Fisher Iris data")
                   MANOVA Canonical Variable plot for Fisher Iris data
           +-+-------+--------+-------+-------+--------+-------+-------+------+
           |  1                       :                                       |
           |                          :                              3        |
       0.75+                          :                          3 3 3        +
           |                          :                        3 3      3     |
        0.7+ 1  1                     :                         3  33         +
           |1      1                  :                     3  3    3         |
       0.65+  1   11  11              :             2     2    3 3            +
 C         |   1  1                   :               2       3  3  33        |
 a      0.6+    11111                 :            2       3   3    3         +
 n         |    111  11 1             :             2        3     3          |
 V     0.55+      11                  :       2  22 22 2 2 3  333 3           +
 a         |      1111                :           2 222     3  33      3      |
 r      0.5+          11              :          222  2      3                +
 _         |       1 11               :     2   2  22 2   3  3    3          3|
 2     0.45+        1 11              :           2 22             3   3      +
           |                          :     2        2    2        3          |
           |                          :       2 22   2   2                    |
        0.4+              1           :       22      2                       +
           |                          :                      3 3              |
       0.35+                          :          2                            +
           +-+-------+--------+-------+-------+--------+-------+-------+------+
           -0.6    -0.4     -0.2      0      0.2      0.4     0.6     0.8
                                        CanVar_1

Cmd> #10.17 Repeated measures designs

Cmd> accuracy <- vector(45,35,60,50,42,56, 53,41,65,48,45,60,\
60,50,75,61,55,77, 40,30,58,25,30,40, 52,37,54,34,37,39,\
57,47,70,51,43,57, 28,25,40,16,22,31, 37,32,47,23,27,29,\
46,41,50,35,37,46) # accuracy[subjects,period,dial]

Cmd> noise <- factor(rep(rep(run(2),rep(3,2)),9))

Cmd> period <- factor(rep(run(3),rep(18,3)))

Cmd> dial <- factor(rep(rep(run(3),rep(6,3)),3))

Cmd> subjinnoise <- factor(rep(run(3),18))

Cmd> anova("accuracy=noise + E(subjinnoise.noise) +
 period*noise + E(subjinnoise.period.noise) +
 dial*noise + E(subjinnoise.dial.noise) +
 dial*period*noise",fstat:T)
Model used is accuracy=noise + E(subjinnoise.noise) +
 period*noise + E(subjinnoise.period.noise) +
 dial*noise + E(subjinnoise.dial.noise) +
 dial*period*noise
                    DF           SS           MS            F      P-value
CONSTANT             1   1.0587e+05   1.0587e+05    169.99349   0.00019973
noise                1       468.17       468.17      0.75174      0.43484
ERROR1               4       2491.1       622.78     21.21097   0.00025662
period               2       3722.3       1861.2     63.38884   1.2413e-05
noise.period         2          333        166.5      5.67077     0.029268
ERROR2               8       234.89       29.361      2.22526      0.13944
dial                 2       2370.3       1185.2     89.82316   3.3037e-06
noise.dial           2       50.333       25.167      1.90737      0.21022
ERROR3               8       105.56       13.194      1.66084      0.18446
period.dial          4       10.667       2.6667      0.33566      0.84992
noise.period.dial    4       11.333       2.8333      0.35664      0.83567
ERROR4              16       127.11       7.9444

Cmd> ms <- SS/DF # compute mean squares

Cmd> df <- DF;termnames <- TERMNAMES

Cmd> J1 <- termnames=="noise.period"

Cmd> J2 <- termnames=="ERROR2"

Cmd> f <- ms[J1]/ms[J2];f
 noise.period
       5.6708

Cmd> 1 - cumF(f,df[J1],df[J2])
(1)     0.029268

Cmd> accuracy1 <- matrix(accuracy,6);noise1 <- factor(noise[run(6)])

Cmd> manova("accuracy1=noise1") # accuracy1 is 6 by 9
Model used is accuracy1=noise1
WARNING: summaries are sequential
NOTE: SS/SP matrices suppressed because of size; use 'manova(,sssp:T)'
                          SS and SP Matrices
                DF
CONSTANT         1
                    Type 'SS[1,,]' to see SS/SP matrix
noise1           1
                    Type 'SS[2,,]' to see SS/SP matrix
ERROR1           4
                    Type 'SS[3,,]' to see SS/SP matrix

Cmd> list(SS,DF,RESIDUALS) # see shapes of side-effect variables
DF              REAL   3    
RESIDUALS       REAL   6     9    
SS              REAL   3     9     9    

Cmd> s <- matrix(SS[3,,],9)/DF[3] # estimated covariance matrix

Cmd> period1 <- rep(run(3),rep(3,3));tmpx <- 1*(period1==run(3)')

Cmd> hconcat(period1,tmpx) # tmpx are dummy variables for period
(1,1)            1            1            0            0
(2,1)            1            1            0            0
(3,1)            1            1            0            0
(4,1)            2            0            1            0
(5,1)            2            0            1            0
(6,1)            2            0            1            0
(7,1)            3            0            0            1
(8,1)            3            0            0            1
(9,1)            3            0            0            1

Cmd> vtot <- tmpx' %*% s %*% tmpx; e <- (1/3)*rep(1,3)%*%rep(1,3)'

Cmd> print(vtot,e)
vtot:
             (1)          (2)          (3)
(1)       1025.8       812.33       524.83
(2)       812.33       690.83       443.08
(3)       524.83       443.08       327.83
e:
(1,1)      0.33333      0.33333      0.33333
(2,1)      0.33333      0.33333      0.33333
(3,1)      0.33333      0.33333      0.33333

Cmd> m <- vtot - e %*% vtot # m is 3 by 3

Cmd> adjustment <- trace(m)^2/(2*trace(m %*% m)); adjustment
(1)       0.6476

Cmd> tmpx1 <- (1/3) * tmpx %*% tmpx' - (1/9)*rep(1,9)%*%rep(1,9)'

Cmd> m1 <- tmpx1 %*% s # (projection matrix) %*% (covariance matrix)

Cmd> trace(m1)^2/(2*trace(m1%*%m1))
(1)       0.6476

Cmd> dfad1 <- adjustment*df[J1] # adjust DF

Cmd> dfad2 <- adjustment*df[J2]

Cmd> 1 - cumF(f,dfad1,dfad2) # adjusted P value
(1)      0.05694

Cmd> dial1 <- rep(run(3),3);tmpy <- 1*(dial1==run(3)')

Cmd> hconcat(dial1,tmpy)
(1,1)            1            1            0            0
(2,1)            2            0            1            0
(3,1)            3            0            0            1
(4,1)            1            1            0            0
(5,1)            2            0            1            0
(6,1)            3            0            0            1
(7,1)            1            1            0            0
(8,1)            2            0            1            0
(9,1)            3            0            0            1

Cmd> vtot <- tmpy %c% s %*% tmpy

Cmd> m <- vtot - e %*% vtot # same e == (1/3)*rep(1,3) %*% rep(1,3)'

Cmd> adjustment <- trace(m)^2/(2*trace(m %*% m)); adjustment
(1)      0.91707

Cmd> tmpy1 <- (1/3) * tmpy %*% tmpy' - (1/9)*rep(1,9)%*%rep(1,9)'

Cmd> m1 <- tmpy1 %*% s

Cmd> trace(m1)^2/(2*trace(m1%*%m1))
(1)      0.91707

Cmd> I9 <- dmat(9,1) # 9 by 9 diagonal matrix with 1's on diagonal

Cmd> tmpxy <- (I9 - tmpx%*%tmpx'/3) %*% (I9 - tmpy%*%tmpy'/3)

Cmd> # or tmpxy <- I9-tmpx1-tmpy1-(1/9)*rep(1,9) %*% rep(1,9)'

Cmd> m1 <- tmpxy %*% s

Cmd> adjustment <- trace(m1)^2/(4*trace(m1%*%m1)); adjustment
(1)      0.51342

Cmd> #10.18  Logistic regression

Cmd> y <- vector(0,9,21,47,60,63) # counts

Cmd> n <- 70 # or n <- rep(70,6)

Cmd> current <- run(0,5)

Cmd> logistic("y=current",n,inc:T,pvals:T)
Model used is y=current
WARNING: summaries are sequential
                DF     Deviance         MDev      P-value
CONSTANT         1      0.95274      0.95274      0.32902
current          1       241.13       241.13            0
ERROR1           4       9.3526       2.3382     0.052865

Cmd> regcoefs() # see Sec. 3.13.1
NOTE: standard errors assume scale parameter is 1
                 Coef       StdErr            t
CONSTANT       -3.301       0.3238      -10.195
current        1.2459      0.11193       11.132

Cmd> sqrt(SS[2])*coefs(2)/abs(coefs(2))
      current
       15.528

Cmd> logistic("y=current-1+1",n,inc:T,silent:T) # make CONSTANT last

Cmd> sqrt(SS[2])*coefs(2)/abs(coefs(2))
     CONSTANT
      -13.736

Cmd> logistic("y=P2(current)",n,inc:T,pvals:T) # quadratic model
Model used is y=P2(current)
WARNING: summaries are sequential
                DF     Deviance         MDev      P-value
CONSTANT         1      0.95274      0.95274      0.32902
{current}        1       241.13       241.13            0
{(current)^2}    1       5.5426       5.5426     0.018559
ERROR1           3         3.81         1.27      0.28273

Cmd> regcoefs()
NOTE: standard errors assume scale parameter is 1
                     Coef       StdErr            t
CONSTANT          -4.4043      0.65897      -6.6836
{current}          2.2404      0.47723       4.6946
{(current)^2     -0.18078     0.079678      -2.2689

Cmd> #10.19 Poisson regression

Cmd> counts <- vector(141,68,17,54,44,11,40,51,19)

Cmd> parentuse <- factor(rep(run(3),rep(3,3))) # 1,1,1,2,2,2,3,3,3

Cmd> studentuse <- factor(rep(run(3),3)) # 1,2,3,1,2,3,1,2,3

Cmd> poisson("counts=studentuse + parentuse",inc:T,pvals:T)
Model used is counts=studentuse + parentuse
WARNING: summaries are sequential
                DF     Deviance         MDev      P-value
CONSTANT         1       2599.8       2599.8            0
studentuse       2       138.96        69.48            0
parentuse        2        57.38        28.69   3.4683e-13
ERROR1           4       22.254       5.5634   0.00017841

Cmd> secoefs("parentuse") # or secoefs(3)
NOTE: standard errors assume scale parameter is 1
component: coefs
(1)      0.48308     -0.24611     -0.23697
component: se
(1)     0.063214     0.074694     0.074508

Cmd> HII # "leverages", see Sec. 3.6
(1)      0.76776      0.68813      0.55984      0.64368      0.52151
(6)      0.32469      0.64474      0.52294       0.3267

Cmd> matrix(WTDRESIDUALS,3) # weighted residuals in the log scale
(1,1)       1.9819     -0.46946      -2.3735
(2,1)      -1.6247      0.64478       1.6869
(3,1)      -1.4061     -0.15101       2.1658

Cmd> predtable() # table of fitted values
(1,1)       119.35       57.562        58.09
(2,1)       82.782       39.926       40.292
(3,1)        23.87       11.512       11.618

Cmd> onedof <- outer(coefs("parentuse"),coefs("studentuse"))

Cmd> #See Sec 3.16 for another use of outer()

Cmd> print(onedof) # display as table
onedof:
(1,1)      0.31807      0.14134     -0.45942
(2,1)     -0.16204    -0.072008      0.23405
(3,1)     -0.15603    -0.069336      0.22537

Cmd> onedof <- vector(onedof) # make a vector

Cmd> poisson("counts=studentuse+parentuse+onedof",inc:T,pvals:T)
Model used is counts=studentuse+parentuse+onedof
WARNING: summaries are sequential
                DF     Deviance         MDev      P-value
CONSTANT         1       2599.8       2599.8            0
studentuse       2       138.96        69.48            0
parentuse        2        57.38        28.69   3.4683e-13
onedof           1       19.238       19.238   1.1539e-05
ERROR1           3       3.0155       1.0052      0.38923

Cmd> ipf("counts=studentuse+parentuse",inc:T)
Model used is counts=studentuse+parentuse
WARNING: summaries are sequential
                DF     Deviance         MDev
CONSTANT         1       2599.8       2599.8
studentuse       2       138.96        69.48
parentuse        2        57.38        28.69
ERROR1           4       22.254       5.5634

Cmd> secoefs("parentuse")
ERROR: standard errors not available after ipf()

Cmd> coefs("parentuse")
(1)      0.48308     -0.24611     -0.23697

Cmd> HII  # this is not correct; it sums to the right total, though
(1)      0.55556      0.55556      0.55556      0.55556      0.55556
(6)      0.55556      0.55556      0.55556      0.55556

Cmd> #10.20 Robust regression

Cmd> adsorption[8] <- 300

Cmd> regress("adsorption=iron+aluminum",pvals:T)
Model used is adsorption=iron+aluminum
                 Coef       StdErr            t      P-Value
CONSTANT       29.819       62.529      0.47689      0.64369
iron         -0.26721      0.53277     -0.50155      0.62684
aluminum       1.3827       1.2795       1.0806      0.30525

N: 13, MSE: 6175.4, DF: 10, R^2: 0.13141
Regression F(2,10): 0.75648, P-value: 0.49439, Durbin-Watson: 2.0567
To see the ANOVA table type 'anova()'

Cmd> robust("adsorption=iron+aluminum",fstats:T)
Model used is adsorption=iron+aluminum
WARNING: summaries are sequential
                DF          SS*          MS*           F*     P-value*
CONSTANT         1        12644        12644    688.97181   1.4839e-10
iron             1       3172.8       3172.8    172.88388   1.2313e-07
aluminum         1       626.13       626.13     34.11716   0.00016359
ERROR1          10       183.52       18.352
* ANOVA is approximate and should be interpreted with caution

Robust estimate of sigma: 4.5468

Cmd> f <- (SS[3]/DF[3])/(SS[4]/DF[4]) # F-statistic

Cmd> vector(f,1-cumF(f,DF[3],DF[4]))# F-statistic & P value
(1)       34.117   0.00016359

Cmd> temp <- secoefs(byterm:F) # get coefs and their stderrs

Cmd> compnames(temp) # coefs(byterm:F) has 2 components
(1) "coefs"
(2) "se"

Cmd> coef <- vector(temp[1]);se <- vector(temp[2]);tstats <- coef/se

Cmd> hconcat(coef,se,tstats) # table of coeffs, std errors, t-stats
(1,1)      -6.1327       3.4087      -1.7991
(2,1)     0.097179     0.029044       3.3459
(3,1)      0.40743     0.069753        5.841

Cmd> WTDRESIDUALS/RESIDUALS
 (1)            1            1            1            1      0.51868
 (6)            1            1     0.012877            1      0.32462
(11)            1            1      0.71154

Cmd> WTDRESIDUALS
 (1)      -1.0916      -1.4297     -0.43231       2.7117      -3.4105
 (6)     -0.16143      -2.7356       3.4105       2.6945      -3.4105
(11)      0.52653    -0.081364       3.4105

Cmd> robust("adsorption=iron+aluminum",fstats:T,marginal:T)
Model used is adsorption=iron+aluminum
WARNING: SS are Type III sums of squares
                DF          SS*          MS*           F*     P-value*
CONSTANT         1       59.402       59.402      3.23674       0.1022
iron             1       205.46       205.46     11.19513    0.0074162
aluminum         1       626.13       626.13     34.11716   0.00016359
ERROR1          10       183.52       18.352
* ANOVA is approximate and should be interpreted with caution

Robust estimate of sigma: 4.5468

Cmd> m <- sum(round(WTDRESIDUALS/RESIDUALS,12) == 1 )

Cmd> n <- length(adsorption); p <- sum(DF[-4]); vector(n,m,p)
(1)           13            9            3

Cmd> mu <- m/n; K <- (1 + (p/n)*(1-mu)/mu)/mu; K
(1)       1.5926

Cmd> fit <- adsorption - RESIDUALS; y1 <- fit + K * WTDRESIDUALS

Cmd> anova("y1=iron+aluminum",marginal:T) # same result as robust()
Model used is y1=iron+aluminum
WARNING: SS are Type III sums of squares
                DF           SS           MS
CONSTANT         1       59.402       59.402
iron             1       205.46       205.46
aluminum         1       626.13       626.13
ERROR1          10       183.52       18.352

Cmd> adsorption[8] <- ? # set case 8 to MISSING

Cmd> regress("adsorption=iron+aluminum",pvals:T)
Model used is adsorption=iron+aluminum
WARNING: cases with missing values deleted
                 Coef       StdErr            t      P-Value
CONSTANT      -6.8757       3.6342      -1.8919     0.091053
iron          0.10788     0.031191       3.4586    0.0071765
aluminum      0.36221      0.07541       4.8032   0.00096936

N: 13, Omitted: 1, MSE: 20.164, DF: 9, R^2: 0.95124
Regression F(2,9): 87.787, P-value: 1.2483e-06, Durbin-Watson: 2.6746
To see the ANOVA table type 'anova()'

Cmd> anova(,marginal:T,fstat:T)
Model used is adsorption=iron+aluminum
WARNING: cases with missing values deleted
WARNING: SS are Type III sums of squares
                DF           SS           MS            F      P-value
CONSTANT         1       72.174       72.174      3.57944     0.091053
iron             1        241.2        241.2     11.96200    0.0071765
aluminum         1       465.19       465.19     23.07090   0.00096936
ERROR1           9       181.47       20.164

Cmd> asciisave("chkpnt10.sav") # make sure all variables created are savable
Workspace asciisaved on file chkpnt10.sav

Cmd> restore("chkpnt10.sav")
Restoring workspace from file chkpnt10.sav
Workspace saved Thu Mar 25 13:41:10 1999

Cmd> print(paste("Total time =",gettime(quiet:T,keep:T) - startTime))
Total time = 6.74

Cmd> 
Normal termination by end of file on input
