Functie voor het berekenen van het Minimale Significante Verschil (in R)
In het kader wordt een functie gegeven om het Minimaal Significante Verschil (Minimum Significant Difference, MSD) te bepalen tussen een controle meting en de meting in een behandeling als de non-parametrische Wilcoxon's toets (=Mann-Whitney-U-toets) gebruikt is als significantie test. De functie is geschreven in R, een vrij beschikbare statistische software omgeving.
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
}
Dit programma wordt door ECOSTAT vrij beschikbaar gesteld. ECOSTAT wijst nadrukkelijk elke vorm van garantie voor de software af. De software wordt ter beschikking gesteld 'zoals deze is' zonder enige waarborg. Als dit programma gebruikt wordt als onderdeel bij het ontwikkelen van een ander software product, wordt verwacht dat daarbij de bijdrage van ECOSTAT en van de auteur van dit programma, N. van der Hoeven, erkend wordt.
Sluiten software informatie
Gratis software en computertips van ECOSTAT
Lijst van door ECOSTAT ontwikkelde vrij beschikbare software en computertips van ECOSTAT
| |
