alpha <- 0.05
n <- length(y)
t <- year - min(year)
fit <- lm(log(y) ~ t)

cgr    <- 100 * expm1(coef(fit)[["t"]])
cgr_ci <- 100 * expm1(confint(fit, "t", level = 1 - alpha))
cgr_p  <- summary(fit)$coefficients["t", "Pr(>|t|)"]

cagr <- 100 * ((y[n] / y[1])^(1 / (year[n] - year[1])) - 1)
aagr <- 100 * mean(diff(y) / head(y, -1))

models <- list(Linear      = lm(y ~ t),
               Quadratic   = lm(y ~ t + I(t^2)),
               Cubic       = lm(y ~ t + I(t^2) + I(t^3)),
               Exponential = lm(log(y) ~ t))

trend <- do.call(rbind, lapply(names(models), function(name) {
  m  <- models[[name]]
  ll <- logLik(m); k <- attr(ll, "df")
  aic <- -2 * as.numeric(ll) + 2 * k
  mean_fit <- fitted(m)
  if (name == "Exponential") {
    aic <- aic + 2 * sum(log(y))
    mean_fit <- exp(mean_fit + mean(residuals(m)^2) / 2)
  }
  data.frame(Model = name,
             AICc  = aic + 2 * k * (k + 1) / (n - k - 1),
             R2    = 1 - sum((y - mean_fit)^2) / sum((y - mean(y))^2))
}))

S <- sum(sapply(1:(n - 1), function(i) sum(sign(y[(i + 1):n] - y[i]))))
ties <- as.numeric(table(y))
varS <- (n * (n - 1) * (2 * n + 5) - sum(ties * (ties - 1) * (2 * ties + 5))) / 18
Z <- (S - sign(S)) / sqrt(varS)
p <- 2 * pnorm(-abs(Z))

pairs  <- combn(n, 2)
slopes <- sort((y[pairs[2, ]] - y[pairs[1, ]]) / (year[pairs[2, ]] - year[pairs[1, ]]))
sen <- median(slopes)
C <- qnorm(1 - alpha / 2) * sqrt(varS)
sen_ci <- slopes[c(ceiling((length(slopes) - C) / 2), floor((length(slopes) + C) / 2) + 1)]

cv    <- 100 * sd(y) / mean(y)
r2adj <- summary(lm(y ~ t))$adj.r.squared
cdvi  <- cv * sqrt(1 - max(min(r2adj, 1), 0))

# area and production: first and last years
yield <- production / area
effects <- c(Area        = yield[1] * diff(area),
             Yield       = area[1] * diff(yield),
             Interaction = diff(area) * diff(yield),
             Total       = diff(production))

U <- sapply(1:(n - 1), function(k) sum(sign(outer(y[1:k], y[(k + 1):n], function(a, b) b - a))))
K <- max(abs(U))
change_after <- year[which.max(abs(U))]
pettitt_p <- min(1, 2 * exp(-6 * K^2 / (n^3 + n^2)))