Uncertainty of fitness ratios

Version: 2012-12-19 20:03:40

From Karen C.:

q.mean = 0.915
q.var = 0.00014
h.lo.mean = 0.837
h.lo.var = 0.000266
h.hi.mean = 0.9897
h.hi.var = 0.000199

Delta method by hand (i.e. using \[ \mbox{CV}^2(q/h) = \mbox{CV}^2(q)+\mbox{CV}^2(h) \\ \mbox{Var}(q/h)/(q^2/h^2) = \mbox{Var}(q)/q^2 + \mbox{Var}(h)/h^2 \\ \mbox{Var}(q/h) = 1/h^2 \mbox{Var}(q) + q^2/h^4 \mbox{Var}(h) \]

(qh.lo.mean <- q.mean/h.lo.mean)
## [1] 1.093
(qh.lo.var <- q.var/(h.lo.mean)^2 + h.lo.var * q.mean^2/h.lo.mean^4)
## [1] 0.0006536
(qh.lo.se <- sqrt(qh.lo.var))
## [1] 0.02557

Using an automated delta-variance calculator from the emdbook package:

library(emdbook)
deltavar(q/h, mean = c(q = q.mean, h = h.lo.mean), Sigma = c(q.var, h.lo.var))
## [1] 0.0006536

Confidence limits:

(qh.lo.ci1 <- qh.lo.mean + c(-1, 1) * 1.96 * qh.lo.se)
## [1] 1.043 1.143

Monte Carlo approach: pick many random normal deviates for \( q \) and for \( h \) and calculate their ratio:

mcsim_N <- 10000
set.seed(101)
q.mc <- rnorm(mcsim_N, q.mean, sqrt(q.var))
h.lo.mc <- rnorm(mcsim_N, h.lo.mean, sqrt(h.lo.var))
qh.lo.mc <- q.mc/h.lo.mc
par(las = 1, bty = "l")
hist(qh.lo.mc, col = "gray", breaks = 50, main = "", freq = FALSE)
curve(dnorm(x, mean = mean(qh.lo.mc), sd = sd(qh.lo.mc)), add = TRUE, col = 2)
abline(v = 1, col = "purple", lwd = 2)
curve(dnorm(x, mean = qh.lo.mean, sd = qh.lo.se), add = TRUE, col = 4)

plot of chunk unnamed-chunk-2

sd(qh.lo.mc)  ## estimated std err from MC
## [1] 0.02536
quantile(qh.lo.mc, c(0.025, 0.975))  ## Monte Carlo confidence intervals
##  2.5% 97.5% 
## 1.045 1.145

(The red line is the normal approximation to the MC distribution: the blue line is the distribution assumed/estimated by the delta method.) This indicates that the delta method is not too bad …

MC for high case:

h.hi.mc <- rnorm(mcsim_N, h.hi.mean, sqrt(h.hi.var))
qh.hi.mc <- q.mc/h.hi.mc

If you want do this for \( q/(1-q) \cdot (1-h)/h \), it will be easiest to use the automated deltavar calculation or the MC approach (I could figure out the equation by hand but I don't want to bother), e.g.

qh2.mean <- q.mean/h.lo.mean * (1 - h.lo.mean)/(1 - q.mean)
qh2.mc <- q.mc/h.lo.mc * (1 - h.lo.mc)/(1 - q.mc)
quantile(qh2.mc, c(0.025, 0.975))
##  2.5% 97.5% 
## 1.461 3.151
qh2.dv <- deltavar(q/h * (1 - h)/(1 - q), mean = c(q = q.mean, h = h.lo.mean), 
    Sigma = c(q.var, h.lo.var))
qh2.se <- sqrt(qh2.dv)
par(las = 1, bty = "l")
hist(qh2.mc, col = "gray", breaks = 50, main = "", freq = FALSE, xlim = c(1, 
    3.5))
curve(dnorm(x, mean = mean(qh2.mc), sd = sd(qh2.mc)), add = TRUE, col = 2)
curve(dnorm(x, mean = qh2.mean, sd = qh2.se), add = TRUE, col = 4)
abline(v = 1, col = "purple", lwd = 2)

plot of chunk f2

library(ggplot2)
d <- data.frame(qh = c(qh.lo.mc, qh.hi.mc), v = rep(c("low", "high"), each = mcsim_N))
theme_set(theme_bw())
ff <- function(x) {
    result <- data.frame(y = mean(x), ymin = mean(x) - 1.96 * sd(x), ymax = mean(x) + 
        1.96 * sd(x))
}
ggplot(d, aes(x = v, y = qh)) + geom_violin(fill = "gray") + geom_hline(yintercept = 1, 
    colour = "red", lwd = 2, alpha = 0.5) + stat_summary(fun.data = ff, colour = "blue", 
    size = 2) + labs(x = "", y = "Relative fitness of type I males")

plot of chunk violins