Skip to content

Instantly share code, notes, and snippets.

View eric-pedersen's full-sized avatar

Eric Pedersen eric-pedersen

View GitHub Profile
@eric-pedersen
eric-pedersen / link-function-demo.Rmd
Last active July 20, 2022 20:06
testing predictive accuracy when changing link versus a nonlinear functional response
@eric-pedersen
eric-pedersen / prof-app.R
Last active February 6, 2021 00:15
Shiny apps for interactively demonstrating linear regression for Galton's height data. prof-app.R is for the instructor to show guessed fits, and line of best fit. student-app.R is to share with students interactively. Data is shared between apps via a mongodb database
#
# This is a Shiny web application. You can run the application by clicking
# the 'Run App' button above.
#
# Find out more about building applications with Shiny here:
#
# http://shiny.rstudio.com/
#
library(shiny)
library(mgcv)
set.seed(3);n <- 400
############################################
## First example: simulated Tweedie model...
############################################
dat <- gamSim(1,n=n,dist="poisson",scale=.2)
dat$y <- rTweedie(exp(dat$f),p=1.3,phi=.5) ## Tweedie response
b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=tw(),
#Create an explicit binomially distributed set of numbers
n = 1000
frac = 0.9
x = rep(c(1,0),times = c(n*frac, n*(1-frac)))
#Fit a Gaussian model and a binomial model to the same data
gauss_mod = glm(x~1,family = gaussian)
binom_mod = glm(x~1, family= binomial)
#Compare AIC
library(mgcv)
#creating data for the smooth matrix and for predictions
inner_range = c(1900, 2000)
outer_range = c(1850, 2050)
dat = data.frame(year = seq(inner_range[1],inner_range[2],length=500))
pred = data.frame(year = seq(outer_range[1],outer_range[2],length=500))
library(mgcv)
#creating data for the smooth matrix and for predictions
inner_range = c(1900, 2000)
outer_range = c(1850, 2050)
dat = data.frame(year = seq(inner_range[1],inner_range[2],length=500))
pred = data.frame(year = seq(outer_range[1],outer_range[2],length=500))
@eric-pedersen
eric-pedersen / extr2.R
Last active May 29, 2020 18:51 — forked from dill/extr.R
Forecasting with B-splines
# extrapolating into the future with B-splines
# based on code in ?smooth.construct.bs.smooth.spec
library(mgcv)
# annual diameter of women’s skirts at the hem 1866 to 1911
# Hipel and McLeod, 1994
skirts <- scan("http://robjhyndman.com/tsdldata/roberts/skirts.dat",skip=5)
skirtseries <- data.frame(year=1866:1911, diam=skirts)
@eric-pedersen
eric-pedersen / extr-varying-starts.R
Last active May 30, 2020 17:38 — forked from dill/extr.R
Forecasting with B-splines
library(mgcv)
# annual diameter of women’s skirts at the hem 1866 to 1911
# Hipel and McLeod, 1994
skirts <- scan("http://robjhyndman.com/tsdldata/roberts/skirts.dat",skip=5)
skirtseries <- data.frame(year=1866:1911, diam=skirts)
# prediction grid
pd <- data.frame(year=seq(1864, 1960, by=1))
@eric-pedersen
eric-pedersen / testing_maxima.R
Created September 19, 2017 19:53
alternate (GAM-based) test for extrema in fitted curves
library(dplyr)
library(mgcv)
#You need the multivariate normal RNG from the MASS package
mvrnorm = MASS::mvrnorm
# Functions to calculate the 1st and 2nd derivatives of a given time series with a given step size
# Uses two-point approximations for both 1st and 2nd derivs.
calc_1st_deriv = function(y,delta) (lead(y,1) - lag(y,1))/(2*delta)
calc_2nd_deriv = function(y,delta) (lead(y,1) + lag(y,1)-2*y)/delta^2
@eric-pedersen
eric-pedersen / gam_cc_midpoint_test.R
Last active October 19, 2016 19:55
I wanted to see how the choice of midpoint in a circular cubic smooth in mgcv affected model fit. When fitting monthly models, we know that the break between month 1 &12 should occur between them (not right at the values). How does the choice of midpoint affect bias/variance of the fits?
library(dplyr)
library(tidyr)
library(ggplot2)
library(mgcv)
library(purrr)
nreps = 100
#Simulate data from a sinusoidal function. Point of equality should be at
#0/12 (or 0.5/12.5, or 1/13, etc.) by design.
sims = crossing(x= 1:12, rep = 1:nreps)%>%
mutate(y= 4*sin(2*pi*x/12)+rnorm(nreps,0,1))