Interval between extremely wet years increasing?

16 07 2021

The other day I was playing around with some Bureau of Meteorology data for my little patch of the Adelaide Hills (free data — how can I resist?), when I discovered an interesting trend.

Living on a little farm with a small vineyard, I’m rather keen on understanding our local weather trends. Being a scientist, I’m also rather inclined to analyse data.

My first question was given the strong warming trend here and everywhere else, plus ample evidence of changing rainfall patterns in Australia (e.g., see here, here, here, here, here), was it drying out, getting wetter, or was the seasonal pattern of rainfall in my area changing?

I first looked to see if there was any long-term trend in total annual rainfall over time. Luckily, the station records nearest my farm go all the way back to 1890:

While the red line might suggest a slight decrease since the late 19th Century, it’s no different to an intercept-only model (evidence ratio = 0.84) — no trend.

Here’s the R code to do that analysis (you can download the data here, or provide your own data in the same format):

## IMPORT MONTHLY PRECIPITATION DATA
dat <- read.table("monthlyprecipdata.csv", header=T, sep=",")

## CALCULATE ANNUAL VECTORS
precip.yr.sum <- xtabs(dat$Monthly.Precipitation.Total..millimetres. ~ dat$Year)
precip.yr.sum <- precip.yr.sum[-length(precip.yr.sum)]
year.vec <- as.numeric(names(precip.yr.sum))

## PLOT
plot(year.vec, as.numeric(precip.yr.sum), type="l", pch=19, xlab="year", ylab="annual precipitation (mm)")
fit.yr <- lm(precip.yr.sum ~ year.vec)
abline(fit.yr, lty=2, lwd=2, col="red")
abline(h=mean(as.numeric(precip.yr.sum)),lty=2, lwd=3)

## TEST FOR TREND
# functions
AICc <- function(...) {
  models <- list(...)
  num.mod <- length(models)
  AICcs <- numeric(num.mod)
  ns <- numeric(num.mod)
  ks <- numeric(num.mod)
  AICc.vec <- rep(0,num.mod)
  for (i in 1:num.mod) {
    if (length(models[[i]]$df.residual) == 0) n <- models[[i]]$dims$N else n <- length(models[[i]]$residuals)
    if (length(models[[i]]$df.residual) == 0) k <- sum(models[[i]]$dims$ncol) else k <- (length(models[[i]]$coeff))+1
    AICcs[i] <- (-2*logLik(models[[i]])) + ((2*k*n)/(n-k-1))
    ns[i] <- n
    ks[i] <- k
    AICc.vec[i] <- AICcs[i]
  }
  return(AICc.vec)
}

delta.AIC <- function(x) x - min(x) ## where x is a vector of AIC
weight.AIC <- function(x) (exp(-0.5*x))/sum(exp(-0.5*x)) ## Where x is a vector of dAIC
ch.dev <- function(x) ((( as.numeric(x$null.deviance) - as.numeric(x$deviance) )/ as.numeric(x$null.deviance))*100) ## % change in deviance, where x is glm object

linreg.ER <- function(x,y) { # where x and y are vectors of the same length; calls AICc, delta.AIC, weight.AIC functions
  fit.full <- lm(y ~ x); fit.null <- lm(y ~ 1)
  AIC.vec <- c(AICc(fit.full),AICc(fit.null))
  dAIC.vec <- delta.AIC(AIC.vec); wAIC.vec <- weight.AIC(dAIC.vec)
  ER <- wAIC.vec[1]/wAIC.vec[2]
  r.sq.adj <- as.numeric(summary(fit.full)[9])
  return(c(ER,r.sq.adj))
}

linreg.ER(year.vec, as.numeric(precip.yr.sum))
Read the rest of this entry »




Influential conservation ecology papers of 2017

27 12 2017

Gannet Shallow Diving 03
As I have done for the last four years (20162015, 2014, 2013), here’s another retrospective list of the top 20 influential conservation papers of 2017 as assessed by experts in F1000 Prime.

Read the rest of this entry »