SNStatComp/rtrim

Wrong alpha level in the overall function

Opened this issue · 0 comments

In the overall() function, alpha = c(0.05, 0.001); however, trend categories below contain p<0.01 level, not p<0.001.

function (x, which = c("imputed", "fitted"), changepoints = numeric(0), 
    bc = FALSE) 
{
    stopifnot(class(x) == "trim")
    which <- match.arg(which)
    if (is.character(changepoints) && changepoints == "model") {
        changepoints <- x$changepoints
    }
    if (all(changepoints %in% x$time.id)) 
        changepoints <- match(changepoints, x$time.id)
    tt_mod <- x$tt_mod
    tt_imp <- x$tt_imp
    var_tt_mod <- x$var_tt_mod
    var_tt_imp <- x$var_tt_imp
    J <- ntime <- x$ntime
    if (length(changepoints) == 0) 
        changepoints <- 1
    .meaning <- function(bhat, berr, df) {
        if (df <= 0) 
            return("Unknown (df<=0)")

        alpha = c(0.05, 0.001) # HERE!!!

        stopifnot(df > 0)
        if (bc) {
            tval <- qnorm((1 - alpha/2))
        }
        else {
            tval <- qt((1 - alpha/2), df)
        }
        blo <- bhat - tval * berr
        bhi <- bhat + tval * berr
        multiplicative <- TRUE
        if (multiplicative) {
            blo <- exp(blo)
            bhi <- exp(bhi)
            if (blo[2] > 1.05) 
                return("Strong increase (p<0.01)")
            if (bhi[2] < 0.95) 
                return("Strong decrease (p<0.01)")
            if (blo[1] > 1.05) 
                return("Strong increase (p<0.05)")
            if (bhi[1] < 0.95) 
                return("Strong decrease (p<0.05)")
            eps = 1e-07
            if (blo[2] > 1 + eps) 
                return("Moderate increase (p<0.01)")
            if (bhi[2] < 1 - eps) 
                return("Moderate decrease (p<0.01)")
            if (blo[1] > 1 + eps) 
                return("Moderate increase (p<0.05)")
            if (bhi[1] < 1 - eps) 
                return("Moderate decrease (p<0.05)")
            if (blo[1] > 0.95 && bhi[1] < 1.05) 
                return("Stable")
            return("Uncertain")
        }
        else {
            if (blo[2] > +0.05) 
                return("Strong increase (p<0.01)")
            if (bhi[2] < -0.05) 
                return("Strong decrease (p<0.01)")
            if (blo[1] > +0.05) 
                return("Strong increase (p<0.05)")
            if (bhi[1] < -0.05) 
                return("Strong decrease (p<0.05)")
            eps = 1e-07
            if (blo[2] > +eps) 
                return("Moderate increase (p<0.01)")
            if (bhi[2] < -eps) 
                return("Moderate decrease (p<0.01)")
            if (blo[1] > +eps) 
                return("Moderate increase (p<0.05)")
            if (bhi[1] < -eps) 
                return("Moderate decrease (p<0.05)")
            if (blo[1] > -0.05 && bhi[1] < 0.05) 
                return("Stable")
            return("Uncertain")
        }
    }
    .compute.overall.slope <- function(tpt, tt, var_tt, src) {
        n <- length(tpt)
        stopifnot(length(tt) == n)
        stopifnot(nrow(var_tt) == n && ncol(var_tt) == n)
        problem = tt < 1e-06
        log_tt = log(tt)
        log_tt[problem] <- 0
        alt_tt <- exp(log_tt)
        X <- cbind(1, tpt)
        y <- matrix(log_tt)
        bhat <- solve(t(X) %*% X) %*% t(X) %*% y
        yhat <- X %*% bhat
        dvtt <- 1/alt_tt
        Om <- diag(dvtt) %*% var_tt %*% diag(dvtt)
        var_beta <- solve(t(X) %*% X) %*% t(X) %*% Om %*% X %*% 
            solve(t(X) %*% X)
        b_err <- sqrt(diag(var_beta))
        df <- n - 2
        t_val <- bhat/b_err
        if (df > 0) 
            p <- 2 * pt(abs(t_val), df, lower.tail = FALSE)
        else p <- c(NA, NA)
        j <- 1:n
        D <- sum((j - mean(j))^2)
        SSR <- b_err[2]^2 * D * (n - 2)
        z <- data.frame(add = bhat, se_add = b_err, mul = exp(bhat), 
            se_mul = exp(bhat) * b_err, p = p, row.names = c("intercept", 
                "slope"))
        z$meaning = c("<none>", .meaning(z$add[2], z$se_add[2], 
            n - 2))
        list(src = src, coef = z, SSR = SSR)
    }
    if (which == "imputed") {
        tt <- tt_imp
        var_tt <- var_tt_imp
        src = "imputed"
    }
    else if (which == "fitted") {
        tt <- tt_mod
        var_tt <- var_tt_mod
        src = "fitted"
    }
    J = length(tt)
    if (length(changepoints) == 0) {
        out <- .compute.overall.slope(1:J, tt, var_tt, src)
        out$type <- "normal"
    }
    else {
        stopifnot(min(changepoints) >= 1)
        stopifnot(max(changepoints) < J)
        stopifnot(all(diff(changepoints) > 0))
        ncp <- length(changepoints)
        cpx <- c(changepoints, J)
        int.collector <- data.frame()
        slp.collector <- data.frame()
        SSR.collector <- numeric(ncp)
        for (i in 1:ncp) {
            idx <- cpx[i]:cpx[i + 1]
            tmp <- .compute.overall.slope(idx, tt[idx], var_tt[idx, 
                idx], src)
            prefix <- data.frame(from = x$time.id[cpx[i]], upto = x$time.id[cpx[i + 
                1]])
            intercept <- tmp$coef[1, ]
            int.collector <- rbind(int.collector, cbind(prefix, 
                intercept))
            slope <- tmp$coef[2, ]
            slp.collector <- rbind(slp.collector, cbind(prefix, 
                slope))
            SSR.collector[i] = tmp$SSR
        }
        out <- list(src = src, slope = slp.collector, intercept = int.collector, 
            SSR = SSR.collector)
    }
    out$J = J
    out$tt = tt
    out$err = sqrt(diag(var_tt))
    out$timept <- x$time.id
    structure(out, class = "trim.overall")
}