Wrong alpha level in the overall function
Opened this issue · 0 comments
Odd-Bird commented
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")
}