Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Wrong alpha level in the overall function #33

Open
Odd-Bird opened this issue Jul 11, 2024 · 0 comments
Open

Wrong alpha level in the overall function #33

Odd-Bird opened this issue Jul 11, 2024 · 0 comments

Comments

@Odd-Bird
Copy link

Odd-Bird commented Jul 11, 2024

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")
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant