[R] Ordinal categorical data with GLM
John Fox
jfox at mcmaster.ca
Tue Aug 13 02:23:25 CEST 2002
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list