#The following MH functions were modified from code originally 
#developed by former student Imran Quraishi

#for rate ratios
mh.test.rate <- function (e1,y1,e2,y2) { # (in terms of person years)
        et <- e1+e2
        yt <- y1+y2
        rr.mh <- sum(e2*y1/yt)/sum(e1*y2/yt)
        se.log.rr <- sqrt(sum(y1*et/yt^2*y2))/(sqrt(rr.mh)*sum(y1*et/(yt*(y1+rr.mh*y2))*y2))
        rr.mh.ul <- exp(log(rr.mh)+1.96*se.log.rr)
        rr.mh.ll <- exp(log(rr.mh)-1.96*se.log.rr)
        e1.est <- et*y1/(y1+rr.mh*y2)
        e2.est <- et*y2*rr.mh/(y1+rr.mh*y2)
        homog.test.chisq <- sum((e1-e1.est)^2/e1.est)+sum((e2-e2.est)^2/e2.est)
                 s<-length(e1) #number of statra
        homog.test.p <- 1-pchisq(homog.test.chisq,s-1)
        rr.test.chisq <- (sum(e2)-sum(y2*et/yt))^2/sum(y1*et/yt^2*y2)
        rr.test.p <- 1-pchisq(rr.test.chisq,1)
        ret.val<-list(rr.mh,se.log.rr,rr.mh.ll,rr.mh.ul,homog.test.chisq,homog.test.p,rr.test.chisq,
          rr.test.p)
        names(ret.val)<-c('rr.mh','se.log.rr','rr.mh.ll', 'rr.mh.ul', 'homog.test',
          'homog.test.p', 'rr.test.chisq', 'rr.test.p')
        ret.val
}

#data from Woodward pg 231 table 5.13, Coronary Events by Housing Status and Age
er<-c(2,24,31,28,28,2) #coronary events, renters 
eo<-c(3,19,25,27,26,4) #coronary events, owners
yr<-c(1107447,3058986,3506530,3756650,2419622,351710) #person-years, renters
yo<-c(1619328,4550166,4857904,4536832,2680843,356394) #person-years, owners
round((er/yr)/(eo/yo),2)
mh.test.r(eo,yo,er,yr)

# For odds ratios:
mh.test.odds <- function (a,b,c,d) {
        n <- a+b+c+d
        or.mh <- sum(a*d/n)/sum(b*c/n)
        P <- (a+d)/n
        R <- a*d/n
        Q <- (b+c)/n
        S <- b*c/n
        se.log.or <- sqrt( sum(P*R)/2/(sum(R))^2 + (sum(P*S)+sum(Q*R))/2/sum(R)/sum(S)
          +  sum(Q*S)/2/(sum(S))^2 )
        or.mh.ll <- exp(log(or.mh)-1.96*se.log.or)
        or.mh.ul <- exp(log(or.mh)+1.96*se.log.or)

        O <- a+c
        F. <- a+b
        Obar <- b+d
        Fbar <- c+d
        E <- O*F./n
        V <- O*Obar/(n^2*(n-1))*F.*Fbar
                 cmh.test.chisq <- (sum(a)-sum(E))^2/sum(V)
        cmh.test.p <- 1-pchisq(cmh.test.chisq,1)
        P <- (or.mh-1)*(F.+O)+n
        Estar.a <- (P+sqrt(P^2-4*or.mh*(or.mh-1)*F.*O))/2/(or.mh-1)
        Estar.b <- F.-Estar.a
        Estar.c <- O-Estar.a
        Estar.d <- n-Estar.a-Estar.b-Estar.c
        for (i in 1:length(a)) {
                if (Estar.a[i]<=0 | Estar.b[i]<=0 | Estar.c[i]<=0 | Estar.d[i]<=0) {
                        Estar.a[i]<-(P[i]-sqrt(P[i]^2-4*or.mh*(or.mh-1)*F.[i]*O[i]))/2/(or.mh-1)
                        Estar.b[i] <- F.[i]-Estar.a[i]
                        Estar.c[i] <- O[i]-Estar.a[i]
                        Estar.d[i] <- n[i]-Estar.a[i]-Estar.b[i]-Estar.c[i]
                }
        }
        Vstar <- 1/(1/Estar.a+1/Estar.b+1/Estar.c+1/Estar.d)
        homog.test.chisq <- sum((a-Estar.a)^2/Vstar)
                 s<-length(a) #number of statra
        homog.test.p <- 1-pchisq(homog.test.chisq,s-1)
        ret.val<-list(or.mh,or.mh.ll,or.mh.ul,cmh.test.chisq,cmh.test.p,homog.test.chisq,
          homog.test.p)
        names(ret.val) <- c('or.mh','or.mh.ll','or.mh.ul','cmh.test.chisq','cmh.test.p',
          'homog.test.chisq','homog.test.p')
        ret.val
}
# For risk ratios:
mh.test.risk <- function (a,b,c,d) {
        n <- a+b+c+d
        O <- a+c
        F. <- a+b
        Obar <- b+d
        Fbar <- c+d
         rr.mh<-sum(a*Fbar/n)/sum(c*F./n)
                 num<-sum((F.*Fbar*O-a*c*n)/n^2)
                 den<-sum(a*Fbar/n)*sum(c*F./n)
                 se.log.rr<-sqrt(num/den)
        rr.mh.ll <- exp(log(rr.mh)-1.96*se.log.rr)
        rr.mh.ul <- exp(log(rr.mh)+1.96*se.log.rr)

        E <- O*F./n
        V <- O*Obar/(n^2*(n-1))*F.*Fbar
                 cmh.test.chisq <- (sum(a)-sum(E))^2/sum(V)
        cmh.test.p <- 1-pchisq(cmh.test.chisq,1)
       
        Estar.a <- F.*O*rr.mh/(Fbar+F.*rr.mh)
        Estar.b <- F.-Estar.a
        Estar.c <- O-Estar.a
        Estar.d <- n-Estar.a-Estar.b-Estar.c
 
        Vstar <- 1/(1/Estar.a+1/Estar.b+1/Estar.c+1/Estar.d)
        homog.test.chisq <- sum((a-Estar.a)^2/Vstar)
                 s<-length(a) #number of statra
        homog.test.p <- 1-pchisq(homog.test.chisq,s-1)
        ret.val<-list(rr.mh,rr.mh.ll,rr.mh.ul,cmh.test.chisq,cmh.test.p,homog.test.chisq,
          homog.test.p)
        names(ret.val) <- c('rr.mh','rr.mh.ll','rr.mh.ul','cmh.test.chisq','cmh.test.p',
          'homog.test.chisq','homog.test.p')
        ret.val
}

#Example 1 from notes, Woodward pg 179, table 4.16, Ex. 4.15, pg 182
mh.test.risk(c(67,8),c(2061,51),c(46,11),c(3454,41)) 
mh.test.odds(c(67,8),c(2061,51),c(46,11),c(3454,41)) 
#Example 2 from notes, Woodward pg 169, table 4.13 pg 169, Ex. 4.9 and Ex. 4.11 
mh.test.risk(c(33,52),c(923,898),c(48,29),c(1722,678))