Skip to content

Instantly share code, notes, and snippets.

@jebyrnes
Created April 8, 2019 19:13
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save jebyrnes/afdc436a0052a9c8069161fb6b8f6718 to your computer and use it in GitHub Desktop.
Save jebyrnes/afdc436a0052a9c8069161fb6b8f6718 to your computer and use it in GitHub Desktop.
Does a gam act like a fixed effect for autocorrelation?
library(MASS)
library(tidyverse)
library(ggplot2)
library(mgcv)
gaussprocess <- function(from = 0, to = 1, K = function(s, t) {min(s, t)},
start = NULL, m = 1000) {
# Simulates a Gaussian process with a given kernel
#
# args:
# from: numeric for the starting location of the sequence
# to: numeric for the ending location of the sequence
# K: a function that corresponds to the kernel (covariance function) of
# the process; must give numeric outputs, and if this won't produce a
# positive semi-definite matrix, it could fail; default is a Wiener
# process
# start: numeric for the starting position of the process; if NULL, could
# be randomly generated with the rest
# m: positive integer for the number of points in the process to simulate
#
# return:
# A data.frame with variables "t" for the time index and "xt" for the value
# of the process
t <- seq(from = from, to = to, length.out = m)
Sigma <- sapply(t, function(s1) {
sapply(t, function(s2) {
K(s1, s2)
})
})
path <- mvrnorm(mu = rep(0, times = m), Sigma = Sigma)
if (!is.null(start)) {
path <- path - path[1] + start # Must always start at "start"
}
return(data.frame(xt = path))
}
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}
### Explore
dat <- as_tibble(gaussprocess(m=100)) %>%
mutate(x1 = rnorm(n(), xt, 0.25),
x2 = rnorm(n(), xt, 0.25),
y = rnorm(n(), x1 - x2, 1),
y_only_x1 = rnorm(n(), x1, 0.5),
time = 1:n()
)
pairs(dat, lower.panel = panel.cor)
ggplot(dat %>% slice(1:100) %>% gather(variable, value, -time),
aes(x = time, y = value, color = variable)) +
geom_line()
### Models - linear
mod_linear <- lm(y ~ x1, data = dat)
mod_linear_correct <- lm(y ~ x1+x2, data = dat)
summary(mod_linear)
summary(mod_linear_correct)
plot(mod_linear, which = 2)
acf(mod_linear$residuals)
#correlation in residuals
mod_cor_resid <- gam(y ~ x1, data=dat, correlation=corExp(form = ~time, nugget = TRUE))
summary(mod_cor_resid)
# with basis function
mod_gp <- gam(y ~ x1 + s(time, bs='gp'), data=dat)
summary(mod_gp)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment