On Wed, 24 Apr 2013, meng wrote:
Hi,Achim:
Can all the analysis you mentioned via loglm be performed via
glm(...,family=poisson)?
Yes.
## transform table back to data.frame
df <- as.data.frame(tab)
## fit models: conditional independence, no-three way interaction,
## and saturated
g1 <- glm(Freq ~ age/(drug + case), data = df, family = poisson)
g2 <- glm(Freq ~ (age + drug + case)^2, data = df, family = poisson)
g3 <- glm(Freq ~ age * drug * case, data = df, family = poisson)
## likelihood-ratio tests (against saturated)
anova(g1, g3, test = "Chisq")
anova(g2, g3, test = "Chisq")
## compare fitted frequencies (which are essentially equal)
all.equal(as.numeric(fitted(g1)),
as.data.frame(as.table(fitted(m1)))$Freq)
all.equal(as.numeric(fitted(g2)),
as.data.frame(as.table(fitted(m2)))$Freq)
The difference is mainly that loglm() has a specialized user interface and
that it uses a different optimizer (iterative proportional fitting rather
than iterative reweighted least squares).
Best,
Z
Many thanks.
At 2013-04-24 19:37:10,"Achim Zeileis" <achim.zeil...@uibk.ac.at> wrote:
>On Wed, 24 Apr 2013, meng wrote:
>
>> Hi all:
>> For stratified count data,how to perform regression analysis?
>>
>> My data:
>> age case oc count
>> 1 1 1 21
>> 1 1 2 26
>> 1 2 1 17
>> 1 2 2 59
>> 2 1 1 18
>> 2 1 2 88
>> 2 2 1 7
>> 2 2 2 95
>>
>> age:
>> 1:<40y
>> 2:>40y
>>
>> case:
>> 1:patient
>> 2:health
>>
>> oc:
>> 1:use drug
>> 2:not use drug
>>
>> My purpose:
>> Anaysis whether case and oc are correlated, and age is a stratified varia
ble.
>>
>> My solution:
>> 1,Mantel-Haenszel test by using function "mantelhaen.test"
>> 2,loglinear regression by using function glm(count~case*oc,family=poisson
).But I don't know how to handle variable "age",which is the stratified vari
able.
>
>Instead of using glm(family = poisson) it is also convenient to use
>loglm() from package MASS for the associated convenience table.
>
>The code below shows how to set up the contingency table, visualize it
>using package vcd, and then fit two models using loglm. The models
>considered are conditional independence of case and drug given age and the
>"no three-way interaction" already suggested by Peter. Both models are
>also accompanied by visualizations of the residuals.
>
>## contingency table with nice labels
>tab <- expand.grid(drug = 1:2, case = 1:2, age = 1:2)
>tab$count <- c(21, 26, 17, 59, 18, 88, 7, 95)
>tab$age <- factor(tab$age, levels = 1:2, labels = c("<40", ">40"))
>tab$case <- factor(tab$case, levels = 1:2, labels = c("patient",
>"healthy"))
>tab$drug <- factor(tab$drug, levels = 1:2, labels = c("yes", "no"))
>tab <- xtabs(count ~ age + drug + case, data = tab)
>
>## visualize case explained by drug given age
>library("vcd")
>mosaic(case ~ drug | age, data = tab,
> split_vertical = c(TRUE, TRUE, FALSE))
>
>## test wheter drug and case are independent given age
>m1 <- loglm(~ age/(drug + case), data = tab)
>m1
>
>## visualize corresponding residuals from independence model
>mosaic(case ~ drug | age, data = tab,
> split_vertical = c(TRUE, TRUE, FALSE),
> residuals_type = "deviance",
> gp = shading_hcl, gp_args = list(interpolate = 1.2))
>mosaic(case ~ drug | age, data = tab,
> split_vertical = c(TRUE, TRUE, FALSE),
> residuals_type = "pearson",
> gp = shading_hcl, gp_args = list(interpolate = 1.2))
>
>## test whether there is no three-way interaction
>## (i.e., dependence of case on drug is the same for both age groups)
>m2 <- loglm(~ (age + drug + case)^2, data = tab)
>m2
>
>## visualization (with default pearson residuals)
>mosaic(case ~ drug | age, data = tab,
> expected = ~ (age + drug + case)^2,
> split_vertical = c(TRUE, TRUE, FALSE),
> gp = shading_hcl, gp_args = list(interpolate = 1.2))
>
>
>> Many thanks for your help.
>>
>> My best.
>> [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> 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.h
tml
>> and provide commented, minimal, self-contained, reproducible code.
>>
______________________________________________
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.
