LOGIT MODEL ANALYSIS OF PNEUMOCONIOSIS DATA
-------------------------------------------

> pneumo.df <- read.table("data/pneumo.dat",header=T)
> pneumo.df
  time normal mild severe
1  5.8     98    0      0
2 15.0     51    2      1
3 21.5     34    6      3
4 27.5     35    5      8
5 33.5     32   10      9
6 39.5     23    7      8
7 46.0     12    6     10
8 51.5      2    2      5
> attach(pneumo.df)
> 
> # Proportional odds model fit
> logexpos <- log(rep(time,3))
> freq <- c(normal,mild,severe)
> response <- as.factor(rep(1:3,each=8))
> pnewmo.df <- data.frame(response,logexpos,freq)
> pnewmo.df
   response logexpos freq
1         1 1.757858   98
2         1 2.708050   51
3         1 3.068053   34
4         1 3.314186   35
5         1 3.511545   32
6         1 3.676301   23
7         1 3.828641   12
8         1 3.941582    2
9         2 1.757858    0
10        2 2.708050    2
11        2 3.068053    6
12        2 3.314186    5
13        2 3.511545   10
14        2 3.676301    7
15        2 3.828641    6
16        2 3.941582    2
17        3 1.757858    0
18        3 2.708050    1
19        3 3.068053    3
20        3 3.314186    8
21        3 3.511545    9
22        3 3.676301    8
23        3 3.828641   10
24        3 3.941582    5

> library(MASS)
> pnewmo.po <- polr(response~logexpos,weights=freq,data=pnewmo.df)
> summary(pnewmo.po)

Coefficients:
     Value Std. Error    t value 
  2.691131   0.392416   6.857852 

Intercepts:
    Value   Std. Error t value
1|2  9.9758  1.3627     7.3205
2|3 10.8919  1.3842     7.8687

Residual Deviance: 404.4111 
AIC: 410.4111 
> 
> # Continuation ratio model fit
> m <- c(normal+mild+severe,mild+severe)
> y <- c(mild+severe,severe)
> x <- as.factor(rep(1:2,each=8))
> lt <- log(rep(time,2))
> data.frame(m,y,x,lt)
    m  y x       lt
1  98  0 1 1.757858
2  54  3 1 2.708050
3  43  9 1 3.068053
4  48 13 1 3.314186
5  51 19 1 3.511545
6  38 15 1 3.676301
7  28 16 1 3.828641
8   9  7 1 3.941582
9   0  0 2 1.757858
10  3  1 2 2.708050
11  9  3 2 3.068053
12 13  8 2 3.314186
13 19  9 2 3.511545
14 15  8 2 3.676301
15 16 10 2 3.828641
16  7  5 2 3.941582
> pneumo.cr <- glm(cbind(y,m-y)~x+lt-1,binomial)
> summary(pneumo.cr)

Coefficients:
   Estimate Std. Error z value Pr(>|z|)    
x1  -8.9455     1.1528  -7.759 8.53e-15 ***
x2  -8.2949     1.2055  -6.881 5.95e-12 ***
lt   2.3902     0.3342   7.152 8.52e-13 ***
---
(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 226.5831  on 15  degrees of freedom
Residual deviance:   8.7503  on 12  degrees of freedom
AIC: 59.449