28

I want to derive the limits for the $100(1-\alpha)\%$ confidence interval for the ratio of two means.
Suppose, $X_1 \sim N(\theta_1, \sigma^2)$ and $X_2 \sim N(\theta_2, \sigma^2)$ being independent, the mean ratio $\Gamma = \theta_1/\theta_2$. I tried to solve: $$\text{Pr}(-z(\alpha/2)) \leq X_1 - \Gamma X_2 / \sigma \sqrt {1 + \gamma^2} \leq z(\alpha/2)) = 1 - \alpha$$ but that equation couldn't be solved for many cases (no roots). Am I doing something wrong? Is there a better approach? Thanks

kjetil b halvorsen
  • 63,378
  • 26
  • 142
  • 467
francogrex
  • 479
  • 1
  • 4
  • 8
  • 1
    The problem is that the ratio of two numbers from two normal distributions follows [Cauchy distribution](http://en.wikipedia.org/wiki/Cauchy_distribution) and thus the variance is undefined. –  Oct 02 '11 at 08:58
  • 7
    @mbq - the Cauchy distribution presents no problems for confidence intervals, as the CDF is the inverse tangent function. Variance need not be defined for CIs to work. And the ratio of two normal RVs with zero mean is Cauchy, but not necessarily two normal RVs with non-zero mean. – probabilityislogic Oct 02 '11 at 10:58
  • 1
    @probabilityislogic Sure, I must stop trying to think on Sunday mornings. –  Oct 02 '11 at 12:13

4 Answers4

33

Fieller's method does what you want -- compute a confidence interval for the quotient of two means, both assumed to be sampled from Gaussian distributions.

  • The original citation is: Fieller EC: The biological standardization of Insulin. Suppl to J R Statist Soc 1940, 7:1-64.

  • The Wikipedia article does a good job of summarizing.

  • I've created an online calculator that does the computation.

  • Here is a page summarizing the math from the first edition of my Intuitive Biostatistics

Harvey Motulsky
  • 14,903
  • 11
  • 51
  • 98
  • It's very good references, I also like that you actually made a calculator for it (+1). As expected though, in your calculator you clearly state that when the confidence interval of the denominator includes zero, it is not possible to compute the CI of the quotient. I think it's the same that happens when I try solving the quadratic equation. suppose variance is 1, mu1 = 0 and mu2=1, N=10000. It's unsolvable. – francogrex Oct 02 '11 at 16:03
  • 2
    thanks for the online calculator Harvey, I'm a typical biologist with insufficient background in statistics and your calculator was exactly what I needed. – Timtico Jul 05 '12 at 08:29
  • Awesome calculator - exactly what I was looking for. Thanks – Alexander May 30 '15 at 12:36
  • @harvey-motulsky the link to the appendix no longer works. I was wondering is the material from this appendix included in the third edition of Intuitive Biostatistics? – Gabriel Southern May 11 '16 at 01:19
  • @GabrielSouthern Thanks for pointing out the link rot. Fixed. – Harvey Motulsky May 11 '16 at 18:10
  • @HarveyMotulsky, thanks the updated link is useful although without the context from the rest of the book I don't completely understand the equations. I am interested in understanding how to write code to calculate the CI (probably similar to your online calculator). Since that seems too much to ask in a comment I've asked it as a separate question instead: http://stats.stackexchange.com/questions/212100/calculating-speedup-confidence-interval – Gabriel Southern May 11 '16 at 20:25
10

R has the package mratios with the function t.test.ratio.

Gemechis Dilba Djira, Mario Hasler, Daniel Gerhard and Frank Schaarschmidt (2011). mratios: Inferences for ratios of coefficients in the general linear model. R package version 1.3.15. http://CRAN.R-project.org/package=mratios

See also http://www.r-project.org/user-2006/Slides/DilbaEtAl.pdf

1

Also if you want to compute Fieller's confidence interval not using mratios (typically because you don't want a simple lm fit but for example a glmer or glmer.nb fit), you can use the following FiellerRatioCI function, with model the output of the model, aname the name of the numerator parameter, bname the name of the denomiator parameter. You can also use directly the FiellerRatioCI_basic function giving, a, b and the covariance matrix between a and b.

FiellerRatioCI <- function (x, ...) { # generic Biomass Equilibrium Level
    UseMethod("FiellerRatioCI", x)
}
FiellerRatioCI_basic <- function(a,b,V,alpha=0.05){
    theta <- a/b
    v11 <- V[1,1]
    v12 <- V[1,2]
    v22 <- V[2,2]

    z <- qnorm(1-alpha/2)
    g <- (z^2)*v22/b^2
    C <- sqrt(v11 - 2*theta*v12 + theta^2 * v22 - g*(v11-v12^2/v22))
    minS <- (1/(1-g))*(theta- g*v12/v22 - z/b * C)
    maxS <- (1/(1-g))*(theta- g*v12/v22 + z/b * C)
    return(c(ratio=theta,min=minS,max=maxS))
}
FiellerRatioCI.glmerMod <- function(model,aname,bname){
    V <- vcov(model)
    a<-as.numeric(unique(coef(model)$culture[aname]))
    b<-as.numeric(unique(coef(model)$culture[bname]))
    return(FiellerRatioCI_basic(a,b,V[c(aname,bname),c(aname,bname)]))
}
FiellerRatioCI.glm <- function(model,aname,bname){
    V <- vcov(model)
    a <- coef(model)[aname]
    b <- coef(model)[bname]
    return(FiellerRatioCI_basic(a,b,V[c(aname,bname),c(aname,bname)]))
}

Example (based on standard glm basic example):

 counts <- c(18,17,15,20,10,20,25,13,12)
 outcome <- gl(3,1,9)
 treatment <- gl(3,3)
 glm.D93 <- glm(counts ~ outcome + treatment, family = poisson())

 FiellerRatioCI(glm.D93,"outcome2","outcome3")
ratio.outcome2            min            max 
      1.550427      -2.226870      17.880574
cmbarbu
  • 161
  • 1
  • 7
0

You can calculate it through:

  • Fieller's method
  • The Taylor method, also called Delta method: it's easier than Fieller's but will fail if the denominator approaches zero.
  • The Hwang–bootstrap method, a bootstrap technique that does not result in unbounded confidence limits.

Here you can find a thorough description and comparison of these methods.