Now consistent output for overall() applied to whole period, or chang…

…epoints

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("<none>", .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)

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