                             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 4 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> #4.  Generalized linear models and robustness (advanced topic.)

Cmd> #4.1 Alternatives to Least Squares

Cmd> #4.2 Generalized Linear Model Commands

Cmd> y <- vector(34,34,37,23,23,38,12,27,29,17,17,29)

Cmd> n <- vector(87,88,82,82,84,80,87,82,80,84,85,84)

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

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

Cmd> #4.2.2 Side effect variables

Cmd> #4.2.3 Keywords

Cmd> #4.2.4 logistic() and probit()

Cmd> logistic("y=a+b",n,increment:T) # split Overall model deviance
Model used is y=a+b
WARNING: summaries are sequential
                DF     Deviance         MDev
CONSTANT         1       135.64       135.64
a                2       19.216       9.6078
b                3       18.694       6.2314
ERROR1           6       9.8484       1.6414

Cmd> SS
     CONSTANT            a            b       ERROR1
       135.64       19.216       18.694       9.8484

Cmd> SS[1] + SS[2] + SS[3] # total deviance for model
     CONSTANT
       173.55

Cmd> logistic("y=a+b",n)
Model used is y=a+b
WARNING: summaries are sequential
                DF     Deviance         MDev
Overall model    6       173.55       28.925
ERROR1           6       9.8484       1.6414

Cmd> # Overall model deviance is SS[1]+SS[2]+SS[3] in preceding

Cmd> SS #  All but last two elements are 0
                          Overall model       ERROR1
            0            0       173.55       9.8484

Cmd> DF # Degrees of freedom for 1st 3 terms are not pooled
(1)            0            0            6            6

Cmd> 1 - cumchi(SS[4],DF[4]) # P value for goodness of fit
(1)      0.13118

Cmd> probit("y=a+b",n)
Model used is y=a+b
WARNING: summaries are sequential
                DF     Deviance         MDev
Overall model    6       173.77       28.962
ERROR1           6       9.6263       1.6044

Cmd> groups <- factor(run(3)); counts <- vector(0, 7, 10)

Cmd> logistic("counts=groups",10,inc:T)
Model used is counts=groups
WARNING: problimit = 1e-08 was hit by logistic() at least once
WARNING: summaries are sequential
                DF     Deviance         MDev
CONSTANT         1      0.53492      0.53492
groups           2       28.837       14.418
ERROR1           0            0    undefined

Cmd> secoefs() # see Sec. 3.13
NOTE: standard errors assume scale parameter is 1
component: CONSTANT
  component: coefs
(1)      0.28243
  component: se
(1)       1490.7
component: groups
  component: coefs
(1)      -19.703      0.56487       19.138
  component: se
(1)         2357       1490.7         2357

Cmd> glmtable() # compute fitted probabilities; see Sec. 3.18
NOTE: standard errors assume scale parameter is 1
component: estimate
(1)   3.6788e-09          0.7            1
component: SEest
(1)   1.1633e-05      0.14491   1.1633e-05

Cmd> #4.2.5 poisson() and ipf()

Cmd> poisson("y=a+b",increment:T)
Model used is y=a+b
WARNING: summaries are sequential
                DF     Deviance         MDev
CONSTANT         1       1485.4       1485.4
a                2       10.621       5.3107
b                3         13.1       4.3666
ERROR1           6       7.4863       1.2477

Cmd> ipf("y=a+b")
Model used is y=a+b
WARNING: summaries are sequential
                DF     Deviance         MDev
Overall model    6       1509.1       251.52
ERROR1           6       7.4863       1.2477

Cmd> #4.2.6 glmfit()

Cmd> glmfit("y=a+b") # link:"identity", distrib:"normal" assumed
Model used is y=a+b
WARNING: summaries are sequential
                DF           SS           MS
CONSTANT         1       8533.3       8533.3
a                2       288.17       144.08
b                3          358       119.33
ERROR1           6        136.5        22.75

Response distribution is normal, link function is identity

Cmd> # This table is identical with table from anova("y=a+b")

Cmd> glmfit("y=a+b",distrib:"poisson",link:"log") #same as poisson()
Model used is y=a+b
WARNING: summaries are sequential
                DF     Deviance         MDev
Overall model    6       1509.1       251.52
ERROR1           6       7.4863       1.2477

Response distribution is poisson, link function is log

Cmd> glmfit("y=a+b",n:n,distrib:"binomial",link:"probit") #probit()
Model used is y=a+b
WARNING: summaries are sequential
                DF     Deviance         MDev
Overall model    6       173.77       28.962
ERROR1           6       9.6263       1.6044

Response distribution is binomial, link function is probit

Cmd> glmfit("y=a+b",distrib:"normal",link:"log")
ERROR: glmfit(Model,link:"log",distrib:"normal") is not implemented

Cmd> #4.3 Robust regression using robust()

Cmd> y1 <- y; y1[12] <- 130 # make modified y with outlier

Cmd> anova("y1=a+b") # fit by least squares
Model used is y1=a+b
                DF           SS           MS
CONSTANT         1        14770        14770
a                2       3318.2       1659.1
b                3       1763.6       587.86
ERROR1           6       5523.2       920.53

Cmd> robust("y1=a+b") # fit robustly, truncation c = .75
Model used is y1=a+b
WARNING: summaries are sequential
                DF          SS*          MS*
CONSTANT         1       9332.6       9332.6
a                2       516.85       258.42
b                3       245.62       81.873
ERROR1           6       371.15       61.858
* ANOVA is approximate and should be interpreted with caution

Robust estimate of sigma: 11.657

Cmd> WTDRESIDUALS/RESIDUALS # only 1 residual was truncated
 (1)            1            1            1            1            1
 (6)            1            1            1            1            1
(11)            1     0.091936

Cmd> modelinfo(weights:T,scale:T,sigmahat:T)
component: scale
(1)        7.865
component: weights
 (1)            1            1            1            1            1
 (6)            1            1            1            1            1
(11)            1    0.0084523
component: sigmahat
(1)       11.657

Cmd> sqrt(SS[4]/DF[4]) # value of scale
       ERROR1
        7.865

Cmd> robust("y1=a+b",trunc:.60) # more severe truncation
Model used is y1=a+b
WARNING: summaries are sequential
                DF          SS*          MS*
CONSTANT         1         9233         9233
a                2       485.06       242.53
b                3        253.6       84.532
ERROR1           6       332.51       55.418
* ANOVA is approximate and should be interpreted with caution

Robust estimate of sigma: 13.078

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

Cmd> restore("chkpnt4.sav")
Restoring workspace from file chkpnt4.sav
Workspace saved Thu Mar 25 13:40:23 1999

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

Cmd> 
Normal termination by end of file on input
