[R] quesion about model g1
Achim Zeileis
Achim.Zeileis at uibk.ac.at
Mon Apr 29 08:33:25 CEST 2013
On Mon, 29 Apr 2013, meng wrote:
> Hello Achim:
> Sorry for another question about the model g1 in the last mail.
>
> As to model g2 and g3:
> g2 <- glm(Freq ~ (age + drug + case)^2, data = df, family = poisson)
> g3 <- glm(Freq ~ age * drug * case, data = df, family = poisson)
> anova(g2, g3, test = "Chisq")
>
> I know clearly that the only difference between g2 and g3 is that g2 has
> no 3-way interaction while g3 has,and anova tests whether this "only
> difference(i.e. 3-way interaction)" is significant or not.
>
> But as to g1 and g3:
> g1 <- glm(Freq ~ age/(drug + case), data = df, family = poisson)
>
> I can't find out the "only difference" between g1 and g3,so I don't know
> what "anova(g1, g3, test = "Chisq")" tests for. Also, what "/" sign
> following age in g1 refers to?
The "/" could be replaced by a "*" here and the fitted values and
corresponding log-likelihood would not change. Only the coefficients
change: "/" induces a nested coding while "*" employs the interaction
coding.
Breaking everything down to main and interaction effects (and ignoring the
particular coding of the coefficients), the three models are
g1: a + d + c + a:d + a:c
g2: a + d + c + a:d + a:c + d:c
g3: a + d + c + a:d + a:c + d:c + a:d:c
with interpretations:
g1: conditional independence of drug and case given age
g2: no three-way interaction (case depends on drug but in the same way for
different levels of age)
g3: saturated model
>
> Many thanks and sorry for many quesions.
>
>
> Best
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> At 2013-04-24 22:22:55,"Achim Zeileis" <Achim.Zeileis at 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.Zeileis at 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 at 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 at 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.
>
More information about the R-help
mailing list