Bug in MUS::MUS.calc.n.conservative (Result is wrong)

DESCRIPTION

@@ -3,14 +3,17 @@ Package: MUS
 Type: Package
 Title: Monetary Unit Sampling and Estimation Methods, Widely Used in
         Auditing
-Version: 0.1.4
-Date: 2017-05-29
-Author: Henning Prömpers
+Version: 0.1.6
+Date: 2019-09-15
+Author: Henning Prömpers, André Guimarães
 Maintainer: Henning Prömpers <henning@proempers.net>
 Description: Sampling and evaluation methods to apply Monetary Unit Sampling (or in older literature Dollar Unit Sampling) during an audit of financial statements.
-Depends: R (>= 3.3.2), stats
+Depends: R (>= 3.4.0), stats
+Suggests: DescTools, pander
 License: GPL (>= 2)
 NeedsCompilation: no
-Packaged: 2017-05-29 11:35:09 UTC; henning
+BugReports: https://github.com/alsguimaraes/MUS
+Packaged: 2019-09-15 11:06:00 UTC; henning
 Repository: CRAN
-Date/Publication: 2017-05-29 12:59:00 UTC
+Date/Publication: 2019-09-15 11:06:00 UTC
+RoxygenNote: 6.1.1

MD5

@@ -1,12 +1,38 @@
-92ef6028f21f3e52cfa959f4b175d688 *DESCRIPTION
-ebc82eeabe85c8e859c9188db3f4e2b8 *NAMESPACE
-28f87c6f5f4acfa082b367640e383b15 *R/MUS.evaluation.R
-131030fc311f00ca7d307d0ab511cf7d *R/MUS.extraction.R
-1d9a0a01a4994d8b473b3f76f4bf9bf4 *R/MUS.planning.R
-b06b9215961c94e2b5c5d23133f578d1 *R/print.MUS.evaluation.result.R
-466cdf8057b5044796def504e2ae97b5 *build/partial.rdb
-8d66b1de524e389d4a5c400f05590d6f *man/MUS-package.Rd
-80a60a86dd5f76a1d99a4739ccdd3bc7 *man/MUS.evaluation.Rd
-8c60a29255dfa4a933cb81811bc64113 *man/MUS.extraction.Rd
-88f5baa485a04c1a4cc3006284fb44dc *man/MUS.planning.Rd
-e8422b46140bf162d6ca7c2c62b18502 *man/print.MUS.evaluation.result.Rd
+a5cc879bbbde2779963efd937f75589b *build/partial.rdb
+5f0bc17891b9e08231936ee3590e79ac *DESCRIPTION
+5f9db751af196e6747c5fd418b66583b *examples/data.csv
+da533b2e22782e7c14d31387afb1bd4d *examples/example.R
+bd1409a941923f8e6e9a83db107ea854 *examples/example.Rmd
+49c905092eb5ae1fc00206fcdf7ac445 *examples/logo.png
+2553531635d707000aa31253b5f274ef *examples/workflow.R
+32631548c047c541ee897b741ad2953e *inst/po/en@quot/LC_MESSAGES/R-MUS.mo
+31f2543e08e9774dfef16d32bf340954 *inst/po/pt_BR/LC_MESSAGES/R-MUS.mo
+b8d9689af48ca28fc3c4c470e293c1c6 *man/MUS-package.Rd
+56682d1d63ca32060c7463c69b094219 *man/MUS.binomial.bound.Rd
+f01f9da0bce2624bff8da2d4696745d4 *man/MUS.calc.n.conservative.Rd
+eeef03ed9e4cb2fbb6cec4bfeea7a3a1 *man/MUS.combine.Rd
+f7d247df7d06a60fda55dd180738419e *man/MUS.combined.high.error.rate.Rd
+be4052a145e687017fca8029da1c7ab4 *man/MUS.evaluation.Rd
+2330e71495709ce7a90a2d5461fb8851 *man/MUS.extend.Rd
+7f2320966a4f6827ccaa61dbf8139ed8 *man/MUS.extraction.Rd
+292c6ed68ffd431f5c12c62be370ab3b *man/MUS.factor.Rd
+3c08c4957450564943a3efc8ff29d955 *man/MUS.moment.bound.Rd
+3d79b9c07c9bc7fef7a0a040b0965ad9 *man/MUS.multinomial.bound.Rd
+0daddbe026e7bfa835d7fea40ee69343 *man/MUS.planning.Rd
+c7e85c735e646069ec888201aaf083b1 *man/print.MUS.evaluation.result.Rd
+4e1716e1ea1a4f37d589ec334a6b543b *man/print.MUS.extraction.result.Rd
+d0f925ae487b32a55e2487062eef2fb7 *man/print.MUS.planning.result.Rd
+1faebf3391ec44409007b4658382bfe9 *NAMESPACE
+8e6549ce43e86d56f79f71b83d51f1b8 *po/R-MUS.pot
+31f2543e08e9774dfef16d32bf340954 *po/R-pt_BR.mo
+a0de8307b4c2e8554329170c9c4c015f *po/R-pt_BR.po
+177c1abba09291d8b502a473533589b4 *R/MUS.bounds.R
+4108a5e86888565eb8245a56e9cc006d *R/MUS.combine.R
+b2a016cc212874e2af2cfce533b8415d *R/MUS.evaluation.R
+5b34a0c94786bb478ae075f297ae97f2 *R/MUS.extend.R
+6d3e017a312482d5a3ff2b04c5ba756c *R/MUS.extraction.R
+4502183f507ae8473a44e9bc5a36cd20 *R/MUS.planning.R
+11a67054d6fa4ad2eea4099c74404b71 *R/print.MUS.evaluation.result.R
+7758ec42e5ebb60844c92cd8fae4bde5 *R/print.MUS.extraction.result.R
+b1214e49ec92ade404dcc845373be121 *R/print.MUS.planning.result.R
+42daef4e6a76933a0ac4ebecf31fd3fc *README.md

NAMESPACE

@@ -1,3 +1,5 @@
 exportPattern("^[[:alpha:]]+")
-importFrom("stats", "phyper", "runif", "uniroot")
+importFrom("stats", "phyper", "runif", "uniroot", "aggregate", "qgamma", "qnorm", "qt", "sd")
 S3method(print, MUS.evaluation.result)
+S3method(print, MUS.planning.result)
+S3method(print, MUS.extraction.result)

R/MUS.bounds.R

@@ -0,0 +1,128 @@
+MUS.moment.bound <- function(x, confidence.level=0.95, as.pct=FALSE, include.high.values=TRUE) {
+    # Dworking & Grimlund, 1984
+    # data = c(rep(0, 96), -.16, .04, .18, .47)
+    if (!class(x)=="MUS.evaluation.result" && !is.vector(x)) stop("x has to be a vector or an object of type MUS.evaluation.result. Use function MUS.evaluate to create such an object.")
+
+    if (class(x)=="MUS.evaluation.result") {
+        data <- c(1 - x$filled.sample[,x$col.name.audit.values] / x$filled.sample[,x$col.name.book.values])
+        if (include.high.values && is.data.frame(x$filled.high.values)) {
+            data <- c(data, (1 - x$filled.high.values[,x$col.name.audit.values] / x$filled.high.values[,x$col.name.book.values]))
+        }
+        confidence.level <- x$confidence.level
+        mult <- ifelse(as.pct, 100, x$Results.Total$Net.most.likely.error[1])
+    } else {
+        data <- x
+        mult <- 0
+    }
+
+    taintings = data[data!=0]
+    n <- length(taintings)
+    N <- length(data)
+
+    hypothetical.tainting = 0.81*(1-0.667*tanh(10*mean(taintings)))*(1+0.667*tanh(n/10))
+    TN <- sapply(1:3, function(j) (hypothetical.tainting^j+sum(taintings^j))/(n+1))
+    RN <- c( ((1+n)/(2+N)), ((1+n)/(2+N)) * ((2+n)/(3+N)), ((1+n)/(2+N)) * ((2+n)/(3+N)) * ((3+n)/(4+N)) )
+    UN <- c(  RN[1]*TN[1], (RN[1]*TN[2] + (N-1)*RN[2]*TN[1]^2)/N, ( RN[1]*TN[3] + 3*(N-1)*RN[2]*TN[1]*TN[2] + (N-1)*(N-2)*RN[3]*TN[1]^3)/N^2 )
+    UC <- c( 0, UN[2]-UN[1]^2, UN[3]-3*UN[1]*UN[2]+2*UN[1]^3 )
+
+    A <- (4*UC[2]^3/UC[3]^2)
+    B <- (0.5*UC[3]/UC[2])
+    G <- (UN[1]-2*UC[2]^2/UC[3])
+    Z <- qnorm(confidence.level)
+    CB <- G + A*B*(1+Z/sqrt(9*A)-1/(9*A))^3
+    ifelse(mult == 0, CB*100, (1 + CB) * (mult))
+}
+
+MUS.binomial.bound <- function(x, scope="qty", as.pct=FALSE, include.high.values=TRUE, confidence.level=0.95) {
+    if (!class(x)=="MUS.evaluation.result" && !is.vector(x)) stop("x has to be a vector or an object of type MUS.evaluation.result. Use function MUS.evaluate to create such an object.")
+    if (class(x)=="MUS.evaluation.result") {
+        if (scope == "value") {
+            if (include.high.values && is.data.frame(x$filled.high.values) ) {
+                misstatement <- sum(x$filled.sample[,x$col.name.book.values] - x$filled.sample[,x$col.name.audit.values]) +
+                    sum(x$filled.high.values[,x$col.name.book.values] - x$filled.high.values[,x$col.name.audit.values])
+                audited.value <- sum(x$filled.sample[,x$col.name.book.values]) + sum(x$filled.high.values[,x$col.name.book.values])
+                book.value <- length(x$filled.sample[,x$col.name.book.values]) + length(x$filled.high.values[,x$col.name.book.values])
+            } else {
+                misstatement <- sum(x$filled.sample[,x$col.name.book.values] - x$filled.sample[,x$col.name.audit.values])
+                audited.value <- sum(x$filled.sample[,x$col.name.book.values])
+                book.value <- length(x$filled.sample[,x$col.name.book.values])
+            }
+            misstatement <- ceiling(misstatement / audited.value * book.value )
+        } else {
+            if (include.high.values && is.data.frame(x$filled.high.values) ) {
+                misstatement <- sum(x$filled.sample[,x$col.name.book.values] != x$filled.sample[,x$col.name.audit.values]) +
+                    sum(x$filled.high.values[,x$col.name.book.values] != x$filled.high.values[,x$col.name.audit.values])
+                book.value <- length(x$filled.sample[,x$col.name.book.values]) + length(x$filled.high.values[,x$col.name.book.values])
+            } else {
+                misstatement <- sum(x$filled.sample[,x$col.name.book.values] != x$filled.sample[,x$col.name.audit.values])
+                book.value <- length(x$filled.sample[,x$col.name.book.values])
+            }
+        }
+        materiality <- x$tolerable.error / x$book.value
+        confidence.level = x$confidence.level
+        mult <- ifelse(as.pct, 100, x$book.value)
+    } else {
+        # x are taintings
+        mult <- 100
+        book.value <- length(x)
+        misstatement <- sum(x)
+    }
+
+    if (requireNamespace("DescTools", quietly = TRUE)) {
+        bc <- DescTools::BinomCI(misstatement, book.value, conf.level=1-(1-confidence.level) * 2, method = "clopper-pearson")
+        ifelse(as.pct, (bc[3] * mult), round(bc[3] * mult))
+    } else {
+        message("MUS.binomial.bound requires package DescTools.")
+	NULL
+    }
+}
+
+MUS.multinomial.bound <- function(x, as.pct=FALSE, include.high.values=TRUE) {
+	if (class(x)!="MUS.evaluation.result") stop("x has to be an object from type MUS.evaluation.result. Use function MUS.evaluate to create such an object.")
+    res <- NA
+    if (requireNamespace("DescTools", quietly = TRUE)) {
+        if (include.high.values && is.data.frame(x$filled.high.values) ) {
+            misstatement <- ceiling(c(1 - x$filled.sample[,x$col.name.audit.values] / x$filled.sample[,x$col.name.book.values],
+                1 - x$filled.high.values[,x$col.name.audit.values] / x$filled.high.values[,x$col.name.book.values]))
+        } else {
+            misstatement <- ceiling(c(1 - x$filled.sample[,x$col.name.audit.values] / x$filled.sample[,x$col.name.book.values]))
+        }
+
+        observed <- aggregate(data.frame(count = misstatement), list(value = misstatement), length)
+        res <- DescTools::MultinomCI(observed$count, conf.level=1-(1-x$confidence.level) * 2, method="sisonglaz")
+        mult <- ifelse(as.pct, 100, x$book.value)
+        ifelse(as.pct, ((1-res[observed$value==0][2])) * mult, round(((1-res[observed$value==0][2])) * mult))
+    } else {
+        message("MUS.multinomial.bound requires package DescTools.")
+	NULL
+    }
+}
+
+MUS.combined.high.error.rate <- function(evaluation, interval.type="one-sided"){
+	filled.sample <- evaluation$filled.sample
+	filled.high.values <- evaluation$filled.high.values
+	col.name.audit.values <- evaluation$col.name.audit.values
+	col.name.book.values <- evaluation$col.name.book.values
+	confidence.level <- evaluation$confidence.level
+
+	ratios <- 1 - filled.sample[,col.name.audit.values]/filled.sample[,col.name.book.values]
+	qty_errors <- sum(ratios!=0)
+	ratios_mean <- mean(ratios)
+	ratios_sd <- sd(ratios)
+    if (is.data.frame(filled.high.values)) {
+        N <- nrow(evaluation$data) - nrow(filled.high.values)
+        Y <- sum(evaluation$data[, col.name.book.values]) - sum(filled.high.values[, col.name.book.values])
+        high.values.error <- sum(filled.high.values[, col.name.book.values]-filled.high.values[, col.name.audit.values])
+    } else {
+        N <- nrow(evaluation$data)
+        Y <- sum(evaluation$data[, col.name.book.values])
+        high.values.error <- 0
+    }
+	R <- ifelse(interval.type == "two-sided",  1 - (1 - confidence.level) / 2, confidence.level)
+	U <- qt(R, qty_errors - 1)
+
+	most.likely.error <- ratios_mean * Y
+	precision <- U * Y * ratios_sd / sqrt(nrow(filled.sample))
+	upper.error.limit <- most.likely.error + precision * sign(most.likely.error) + high.values.error
+	upper.error.limit
+}

R/MUS.combine.R

@@ -0,0 +1,131 @@
+MUS.combine <- function(object.list) {
+    if (!is.list(object.list) || length(object.list)<1) {
+        stop("object.list must be a list with one or more MUS.evaluation.result objects.")
+    }
+    s <- 1
+    x <- object.list[[s]]
+    if (length(object.list)>1) {
+        x$High.value.threshold <- "-"
+        x$Strata <- length(object.list)
+        x$qty.rejected <- ifelse(x$acceptable, 0, 1)
+        x$qty.accepted <- ifelse(x$acceptable, 1, 0)
+        for (s in 2:length(object.list)) {
+            y <- object.list[[s]]
+            x$data <- rbind(x$data, y$data)
+            x$sample <- rbind(x$sample, y$sample)
+            if ("filled.sample" %in% names(x)) {
+                x$filled.sample <- rbind(x$filled.sample, y$filled.sample)
+            }
+            if ("high.values" %in% names(x) && class(y$high.values)=="data.frame") {
+                if (class(x$high.values)=="data.frame") {
+                    x$high.values <- rbind(x$high.values, y$high.values)
+                } else {
+                    x$high.values <- y$high.values
+                }
+            }
+            if ("filled.high.values" %in% names(x) && class(y$filled.high.values)=="data.frame") {
+                if (class(x$filled.high.values)=="data.frame") {
+                    x$filled.high.values <- rbind(x$filled.high.values, y$filled.high.values)
+                } else {
+                    x$filled.high.values <- y$filled.high.values
+                }
+            }
+            if ("sample.population" %in% names(x)) {
+                x$sample.population <- rbind(x$sample.population, y$sample.population)
+            }
+            if ("Results.Total" %in% names(x)) {
+                for (j in 1:2) {
+                    x$Results.Total$Number.of.Errors[j] <-
+                        x$Results.Total$Number.of.Errors[j] +
+                        y$Results.Total$Number.of.Errors[j]
+
+                    x$Results.Total$Net.upper.error.limit[j] <-
+                        x$Results.Total$Net.upper.error.limit[j] +
+                        y$Results.Total$Net.upper.error.limit[j]
+
+                    x$Results.Total$Gross.upper.error.limit[j] <-
+                        x$Results.Total$Gross.upper.error.limit[j] +
+                        y$Results.Total$Gross.upper.error.limit[j]
+
+                    x$Results.Sample$Number.of.Errors[j] <-
+                        x$Results.Sample$Number.of.Errors[j] +
+                        y$Results.Sample$Number.of.Errors[j]
+
+                    x$Results.Sample$Net.upper.error.limit[j] <-
+                        x$Results.Sample$Net.upper.error.limit[j] +
+                        y$Results.Sample$Net.upper.error.limit[j]
+
+                    x$Results.Sample$Gross.upper.error.limit[j] <-
+                        x$Results.Sample$Gross.upper.error.limit[j] +
+                        y$Results.Sample$Gross.upper.error.limit[j]
+
+                    x$Results.High.values$Number.of.Errors[j] <-
+                        x$Results.High.values$Number.of.Errors[j] +
+                        y$Results.High.values$Number.of.Errors[j]
+
+                    x$Results.High.values$Gross.Value.of.Errors[j] <-
+                        x$Results.High.values$Gross.Value.of.Errors[j] +
+                        y$Results.High.values$Gross.Value.of.Errors[j]
+
+                    x$Results.Sample$Precision.Gap.widening[j] <-
+                        x$Results.Sample$Precision.Gap.widening[j] +
+                        y$Results.Sample$Precision.Gap.widening[j]
+
+                    x$Results.Sample$Total.Precision[j] <-
+                        x$Results.Sample$Total.Precision[j] +
+                        y$Results.Sample$Total.Precision[j]
+
+                    x$Results.Total$Net.most.likely.error[j] <-
+                        x$Results.Total$Net.most.likely.error[j] +
+                        y$Results.Total$Net.most.likely.error[j]
+                }
+
+                x$Results.Total$Total.number.of.items.examined <-
+                    x$Results.Total$Total.number.of.items.examined +
+                    y$Results.Total$Total.number.of.items.examined
+
+                x$Results.Sample$Sample.Size <-
+                    x$Results.Sample$Sample.Size +
+                    y$Results.Sample$Sample.Size
+
+                x$Results.Sample$Basic.Precision <-
+                    x$Results.Sample$Basic.Precision +
+                    y$Results.Sample$Basic.Precision
+
+                x$Results.High.values$Number.of.high.value.items <-
+                    x$Results.High.values$Number.of.high.value.items +
+                    y$Results.High.values$Number.of.high.value.items
+
+                x$Results.High.values$Net.Value.of.Errors <-
+                    x$Results.High.values$Net.Value.of.Errors +
+                    y$Results.High.values$Net.Value.of.Errors
+
+                x$UEL.low.error.rate <- x$UEL.low.error.rate + y$UEL.low.error.rate
+                x$UEL.high.error.rate <- x$UEL.high.error.rate + y$UEL.high.error.rate
+
+                if ("acceptable.low.error.rate" %in% names(x)) {
+                    x$acceptable.low.error.rate <-
+                        ifelse(y$acceptable.low.error.rate, x$acceptable.low.error.rate, y$acceptable.low.error.rate)
+                }
+                if ("acceptable.high.error.rate" %in% names(x)) {
+                    x$acceptable.high.error.rate <-
+                        ifelse(y$acceptable.high.error.rate, x$acceptable.high.error.rate, y$acceptable.high.error.rate)
+                }
+                if ("acceptable" %in% names(x)) {
+                    x$acceptable <- ifelse(y$acceptable, x$acceptable, y$acceptable)
+                    x$qty.rejected <- x$qty.rejected + ifelse(y$acceptable, 0, 1)
+                    x$qty.accepted <- x$qty.accepted + ifelse(y$acceptable, 1, 0)
+                }
+
+            }
+
+            x$book.value <- x$book.value + y$book.value
+            x$tolerable.error <- x$tolerable.error + y$tolerable.error
+            x$expected.error <- x$expected.error + y$expected.error
+            x$n <- x$n + y$n
+            x$n.min <- x$n.min + y$n.min
+        }
+        x$combined <- TRUE
+    }
+    x
+}