walker: Bayesian regression with random walk coefficients

Description Usage Arguments Details Value Note See Also Examples

View source: R/walker.R

Description

Function walker performs Bayesian inference of a linear regression model with time-varying, random walk regression coefficients, i.e. ordinary regression model where instead of constant coefficients the coefficients follow first or second order random walks. All Markov chain Monte Carlo computations are done using Hamiltonian Monte Carlo provided by Stan, using a state space representation of the model in order to marginalise over the coefficients for efficient sampling.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
walker(
  formula,
  data,
  sigma_y_prior = c(2, 0.01),
  beta,
  init,
  chains,
  return_x_reg = FALSE,
  gamma_y = NULL,
  return_data = TRUE,
  ...
)

Arguments

formula

An object of class {formula} with additional terms rw1 and/or rw2 e.g. y ~ x1 + rw1(~ -1 + x2). See details.

data

An optional data.frame or object coercible to such, as in {lm}.

sigma_y_prior

A vector of length two, defining the a Gamma prior for the observation level standard deviation with first element corresponding to the shape parameter and second to rate parameter. Default is Gamma(2, 0.0001). Not used in walker_glm.

beta

A length vector of length two which defines the prior mean and standard deviation of the Gaussian prior for time-invariant coefficients

init

Initial value specification, see sampling. Note that compared to default in rstan, here the default is a to sample from the priors.

chains

Number of Markov chains. Default is 4.

return_x_reg

If TRUE, does not perform sampling, but instead returns the matrix of predictors after processing the formula.

gamma_y

An optional vector defining a damping of the observational level noise. More specifically, σ_t = gamma_t * σ_y.

return_data

if TRUE, returns data input to sampling. This is needed for lfo.

...

Further arguments to sampling.

Details

The rw1 and rw2 functions used in the formula define new formulas for the first and second order random walks. In addition, these functions need to be supplied with priors for initial coefficients and the standard deviations. For second order random walk model, these sigma priors correspond to the standard deviation of slope disturbances. For rw2, also a prior for the initial slope nu needs to be defined. See examples.

Value

A list containing the stanfit object, observations y, and covariates xreg and xreg_new.

Note

Beware of overfitting and identifiability issues. In particular, be careful in not defining multiple intercept terms (only one should be present). By default rw1 and rw2 calls add their own time-varying intercepts, so you should use 0 or -1 to remove some of them (or the time-invariant intercept in the fixed-part of the formula).

See Also

walker_glm for non-Gaussian models.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
set.seed(1)
x <- rnorm(10)
y <- x + rnorm(10)

# different intercept definitions:

# both fixed intercept and time-varying level,
# can be unidentifiable without strong priors:
fit1 <- walker(y ~ rw1(~ x, beta = c(0, 1)), 
  beta = c(0, 1), chains = 1, iter = 1000) 
## Not run: 
# only time-varying level, using 0 or -1 removes intercept:
fit2 <- walker(y ~ 0 + rw1(~ x, beta = c(0, 1)), chains = 1, iter = 1000)

# time-varying level, no covariates:
fit3 <- walker(y ~ 0 + rw1(~ 1, beta = c(0, 1)), chains = 1, iter = 1000)

# fixed intercept no time-varying level:
fit4 <- walker(y ~ rw1(~ 0 + x, beta = c(0, 1)), 
  beta = c(0, 1), chains = 1, iter = 1000) 

# only time-varying effect of x:
fit5 <- walker(y ~ 0 + rw1(~ 0 + x, beta = c(0, 1)), chains = 1, iter = 1000) 

## End(Not run)

## Not run: 

rw1_fit <- walker(Nile ~ -1 + 
  rw1(~ 1, 
    beta = c(1000, 100), 
    sigma = c(2, 0.001)), 
  sigma_y_prior = c(2, 0.005), 
  iter = 2000, chains = 1)
  
rw2_fit <- walker(Nile ~ -1 + 
  rw2(~ 1,
    beta = c(1000, 100), 
    sigma = c(2, 0.001), 
    nu = c(0, 100)), 
  sigma_y_prior = c(2, 0.005), 
  iter = 2000, chains = 1)
  
g_y <- geom_point(data = data.frame(y = Nile, x = time(Nile)), 
  aes(x, y, alpha = 0.5), inherit.aes = FALSE) 
g_rw1 <- plot_coefs(rw1_fit) + g_y
g_rw2 <- plot_coefs(rw2_fit) + g_y
if(require("gridExtra")) {
  grid.arrange(g_rw1, g_rw2, ncol=2, top = "RW1 (left) versus RW2 (right)")
} else {
  g_rw1
  g_rw2
}

y <- window(log10(UKgas), end = time(UKgas)[100])
n <- 100
cos_t <- cos(2 * pi * 1:n / 4)
sin_t <- sin(2 * pi * 1:n / 4)
dat <- data.frame(y, cos_t, sin_t)
fit <- walker(y ~ -1 + 
  rw1(~ cos_t + sin_t, beta = c(0, 10), sigma = c(2, 1)), 
  sigma_y_prior = c(2, 10), data = dat, chains = 1, iter = 2000)
print(fit$stanfit, pars = c("sigma_y", "sigma_rw1"))

plot_coefs(fit)
# posterior predictive check:
pp_check(fit)

newdata <- data.frame(
  cos_t = cos(2 * pi * 101:108 / 4), 
  sin_t = sin(2 * pi * 101:108 / 4))
pred <- predict(fit, newdata)
plot_predict(pred)

# example on scalability
set.seed(1)
n <- 2^12
beta1 <- cumsum(c(0.5, rnorm(n - 1, 0, sd = 0.05)))
beta2 <- cumsum(c(-1, rnorm(n - 1, 0, sd = 0.15)))
x1 <- rnorm(n, mean = 2)
x2 <- cos(1:n)
rw <- cumsum(rnorm(n, 0, 0.5))
signal <- rw + beta1 * x1 + beta2 * x2
y <- rnorm(n, signal, 0.5)

d <- data.frame(y, x1, x2)

n <- 2^(6:12)
times <- numeric(length(n))
for(i in seq_along(n)) {
  times[i] <- sum(get_elapsed_time(
    walker(y ~ 0 + rw1(~ x1 + x2, 
      beta = c(0, 10)), 
      data = d[1:n[i],],
      chains = 1, seed = 1, refresh = 0)$stanfit))
}
plot(log2(n), log2(times))

## End(Not run)

walker documentation built on Sept. 11, 2021, 9:07 a.m.