DGB: DGB

# DGB
## Ratio of Medians from Paired Samples

ci.ratio.median.ps <- function(alpha, y1, y2) {
 # Computes confidence interval for a ratio of
 # population medians in a 2-level within-subjects design
 # Arguments:
 #   alpha: alpha level for 1-alpha confidence
 #   y1:    vector of scores for level 1
 #   y2:    vector of scores for level 2  
 # Values:
 #   medians, median ratio, lower limit, upper limit
 z <- qnorm(1 - alpha/2)
 n <- length(y1)
 y1 <- sort(y1)
 y2 <- sort(y2)
 median1 <- median(y1)
 median2 <- median(y2)
 y1 <- log(y1 + 1)
 y2 <- log(y2 + 1)
 a <- round((n + 1)/2 - sqrt(n))
 if (a < 1) {a1 = 1}
 p <- pbinom(a - 1, size = n, prob = .5)
 z0 <- qnorm(1 - p)
 L1 <- log(exp(y1[a]) - 1)
 U1 <- log(exp(y1[n - a + 1]) - 1)
 se1 <- (U1 - L1)/(2*z0)
 L2 <- log(exp(y2[a]) - 1)
 U2 <- log(exp(y2[n - a + 1]) - 1)
 se2 <- (U2 - L2)/(2*z0)
 a1 <- (y1 < log(median1))
 a2 <- (y2 < log(median2))
 a3 <- a1 + a2
 a4 <- sum(a3 == 2)
 if (n/2 == trunc(n/2)) {
   p00 <- (sum(a4) + .25)/(n + 1)
 } else {
   p00 <- (sum(a4) + .25)/n 
 }
 cov <- (4*p00 - 1)*se1*se2
 diff <- log(median1) - log(median2)
 se <- sqrt(se1^2 + se2^2 - 2*cov)
 L <- exp(diff - z*se)
 U <- exp(diff + z*se)
 out <- t(c(median1, median2, exp(diff), L, U))
 colnames(out) <- c("Median1", "Median2", "Median1/Median2", "LL", "UL")
 return(out)
}