Gavin Simpson 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ## 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
December 6, 2021 17:05
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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ################################################ | |
| ## 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| `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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ## 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| `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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| 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
July 21, 2017 13:29
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)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ## Compare Brooks et al glmmTMB paper with mgcv | |
| ## Packages | |
| library("glmmTMB") | |
| library("mgcv") | |
| library("ggplot2") | |
| theme_set(theme_bw()) | |
| library("ggstance") | |
| ## Salamander |
OlderNewer