Skip to content

Instantly share code, notes, and snippets.

View gavinsimpson's full-sized avatar

Gavin Simpson gavinsimpson

View GitHub Profile
@gavinsimpson
gavinsimpson / load-process-cet-monthly.R
Created May 12, 2014 19:24
R code to load and process the monthly data set from the Central England Temperature time series
## read in the CET data
CET <- url("http://www.metoffice.gov.uk/hadobs/hadcet/cetml1659on.dat")
cet <- read.table(CET, sep = "", skip = 6, header = TRUE,
fill = TRUE, na.string = c(-99.99, -99.9))
names(cet) <- c(month.abb, "Annual") ## fix up df names
## remove last row of incomplete data
cet <- cet[-nrow(cet), ]
## get rid of the annual too - store for plotting
rn <- as.numeric(rownames(cet))
@gavinsimpson
gavinsimpson / Deriv.R
Last active May 25, 2024 08:50
R functions to compute and plot the first derivative of a spline term in a GAM(M) fitted by `gam()` or `gamm()` in package mgcv
################################################
## Functions for derivatives of GAM(M) models ##
################################################
Deriv <- function(mod, n = 200, eps = 1e-7, newdata, term) {
if(inherits(mod, "gamm"))
mod <- mod$gam
m.terms <- attr(terms(mod), "term.labels")
if(missing(newdata)) {
newD <- sapply(model.frame(mod)[, m.terms, drop = FALSE],
function(x) seq(min(x), max(x), length = n))
@gavinsimpson
gavinsimpson / derivSimulCI.R
Last active March 21, 2019 08:23
First derivatives and simultaneous confidence intervals for GAM spline terms
`derivSimulCI` <- function(mod, n = 200, eps = 1e-7, newdata, term,
samples = 10000) {
stopifnot(require("MASS"))
if(inherits(mod, "gamm"))
mod <- mod$gam
m.terms <- attr(terms(mod), "term.labels")
if(missing(newdata)) {
newD <- sapply(model.frame(mod)[, m.terms, drop = FALSE],
function(x) seq(min(x), max(x) - (2*eps), length = n))
names(newD) <- m.terms
@gavinsimpson
gavinsimpson / genURLS.R
Last active November 6, 2018 16:13
Functions to automate downloading data from the Government of Canada's Climate/Weather website
genURLS <- function(id, start, end, timeframe = c("hourly", "daily", "monthly")) {
years <- seq(start, end, by = 1)
nyears <- length(years)
timeframe <- match.arg(timeframe)
if (isTRUE(all.equal(timeframe, "hourly"))) {
years <- rep(years, each = 12)
months <- rep(1:12, times = nyears)
ids <- rep(id, nyears * 12)
} else if (isTRUE(all.equal(timeframe, "daily"))) {
months <- 1 # this is essentially arbitrary & ignored if daily
@gavinsimpson
gavinsimpson / loadCET
Last active November 21, 2015 23:02
A simple function to load and process the CET data into a nice format for modelling
## read in the CET data
`loadCET` <- function() {
CET <- url("http://www.metoffice.gov.uk/hadobs/hadcet/cetml1659on.dat")
on.exit(close(CET))
cet <- read.table(CET, sep = "", skip = 6, header = TRUE,
fill = TRUE, na.string = c(-99.99, -99.9))
names(cet) <- c(month.abb, "Annual")
## remove last row of incomplete data
cet <- cet[-nrow(cet), ] # FIXME: this removes the last row regardless
## get rid of the annual too - store for plotting
@gavinsimpson
gavinsimpson / simulate.gamm.R
Last active October 26, 2018 08:46
S3 method for simulate() for "gamm" objects from package mgcv
`simulate.gamm` <- function(object, nsim = 1, seed = NULL, newdata,
freq = FALSE, unconditional = FALSE, ...) {
if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
runif(1)
if (is.null(seed))
RNGstate <- get(".Random.seed", envir = .GlobalEnv)
else {
R.seed <- get(".Random.seed", envir = .GlobalEnv)
set.seed(seed)
RNGstate <- structure(seed, kind = as.list(RNGkind()))
Verifying that +ucfagls is my blockchain ID. https://onename.com/ucfagls
@gavinsimpson
gavinsimpson / asciify.R
Last active September 7, 2016 05:32
Create a MySQL-like ASCII table of data
asciify <- function(df, pad = 1, ...) {
## error checking
stopifnot(is.data.frame(df))
## internal functions
SepLine <- function(n, pad = 1) {
tmp <- lapply(n, function(x, pad) paste(rep("-", x + (2* pad)),
collapse = ""),
pad = pad)
paste0("+", paste(tmp, collapse = "+"), "+")
}
@gavinsimpson
gavinsimpson / modelled-nuuk-rainfall.png
Last active April 28, 2021 15:43
R code to download, extract, and fit a Tweedie GAM to monthly rainfall total time series from Nuuk, Greenland, using mgcv
modelled-nuuk-rainfall.png
@gavinsimpson
gavinsimpson / compare-glmmTMB-with-mgcv.R
Last active April 30, 2024 20:09
A quick R script I knocked up to compare the glmmTMB and mgcv packages for fitting zero-inflated GLMMs to the Salamander and Owls data sets from Brooks et al (2017)
## Compare Brooks et al glmmTMB paper with mgcv
## Packages
library("glmmTMB")
library("mgcv")
library("ggplot2")
theme_set(theme_bw())
library("ggstance")
## Salamander