[R] test Breslow-Day for svytable??

[David Winsemius]

On Aug 31, 2012, at 7:20 AM, Diana Marcela Martinez Ruiz wrote:

> Hi all,
> 
> I want to know how to perform the test Breslow-Day test for homogeneity of 
> odds ratios (OR) stratified for svytable. This test is obtained with the following code:
> 
> epi.2by2 (dat = daty, method = "case.control" conf.level = 0.95,

missing comma here ...........................^

>    units = 100, homogeneity = "breslow.day", verbose = TRUE)
> 
> where "daty" is the object type table  svytable consider it, but when I run the code
> does not throw the homogeneity test.

You are asked in the Posting guide to copy all errors and warnings when asking about unexpected behavior. When I run epi.2y2 on the output of a syvtable object I get no errors, but I do get warnings which I think are due to non-integer entries in the weighted table. I also get from a svytable() usingits first example on the help page an object that is NOT a set of 2 x 2 tables in an array of the structure as expected by epi.2by2(). The fact that epi.2by2() will report numbers with labels for a 2 x 3 table means that its error checking is weak.

This is the output of str(dat) from one of the example on epi.2by2's help page:

> str(dat)
 table [1:2, 1:2, 1:3] 41 13 6 53 66 37 25 83 23 37 ...
 - attr(*, "dimnames")=List of 3
  ..$ Exposure: chr [1:2] "+" "-"
  ..$ Disease : chr [1:2] "+" "-"
  ..$ Strata  : chr [1:3] "20-29 yrs" "30-39 yrs" "40+ yrs"

Notice that is is a 2 x 2 x n array. (Caveat:: from here on out I am simply reading the help pages and using str() to look at the objects created to get an idea regarding success or failure. I am not an experienced user of either package.)  I doubt that  what you got from svytable is a 2 x 2 table. As another example you can build a 2 x 2 x n table from the built-in dataset: UCBAdmissions 

DF <- as.data.frame(UCBAdmissions)
## Now 'DF' is a data frame with a grid of the factors and the counts
## in variable 'Freq'.
dat2 <- xtabs(Freq ~ Gender + Admit+Dept, DF)
epiR::epi.2by2(dat = dat2, method = "case.control", conf.level = 0.95, 
     units = 100, homogeneity = "breslow.day", verbose = TRUE)$OR.homog
#-------------
  test.statistic df    p.value
1       18.82551  5 0.00207139

Using svydesign and svytable I _think_ this is how one would go about constructing a 2 x 2 table:

tbl2<-svydesign(  ~ Gender + Admit+Dept, weights=~Freq, data=DF)
  summary(dclus1)
(tbl2by2 <- svytable(~ Gender + Admit+Dept, tbl2))
 epiR::epi.2by2(dat = tbl, method = "case.control", conf.level = 0.95, 
     units = 100, homogeneity = "breslow.day", verbose = TRUE)$OR.homog
#-----------------------
  test.statistic df    p.value
1       18.82551  5 0.00207139

(At least I got internal consistency. I see you copied Thomas Lumley, which is a good idea. I'll be happy to get corrected on any point. I'm adding the maintainer of epiR to the recipients.)

-- 
David Winsemius, MD
Alameda, CA, USA