Nothing
#' Sign test for censored data
#'
#' @description Performs a nonparametric sign test of whether the median difference between two columns of paired censored data equals 0. Uses the Fong adjustment for pairs of equal values.
#' @param xd The first column of data values plus detection limits
#' @param xc The column of censoring indicators, where 1 (or `TRUE`) indicates a detection limit in the xd column, and 0 (or `FALSE`) indicates a detected value in xd.
#' @param yd The second column of data values plus detection limits.
#' @param yc The column of censoring indicators, where 1 (or `TRUE`) indicates a detection limit in the `yd` column, and 0 (or `FALSE`) indicates a detected value in `yd`.
#' @param alternative The usual notation for the alternate hypothesis. Default is `“two.sided”`. Options are `“greater”` or `“less”`.
#' @param printstat Logical `TRUE`/`FALSE` option of whether to print the resulting statistics in the console window, or not. Default is `TRUE.`
#'
#' @importFrom stats dbinom pbinom binom.test
#' @export
#' @return Returns the number of `xd` and `yd` values greater than, less than and tied. Corrected and uncorrected `p-value` for ties also displayed.
#' @references
#' Fong, D.Y.T., Kwan, C.W., Lam, K.F., Lam, K.S.L., 2003. Use of the Sign Test for the Median in the Presence of Ties. The American Statistician 57, 237–240.
#'
#' Helsel, D.R., 2011. Statistics for censored environmental data using Minitab and R, 2nd ed. John Wiley & Sons, USA, N.J.
#'#'
#' @examples
#'
#' data(atrazine)
#'
#' cen_signtest(atrazine$June,atrazine$JuneCen,atrazine$Sept,atrazine$SeptCen)
cen_signtest <- function(xd, xc, yd, yc, alternative="two.sided",printstat=TRUE) {
xname <- deparse(substitute(xd))
yname <- deparse(substitute(yd))
nonas <- na.omit(data.frame(xd, xc, yd, yc))
if(length(nonas[,1]) != length(nonas[,3]))
stop("Lengths of x and y must be the same for paired data.")
if(alternative == "two.sided")
{txt3 <- paste(" alternative hypothesis: true median difference is not = 0"); dir <- 0}
else if(alternative == "less")
{txt3 <- paste(" alternative hypothesis: true median difference < 0"); dir <- (-1)}
else if (alternative == "greater")
{txt3 <- paste(" alternative hypothesis: true median difference > 0"); dir <- 1}
else {stop('Invalid choice for alternative, must match "two.sided", "greater", or "less"')}
N = length(nonas[,1])
for(i in seq(N)) {
if(nonas[i,2] && nonas[i,4]) { # Both censored, make same
xyMax <- max(nonas[i,1], nonas[i,3])
nonas[i,1] <- nonas[i,3] <- xyMax
}
else if(nonas[i,2] && nonas[i,1] > nonas[i,3]) { # x censored > y
nonas[i,3] <- nonas[i,1]
nonas[i,4] <- TRUE
}
else if(nonas[i,4] && nonas[i,3] > nonas[i,1]) { # y censored > x
nonas[i,1] <- nonas[i,3]
nonas[i,2] <- TRUE
}
} # Done, no need to check uncensored observations
# setting <s to slightly below detected obs at that value
nonas[,1] <- nonas[,1] - nonas[,2] * (0.001) * nonas[,1]
nonas[,3] <- nonas[,3] - nonas[,4] * (0.001) * nonas[,3]
x <- nonas[,1]
names(x) <- xname
y <- nonas[,3]
names(y) <- yname
# test.dat <- data.frame(x, nonas[2], y, nonas[4])
# print(test.dat)
x.gt.y <- sum(x > y)
x.lt.y <- sum (y > x)
ties <- length(nonas[,1]) - x.gt.y - x.lt.y
# result.counts <- data.frame(x.gt.y, x.lt.y, ties)
# print(result.counts)
s.binom <- binom.test(x.gt.y, length(nonas[,1])-ties, p=0.5, alternative = alternative)
# LSA version in the coin package. Not currently used.
# s.out <- sign_test(x~y, alternative = alternative)
# txt.s <- paste("n =", N, " Z=", round(statistic(s.out), 4), " p-value =", round(pvalue(s.out), 4))
# print(txt.s)
# Fong's Modified Sign Test correction for ties
dir <- dir*sign(x.gt.y - x.lt.y) #dir=1: alt in same direction as data
num.Fong <- 1-pbinom(max(x.gt.y, x.lt.y)-1, N, prob=0.5)
denom.Fong <- 1-pbinom(floor((N-ties+1)/2)-1, N, prob=0.5)
p.Fong <- num.Fong/denom.Fong
if (dir== 1) {num.Fong <- 1-pbinom((max(x.gt.y, x.lt.y)-1), N, prob=0.5)
denom.Fong <- 2*(1-pbinom(floor((N-ties+1)/2)-1, N, prob=0.5))
p.Fong <- num.Fong/denom.Fong}
if (dir== -1) p.Fong <- 1-p.Fong/2+dbinom(min(x.gt.y, x.lt.y), N-ties, prob=0.5)
#{num.Fong <- pbinom((max(x.gt.y, x.lt.y)), N, prob=0.5)
# denom.Fong <- 2*(1-pbinom(floor((N-ties+1)/2)-1, N, prob=0.5))
# p.Fong <- num.Fong/denom.Fong}
if(printstat==TRUE){
txt <- paste("Censored sign test for median(x:", xname, " - ", "y:", yname, ") equals 0", sep = "")
txt2 <- paste(" n =", N, " n+ =", x.gt.y, " n- =", x.lt.y, " ties:", ties, "\n")
cat(txt, "\n",txt3, "\n", txt2, "\n", " No correction for ties:", " p-value =", signif(s.binom$p.value, 4), "\n")
if (ties != 0) cat("Fong correction for ties:", " p-value =", signif(p.Fong, 4))
}
x=list(n=N,nplus=x.gt.y,nminu=x.lt.y,ties=ties,
p.value=signif(s.binom$p.value, 4),
p.value.corr=if(ties!=0){signif(p.Fong,4)}else{NA})
invisible(x)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.