Skip to content

Commit

Permalink
some edits to vignettes
Browse files Browse the repository at this point in the history
  • Loading branch information
santikka committed Mar 30, 2023
1 parent d564e15 commit 310f14d
Show file tree
Hide file tree
Showing 2 changed files with 26 additions and 21 deletions.
27 changes: 13 additions & 14 deletions vignettes/dynamite_priors.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -37,53 +37,52 @@ set.seed(0)

# Default prior distributions in dynamite

The default priors in `dynamite` or chosen to be relatively uninformative (i.e. weakly informative), such that they provide computational stability by slightly regularizing the posterior. The motivation behind the priors is thus similar to other popular Stan-based R packages such as `brms`^[https://CRAN.R-project.org/package=brms] and `rstanarm`^[https://CRAN.R-project.org/package=rstanarm].
The default priors in `dynamite` are chosen to be relatively uninformative (i.e., weakly informative) such that they provide computational stability by slightly regularizing the posterior. The motivation behind the default priors is thus similar to other popular Stan-based R packages such as `brms`^[https://CRAN.R-project.org/package=brms] and `rstanarm`^[https://CRAN.R-project.org/package=rstanarm].

Define $\sigma_x=\max(1, SD(x))$, where $SD(x)$ is the standard deviation of the predictor variable $x$ over groups and non-fixed^[See section [Lagged responses and predictors](https://docs.ropensci.org/dynamite/articles/dynamite.html#sec:lags) in the main vignette.] time points. Define also
Define $\sigma_x=\max(1, \text{SD}(x))$, where $\mbox{SD}(x)$ is the standard deviation of the predictor variable $x$ over groups and non-fixed time points (see section "Lagged responses and predictors" in the main vignette for more information: `vignette("dynamite", package = "dynamite")`). Define also
\[
\sigma_y = \begin{cases}
\max(1, SD(y)), &\text{if family is gaussian or student}\\
\max(1, \text{SD}(y)), &\text{if family is gaussian or student}\\
1, &\text{otherwise}
\end{cases},
\]
where $SD(y)$ is the standard deviation of the response variable as over groups and non-fixed time points.
where $SD(y)$ is the standard deviation of the response variable as over groups and non-fixed time points.

## Regression coefficients

We define the default prior for a time-invariant regression coefficient $\beta$ as well for the first time-varying coefficient $\delta_1$ as zero-mean normal distribution with standard deviation of
We define the default prior for a time-invariant regression coefficients $\beta$ as well for the first time-varying coefficient $\delta_1$ as zero-mean normal distribution with standard deviation of
$2\sigma_y/\sigma_x$. The two maximums used in definitions of $\sigma_y$ and $\sigma_x$ ensure that the prior standard deviation is at least 2.

## Intercept

As the posterior correlations between the intercept $\alpha$ term and the regression coefficients $\beta$ and $\delta$ can be cause computational inefficiencies with gradient based sampling algorithms of Stan, sampling of the intercept is performed indirectly via parameter $a$, so that the intercept $\alpha$ (or $\alpha_1$ in case of time-varying intercept) is constructed as
As the posterior correlations between the intercept $\alpha$ and the regression coefficients $\beta$ and $\delta$ can be cause computational inefficiencies with gradient-based sampling algorithms of Stan, sampling of the intercept is performed indirectly via parameter $a$, so that the intercept $\alpha$ (or $\alpha_1$ in case of time-varying intercept) is constructed as
\[
\alpha = a - \bar X^\beta_1\beta - \bar X_1^\delta\delta_1,
\]
where $\bar X^\beta$ and $\bar X^\delta$ are the means of the corresponding predictors at first (non-fixed) time point. The prior is then defined for $a$ as
\[
a \sim \mbox{N}(\bar y, 2\sigma_y),
\]
where $\bar y$ is mean of the response variable values at first time point after applying the link function to them (except in case of categorical and multinomial case where $\bar y$ is set to zero).
where $\bar y$ is mean of the response variable values at first time point after applying the link function (except in case of categorical and multinomial response where $\bar y$ is set to zero).

## Standard deviation parameters

The prior for the standard deviation parameter $\tau$ of random walk prior of the spline coefficients is half-normal with the standard deviation of $2\sigma_y/\sigma_x$ (with $\sigma_x=1$ in case of varying intercept). Same prior distribution is also used for the the standard deviations of the random effects. Note that this prior can be too uninformative in some cases especially for the $\tau$ parameters.
The prior for the standard deviation parameter $\tau$ of the random walk prior of the spline coefficients is half-normal with the standard deviation of $2\sigma_y/\sigma_x$ (with $\sigma_x=1$ in case of a time-varying intercept). The same prior distribution is also used for the the standard deviations of the random effects. Note that this prior can be too uninformative in some cases especially for the $\tau$ parameters.

## Correlation matrices

The correlation matrix of the random effects uses LKJ(1) prior using the Cholesky parameterization, see Stan documentation of `lkj_cholesky_corr` for details^[https://mc-stan.org/docs/functions-reference/cholesky-lkj-correlation-distribution.html].
The correlation matrix of the random effects uses a $\mbox{LKJ}(1)$ prior with the Cholesky parameterization, see Stan documentation of `lkj_cholesky_corr` for details^[https://mc-stan.org/docs/functions-reference/cholesky-lkj-correlation-distribution.html].

## Family-specific parameters

For standard deviation parameter of gaussian and student responses, we use exponential prior with rate parameter $\frac{1}{2\max(1,\sigma_y)}$. The degrees-of-freedom parameter of the student's $t$-distribution has a $\mbox{Gamma}(2, 0.1)$ prior, whereas other family specific parameters $\phi \sim \mbox{Exponential}(1)$.
For standard deviation parameter of gaussian and student's $t$ responses, we use exponential prior with a rate parameter $\frac{1}{2\max(1,\sigma_y)}$. The degrees-of-freedom parameter of the student's $t$-distribution has a $\mbox{Gamma}(2, 0.1)$ prior, whereas for other family specific parameters we set $\phi \sim \mbox{Exponential}(1)$.

# Defining priors of dynamite models

While `dynamite()` can be used with default prior choices, we recommend to check whether the defaults make sense for your particular problem and modify them accordingly based on the domain knowledge.
While `dynamite()` can be used with the default prior choices, we recommend to check whether the defaults make sense for your particular problem and to odify them accordingly based on the domain knowledge.

Any univariate unbounded continuous distributions supported by Stan can be used as a prior for univariate model parameters (in case the parameter is constrained to be positive, the distribution is automatically truncated). Naturally, any univariate distribution bounded to the positive real line can be used as a prior for parameters constrained to be positive (e.g. standard deviation parameter).
See Stan function reference at \url{https://mc-stan.org/users/documentation/} for possible distribution.
For custom priors, you should first get the default priors with `get_priors()` function, and then modify the `priors` column of the obtained data frame before supplying it to the `dynamite()` function:
Any univariate unbounded continuous distributions supported by Stan can be used as a prior for univariate model parameters (in case the parameter is constrained to be positive, the distribution is automatically truncated). Naturally, any univariate distribution bounded to the positive real line can be used as a prior for parameters constrained to be positive (e.g., a standard deviation parameter).
See the Stan function reference at \url{https://mc-stan.org/users/documentation/} for possible distributions. For custom priors, you should first get the default priors with the `get_priors()` function, and then modify the `priors` column of the obtained data frame before supplying it to the `dynamite()` function:

```{r}
f <- obs(y ~ -1 + z + varying(~ x + lag(y)) + random(~1 + z), "gaussian") +
Expand Down
20 changes: 13 additions & 7 deletions vignettes/dynamite_simulation.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ suppressPackageStartupMessages(library("ggplot2"))

Dynamic multivariate panel models (DMPM) are conceptually challenging due to their multichannel and multivariate nature. These models can also contain a large number of various parameters, especially when the model contains time-varying effects or group-level random effects. This poses a challenge for users that wish to generate data according to such models, as manual specification of the parameters becomes cumbersome.

While the main purpose of `dynamite` is to fit dynamic multivariate panel models to data and to obtain predictions, the package can fortunately also be used to generate data according to a specified model. The used need only specify the model formula, which along with data on the fixed covariates can be used the obtain the dimensions of all parameters in the model via the `get_parameter_dims()` function. As the name suggests, this function will return the names and dimensions of all parameters in the model which need to be specified for data simulation. However, in practice it is not necessary to define values for all parameters, as Stan will automatically generate random initial values for any missing parameters.
While the main purpose of `dynamite` is to fit dynamic multivariate panel models to data and to obtain predictions, the package can fortunately also be used to generate data according to a specified model. The used need only specify the model formula, which along with data on the fixed covariates can be used the obtain the dimensions of all parameters in the model via the `get_parameter_dims()` function. As the name suggests, this function will return the names and dimensions of all parameters in the model which need to be specified for data simulation. However, in practice it is not necessary to define values for all parameters, as Stan will automatically generate random initial values for any missing parameters. See the documentation of `as.data.frame.dynamitefit` for more details on the various parameters.

After the user has specified the parameter values, they should be supplied to `dynamite()` as a list via the argument `init`, which will be used by Stan to initialize the parameters in the underlying Stan model. Furthermore, the number of iterations should be set to `1` and the simulation algorithm to `"Fixed_param"` so that no posterior sampling is carried out. This will result in a `dynamitefit` object that can be subsequently used in `predict` to obtain the desired simulated data.

Expand All @@ -38,7 +38,6 @@ library(dynamite)
set.seed(1)
n_id <- 100L
n_time <- 20L
d <- data.frame(
y = sample(factor(c("a", "b", "c")), size = n_id, replace = TRUE),
x = sample(factor(c("A", "B", "C")), size = n_id, replace = TRUE),
Expand Down Expand Up @@ -72,7 +71,7 @@ get_parameter_dims(x = f, data = d, time = "time", group = "id")
# This is actually computed
get_parameter_dims(categorical_example_fit)
```
Now we have the required information to specify the parameters of the model. The actual values to be chosen for the parameters naturally depends on the scenario and is up to the user. We set the following values for the simulation as the list `init`.
In other words, all `beta` parameters are vectors of length 5, and the `a` parameters are scalars. The `a` type parameters are centered versions of the intercepts `alpha` at the first time index (see the package vignette on default priors for more information: `vignette("dynamite_priors", package = "dynamite")`). Now we have the required information to specify the values for the parameters of the model. The actual values to be chosen naturally depends on the scenario and is up to the user. We set the following values for the simulation as the list `init`.
```{r categparamvals, eval=FALSE, echo=TRUE}
init <- list(
beta_x_B = c(2, 0.8, 0.2, 0, 0),
Expand Down Expand Up @@ -107,14 +106,13 @@ categorical_example <- predict(fit, type = "response") |>

# Model with time-varying effects

As mentioned earlier, simulation from a model with time-varying effects is more challenging due to the inclusion of p-splines. We consider a single-channel model with a gaussian response that has a time-varying effect on a single covariate.
As mentioned earlier, simulation from a model with time-varying effects is more challenging due to the inclusion of p-splines which means that we must also define initial values for the spline coefficients. We consider a single-channel model with a gaussian response that has a time-varying effect on a single covariate.

```{r gaussiainit}
library(dynamite)
set.seed(1)
n_id <- 100L
n_time <- 20L
d <- data.frame(
y = rnorm(n_id, 1, 0.5),
time = 1,
Expand All @@ -139,10 +137,18 @@ f <- obs(y ~ 1 + varying(~ -1 + x + lag(y)), family = "gaussian") +
splines(df = 10)
```
Again, we apply the `get_parameter_dims()` function to get the model parameters and their dimensions.
```{r gaussiandims}
```{r gaussiandimsecho, eval=TRUE, echo=FALSE}
get_parameter_dims(x = f, data = d, time = "time", group = "id")
```
Here, `omega_y` defines the spline coefficients for the two time-varying effects and `tau_y` defines the standard deviations of the random walk priors. We choose the following values for the model parameters:
```{r gaussiandimseval, eval=TRUE, echo=FALSE}
list(
omega_y = c(2L, 10L),
tau_y = 2L,
a_y = 1L,
sigma_y = 1L
)
```
Here, `omega_y` defines the spline coefficients for the time-varying effects. This parameter is a 2 by 10 matrix because we have 2 time-varying effects, `x` and `lag(y)`, and the degrees of freedom of the splines is 10. The parameter `tau_y` defines the standard deviations of the random walk priors for the two time-varying effects. We choose the following values for the model parameters:
```{r gaussianparamvals}
init <- list(
omega_y = rbind(
Expand Down

0 comments on commit 310f14d

Please sign in to comment.