R Code for finding the lower/upper confidence limit of
a likelihood ratio test (wilks) based confidence interval.

The data are two samples of binomial type data (success/failure)
The parameter is the ratio of two success probabilities.

There are basically two pieces of computations to do: the code to do
likelihood ratio test, which should return the -2LLR; and the code to find
the upper and lower confidence limit of the parameter being tested.

We can either write two functions, one for the likelihood ratio test, and one to find
upper and lower confidence limit. Or we can write one function to do both.


The following function is [one function] specifically for the situation for data n1, x1, n2, x2
(where n1 is sample size for sample one, and x1 is the observed number of success in sample one, etc)
We are finding the confidence interval for p2/p1.  [not p1/p2]


==================== COPY and PASTE into R ==========================================================

BinomRatio <- function(n1, x1, n2, x2, step=0.05, initStep=0, level=3.84) {
#######################
    if( x1 >= n1 ) stop("sample one is 100% success?")
    if( x2 >= n2 ) stop("sample two is 100% success?")
    if( x1 <= 0 ) stop("x1 must > 0")
    if( x2 <= 0 ) stop("x2 must > 0")
MLE <- (x2*n1)/(x1*n2)   #### recall we are to find CI for p2/p1, not p1/p2.

Theta <- function(theta, n1=n1, x1=x1, n2=n2, x2=x2)
  {
       llr <- function(const, n1, x1, n2, x2, theta) {
              npllik1 <- -2*(x1*log(const) + (n1-x1)*log(1-const) -
              x1*log(x1/n1) - (n1-x1)*log(1-x1/n1))
              npllik2 <-
              -2*(x2*log(const*theta)+(n2-x2)*log(1-const*theta)-x2*log(x2/n2)-(n2-x2)*log(1-x2/n2))
        return(npllik1 + npllik2)
        }
       upBD <- min( 0.999999, 1/theta -0.000001)
       temp <- optimize(f = llr, lower = 0.000001, upper = upBD, n1 = n1, x1 = x1, n2 = n2, x2 = x2,theta = theta)

    cstar <- temp$minimum
    val <- temp$objective
    pvalue <- 1 - pchisq( val, df=1)
    list(`-2LLR` = val, cstar = cstar, Pval=pvalue)
  }
value <- 0
step1 <- step
Lbeta <- MLE - initStep
for( i in 1:8 ) {
while(value < level) {
Lbeta <- Lbeta - step1
value <- Theta(theta = Lbeta, n1=n1, x1=x1, n2=n2, x2=x2)$'-2LLR'
}
Lbeta <- Lbeta + step1
step1 <- step1/10
value <- Theta( theta =Lbeta, n1=n1, x1=x1, n2=n2, x2=x2)$'-2LLR'
}

value1 <- value
value <- 0

Ubeta <- MLE + initStep
for( i in 1:8 ) {
   while(value < level) {
    Ubeta <- Ubeta + step
    value <- Theta(theta=Ubeta, n1=n1, x1=x1, n2=n2, x2=x2 )$'-2LLR'
   }
   Ubeta <- Ubeta - step
   step <- step/10
   value <- Theta(theta=Ubeta, n1=n1, x1=x1, n2=n2, x2=x2)$'-2LLR'
}
return( list(Low=Lbeta, Up=Ubeta, FstepL=step1, FstepU=step,
Lvalue=value1, Uvalue=value) )
}




===================COPY AND PASTE above into R==================================================

Here are two examples how it is run and getting results:

> BinomRatio(n1=100, x1=50, n2=120, x2=65)
$Low
[1] 0.8406968

$Up
[1] 1.409807

$FstepL
[1] 5e-10

$FstepU
[1] 5e-10

$Lvalue
[1] 3.84

$Uvalue
[1] 3.84

> BinomRatio(n1=100, x1=50, n2=120, x2=65, level=2.7055)
$Low
[1] 0.8756334

$Up
[1] 1.34966

$FstepL
[1] 5e-10

$FstepU
[1] 5e-10

$Lvalue
[1] 2.7055

$Uvalue
[1] 2.7055


The nice feature of the Wilks confidence interval is that when you switch the label of sample 1 and sample 2
then the new confidence interval is just the reciprocal of the original confidence interval---someting we shall call
"invariance property".

The Wald confidence interval do not have "invariance property". There for the confidence interval for p2/p1 and
for p1/p2 are NOT just a reciprocal of each other.


So, try

BinomRatio(n1=120, x1=65, n2=100, x2=50)

you should get a confidence interval that is the reciprocal of the one in above.