Skip to content

Wrong alpha level in the overall function #33

@Odd-Bird

Description

@Odd-Bird

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")
}

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions