R help - Apr 2013 - Regression on stratified count data

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.
 
Many thanks for your help.

My best.
	[[alternative HTML version deleted]]

On Apr 24, 2013, at 06:15 , 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.
The canonical way is to fit the model without 2nd order interaction: 

count~case*oc*age-case:oc:case  . 

(It may take the back of an envelope or two to realize that this is equivalent
to the common OR assumption of the MH test.)

Alternatively, use logistic regression 

glm(case ~ oc + age, family=binomial, weight=count, data=dd)

(NB: it is important that case is a factor here!)

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com

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 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.
>