See comments in line.
On Wed, 14 Apr 2010, Sachi Ito wrote:
Hi all,
I've been running loglinear models for three-way tables: one of the
variables having three levels, and the other two having two levels each.
An example looks like below:
yes.no <- c("Yes","No")
switch <- c("On","Off")
att <- c("BB","AA","CC")
L <- gl(2,1,12,yes.no)
T <- gl(2,2,12,switch)
A <- gl(3,4,12,att)
n <- c(1136,4998,25,339,305,2752,31,692,251,1677,17,557)
d.table <- data.frame(A,T,L,n)
d.table
A T L n
1 BB On Yes 1136
2 BB On No 4998
3 BB Off Yes 25
4 BB Off No 339
5 AA On Yes 305
6 AA On No 2752
7 AA Off Yes 31
8 AA Off No 692
9 CC On Yes 251
10 CC On No 1677
11 CC Off Yes 17
12 CC Off No 557
First, I run the independence model and found a poor fit:
library(MASS)
loglm(n~A+T+L)
Call:
loglm(formula = n ~ A + T + L)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 1001.431 7 0
Pearson 1006.287 7 0
Thus, I went on and run the two-way association model and found a good fit:
loglm(n~A:T+A:L+T:L)
Call:
loglm(formula = n ~ A:T + A:L + T:L)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 4.827261 2 0.08948981
Pearson 4.680124 2 0.09632168
I compared the independence model (Model1), two-way association model (Model
2), and three-way interaction model (Saturated), and found that the two-way
association model was the most parsimonious one:
ind <- loglm(n~A+T+L)
twoway <- loglm(n~A:T+A:L+T:L)
anova(ind,twoway)
LR tests for hierarchical log-linear models
Model 1:
n ~ T + A + L
Model 2:
n ~ A:L + A:T + T:L
Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
Model 1 1001.430955 7
Model 2 4.827261 2 996.603694 5 0.00000
Saturated 0.000000 0 4.827261 2 0.08949
By running a Chi-square test, I found that all of the three two-way
associations are significant.
drop1(twoway,test="Chisq")
Single term deletions
Model:
n ~ A:T + A:L + T:L
Df AIC LRT Pr(Chi)
<none> 24.83
A:T 2 645.91 625.08 < 2.2e-16 ***
A:L 2 152.93 132.10 < 2.2e-16 ***
T:L 1 143.60 120.77 < 2.2e-16 ***
---
Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
Then, I got the coefficients:
coef(twoway)
$`(Intercept)`
[1] 5.866527
$A
BB AA CC
0.27277069 -0.01475892 -0.25801177
$T
On Off
1.156143 -1.156143
$L
Yes No
-1.225228 1.225228
$A.T
T
A On Off
BB 0.4809533 -0.4809533
AA -0.1783651 0.1783651
CC -0.3025882 0.3025882
$A.L
L
A Yes No
BB 0.19166429 -0.19166429
AA -0.15632604 0.15632604
CC -0.03533825 0.03533825
$T.L
L
T Yes No
On 0.2933774 -0.2933774
Off -0.2933774 0.2933774
I, then, hand-calculated odds ratio for A x T, A x L, and T x L.
T x L:
*?TL *=* e4(.293) *= 3.23
A x L:
*?AL(BB vs. AA) *= *e 2(.19166) + 2(.1563) = *2.01
*?AL(BB vs. CC) *= *e 2(.19166) + 2(.03533) = *1.57
A x T:
*?AT(BB vs. AA) *= *e 2(.48095) + 2(.17837) = 3.74*
*
*
*?AT(BB vs. CC) = e 2(.48095) + 2(.30259) = 4.79 *
Now, I'd like to know if BB and AA (or BB and CC) are significantly
different from each other (i.e., the odds of BB to be 2.01 times larger than
AA is significant) and the differences between BB and CC are significant
(i.e., the odds of BB to be 1.6 times bigger is significant) etc.
It will be easier to answer this if you formulate the problem as a
surrogate Poisson regression via glm().
The relationship between the surrogate Poisson and the loglinear model is
well known. See V&R or McCullagh and Nelder, for example.
Then the hypotheses about odds ratios become hypotheses about
coefficients, which you can test via summary(), or linear combinations of
coefficients, which you can test with the pieces that vcov() and coef()
provide.
HTH,
Chuck
I'd really appreciate if someone can answer this question!
Thank you,
Sachi
[[alternative HTML version deleted]]
Charles C. Berry (858) 534-2098
Dept of Family/Preventive Medicine
E mailto:cbe...@tajo.ucsd.edu UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
