CUSUM function in qcc package

I would like to know the method of calculation for the positive and negative cusum's which is used in the CUSUM function part of the qcc package.
By usuing the method for tabular CUSUM provided in the book "Introduction to statistical quality control 6th edition" by Douglas C.Montgomery which consistes of using the following formulas:


But was unable to obtain the same results as the ones provided by the cusum function.

When using the cusum (qcc)

cusum(data, center = 0.0122, std.dev = 0.006, decision.interval = 4, se.shift = 0.5)

 Samples

Group [,1]
1 0.030000000
2 0.014277869
3 0.011239112
4 0.010664081
5 0.008917331
6 0.010301526
7 0.009042467
8 0.010976948
9 0.009000000
10 0.010000000
11 0.010300000

cusumchart[["statistics"]]
1 2 3 4 5 6 7 8 9 10 11
0.030000000 0.014277869 0.011239112 0.010664081 0.008917331 0.010301526 0.009042467 0.010976948 0.009000000 0.010000000 0.010300000
cusumchart[["pos"]]
[1] 2.7166667 2.8129782 2.4028302 1.8968438 1.0997323 0.5333201 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
cusumchart[["neg"]]
[1] 0.000000000 0.000000000 0.000000000 -0.005986428 -0.303097912 -0.369510173 -0.645765623 -0.599607555 -0.882940888 -0.999607555 -1.066274222
'''

You can get the code for the function by typing cusum at the console prompt. Here is what I got.

> cusum
function (data, sizes, center, std.dev, head.start = 0, decision.interval = 5, 
    se.shift = 1, data.name, labels, newdata, newsizes, newlabels, 
    plot = TRUE, ...) 
{
    call <- match.call()
    if (missing(data)) 
        stop("'data' argument is not specified")
    if (missing(data.name)) 
        data.name <- deparse(substitute(data))
    data <- data.matrix(data)
    if (missing(sizes)) {
        sizes <- apply(data, 1, function(x) sum(!is.na(x)))
    }
    else {
        if (length(sizes) == 1) 
            sizes <- rep(sizes, nrow(data))
        else if (length(sizes) != nrow(data)) 
            stop("sizes length doesn't match with data")
    }
    if (decision.interval <= 0) 
        stop("decision.interval must be positive")
    if (head.start < 0 || head.start >= decision.interval) 
        stop("head.start must be non-negative and less than decision.interval")
    type <- if (any(sizes == 1)) 
        "xbar.one"
    else "xbar"
    if (ncol(data) == 1 & any(sizes > 1) & missing(std.dev)) 
        stop("sizes larger than 1 but data appears to be single samples. In this case you must provide also the std.dev")
    if (missing(labels)) {
        if (is.null(rownames(data))) 
            labels <- 1:nrow(data)
        else labels <- rownames(data)
    }
    stats <- paste("stats.", type, sep = "")
    if (!exists(stats, mode = "function")) 
        stop(paste("function", stats, "is not defined"))
    stats <- do.call(stats, list(data, sizes))
    statistics <- stats$statistics
    if (missing(center)) 
        center <- stats$center
    sd <- paste("sd.", type, sep = "")
    if (!exists(sd, mode = "function")) 
        stop(paste("function", sd, "is not defined!"))
    if (missing(std.dev)) {
        std.dev <- switch(type, xbar = {
            if (any(sizes > 25)) "RMSDF" else "UWAVE-R"
        }, NULL)
        std.dev <- do.call(sd, list(data, sizes, std.dev))
    }
    else {
        if (is.character(std.dev)) {
            std.dev <- do.call(sd, list(data, sizes, std.dev))
        }
        else {
            if (!is.numeric(std.dev)) 
                stop("if provided the argument 'std.dev' must be a method available or a numerical value. See help(qcc).")
        }
    }
    names(statistics) <- rownames(data) <- labels
    names(dimnames(data)) <- list("Group", "Samples")
    object <- list(call = call, type = "cusum", data.name = data.name, 
        data = data, statistics = statistics, sizes = sizes, 
        center = center, std.dev = std.dev)
    if (!missing(newdata)) {
        newdata.name <- deparse(substitute(newdata))
        newdata <- data.matrix(newdata)
        if (missing(newsizes)) {
            newsizes <- apply(newdata, 1, function(x) sum(!is.na(x)))
        }
        else {
            if (length(newsizes) == 1) 
                newsizes <- rep(newsizes, nrow(newdata))
            else if (length(newsizes) != nrow(newdata)) 
                stop("newsizes length doesn't match with newdata")
        }
        stats <- paste("stats.", type, sep = "")
        if (!exists(stats, mode = "function")) 
            stop(paste("function", stats, "is not defined"))
        newstats <- do.call(stats, list(newdata, newsizes))$statistics
        if (missing(newlabels)) {
            if (is.null(rownames(newdata))) {
                start <- length(statistics)
                newlabels <- seq(start + 1, start + length(newstats))
            }
            else {
                newlabels <- rownames(newdata)
            }
        }
        names(newstats) <- newlabels
        object$newstats <- newstats
        object$newdata <- newdata
        object$newsizes <- newsizes
        object$newdata.name <- newdata.name
        statistics <- c(statistics, newstats)
        sizes <- c(sizes, newsizes)
    }
    n <- length(statistics)
    z <- (statistics - center)/(std.dev/sqrt(sizes))
    ldb <- -decision.interval
    udb <- decision.interval
    z.f <- z - se.shift/2
    cusum.pos <- rep(NA, n)
    cusum.pos[1] <- max(0, head.start + z.f[1])
    for (i in 2:n) cusum.pos[i] <- max(0, cusum.pos[i - 1] + 
        z.f[i])
    z.f <- z + se.shift/2
    cusum.neg <- rep(NA, n)
    cusum.neg[1] <- max(0, head.start - z.f[1])
    for (i in 2:n) cusum.neg[i] <- max(0, cusum.neg[i - 1] - 
        z.f[i])
    cusum.neg <- -cusum.neg
    violations <- list(lower = which(cusum.neg < ldb), upper = which(cusum.pos > 
        udb))
    object$type <- "cusum"
    object$pos <- cusum.pos
    object$neg <- cusum.neg
    object$head.start <- head.start
    object$decision.interval <- decision.interval
    object$se.shift <- se.shift
    object$violations <- violations
    class(object) <- "cusum.qcc"
    if (plot) 
        plot(object, ...)
    return(object)
}
<bytecode: 0x000001caa0dfec48>
<environment: namespace:qcc>

This topic was automatically closed 42 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.