Hi,Achim:
Can all the analysis you mentioned via loglm be performed via glm(...,family=poisson)? 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 variable. >> >> 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 variable. > >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.html >> 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.