diff --git a/pkg/R/trim_overall.R b/pkg/R/trim_overall.R index d888b3f..c5d7e02 100644 --- a/pkg/R/trim_overall.R +++ b/pkg/R/trim_overall.R @@ -20,8 +20,8 @@ #' specifying the \code{trim} model (See also the example below). #' #' -#' @return a list of class \code{trim.overall} containing overall slope -#' coefficients (\code{coef}), augmented wih p-values and an interpretation). +#' @return a list of class \code{trim.overall} containing, a.o., overall slope +#' coefficients (\code{slope}), augmented wih p-values and an interpretation). #' @export #' #' @family analyses @@ -36,7 +36,7 @@ #' # Overall is a list, you can get information out if it using the $ syntax, #' # for example #' L <- overall(z) -#' L$coef +#' L$slope #' #' # Obtain the slope from changepoint to changepoint #' z <- trim(count ~ site + time, data=skylark, model=2,changepoints=c(1,4,6)) @@ -59,7 +59,10 @@ overall <- function(x, which=c("imputed","fitted"), changepoints=numeric(0)) { tt_imp <- x$tt_imp var_tt_mod <- x$var_tt_mod var_tt_imp <- x$var_tt_imp - ntime <- x$ntime + J <- ntime <- x$ntime + + # Set changepoints in case none are given + if (length(changepoints)==0) changepoints <- 1 # if (base>0) { # use index instead # browser() @@ -87,23 +90,24 @@ overall <- function(x, which=c("imputed","fitted"), changepoints=numeric(0)) { tval <- qt((1-alpha/2), df) blo <- bhat - tval * berr bhi <- bhat + tval * berr + + # First priority: evidece for a strong trend? if (blo[2] > 1.05) return("Strong increase (p<0.001)") if (bhi[2] < 0.95) return("Strong decrease (p<0.001)") if (blo[1] > 1.05) return("Strong increase (p<0.05)") if (bhi[1] < 0.95) return("Strong decrease (p<0.05)") + # Second prority: evidence for a moderate trend? eps = 1e-7 # required to get a correct interpretation for slope=0.0 (Stable) if (blo[2] > 1.0+eps) return("Moderate increase (p<0.001)") if (bhi[2] < 1.0-eps) return("Moderate decrease (p<0.001)") if (blo[1] > 1.0+eps) return("Moderate increase (p<0.05)") if (bhi[1] < 1.0-eps) return("Moderate decrease (p<0.05)") - # TO discuss: order of evaluation matters. What's more important? - # - strong increase p<0.05 or - # - moderate increase p<0.001 - + # Third priority: evidency for stability? if (blo[1]>0.95 && bhi[1]<1.05) return("Stable") + # Leftover category: uncertain return("Uncertain") } @@ -157,8 +161,8 @@ overall <- function(x, which=c("imputed","fitted"), changepoints=numeric(0)) { row.names = c("intercept","slope") ) - z$t = z$mul / z$se_mul - z$p = if (df>0) 2 * pt(abs(z$t), df, lower.tail=FALSE) + tval = z$mul / z$se_mul + z$p = if (df>0) 2 * pt(abs(tval), df, lower.tail=FALSE) else NA z$meaning = c("", .meaning(z$mul[2], z$se_mul[2], n-2)) list(src=src, coef=z, SSR=SSR) @@ -186,22 +190,20 @@ overall <- function(x, which=c("imputed","fitted"), changepoints=numeric(0)) { stopifnot(all(diff(changepoints)>0)) ncp <- length(changepoints) cpx <- c(changepoints, J) # Extend list of overall changepoints with final year - coef.collector = data.frame() - int.collector = data.frame() # here go the intercepts - SSR.collector = numeric(ncp) + int.collector <- data.frame() # Here go the intercepts + slp.collector <- data.frame() # Here go the slopes + SSR.collector <- numeric(ncp) # Here go the SSR info for (i in 1:ncp) { idx <- cpx[i] : cpx[i+1] tmp <- .compute.overall.slope(idx, tt[idx], var_tt[idx,idx], src) - slope <- tmp$coef[2,] # slope is on second row of output dataframe prefix <- data.frame(from=x$time.id[cpx[i]], upto=x$time.id[cpx[i+1]]) - coef.collector <- rbind(coef.collector, cbind(prefix, slope)) - intercept <- tmp$coef[1,] # slope is on second row of output dataframe + intercept <- tmp$coef[1,] # Intercept is on first row of output dataframe int.collector <- rbind(int.collector, cbind(prefix, intercept)) + slope <- tmp$coef[2,] # Slope is on second row of output dataframe + slp.collector <- rbind(slp.collector, cbind(prefix, slope)) SSR.collector[i] = tmp$SSR } - out <- list(src=src - , coef=coef.collector, intercept = int.collector, SSR=SSR.collector - , type="changept") #store 'n' mark + out <- list(src=src, slope=slp.collector, intercept=int.collector, SSR=SSR.collector) } out$J = J out$tt = tt @@ -220,24 +222,7 @@ overall <- function(x, which=c("imputed","fitted"), changepoints=numeric(0)) { #' @export #' @keywords internal print.trim.overall <- function(x,...) { - if (x$type=="normal") { - # Compute 95\% confidence interval of multiplicative slope - bhat <- x$coef[[3]][2] # multiplicative trend (i.e., not log-transformed) - berr <- x$coef[[4]][2] # corresponding standard error - level <- 0.95 - alpha <- 1 - level - df <- x$J-2 - tval <- qt((1-alpha/2), df) - blo <- bhat - tval * berr - bhi <- bhat + tval * berr - # Build an informative string - info <- sprintf("conf.int (mul; level=%d%%)=[%f, %f]", level*100, blo, bhi) - printf("Overall slope (%s): %s\n", x$src, info) - print(x$coef, row.names=TRUE) - } else if (x$type=="changept") { - printf("Overall slope (%s) for changepoints\n", x$src) - print(x$coef, row.names=FALSE) - } else stop("Can't happen") + print(x$slope, row.names=FALSE) } #--------------------------------------------------------------------- Plot ---- @@ -281,6 +266,7 @@ plot.trim.overall <- function(x, imputed=TRUE, ...) { trend.line <- NULL conf.band <- NULL + X$type <- "changept" # Hack for merging overall/changepts if (X$type=="normal") { # Trend line a <- X$coef[[1]][1] # intercept @@ -304,14 +290,14 @@ plot.trim.overall <- function(x, imputed=TRUE, ...) { yconf <- c(ylo, rev(yhi)) conf.band <- cbind(xconf, yconf) } else if (X$type=="changept") { - nsegment = nrow(X$coef) + nsegment = nrow(X$slope) for (i in 1:nsegment) { # Trend line - a <- X$intercept[i,3] # intercept - b <- X$coef[i,3] # slope - from <- which(tpt==X$coef[i,1]) # convert year -> time - upto <- which(tpt==X$coef[i,2]) + a <- X$intercept[i,3] + b <- X$slope[i,3] + from <- which(tpt==X$slope[i,1]) # convert year -> time + upto <- which(tpt==X$slope[i,2]) delta = (upto-from)*10 x <- seq(from, upto, length.out=delta) # continue timepoint 1..J ytrend <- exp(a + b*x) diff --git a/pkg/tests/testthat/test_trim.R b/pkg/tests/testthat/test_trim.R index 02faa6a..53d0ba8 100644 --- a/pkg/tests/testthat/test_trim.R +++ b/pkg/tests/testthat/test_trim.R @@ -32,9 +32,9 @@ trimtest <- function(m, to, tc, vcv=NULL) { # Overall slope tgt <- get_overal_imputed_slope(to) - out <- overall(m,"imputed")$coef[2,] # Slope coefs are in second row + out <- overall(m,"imputed")$slope[1,] for (i in seq_len(ncol(tgt))) { - expect_true(max(abs(out[,i] - tgt[,i])) < 1e-4, + expect_true(max(abs(out[,i+2] - tgt[,i])) < 1e-4, info = sprintf("Overall slope column %d", i)) } @@ -43,7 +43,7 @@ trimtest <- function(m, to, tc, vcv=NULL) { tgt <- get_overal_cp_imputed_slope(to) out <- overall(m, "imputed", tc$overallchangepoints) for (i in seq_len(ncol(tgt))) { - expect_true(max(abs(out$coef[,i]-tgt[,i]), na.rm=TRUE) < 1e-4, + expect_true(max(abs(out$slope[,i]-tgt[,i]), na.rm=TRUE) < 1e-4, info=sprintf("Overall slope (changepts) column %d",i)) } }