Understand.
Many thanks for your detailed explaination.
At 2013-04-24 22:22:55,"Achim Zeileis" <achim.zeil...@uibk.ac.at> wrote:
>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.
>> >>
>>
>>
>>
>>
[[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.html
and provide commented, minimal, self-contained, reproducible code.
