[R] test Breslow-Day for svytable??

[Thomas Lumley]

On Sat, Sep 1, 2012 at 4:27 AM, David Winsemius <dwinsemius at comcast.net> wrote:
>
> 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.)
>

Yes, that will give internal consistency from a data structure point
of view.  It won't give a valid test in real examples, though --
epi.2by2 doesn't know about complex sampling, and what you're passing
it is just an estimate of the population 2x2xK table.

What would work, though it's not quite the same as the Breslow-Day
test, is to use svyloglin() and do a Rao-Scott test comparing the
model with all two-way interactions ~(Gender+Dept+Admit)^2 to the
saturated model ~Gender*Dept*Admit.

    -thomas


-- 
Thomas Lumley
Professor of Biostatistics
University of Auckland