# Mentel-Haneszel Method with Odds Ratio Estimation

# Copy and paste this R code into R Commander Script Window
# Update the count data and labels and submit it for result 

library("lawstat")

boymatrix <- matrix(c(20,100,15,150),nrow=2,byrow=TRUE,   # enter data counts
            dimnames = list("Sleep"  = c("Low","High"),   # label column variable categories
                            "Result" = c("Fail","Pass"))) # label row variable categories
print(boymatrix)


girlmatrix <- matrix(c(30,100,25,200),nrow=2,byrow=TRUE,  # enter data counts
            dimnames = list("Sleep"  = c("Low","High"),   # label column variable categories
                            "Result" = c("Fail","Pass"))) # label row variable categories
print(girlmatrix)


myarray <- array(c(boymatrix,girlmatrix),dim=c(2,2,2))    # prepare data for M-H method
cmh.test(myarray)                                         # run Mantel-Haneszel Method


#Tarones Test for Homogeneity of Odds Ratios (need R version 2.12.2 or higher)
library("metafor")

calcTaronesTest <- function(mylist,referencerow=2)
 {
    numstrata <- length(mylist)
    # make an array "ntrt" of the number of people in the exposed group, in each stratum
    # make an array "nctrl" of the number of people in the unexposed group, in each stratum
    # make an array "ptrt" of the number of people in the exposed group that have the disease,
    # in each stratum
    # make an array "pctrl" of the number of people in the unexposed group that have the disease,
    # in each stratum
    # make an array "htrt" of the number of people in the exposed group that don't have the
    # disease, in each stratum
    # make an array "hctrl" of the number of people in the unexposed group that don't have the
    # disease, in each stratum
    ntrt <- vector()
    nctrl <- vector()
    ptrt <- vector()
    pctrl <- vector()
    htrt <- vector()
    hctrl <- vector()
    if (referencerow == 1) { nonreferencerow <- 2 }
    else                   { nonreferencerow <- 1 }
    for (i in 1:numstrata)
    {
       mymatrix <- mylist[[i]]
       DiseaseUnexposed <- mymatrix[referencerow,1]
       ControlUnexposed <- mymatrix[referencerow,2]
       totUnexposed <- DiseaseUnexposed + ControlUnexposed
       nctrl[i] <- totUnexposed
       pctrl[i] <- DiseaseUnexposed
       hctrl[i] <- ControlUnexposed
       DiseaseExposed <- mymatrix[nonreferencerow,1]
       ControlExposed <- mymatrix[nonreferencerow,2]
       totExposed <- DiseaseExposed + ControlExposed
       ntrt[i] <- totExposed
       ptrt[i] <- DiseaseExposed
       htrt[i] <- ControlExposed
    }
    # calculate Tarone's test of homogeneity, using the rma.mh function from the
    # "metafor" package
    tarone <- rma.mh(ptrt, htrt, pctrl, hctrl, ntrt, nctrl)
    pvalue <- tarone$TAp
    print(paste("Pvalue for Tarone's test =", pvalue))
}

mylist <- list(boymatrix,girlmatrix) # list all data matrices
calcTaronesTest(mylist)