R/ratio.median.is.R

Defines functions ci.ratio.median.is

Documented in ci.ratio.median.is

# DGB
## Ratio of Medians from Two (Independent) Samples

ci.ratio.median.is <- function(alpha, y1, y2) {
 # Computes confidence interval for a ratio of
 # population medians of ratio scale measurements
 # in a 2-group design.
 # Arguments:
 #   alpha: alpha level for 1-alpha confidence
 #   y1:    vector of scores for group 1
 #   y2:    vector of scores for group 2  
 # Values:
 #   medians, median ratio, lower limit, upper limit
 z <- qnorm(1 - alpha/2)
 n1 <- length(y1)
 y1 <- sort(y1)
 n2 <- length(y2)
 y2 <- sort(y2)
 median1 <- median(y1)
 median2 <- median(y2)
 y1 <- log(y1 + 1)
 y2 <- log(y2 + 1)
 a1 <- round((n1 + 1)/2 - sqrt(n1))
 if (a1 < 1) {a1 = 1}
 L1 <- log(exp(y1[a1]) - 1)
 U1 <- log(exp(y1[n1 - a1 + 1]) - 1)
 p <- pbinom(a1 - 1, size = n1, prob = .5)
 z0 <- qnorm(1 - p)
 se1 <- (U1 - L1)/(2*z0)
 a2 <- round((n2 + 1)/2 - sqrt(n2))
 if (a2 < 1) {a2 = 1}
 L2 <- log(exp(y2[a2]) - 1)
 U2 <- log(exp(y2[n2 - a2 + 1]) - 1)
 p <- pbinom(a2 - 1, size = n2, prob = .5)
 z0 <- qnorm(1 - p)
 se2 <- (U2 - L2)/(2*z0)
 se <- sqrt(se1^2 + se2^2)
 diff <- log(median1) - log(median2)
 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)
}
cwendorf/dgb documentation built on May 3, 2022, 9:35 p.m.