Function to calculate the Minimum Significant Difference (in R)
In the box, a function is given to calculate the Minimum Significant Difference (MSD) between a control and a treatment if the non-parametric Wilcoxon's test (=Mann-Whitney-U-test) is used as significance test. The function is written in R, a free statistical software environment.
Kopieer de tekst en plak hem in een R-werkruimte.
MSD.wilcox <- function(x, y, alpha = 0.05, exp.eff="decrease", shift="lin")
{
#######################################################################
#
#
# Function to calculate the Minimum Significant Difference #
# if the significance is tested with the one-sided Wilcoxon test #
# (= the Mann-Whitney-U test)
#
# Method described in N. van der Hoeven, 2008
#
#
Ecotox. Environm. Safety 70: p 61-66 #
#
#
# Author: N. van der Hoeven, ECOSTAT
#
# x: Vector with control data
#
# y: vector with data at NOEC
#
# alpha: significance level of the test, default is 0.05
#
# exp.eff: the direction of the expected effect.
#
# If a decrease is expected, string starting with a d (or D) #
# If an increase is expected, string starting with an i (or I) #
# Default: decreasing
#
# shift: The alternative is either a linear shift or a
#
# proportional shift.
#
# Linear shift: string starting with an l (or L) #
# Proportional shift: string staring with a p (or P) #
# Default: a linear shift
#
#
#
#######################################################################
n <- length(x)
m <- length(y)
nm.u <- length(xy.u <- unique(xy <- c(x,y)))
exp.eff <- substr(exp.eff,1,1) ; shift <- substr(shift,1,1)
if ((exp.eff== "d" ) | (exp.eff == "D" )) {exp.eff <- T}
else {if ((exp.eff== "i" ) | (exp.eff == "I" )) {exp.eff <- F }
else {cat( "Incorrect specification for the direction of the alternative.\n")
cat( "As direction of the alternative, a decrease is assumed.\n\n" )
exp.eff <- T
}
}
if ((shift== "l" ) | (shift == "L" )) {shift <- T}
else {if ((shift== "p" ) | (shift == "P" )) {shift <- F}
else {cat( "Incorrect specification for linear of proportional shift.\n")
cat( "A linear shift is assumed.\n\n" )
shift <- T
}
}
med.x <- median(x) ; med.y <- median(y)
if (shift)
{x1 <- x - med.x ; x1 <- matrix(x1, nrow=n, ncol=m)
y1 <- y - med.y ; y1 <- t(matrix(y1, nrow=m, ncol=n))
v <- as.vector(x1-y1)
}
else
{x1 <- x / med.x ; x1 <- matrix(x1, nrow=n, ncol=m)
y1 <- y / med.y ; y1 <- t(matrix(y1, nrow=m, ncol=n))
v <- as.vector(x1/y1)
}
if (exp.eff)
{ v <- v[order(v, decreasing=FALSE)]
alt <- "g"
meth.e <- "treatment is expected to decrease variable"
}
else
{ v <- v[order(v, decreasing=TRUE)]
alt <- "l"
meth.e <- "treatment is expected to increase variable"
}
if ((n<=50)&(m<=50))
{
# Exact calculation of MSD, not corrected for ties
Qe <- qwilcox(alpha*1.00001,n,m)-1
et <- TRUE
}
else
{
# To many observations for exact calculation of MSD
et <- FALSE
Qe <- NA
}
# Normal approximation of MSD
za <- qnorm(1-alpha)
Qa <- floor((m*n/2)-za*(m*n*(m+n+1)/12)^0.5)
tn <- T <- 1
if (nm.u == (n+m))
{# No ties
tie <- FALSE
Qa.t <- Qa
}
else
{
# Calculation of rang number Q corrected for ties
tie <- TRUE
tn <- summary(as.factor(xy))
T <- sum((tn^3)-tn)
Qa.t <- (m*n/2)-za*((m*n)*(((n+m)^3)-(n+m+T))/((m+n)*(n+m-1)*12))^0.5
Qa.t <- floor(Qa.t)
}
p.val <- wilcox.test(x,y,alternative=alt)
uu <- rep(NA,3)
if (et & (Qe>0)) {uu[1] <- 0.5*(v[Qe]+v[Qe+1]) }
if (Qa>0) {uu[2] <- 0.5*(v[Qa]+v[Qa+1]) }
if ((Qa.t>0)&tie) {uu[3] <- 0.5*(v[Qa.t]+v[Qa.t+1]) }
naam.uu <- c( "exact" , "normal appr." , "tie corr. appr." )
attr(uu, "names" )<- naam.uu
if (shift)
{wil <- (med.y-med.x)
med.MSD <- med.x + uu
meth.s <- "linear shift"
}
else
{wil <- (med.y/med.x)
med.MSD <- med.x * uu
meth.s <- "proportional shift"
}
wil <- c(wil,p.val[[3]][1])
naam.wil <- c( "obs. diff." , "p-value" )
attr(wil, "names" )<- naam.wil
med <- c(med.x,med.y)
naam.med <- c( "control median" , "treatment median" )
attr(med, "names" )<- naam.med
meth <- c(meth.e,meth.s)
out<- list(meth,uu,med,wil,med.MSD)
naam.out <- c( "Alternative and shift tested" , "MSD" )
naam.out <- c(naam.out, "observed medians" , "observed differences" )
naam.out <- c(naam.out, "median at MSD" )
attr(out, "names" )<- naam.out
out
}
This software is provided by ECOSTAT for full, free and open release. It is understood by the recipient/user that ECOSTAT assumes no liability for any errors contained in the code. Although this software is released without conditions or restrictions in its use, it is expected that appropriate credit be given to its author, N. van der Hoeven, and ECOSTAT should the software be included by the recipient as an element in other product development.
Close software information
Free software and computer tips by ECOSTAT
List of the free available software developed by ECOSTAT and some computer tips from ECOSTAT
| |
