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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._