Analyte options
For analytes under inspection (i.e. which are selected in Tab.C3), several analyte specific parameters can be set and will be stored also in the backup file.
Pooling
pooling means that calculations in the material table (mean, sd and uncertainty columns) are not based on the lab means, but rather on all measured values. This is justified when the between lab variance is smaller than the within lab variance (check if the ANOVA P-value is insignificant in Tab.C2, if so pooling might be allowed).
Selecting this option will also affect n in Tab.C3 being either the number of included labs or the number of finite data points from these labs.
Filtering
Individual Samples can be filtered by ID (select show IDs in Fig.C1 to identify outlying IDs). These samples will not be used in any calculations. They are removed from downstream processing right after import. However, information of this filtering step is preserved and included in reports.
Labs can be excluded by Lab-ID. Please note, that this is done after statistical testing for mean and variance outliers (in Tab.C1 and Tab.C2). Filtered labs are depicted in Fig.C1 (using grey color), but excluded from calculating the certified value of the analyte in Tab.C3.
Rounding
Note! Internally, all calculations are performed on non-rounded values (as imported). Only the visual output to the user is rounded to improve readability. Because analytes within an RM might be measured at very different scales it is possible to define analyte specific precision values for rounding. The Precision (Tables) value is used for most values displayed in tables.
The only exception is Tab.C3. Here, Precision (Cert. Val.) is used to round \(\mu_c\) and \(U_{abs}\) to allow to comply with DIN/ISO regulations on rounding and columns containing relative uncertainties are rounded to fixed 4 digit precision.
The user can decide the rounding precision for \(\mu_c\) and \(U_{abs}\) and even specify negative numbers for Precision (Cert. Val.) to indicate a desired rounding digit left of the decimal separator (i.e. \(-1\) to round an uncertainty of 13.1 up to 20). If the uncertainty is even larger, the user should consider to upload data using a different (more appropriate) unit and, for instance, switch from mg/L to g/L.
DIN1333
eCerto suggests the rounding precision for \(\mu_c\) and \(U_{abs}\) according to DIN 1333. The appropriate value is displayed right over the input box of Precision (Cert. Val.). It's background is colored green if it is in accordance with the user selected value and red otherwise.
The significant digit is determined based on the first non-zero digit in \(U_{abs}\). The position will be shown above the analyte parameter Precision (Cert. Val.) upon selection of an analyte in Tab.C3. Once found, the respective digit will be rounded up for \(\mu_c\) if the consecutive digit is \(\ge 5\) and down otherwise. For \(U_{abs}\) the rounding procedure is a bit more complicated. If the respective digit is 1 or 2, than the consecutive digit will be rounded up. If the respective digit is \(\ge 3\), than it will be rounded up itself.
Laboratory statistics to identify potential outliers for lab means and lab variances
For an analyte selected in Tab.C3 various statistical tests regarding lab means and variance will be performed according to ISO 17025 and outlying values will be indicated if observed at the .05 and .01 level respectively. The displayed values n.s. and . indicate either that the test probability was >0.05 or that this lab was not tested because it did not contain an extreme value.
Note! eCerto performs outlier tests at the lower and upper end of the value distribution independently (i.e. two one sided tests). This might lead to confusion as P-values for both independent tests are displayed within the same column of Tab.C1. However, it allows to keep the output compact. Please always use results of Tab.C1 in parallel with Fig.C1 to decide on the removal of Labs.
Except for Cochran, which tests for outliers with respect to variance, all other columns indicate potential outliers regarding lab means.
To compute the statistical tests eCerto uses functions from different packages available for R. Details regarding the conducted statistical test (implementation, parameters) can be found using the following links:
- Scheffe Note! Similar letters indicate similar group and hence no difference. The test is performed at two \(\alpha\) levels (.01 and .05) and re-implemented following the code provided in the
agricolae package.
- Dixon Note! The Dixon-Test was originally developed for small sample sizes of \(3 \le n \le 30\). We use tabulated critical values for Dixon up to \(n=100\) obtained from different tables in this website. Different tables are used according to \(n\).
- Grubbs Note! Test is performed as single (only the most extreme value at each end of the distribution is tested) as well as double version (the two most extreme values are tested against the rest). In accordance with ISO 17025 the double Grubbs Test is only conducted in case that the single Grubbs Test is not significant. Critical values for double Grubbs are obtained from the provided link for \(4 \le n \le 30\) and are estimated for \(31 \le n \le 100\) based on this paper.
- Cochran Note! The Cochran test is performed consecutively, i.e. the Lab showing the highest variance is tested first. If significant, it is removed and the procedure is repeated.
excl indicates that the sd of a lab was too low and the lab was removed from the testing procedure.
Note! Most tests require at least 3 data points (Labs) and finite differences in Lab means and Lab variances. If these conditions are not fulfilled Error might be reported instead of a P-value.
Instead of the qualitative significance levels reached in the test, the options panel to the right of Tab.C1 offers the possibility to display the calculated P-values directly for most of the tests. Further, the calculated test statistic values can be displayed. The user can also select to show the critical values at \(\alpha=0.05\) and \(\alpha=0.01\). Depending on the test these are calculated based on the underlying distributions or obtained from tabulated values.
If any lab is excluded after using the filter section of the analyte panel on top of the page, all tests in Tab.C1 can be conducted without filtered labs selecting the option Exclude filtered Labs. This feature allows to confirm whether or not further outliers are present in the data set. Please use the possibility to remove outliers sequentially with caution, i.e. refrain from over exclusion of initially masked outliers.
Statistical test details
Tests used for outlier detection often rely on tabulated critical values. Decisions to reject or not reject a potential outlier can be dependent on the source of the critical values.
Where ever possible, we tried to omit using fixed tables and rather implemented critical value calculation ab initio. However, to be transparent, we detail the mathematical approach for each test in the following.
Note! Several test functions were implemented based on the functions in the R package outliers by Lukasz Komsta. Modifications were made to allow testing up to n=100 laboratories (was restricted to n=30 in package outliers).
Cochran
The critical value is calculated internally using formula eCerto:::qcochran
## function(p, n, k) {
## f <- stats::qf(p / k, (n - 1) * (k - 1), n - 1)
## c <- 1 / (1 + (k - 1) * f)
## return(c)
## }
## <bytecode: 0x55f7561b7d38>
## <environment: namespace:eCerto>
where parameters p, n and k define the desired alpha level, the number of replicates per group and the number of groups respectively. stats::qf() is the quantile function of the f-distribution as implemented in the R stats package. Likewise, stats::pf() is used to compute the Cochran P-value via function eCerto:::pcochran
## function(q, n, k) {
## f <- (1 / q - 1) / (k - 1)
## p <- 1 - stats::pf(f, (n - 1) * (k - 1), n - 1) * k
## p[p < 0] <- 0
## p[p > 1] <- 1
## return(p)
## }
## <bytecode: 0x55f755e87278>
## <environment: namespace:eCerto>
For comparison with critical values tabulated elsewhere you can generate a respective table in R. Show code
ps <- c(0.01, 0.05)
ns <- c(3:9, seq(10,20,5))
ks <- c(2:14, seq(15,30,5))
out <- lapply(ps, function(p) {
sapply(ns, function(n) {
sapply(ks, function(k) {
eCerto:::qcochran(p=p, n=n, k=k)
})
})
})
names(out) <- paste0("alpha=", ps)
lapply(out, function(x) {
colnames(x) <- paste0("n=", ns)
rownames(x) <- paste0("k=", ks)
round(x, 4)
})
Show Cochran critical values table
## $`alpha=0.01`
## n=3 n=4 n=5 n=6 n=7 n=8 n=9 n=10 n=15 n=20
## k=2 0.9950 0.9794 0.9586 0.9373 0.9172 0.8988 0.8823 0.8674 0.8113 0.7744
## k=3 0.9423 0.8832 0.8335 0.7933 0.7606 0.7335 0.7107 0.6912 0.6241 0.5838
## k=4 0.8643 0.7814 0.7212 0.6761 0.6410 0.6129 0.5897 0.5702 0.5054 0.4678
## k=5 0.7885 0.6957 0.6329 0.5875 0.5531 0.5259 0.5038 0.4853 0.4251 0.3907
## k=6 0.7218 0.6258 0.5635 0.5195 0.4866 0.4609 0.4401 0.4229 0.3672 0.3360
## k=7 0.6644 0.5685 0.5080 0.4659 0.4347 0.4105 0.3911 0.3751 0.3237 0.2950
## k=8 0.6152 0.5210 0.4627 0.4227 0.3932 0.3705 0.3523 0.3373 0.2896 0.2632
## k=9 0.5727 0.4810 0.4251 0.3870 0.3592 0.3378 0.3207 0.3067 0.2622 0.2377
## k=10 0.5358 0.4469 0.3934 0.3572 0.3308 0.3106 0.2945 0.2814 0.2397 0.2169
## k=11 0.5036 0.4175 0.3663 0.3318 0.3068 0.2876 0.2725 0.2601 0.2209 0.1995
## k=12 0.4751 0.3919 0.3428 0.3099 0.2861 0.2680 0.2536 0.2419 0.2049 0.1848
## k=13 0.4498 0.3695 0.3223 0.2909 0.2682 0.2509 0.2373 0.2261 0.1911 0.1721
## k=14 0.4272 0.3495 0.3042 0.2741 0.2525 0.2360 0.2230 0.2124 0.1792 0.1612
## k=15 0.4069 0.3318 0.2882 0.2593 0.2386 0.2228 0.2104 0.2003 0.1687 0.1516
## k=20 0.3297 0.2654 0.2288 0.2048 0.1877 0.1748 0.1646 0.1564 0.1308 0.1170
## k=25 0.2782 0.2220 0.1904 0.1699 0.1553 0.1443 0.1357 0.1287 0.1071 0.0956
## k=30 0.2412 0.1914 0.1635 0.1455 0.1327 0.1231 0.1156 0.1096 0.0909 0.0809
##
## $`alpha=0.05`
## n=3 n=4 n=5 n=6 n=7 n=8 n=9 n=10 n=15 n=20
## k=2 0.9750 0.9392 0.9057 0.8772 0.8534 0.8332 0.8159 0.8010 0.7487 0.7164
## k=3 0.8709 0.7977 0.7457 0.7070 0.6770 0.6531 0.6333 0.6167 0.5613 0.5289
## k=4 0.7679 0.6839 0.6287 0.5894 0.5598 0.5365 0.5175 0.5018 0.4500 0.4205
## k=5 0.6838 0.5981 0.5440 0.5063 0.4783 0.4564 0.4387 0.4241 0.3767 0.3500
## k=6 0.6161 0.5321 0.4803 0.4447 0.4184 0.3980 0.3817 0.3682 0.3247 0.3004
## k=7 0.5612 0.4800 0.4307 0.3972 0.3726 0.3536 0.3384 0.3259 0.2858 0.2635
## k=8 0.5157 0.4377 0.3910 0.3594 0.3362 0.3185 0.3043 0.2927 0.2556 0.2350
## k=9 0.4775 0.4027 0.3584 0.3285 0.3067 0.2901 0.2768 0.2659 0.2313 0.2122
## k=10 0.4450 0.3733 0.3311 0.3028 0.2823 0.2666 0.2541 0.2439 0.2115 0.1936
## k=11 0.4169 0.3482 0.3080 0.2811 0.2616 0.2468 0.2350 0.2254 0.1949 0.1782
## k=12 0.3924 0.3264 0.2880 0.2624 0.2440 0.2299 0.2187 0.2096 0.1808 0.1650
## k=13 0.3709 0.3074 0.2707 0.2463 0.2286 0.2152 0.2046 0.1960 0.1687 0.1538
## k=14 0.3517 0.2907 0.2554 0.2321 0.2153 0.2025 0.1924 0.1841 0.1582 0.1440
## k=15 0.3346 0.2758 0.2419 0.2195 0.2034 0.1912 0.1815 0.1737 0.1489 0.1355
## k=20 0.2705 0.2205 0.1921 0.1735 0.1602 0.1502 0.1422 0.1358 0.1157 0.1048
## k=25 0.2281 0.1846 0.1601 0.1441 0.1328 0.1242 0.1174 0.1120 0.0949 0.0857
## k=30 0.1979 0.1593 0.1377 0.1236 0.1137 0.1061 0.1003 0.0955 0.0806 0.0727
The test is originally based on:
"Snedecor, G.W., Cochran, W.G. (1980). Statistical Methods (seventh edition). Iowa State University Press, Ames, Iowa".
Dixon
The Dixon test statistic is calculated depending on the number of labs (n) using the function eCerto:::dixon.test
## function(x, opposite = FALSE, two.sided = FALSE) {
## x <- sort(x[stats::complete.cases(x)])
## n <- length(x)
## if (n <= 7) {
## type <- "10"
## } else if (n > 7 & n <= 10) {
## type <- "11"
## } else if (n > 10 & n <= 13) {
## type <- "21"
## } else {
## type <- "22"
## }
## if (xor(((x[n] - mean(x)) < (mean(x) - x[1])), opposite)) {
## alt <- paste("lowest value", x[1], "is an outlier")
## q <- switch(type,
## "10" = (x[2] - x[1]) / (x[n] - x[1]),
## "11" = (x[2] - x[1]) / (x[n - 1] - x[1]),
## "21" = (x[3] - x[1]) / (x[n - 1] - x[1]),
## (x[3] - x[1]) / (x[n - 2] - x[1])
## )
## } else {
## alt <- paste("highest value", x[n], "is an outlier")
## q <- switch(type,
## "10" = (x[n] - x[n - 1]) / (x[n] - x[1]),
## "11" = (x[n] - x[n - 1]) / (x[n] - x[2]),
## "21" = (x[n] - x[n - 2]) / (x[n] - x[2]),
## (x[n] - x[n - 2]) / (x[n] - x[3])
## )
## }
## pval <- pdixon(q, n)
## if (two.sided) {
## pval <- 2 * pval
## if (pval > 1) pval <- 2 - pval
## }
## out <- list(
## "statistic" = c(Q = q),
## "alternative" = alt,
## "p.value" = pval,
## "method" = "Dixon test for outliers",
## "data.name" = "eCerto internal"
## )
## class(out) <- "htest"
## return(out)
## }
## <bytecode: 0x55f75297ea30>
## <environment: namespace:eCerto>
In eCerto the one-sided version of the Dixon test is applied consecutively at the lower and upper end. Critical values are stored in tabulated form internally (eCerto::cvals_Dixon) and have been combined in a single table from this source. Show Dixon critical values table
## 0.001 0.002 0.005 0.01 0.02 0.05 0.1 0.2
## 3 0.999 0.998 0.994 0.988 0.976 0.941 0.886 0.782
## 4 0.964 0.949 0.921 0.889 0.847 0.766 0.679 0.561
## 5 0.895 0.869 0.824 0.782 0.729 0.643 0.559 0.452
## 6 0.822 0.792 0.744 0.698 0.646 0.563 0.484 0.387
## 7 0.763 0.731 0.681 0.636 0.587 0.507 0.433 0.344
## 8 0.799 0.769 0.724 0.682 0.633 0.554 0.480 0.386
## 9 0.750 0.720 0.675 0.634 0.586 0.512 0.441 0.352
## 10 0.713 0.683 0.637 0.597 0.551 0.477 0.409 0.325
## 11 0.770 0.746 0.708 0.674 0.636 0.575 0.518 0.445
## 12 0.739 0.714 0.676 0.643 0.605 0.546 0.489 0.420
## 13 0.713 0.687 0.649 0.617 0.580 0.522 0.467 0.399
## 14 0.732 0.708 0.672 0.640 0.603 0.546 0.491 0.422
## 15 0.708 0.685 0.648 0.617 0.582 0.524 0.470 0.403
## 16 0.691 0.667 0.630 0.598 0.562 0.505 0.453 0.386
## 17 0.671 0.647 0.611 0.580 0.545 0.489 0.437 0.373
## 18 0.652 0.628 0.594 0.564 0.529 0.475 0.424 0.361
## 19 0.640 0.617 0.581 0.551 0.517 0.462 0.412 0.349
## 20 0.627 0.604 0.568 0.538 0.503 0.450 0.401 0.339
## 25 0.574 0.550 0.517 0.489 0.457 0.406 0.359 0.302
## 30 0.539 0.517 0.484 0.456 0.425 0.376 0.332 0.278
## 35 0.511 0.490 0.459 0.431 0.400 0.354 0.311 0.260
## 40 0.490 0.469 0.438 0.412 0.382 0.337 0.295 0.246
## 45 0.475 0.454 0.423 0.397 0.368 0.323 0.283 0.234
## 50 0.460 0.439 0.410 0.384 0.355 0.312 0.272 0.226
## 60 0.437 0.417 0.388 0.363 0.336 0.294 0.256 0.211
## 70 0.422 0.403 0.374 0.349 0.321 0.280 0.244 0.201
## 80 0.408 0.389 0.360 0.337 0.310 0.270 0.234 0.192
## 90 0.397 0.377 0.350 0.326 0.300 0.261 0.226 0.185
## 100 0.387 0.368 0.341 0.317 0.292 0.253 0.219 0.179
Values for non-tabulated combinations of p and n are derived by interpolation using function eCerto:::qtable
## function(p, probs, quants) {
## quants <- quants[order(probs)]
## probs <- sort(probs)
## res <- vector()
## for (n in 1:length(p)) {
## pp <- p[n]
## if (pp <= probs[1]) {
## q0 <- quants[c(1, 2)]
## p0 <- probs[c(1, 2)]
## fit <- stats::lm(q0 ~ p0)
## } else if (pp >= probs[length(probs)]) {
## q0 <- quants[c(length(quants) - 1, length(quants))]
## p0 <- probs[c(length(probs) - 1, length(probs))]
## fit <- stats::lm(q0 ~ p0)
## } else {
## x0 <- which(abs(pp - probs) == min(abs(pp - probs)))
## x1 <- which(abs(pp - probs) == sort(abs(pp - probs))[2])
## x <- min(c(x0, x1))
## if (x == 1) {
## x <- 2
## }
## if (x > length(probs) - 2) {
## x <- length(probs) - 2
## }
## i <- c(x - 1, x, x + 1, x + 2)
## q0 <- quants[i]
## p0 <- probs[i]
## fit <- stats::lm(q0 ~ poly(p0, 3))
## }
## res <- c(res, stats::predict(fit, newdata = list(p0 = pp)))
## }
## return(res)
## }
## <bytecode: 0x55f757838000>
## <environment: namespace:eCerto>
The test is originally based on:
"Dixon, W.J. (1950). Analysis of extreme values. Ann. Math. Stat. 21, 4, 488-506"
"Dean, R.B.; Dixon, W.J. (1951). Simplified statistics for small numbers of observations. Anal.Chem. 23, 636-638"
"Dixon, W.J. (1953). Processing data for outliers. J. Biometrics. 9, 74-89"
Grubbs
The Grubbs test statistic is calculated either for a single or a double outlier at both ends of the distribution consecutively using eCerto:::grubbs.test
## function(x, type = 10, tail = c("lower", "upper")) {
## tail <- match.arg(tail)
## if (!any(type == c(10, 20))) {
## message("'grubbs.test' is only implemented for type = 10 or 20. Using type = 20 (double Grubbs).")
## }
## x <- sort(x[stats::complete.cases(x)])
## n <- length(x)
## if (type == 10) {
## # single Grubbs
## if (tail == "lower") {
## alt <- paste("lowest value", x[1], "is an outlier")
## o <- x[1]
## d <- x[2:n]
## } else {
## alt <- paste("highest value", x[n], "is an outlier")
## o <- x[n]
## d <- x[1:(n - 1)]
## }
## g <- abs(o - mean(x)) / stats::sd(x)
## u <- stats::var(d) / stats::var(x) * (n - 2) / (n - 1)
## pval <- 1 - pgrubbs(g, n, type = 10)
## method <- "Grubbs test for one outlier"
## } else {
## # double Grubbs
## if (tail == "lower") {
## alt <- paste("lowest values", x[1], "and", x[2], "are outliers")
## u <- stats::var(x[3:n]) / stats::var(x) * (n - 3) / (n - 1)
## } else {
## alt <- paste("highest values", x[n - 1], "and", x[n], "are outliers")
## u <- stats::var(x[1:(n - 2)]) / stats::var(x) * (n - 3) / (n - 1)
## }
## g <- NULL
## pval <- pgrubbs(u, n, type = 20)
## method <- "Grubbs test for two outliers"
## }
## out <- list(
## "statistic" = c(G = g, U = u),
## "alternative" = alt,
## "p.value" = pval,
## "method" = method,
## "data.name" = "eCerto internal"
## )
## class(out) <- "htest"
## return(out)
## }
## <bytecode: 0x55f757a15120>
## <environment: namespace:eCerto>
Critical values for the single Grubbs test are calculated internally using eCerto:::qgrubbs
## function(p, n, type = 10, rev = FALSE) {
## if (type == 10) {
## if (!rev) {
## t2 <- stats::qt((1 - p) / n, n - 2)^2
## return(((n - 1) / sqrt(n)) * sqrt(t2 / (n - 2 + t2)))
## } else {
## s <- (p^2 * n * (2 - n)) / (p^2 * n - (n - 1)^2)
## t <- sqrt(s)
## if (is.nan(t)) {
## res <- 0
## } else {
## res <- n * (1 - stats::pt(t, n - 2))
## res[res > 1] <- 1
## }
## return(1 - res)
## }
## } else {
## if (n > 30) warning("[qgrubbs] critical value is estimated for n>30")
## gtwo <- eCerto::cvals_Grubbs2
## pp <- as.numeric(colnames(gtwo))
## if (!rev) res <- qtable(p, pp, gtwo[n - 3, ]) else res <- qtable(p, gtwo[n - 3, ], pp)
## res[res < 0] <- 0
## res[res > 1] <- 1
## return(unname(res))
## }
## }
## <bytecode: 0x55f757d4ca10>
## <environment: namespace:eCerto>
For comparison with critical values tabulated elsewhere you can generate a respective table in R. Show code
ps <- c(0.01, 0.025, 0.05, 0.1)
ns <- c(3:10, seq(20,100,40))
out <- sapply(ps, function(p) {
sapply(ns, function(n) {
eCerto:::qgrubbs(p=p, n=n)
})
})
colnames(out) <- paste0("a=", ps)
rownames(out) <- paste0("n=", ns)
round(out, 4)
Show Grubbs (single) critical values table
## a=0.01 a=0.025 a=0.05 a=0.1
## n=3 0.5878 0.6033 0.6289 0.6787
## n=4 0.7575 0.7688 0.7875 0.8250
## n=5 0.8863 0.8961 0.9123 0.9452
## n=6 0.9892 0.9981 1.0131 1.0435
## n=7 1.0744 1.0828 1.0969 1.1258
## n=8 1.1469 1.1549 1.1685 1.1962
## n=9 1.2098 1.2176 1.2307 1.2576
## n=10 1.2652 1.2728 1.2856 1.3118
## n=20 1.6118 1.6184 1.6297 1.6529
## n=60 2.0999 2.1057 2.1155 2.1358
## n=100 2.3040 2.3095 2.3188 2.3381
Critical values for the double Grubbs test are stored in tabulated form internally (eCerto::cvals_Grubbs2) and P-values are derived by comparison. Show Grubbs (double) critical values table
## 0.01 0.025 0.05 0.1
## 4 0.0000 0.0002 0.0008 0.0031
## 5 0.0035 0.0090 0.0183 0.0376
## 6 0.0186 0.0349 0.0565 0.0921
## 7 0.0440 0.0708 0.1020 0.1479
## 8 0.0750 0.1101 0.1478 0.1994
## 9 0.1082 0.1492 0.1909 0.2454
## 10 0.1415 0.1865 0.2305 0.2863
## 20 0.3909 0.4391 0.4804 0.5269
## 60 0.7089 0.7346 0.7553 0.7775
## 100 0.8018 0.8190 0.8326 0.8470
Note! Critical values for double Grubbs are obtained from the outliers package for 4≤n≤30 and are estimated for 31≤n≤100 based on the approximation described in this paper.
The test is originally based on:
"Grubbs, F.E. (1950). Sample Criteria for testing outlying observations. Ann. Math. Stat. 21, 1, 27-58".
Scheffe
The Scheffe multiple comparison test is a re-implementation of the respective function published in the agricolae package. It is calculated on the result of stats::lm(). eCerto:::scheffe.test
## function(y, trt, alpha = 0.05) {
## name.y <- paste(deparse(substitute(y)))
## name.t <- paste(deparse(substitute(trt)))
## A <- y$model
## DFerror <- stats::df.residual(y)
## MSerror <- stats::deviance(y) / DFerror
## y <- A[, 1]
## ipch <- pmatch(trt, names(A))
## nipch <- length(ipch)
## for (i in 1:nipch) {
## if (is.na(ipch[i])) {
## return(trt)
## }
## }
## name.t <- names(A)[ipch][1]
## trt <- A[, ipch]
## if (nipch > 1) {
## trt <- A[, ipch[1]]
## for (i in 2:nipch) {
## name.t <- paste(name.t, names(A)[ipch][i], sep = ":")
## trt <- paste(trt, A[, ipch[i]], sep = ":")
## }
## }
## name.y <- names(A)[1]
## df <- subset(data.frame("value" = y, "Lab" = trt), is.na(y) == FALSE)
## # JL
## means <- plyr::ldply(split(df$value, df$Lab), function(x) {
## data.frame(
## "mean" = mean(x, na.rm = T),
## "sd" = stats::sd(x, na.rm = T),
## "n" = sum(is.finite(x)),
## stringsAsFactors = FALSE
## )
## }, .id = "Lab")
## ntr <- nrow(means)
## Fprob <- stats::qf(1 - alpha, ntr - 1, DFerror)
## Tprob <- sqrt(Fprob * (ntr - 1))
## nr <- unique(means[, 4])
## scheffe <- Tprob * sqrt(2 * MSerror / nr)
## statistics <- data.frame(
## "MSerror" = MSerror,
## "Df" = DFerror,
## "F" = Fprob,
## "Mean" = mean(df[, 1]),
## "CV" = sqrt(MSerror) * 100 / mean(df[, 1])
## )
## if (length(nr) == 1) {
## statistics <- data.frame(statistics, "Scheffe" = Tprob, "CriticalDifference" = scheffe)
## }
## comb <- utils::combn(ntr, 2)
## nn <- ncol(comb)
## dif <- rep(0, nn)
## LCL <- dif
## UCL <- dif
## pval <- rep(0, nn)
## for (k in 1:nn) {
## i <- comb[1, k]
## j <- comb[2, k]
## dif[k] <- means[i, 2] - means[j, 2]
## sdtdif <- sqrt(MSerror * (1 / means[i, 4] + 1 / means[j, 4]))
## pval[k] <- round(1 - stats::pf(abs(dif[k])^2 / ((ntr - 1) * sdtdif^2), ntr - 1, DFerror), 4)
## LCL[k] <- dif[k] - Tprob * sdtdif
## UCL[k] <- dif[k] + Tprob * sdtdif
## }
## pmat <- matrix(1, ncol = ntr, nrow = ntr)
## k <- 0
## for (i in 1:(ntr - 1)) {
## for (j in (i + 1):ntr) {
## k <- k + 1
## pmat[i, j] <- pval[k]
## pmat[j, i] <- pval[k]
## }
## }
## groups <- orderPvalue(means, alpha, pmat)
## names(groups)[1] <- name.y
## parameters <- data.frame(test = "Scheffe", name.t = name.t, ntr = ntr, alpha = alpha)
## rownames(parameters) <- " "
## rownames(statistics) <- " "
## rownames(means) <- means[, 1]
## means <- means[, -1]
## output <- list(
## statistics = statistics,
## parameters = parameters,
## means = means,
## comparison = NULL,
## groups = groups
## )
## class(output) <- "group"
## return(output)
## }
## <bytecode: 0x55f7586215d8>
## <environment: namespace:eCerto>
The test is originally based on:
"Steel, R.; Torri, J.; Dickey, D. (1997) Principles and Procedures of Statistics: A Biometrical Approach. pp189"