[R] Ordinal categorical data with GLM

Duncan Mackay Duncan.Mackay at flinders.edu.au
Tue Aug 13 05:31:08 CEST 2002


Hello,
you may also be interested in the following link (from Dr. Agresti's home
page) that provides a link to a document containing S-plus scripts for most
of the examples in the 1990 edition of "Categorical Data Analysis" by
Agresti.

http://math.cl.uh.edu/~thompsonla/5537/Splusdiscrete.PDF


Cheers,
Duncan
*******************************************
Dr. Duncan Mackay
Biology
Flinders University
GPO Box 2100
Adelaide
S.A.    5001
AUSTRALIA

Ph (08) 8201 2627    FAX (08) 8201 3015

http://www.bio.flinders.edu.au/dam/damres.htm


-----Original Message-----
From: owner-r-help at stat.math.ethz.ch
[mailto:owner-r-help at stat.math.ethz.ch]On Behalf Of John Fox
Sent: Tuesday, 13 August 2002 9:53 AM
To: Andrew Criswell
Cc: r-help at stat.math.ethz.ch
Subject: Re: [R] Ordinal categorical data with GLM


Dear Andrew,

I don't believe that this question has been answered yet.

It's simple to fit the linear-by-linear association (uniform association)
model, since it is essentially a Poisson regression of the cell counts on
factors for the rows and columns plus the products of the row and column
indices.

Adapting your code:

     > df <- data.frame(Freq, X = factor(X), Y = factor(Y), x=X, y=Y)
     > summary(glm(Freq ~ X + Y + I(x*y), data = df, family = poisson))

     Call:
     glm(formula = Freq ~ X + Y + I(x * y), family = poisson, data = df)

     Deviance Residuals:
         Min         1Q     Median         3Q        Max
     -1.034350  -0.175387   0.005087   0.270322   0.579845

     Coefficients:
                 Estimate Std. Error z value Pr(>|z|)
     (Intercept)  3.47292    0.10251  33.878  < 2e-16
     X1           0.29493    0.13241   2.227  0.02592
     X2           0.39634    0.06927   5.722 1.06e-08
     X3          -0.05812    0.06840  -0.850  0.39549
     Y1          -0.80530    0.11878  -6.780 1.21e-11
     Y2          -0.38583    0.08538  -4.519 6.21e-06
     Y3           0.54864    0.06200   8.849  < 2e-16
     I(x * y)     0.11194    0.03641   3.075  0.00211

     (Dispersion parameter for poisson family taken to be 1)

         Null deviance: 445.7627  on 15  degrees of freedom
     Residual deviance:   2.3859  on  8  degrees of freedom
     AIC: 107.42

     Number of Fisher Scoring iterations: 3

     > round(data.3 - t(matrix(residuals(mod, type='response'),4,4)),1)
                 Satisfaction
     Income          Very dissatisfied Little dissatisfied Moderately
satisfied Very satisfied
     <
6,000                    19.3                29.4                 74.9
      82.3
     6,000-15,000               21.4                36.4
103.7          127.4
     15,000-25,000              13.6                25.9
82.4          113.2
     >
25,000                    7.6                16.3                 58.0
      89.1


which is the result given in Agesti.

I hope that this helps,
  John

At 11:31 AM 8/11/2002 +0700, Andrew Criswell wrote:
>Hello All:
>
>I am looking for you help.
>
>I am trying to replicate the results of an example found in Alan Agresti's
>"Categorical Data Analysis" on pages 267-269. The example is one of a 2 x 2
>cross-classification table of ordinal counts: job satisfaction and income.
>
>I am able to get Agresti's results for the independence model (G^2 = 12.03
>with df = 9) assuming as he does that the data is nominal, but I'm unable
to
>derive his model of uniform association (linear-by-linear association, p.
>263-269) for which he gets a value of G^2 = 2.39 with df = 8.
>
>The observed data is represented by table 8.2 on page 268 and as follows:
>
>Freq <- c(20, 24,  80,  82, 22, 38, 104, 125, 13, 28,  81, 113, 7, 18,  54,
>92)
>
>data.3 <- t(matrix(Freq, nrow = 4))
>
>list.3 <- list(Income = c("< 6,000", "6,000-15,000", "15,000-25,000", ">
>25,000"),
>                   Satisfaction = c("Very dissatisfied", "Little
>dissatisfied", "Moderately satisfied", "Very satisfied"))
>
>dimnames(data.3) <- list.3
>
>ftable(data.3)
>
>I am able to obtain Agresti's results for the independence model which
>assumes the data is nominal, not ordinal, using either glm() or loglm().
>
>library(MASS)
>options(contrasts=c("contr.sum", "contr.poly"))
>
>X <- as.integer(gl(4, 4, 16)) - 1
>Y <- as.integer(gl(4, 1, 16)) - 1
>
>data.2 <- data.frame(Freq, X = factor(X), Y = factor(Y))
>
>summary(fm3 <- glm(Freq ~ X + Y, data = data.2, family = poisson()))
>dummy.coef(fm3)
>
>fm4 <- loglm(Freq ~ X + Y, data = data.2, param = T, fit = T)
>fm4;  fm4$param
>
>My question is this: can glm() or some other function be used in the manner
>Agresti employed for ordinal count data?
>
>Thank you,
>ANDREW
>
>Andrew Criswell
>Professor of Finance
>Graduate School
>Bangkok University
>
>
>-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
.-.-
>r-help 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-help-request at stat.math.ethz.ch
>_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
._._

-----------------------------------------------------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario, Canada L8S 4M4
email: jfox at mcmaster.ca
phone: 905-525-9140x23604
web: www.socsci.mcmaster.ca/jfox
-----------------------------------------------------

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
-.-
r-help 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-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
_._


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help 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-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._



More information about the R-help mailing list