-
Notifications
You must be signed in to change notification settings - Fork 14
/
README.Rmd
280 lines (225 loc) · 11.2 KB
/
README.Rmd
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
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
---
output: github_document
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "man/figures/README-",
out.width = "100%",
cache = TRUE
)
```
```{r srr-tags, eval = FALSE, echo = FALSE}
#' @srrstats {G1.2} Contains project status badge.
#' @srrstats {G1.4,G1.4a} Package uses roxygen2 for documentation.
#' @srrstats {G2.0, G2.0a, G2.1, G2.1a, G2.2, G2.4, G2.4a, G2.4b, G2.4c, G2.6}
#' Input types and shapes are tested and checked with autotest and converted
#' explicitly when necessary.
#'
#' @srrstats {G2.3, G2.3a, G2.3b} match.arg and tolower are used where
#' applicable.
#' @srrstats {G1.0, G1.3, G1.4, G1.4a, G1.5, G1.6} General
#' documentation, addressed by the vignettes and the corresponding R
#' Journal paper.
#' @srrstats {G1.1} This is the first software to implement the IS-MCMC by
#' Vihola, Helske, and Franks (2020) and first R package to implement delayed
#' acceptance pseudo-marginal MCMC for state space models. The IS-MCMC method
#' is also available in [walker](github.com/helske/walker) package for a
#' limited class of time-varying GLMss (a small subset of the models
#' supported by this package). Some of the functionality for exponential family
#' state space models is also available in [KFAS](github.com/helske/KFAS), and
#' those models can be converted easily to bssm format for Bayesian analysis.
#' @srrstats {G2.4, G2.4a, G2.4b, G2.4c, G2.6} Explicit conversions are used
#' where necessary.
#'
#' @srrstats {G2.14, G2.14a, G2.14b, G2.14c, G2.15, G2.16} Missing observations
#' (y) are handled automatically as per SSM theory, whereas missing values are
#' not allowed elsewhere. Inputing or ignoring them does not make sense in time
#' series context.
#'
#' @srrstats {G3.0} No floating point equality comparisons are made.
#'
#' @srrstats {G5.4, G5.4a, G5.4b, G5.4c, G5.5, G5.6, G5.6a, G5.6b, G5.7} and
#' @srrstats {BS4.0, BS4.1} The algorithms work as defined per Vihola, Helske,
#' Franks (2020) (all simulations were implemented with the bssm package) and
#' Helske and Vihola (2021). Full replication of the results would take
#' days/weeks (but see also bsm_ng, negbin_series and several testthat tests).
#'
#' @srrstats {G5.8, G5.8a, G5.8b, G5.8c, G5.8d} Tested with autotest and the
#' testthat tests.
#' @srrstats {G5.9, G5.9a, G5.9b} Tested with autotest and the testthat tests.
#'
#' @srrstats {BS1.0, BS1.1, BS1.2, BS1.2a, BS1.2b, BS1.3b} Addressed in the
#' models.R, run_mcmc.R, in vignettes and in the R Journal paper.
#'
#' @srrstats {BS2.1, BS2.1a, BS2.6} Tested and demonstrated by autotest and
#' package examples/tests.
#' @srrstats {BS7.4, BS7.4a} The scales do not matter (in terms of runtime)
#' in random walk Metropolis nor in particle filters, as long as numerical
#' issues are not encountered
```
# bssm
<!-- badges: start -->
[](https://www.repostatus.org/#active)
[](https://github.com/ropensci/software-review/issues/489)
[](https://github.com/helske/bssm/actions)
[](https://app.codecov.io/gh/helske/bssm)
[]( https://CRAN.R-project.org/package=bssm)
[](https://cranlogs.r-pkg.org/badges/bssm)
<!-- badges: end -->
The `bssm` R package provides efficient methods for Bayesian inference of state
space models via particle Markov chain Monte Carlo and importance sampling type
weighted MCMC.
Currently Gaussian, Poisson, binomial, negative binomial, and Gamma observation
densities with linear-Gaussian state dynamics, as well as general non-linear
Gaussian models and discretely observed latent diffusion processes are
supported.
For details, see
* [The bssm paper on The R Journal](https://journal.r-project.org/archive/2021/RJ-2021-103/index.html),
* [Package vignettes at CRAN](https://CRAN.R-project.org/package=bssm)
* Paper on [Importance sampling type estimators based on approximate marginal Markov chain Monte Carlo](https://onlinelibrary.wiley.com/doi/abs/10.1111/sjos.12492)
There are also couple posters and a talk related to IS-correction methodology and bssm package:
* [UseR!2021 talk slides](https://jounihelske.netlify.app/talk/user2021/)
* [SMC 2017 workshop: Accelerating MCMC with an approximation ](http://users.jyu.fi/~jovetale/posters/SMC2017)
* [UseR!2017: Bayesian non-Gaussian state space models in R](http://users.jyu.fi/~jovetale/posters/user2017.pdf)
The `bssm` package was originally developed with the support of Academy of Finland grants 284513, 312605, 311877, and 331817. Current development is focused on increased usability. For recent changes, see NEWS file.
### Citing the package
If you use the `bssm` package in publications, please cite the corresponding R Journal paper:
Jouni Helske and Matti Vihola (2021). "bssm: Bayesian Inference of Non-linear and Non-Gaussian State Space Models in R." The R Journal (2021) 13:2, pages 578-589. https://journal.r-project.org/archive/2021/RJ-2021-103/index.html
## Installation
You can install the released version of bssm from [CRAN](https://CRAN.R-project.org) with:
```{r, eval=FALSE}
install.packages("bssm")
```
And the development version from [GitHub](https://github.com/) with:
```{r, eval=FALSE}
# install.packages("devtools")
devtools::install_github("helske/bssm")
```
Or from R-universe with
```{r, eval = FALSE}
install.packages("bssm", repos = "https://helske.r-universe.dev")
```
## Example
Consider the daily air quality measurements in New Your from May to September 1973, available in the `datasets` package. Let's try to predict the missing ozone levels by simple linear-Gaussian local linear trend model with temperature and wind as explanatory variables (missing response variables are handled naturally in the state space modelling framework, however no missing values in covariates are normally allowed);
```{r example}
library("bssm")
library("dplyr")
library("ggplot2")
set.seed(1)
data("airquality", package = "datasets")
# Covariates as matrix. For complex cases, check out as_bssm function
xreg <- airquality |> select(Wind, Temp) |> as.matrix()
model <- bsm_lg(airquality$Ozone,
xreg = xreg,
# Define priors for hyperparameters (i.e. not the states), see ?bssm_prior
# Initial value followed by parameters of the prior distribution
beta = normal_prior(rep(0, ncol(xreg)), 0, 1),
sd_y = gamma_prior(1, 2, 0.01),
sd_level = gamma_prior(1, 2, 0.01),
sd_slope = gamma_prior(1, 2, 0.01))
fit <- run_mcmc(model, iter = 20000, burnin = 5000)
fit
obs <- data.frame(Time = 1:nrow(airquality),
Ozone = airquality$Ozone) |> filter(!is.na(Ozone))
pred <- fitted(fit, model)
pred |>
ggplot(aes(x = Time, y = Mean)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`),
alpha = 0.5, fill = "steelblue") +
geom_line() +
geom_point(data = obs,
aes(x = Time, y = Ozone), colour = "Tomato") +
theme_bw()
```
Same model but now assuming observations are from Gamma distribution:
```{r gamma-example}
model2 <- bsm_ng(airquality$Ozone,
xreg = xreg,
beta = normal(rep(0, ncol(xreg)), 0, 1),
distribution = "gamma",
phi = gamma_prior(1, 2, 0.01),
sd_level = gamma_prior(1, 2, 0.1),
sd_slope = gamma_prior(1, 2, 0.1))
fit2 <- run_mcmc(model2, iter = 20000, burnin = 5000, particles = 10)
fit2
```
Comparison:
```{r compare}
pred2 <- fitted(fit2, model2)
bind_rows(list(Gaussian = pred, Gamma = pred2), .id = "Model") |>
ggplot(aes(x = Time, y = Mean)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`, fill = Model),
alpha = 0.25) +
geom_line(aes(colour = Model)) +
geom_point(data = obs,
aes(x = Time, y = Ozone)) +
theme_bw()
```
Now let's assume that we also want to use the solar radiation variable as predictor for ozone. As it contains few missing values, we cannot use it directly. As the number of missing time points is very small, simple imputation would likely be acceptable, but let's consider more another approach. For simplicity, the slope terms of the previous models are now omitted, and we focus on the Gaussian case. Let $\mu_t$ be the true solar radiation at time $t$. Now for ozone $O_t$ we assume following model:
$O_t = D_t + \alpha_t + \beta_S \mu_t + \sigma_\epsilon \epsilon_t$\
$\alpha_{t+1} = \alpha_t + \sigma_\eta\eta_t$\
$\alpha_1 \sim N(0, 100^2\textrm{I})$,\
wheere $D_t = \beta X_t$ contains regression terms related to wind and temperature, $\alpha_t$ is the time varying intercept term, and $\beta_S$ is the effect of solar radiation $\mu_t$.
Now for the observed solar radiation $S_t$ we assume
$S_t = \mu_t$\
$\mu_{t+1} = \mu_t + \sigma_\xi\xi_t,$\
$\mu_1 \sim N(0, 100^2)$,\
i.e. we assume as simple random walk for the $\mu$ which we observe without error or not at all (there is no error term in the observation equation $S_t=\mu_t$).
We combine these two models as a bivariate Gaussian model with `ssm_mlg`:
```{r missing-values}
# predictors (not including solar radiation) for ozone
xreg <- airquality |> select(Wind, Temp) |> as.matrix()
# Function which outputs new model components given the parameter vector theta
update_fn <- function(theta) {
D <- rbind(t(xreg %*% theta[1:2]), 1)
Z <- matrix(c(1, 0, theta[3], 1), 2, 2)
R <- diag(exp(theta[4:5]))
H <- diag(c(exp(theta[6]), 0))
# add third dimension so we have p x n x 1, p x m x 1, p x p x 1 arrays
dim(Z)[3] <- dim(R)[3] <- dim(H)[3] <- 1
list(D = D, Z = Z, R = R, H = H)
}
# Function for log-prior density
prior_fn <- function(theta) {
sum(dnorm(theta[1:3], 0, 10, log = TRUE)) +
sum(dgamma(exp(theta[4:6]), 2, 0.01, log = TRUE)) +
sum(theta[4:6]) # log-jacobian
}
init_theta <- c(0, 0, 0, log(50), log(5), log(20))
comps <- update_fn(init_theta)
model <- ssm_mlg(y = cbind(Ozone = airquality$Ozone, Solar = airquality$Solar.R),
Z = comps$Z, D = comps$D, H = comps$H, T = diag(2), R = comps$R,
a1 = rep(0, 2), P1 = diag(100, 2), init_theta = init_theta,
state_names = c("alpha", "mu"), update_fn = update_fn,
prior_fn = prior_fn)
fit <- run_mcmc(model, iter = 60000, burnin = 10000)
fit
```
Draw predictions:
```{r bivariate-fig}
pred <- fitted(fit, model)
obs <- data.frame(Time = 1:nrow(airquality),
Solar = airquality$Solar.R) |> filter(!is.na(Solar))
pred |> filter(Variable == "Solar") |>
ggplot(aes(x = Time, y = Mean)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`),
alpha = 0.5, fill = "steelblue") +
geom_line() +
geom_point(data = obs,
aes(x = Time, y = Solar), colour = "Tomato") +
theme_bw()
obs <- data.frame(Time = 1:nrow(airquality),
Ozone = airquality$Ozone) |> filter(!is.na(Ozone))
pred |> filter(Variable == "Ozone") |>
ggplot(aes(x = Time, y = Mean)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`),
alpha = 0.5, fill = "steelblue") +
geom_line() +
geom_point(data = obs,
aes(x = Time, y = Ozone), colour = "Tomato") +
theme_bw()
```
See more examples in the paper, vignettes, and in the docs.