Created
April 8, 2019 19:13
-
-
Save jebyrnes/afdc436a0052a9c8069161fb6b8f6718 to your computer and use it in GitHub Desktop.
Does a gam act like a fixed effect for autocorrelation?
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
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