gavinsimpson’s gists

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 / 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 / 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 / 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 / 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 / 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()))

gavinsimpson / onename

Created June 8, 2016 22:55

Verifying that +ucfagls is my blockchain ID. https://onename.com/ucfagls

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 / 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

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

Gavin Simpson gavinsimpson