From james.lindsey@luc.ac.be  Mon Dec 18 08:27:58 2000
From: james.lindsey@luc.ac.be (james.lindsey@luc.ac.be)
Date: Mon, 18 Dec 2000 09:27:58 +0100 (MET)
Subject: [Rd] problems with glm (PR#771)
Message-ID: <200012180827.JAA16557@pubhealth.ku.dk>

R1.2.0 with Linux RH5.2

I do not believe that the problems below are new to 1.2 but I only
cover this sort of thing once a year in my course and some of that
happened to be last Friday so too late to report for 1.2. I see that
one or two things that I was going to report have been corrected.
I like the fact that interactions now show : instead of .

Here is some output with comments inserted.


R : Copyright 2000, The R Development Core Team
Version 1.2.0  (2000-12-15)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type `license()' or `licence()' for distribution details.

R is a collaborative project with many contributors.
Type `contributors()' for more information.

Type `demo()' for some demos, `help()' for on-line help, or
`help.start()' for a HTML browser interface to help.
Type `q()' to quit R.

Notice contrasts for below:

> options(contrasts=c("contr.treatment","contr.poly"))
> data <- read.table("tmp.dat",skip=6,head=T)
> data
   Freq Reg66 Reg71
1   118    CC    CC
2    14    LY    CC
3     8    WM    CC
4    12    GL    CC
5    12    CC    LY
6  2127    LY    LY
7    69    WM    LY
8   110    GL    LY
9     7    CC    WM
10   86    LY    WM
11 2548    WM    WM
12   88    GL    WM
13   23    CC    GL
14  130    LY    GL
15  107    WM    GL
16 7712    GL    GL
> attach(data)
> wt <- codes(Reg66)!=codes(Reg71)
> print(z <- glm(Freq~Reg66+Reg71,family=poisson))

Call:  glm(formula = Freq ~ Reg66 + Reg71, family = poisson) 

Coefficients:
(Intercept)      Reg66GL      Reg66LY      Reg66WM      Reg71GL      Reg71LY  
     0.6134       3.9022       2.6900       2.8376       3.9597       2.7245  
    Reg71WM  
     2.8877  

Degrees of Freedom: 15 Total (i.e. Null);  9 Residual
Null Deviance:	    40740 
Residual Deviance: 19880 	AIC: 20000 
> print(z2 <- glm(Freq~Reg66+Reg71,family=poisson,weight=wt))

Call:  glm(formula = Freq ~ Reg66 + Reg71, family = poisson, weights = wt) 

Coefficients:
(Intercept)      Reg66GL      Reg66LY      Reg66WM      Reg71GL      Reg71LY  
     0.4615       2.1242       2.0096       1.7236       2.4555       2.0847  
    Reg71WM  
     1.9140  

Degrees of Freedom: 11 Total (i.e. Null);  5 Residual
Null Deviance:	    486.2 
Residual Deviance: 4.367 	AIC: 82.7 
> z2$fit
          1           2           3           4           5           6 
 118.100000   11.836037    8.891541   13.272421   12.758205 2127.100000 
          7           8           9          10          11          12 
  71.505458  106.736339   10.756661   80.252100 2548.100000   89.991240 
         13          14          15          16 
  18.485138  137.911862  103.603001 7712.100000 

First problem: with diagonal values weighted out in a transition
matrix (mover-stayer model), glm wrongly estimates the diagonals to be
the observed plus 0.1. This was correct in 90.1 a year ago but already
wrong in 1.0.1. I don't have any versions in between installed.

> # correct values from GLIM4
> #   unit   observed    fitted   residual
> #     (1)       118     1.586      0.000
> #      2         14    11.836      0.629
> #      3          8     8.892     -0.299
> #      4         12    13.272     -0.349
> #      5         12    12.758     -0.212
> #     (6)      2127    95.185      0.000
> #      7         69    71.505     -0.296
> #      8        110   106.736      0.316
> #      9          7    10.757     -1.145
> #     10         86    80.252      0.642
> #    (11)      2548    60.287      0.000
> #     12         88    89.991     -0.210
> #     13         23    18.485      1.050
> #     14        130   137.912     -0.674
> #     15        107   103.603      0.334
> #    (16)      7712   154.648      0.000
> # and from R90.1 (wrong in R1.0.1)
> #> z2$fit
> #         1          2          3          4          5          6          7 
> #  1.586454  11.836037   8.891541  13.272421  12.758205  95.184992  71.505458 
> #         8          9         10         11         12         13         14 
> #106.736339  10.756661  80.252100  60.287479  89.991240  18.485138 137.911862 
> #        15         16 
> #103.603001 154.648406 
> 

Second problem: the man page for factors states that the levels are
sorted, but gl does not do this. It keeps them in the order of the
labels. This is inconsistent with everything else, confusing for
students, and liable to induce errors in data analysis. In the example
below, it is not too important because the baseline category is the
same in the two cases but that will not generally be the case.

As an addition comment, I would strongly recommend that, if labels can
be coerced to numeric, they be ordered, not sorted. This is because
the order is very confusing if there are more than 9: for example, the
sorted 10 comes near the beginning and not at the end.

> reg66 <- gl(4,1,16,labels=c("CC","LY","WM","GL"))
> z

Call:  glm(formula = Freq ~ Reg66 + Reg71, family = poisson) 

Coefficients:
(Intercept)      Reg66GL      Reg66LY      Reg66WM      Reg71GL      Reg71LY  
     0.6134       3.9022       2.6900       2.8376       3.9597       2.7245  
    Reg71WM  
     2.8877  

Degrees of Freedom: 15 Total (i.e. Null);  9 Residual
Null Deviance:	    40740 
Residual Deviance: 19880 	AIC: 20000 
> print(za <- glm(Freq~reg66+Reg71,family=poisson))

Call:  glm(formula = Freq ~ reg66 + Reg71, family = poisson) 

Coefficients:
(Intercept)      reg66LY      reg66WM      reg66GL      Reg71GL      Reg71LY  
     0.6134       2.6900       2.8376       3.9022       3.9597       2.7245  
    Reg71WM  
     2.8877  

Degrees of Freedom: 15 Total (i.e. Null);  9 Residual
Null Deviance:	    40740 
Residual Deviance: 19880 	AIC: 20000 

Third problem: if mean value contrasts are used, the level names are
not attached to the variable name, but simply the number of the
category. This rightly drives students mad and is simply
uninterpretable if the levels have been sorted. It takes a great deal
of extreme care to try to find out to what the numbers refer.

> options(contrasts=c("contr.sum","contr.poly"))
> 
> z

Call:  glm(formula = Freq ~ Reg66 + Reg71, family = poisson) 

Coefficients:
(Intercept)      Reg66GL      Reg66LY      Reg66WM      Reg71GL      Reg71LY  
     0.6134       3.9022       2.6900       2.8376       3.9597       2.7245  
    Reg71WM  
     2.8877  

Degrees of Freedom: 15 Total (i.e. Null);  9 Residual
Null Deviance:	    40740 
Residual Deviance: 19880 	AIC: 20000 
> print(zb <- glm(Freq~Reg66+Reg71,family=poisson))

Call:  glm(formula = Freq ~ Reg66 + Reg71, family = poisson) 

Coefficients:
(Intercept)       Reg661       Reg662       Reg663       Reg711       Reg712  
     5.3638      -2.3575       1.5448       0.3325      -2.3930       1.5667  
     Reg713  
     0.3315  

Degrees of Freedom: 15 Total (i.e. Null);  9 Residual
Null Deviance:	    40740 
Residual Deviance: 19880 	AIC: 20000 
> 

Additional comment: I recommend that the aic functions for poisson and
binomial have the explicit calculation of the log density replaced by
the corresponding d function with the log option. This will avoid the
production of NAs in certain cases of zeroes in tables.
  Jim

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._


[Rd] problems with glm (PR#771)
james.lindsey@luc.ac.be