R/median.ps.R

In cwendorf/dgb: DGB

# DGB
## Median Difference from Paired Samples

ci.median.ps <- function(alpha, y1, y2) {
 # Computes confidence interval for a difference 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:
 #   sample medians, SE of difference, lower limit, upper limit
 z <- qnorm(1 - alpha/2)
 n <- length(y1)
 y1 <- sort(y1)
 y2 <- sort(y2)
 median1 <- median(y1)
 median2 <- median(y2)
 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 <- y1[a]
 U1 <- y1[n - a + 1]
 se1 <- (U1 - L1)/(2*z0)
 L2 <- y2[a]
 U2 <- y2[n - a + 1]
 se2 <- (U2 - L2)/(2*z0)
 a1 <- (y1 < median1)
 a2 <- (y2 < 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 <- median1 - median2
 se <- sqrt(se1^2 + se2^2 - 2*cov)
 L <- diff - z*se
 U <- diff + z*se
 out <- t(c(median1, median2, diff, se, L, U))
 colnames(out) <- c("Median1", "Median2", "Median1-Median2", "SE", "LL", "UL")
 return(out)
}