Step by step guide for creating a package that depends on RStan

Stefan Siegert, Jonah Gabry, Martin Lysy, and Ben Goodrich

2020-07-06

Introduction

In this vignette we will walk through the steps necessary for creating an R package that depends on Stan by creating a package with one function that fits a simple linear regression. Before continuing, we recommend that you first read the other vignette Guidelines for Developers of R Packages Interfacing with Stan.

Creating the package skeleton

The rstantools package offers two methods for adding Stan functionality to R packages:

Here we will use rstan_create_package() to initialize a bare-bones package directory. The name of our demo package will be rstanlm; it will fit a simple linear regression model using Stan.

library("rstantools")
rstan_create_package(path = 'rstanlm')
This is rstantools version 2.1.1
Creating package skeleton for package: rstanlm
Package: rstanlm
Title: What the Package Does (One Line, Title Case)
Version: 0.0.0.9000
[email protected] (parsed):
    * First Last <[email protected]> [aut, cre] (<https://orcid.org/YOUR-ORCID-ID>)
Description: What the package does (one paragraph).
License: `use_mit_license()`, `use_gpl3_license()` or friends to pick a
    license
Encoding: UTF-8
LazyData: true
Roxygen: list(markdown = TRUE)
RoxygenNote: 7.1.1
Creating inst/stan/include directory ...
Creating inst/include directory ...
Creating src directory ...
Updating DESCRIPTION ...
Adding 'configure' files ...

Next, add the following lines (e.g., via <package-name>-package.R if using roxygen) to your NAMESPACE:

import(Rcpp)
import(methods)
importFrom(rstan, sampling)
useDynLib(rstanlm, .registration = TRUE)
Done.
Adding rstanlm-package.R file ...
Adding .travis.yml file ...
Adding .gitignore file ...
Adding .Rbuildignore file ...
Configuring Stan compile and module export instructions ...
Further Stan-specific steps are described in 'rstanlm/Read-and-delete-me'.

If we had existing .stan files to include with the package we could use the optional stan_files argument to rstan_create_package() to include them. Another option, which we’ll use below, is to add the Stan files once the basic structure of the package is in place.

We can now set the new working directory to the new package directory and view the contents. (Note: if using RStudio then by default the newly created project for the package will be opened in a new session and you will not need the call to setwd().)

setwd("rstanlm")
list.files(all.files = TRUE)
 [1] "."                  ".."                 ".Rbuildignore"     
 [4] ".gitignore"         ".travis.yml"        "DESCRIPTION"       
 [7] "NAMESPACE"          "R"                  "Read-and-delete-me"
[10] "configure"          "configure.win"      "inst"              
[13] "src"               
file.show("DESCRIPTION")
Package: rstanlm
Title: What the Package Does (One Line, Title Case)
Version: 0.0.0.9000
[email protected]: 
    person(given = "First",
           family = "Last",
           role = c("aut", "cre"),
           email = "[email protected]",
           comment = c(ORCID = "YOUR-ORCID-ID"))
Description: What the package does (one paragraph).
License: `use_mit_license()`, `use_gpl3_license()` or friends to pick a
    license
Encoding: UTF-8
LazyData: true
Roxygen: list(markdown = TRUE)
RoxygenNote: 7.1.1
Biarch: true
Depends: 
    R (>= 3.4.0)
Imports: 
    Rcpp (>= 0.12.0),
    RcppParallel (>= 5.0.1),
    methods,
    rstan (>= 2.18.1),
    rstantools (>= 2.1.1)
LinkingTo: 
    BH (>= 1.66.0),
    Rcpp (>= 0.12.0),
    RcppEigen (>= 0.3.3.3.0),
    RcppParallel (>= 5.0.1),
    StanHeaders (>= 2.18.0),
    rstan (>= 2.18.1)
SystemRequirements: GNU make

Some of the sections in the DESCRIPTION file need to be edited by hand (e.g., Title, Author, Maintainer, and Description, but these also can be set with the fields argument to rstan_create_package()). However, rstan_create_package() has added the necessary packages and versions to Depends, Imports, and LinkingTo to enable Stan functionality.

Read-and-delete-me file

Before deleting the Read-and-delete-me file in the new package directory make sure to read it because it contains some important instructions about customizing your package:

file.show("Read-and-delete-me")
Stan-specific notes:

* All '.stan' files containing stanmodel definitions must be placed in 'inst/stan'.
* Additional files to be included by stanmodel definition files
  (via e.g., #include "mylib.stan") must be placed in any subfolder of 'inst/stan'.
* Additional C++ files needed by any '.stan' file must be placed in 'inst/include',
  and can only interact with the Stan C++ library via '#include' directives
  placed in the file 'inst/include/stan_meta_header.hpp'.
* The precompiled stanmodel objects will appear in a named list called 'stanmodels',
  and you can call them with e.g., 'rstan::sampling(stanmodels$foo, ...)'

You can move this file out of the directory, delete it, or list it in the .Rbuildignore file if you want to keep it in the directory.

file.remove('Read-and-delete-me')
[1] TRUE

Stan files

Our package will call rstan’s sampling() method to use MCMC to fit a simple linear regression model for an outcome variable y with a single predictor x. After writing the necessary Stan program, the file should be saved with a .stan extension in the inst/stan subdirectory. We’ll save the following program to inst/stan/lm.stan:

// Save this file as inst/stan/lm.stan
data {
  int<lower=1> N;
  vector[N] x;
  vector[N] y;
}
parameters {
  real intercept;
  real beta;
  real<lower=0> sigma;
}
model {
  // ... priors, etc.

  y ~ normal(intercept + beta * x, sigma);
}

The inst/stan subdirectory can contain additional Stan programs if required by your package. During installation, all Stan programs will be compiled and saved in the list stanmodels that can then be used by R function in the package. The rule is that the Stan program compiled from the model code in inst/stan/foo.stan is stored as list element stanmodels$foo. Thus, the filename of the Stan program in the inst/stan directory should not contain spaces or dashes and nor should it start with a number or utilize non-ASCII characters.

R files

We next create the file R/lm_stan.R where we define the function lm_stan() in which our compiled Stan model is being used. Setting the rstan_create_package() argument roxygen = TRUE (the default value) enables roxygen2 documentation for the package functions. The following comment block placed in lm_stan.R ensures that the function has a help file and that it is added to the package NAMESPACE:

# Save this file as `R/lm_stan.R`

#' Bayesian linear regression with Stan
#'
#' @export
#' @param x Numeric vector of input values.
#' @param y Numberic vector of output values.
#' @param ... Arguments passed to `rstan::sampling` (e.g. iter, chains).
#' @return An object of class `stanfit` returned by `rstan::sampling`
#'
lm_stan <- function(x, y, ...) {
  standata <- list(x = x, y = y, N = length(y))
  out <- rstan::sampling(stanmodels$lm, data = standata, ...)
  return(out)
}

When roxygen2 documentation is enabled, a top-level package file R/rstanlm-package.R is created by rstan_create_package() to import necessary functions for other packages and to set up the package for compiling Stan C++ code:

file.show(file.path("R", "rstanlm-package.R"))
#' The 'rstanlm' package.
#'
#' @description A DESCRIPTION OF THE PACKAGE
#'
#' @docType package
#' @name rstanlm-package
#' @aliases rstanlm
#' @useDynLib rstanlm, .registration = TRUE
#' @import methods
#' @import Rcpp
#' @importFrom rstan sampling
#'
#' @references
#' Stan Development Team (2020). RStan: the R interface to Stan. R package version 2.19.3. https://mc-stan.org
#'
NULL

The #' @description section can be manually edited to provided specific information about the package.

Documentation

With roxygen documentation enabled, we need to generate the documentation for lm_stan and update the NAMESPACE so the function is exported, i.e., available to users when the package is installed. This can be done with the function roxygen2::roxygenize().

example(source) # defines the sourceDir() function
try(roxygen2::roxygenize(load_code = sourceDir), silent = TRUE)
pkgbuild::compile_dll()
roxygen2::roxygenize()
Writing NAMESPACE
Re-compiling rstanlm
  
─  installing *source* package ‘rstanlm’ ...

  
   ** using staged installation

  
   ** libs

  
   

  
   

  
   Error in get(name, envir = asNamespace(pkg), inherits = FALSE) : 
     object 'CxxFlags' not found
   Calls: ::: -> get
   Execution halted

  
   clang++ -arch x86_64 -ftemplate-depth-256 -std=gnu++14 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"../inst/include" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include/src" -DBOOST_DISABLE_ASSERTS -DEIGEN_NO_DEBUG -DBOOST_MATH_OVERFLOW_ERROR_POLICY=errno_on_error -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/Rcpp/include" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppParallel/include" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/rstan/include" -isysroot /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk -I/usr/local/include   -fPIC  -O3 -mtune=native -march=native -Wno-unused-variable -Wno-unused-function  -Wno-macro-redefined -Wno-unknown-pragmas -Wno-unneeded-internal-declaration -c RcppExports.cpp -o RcppExports.o

  
   

  
   

  
   Error in get(name, envir = asNamespace(pkg), inherits = FALSE) : 
     object 'CxxFlags' not found
   Calls: ::: -> get
   Execution halted

  
   clang++ -arch x86_64 -ftemplate-depth-256 -std=gnu++14 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"../inst/include" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include/src" -DBOOST_DISABLE_ASSERTS -DEIGEN_NO_DEBUG -DBOOST_MATH_OVERFLOW_ERROR_POLICY=errno_on_error -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/Rcpp/include" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppParallel/include" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/rstan/include" -isysroot /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk -I/usr/local/include   -fPIC  -O3 -mtune=native -march=native -Wno-unused-variable -Wno-unused-function  -Wno-macro-redefined -Wno-unknown-pragmas -Wno-unneeded-internal-declaration -c stanExports_lm.cc -o stanExports_lm.o

  
   Error in get(name, envir = asNamespace(pkg), inherits = FALSE) : 
     object 'LdFlags' not found
   Calls: ::: -> get
   Execution halted

  
   Error in get(name, envir = asNamespace(pkg), inherits = FALSE) : 
     object 'LdFlags' not found
   Calls: ::: -> get
   Execution halted

  
   clang++ -arch x86_64 -ftemplate-depth-256 -std=gnu++14 -dynamiclib -Wl,-headerpad_max_install_names -undefined dynamic_lookup -single_module -multiply_defined suppress -L/Library/Frameworks/R.framework/Resources/lib -L/usr/local/lib -o rstanlm.so RcppExports.o stanExports_lm.o -F/Library/Frameworks/R.framework/.. -framework R -Wl,-framework -Wl,CoreFoundation

  
   installing to /private/var/folders/h6/14xy_35x4wd2tz542dn0qhtc0000gn/T/RtmpnbnCxa/devtools_install_15b877172b32f/00LOCK-rstanlm/00new/rstanlm/libs

  
   ** checking absolute paths in shared objects and dynamic libraries

  
─  DONE (rstanlm)
Loading rstanlm
Writing NAMESPACE
Writing lm_stan.Rd
Writing rstanlm-package.Rd

Install and use

Finally, the package is ready to be installed:

# using ../rstanlm because already inside the rstanlm directory
install.packages("../rstanlm", repos = NULL, type = "source")
Installing package into '/private/var/folders/h6/14xy_35x4wd2tz542dn0qhtc0000gn/T/Rtmpgs30S2/Rinst15b611740472b'
(as 'lib' is unspecified)

It is also possible to use devtools::install(quick=FALSE) to install the package. The argument quick=FALSE is necessary if you want to recompile the Stan models. Going forward, if you only make a change to the R code or the documentation, you can set quick=TRUE to speed up the process, or use devtools::load_all().

After installation, the package can be loaded and used like any other R package:

library("rstanlm")
fit <- lm_stan(y = rnorm(10), x = rnorm(10), 
               # arguments passed to sampling
               iter = 2000, refresh = 500)

SAMPLING FOR MODEL 'lm' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 2.4e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.24 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  500 / 2000 [ 25%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1500 / 2000 [ 75%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.024495 seconds (Warm-up)
Chain 1:                0.018478 seconds (Sampling)
Chain 1:                0.042973 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'lm' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 3e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  500 / 2000 [ 25%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1500 / 2000 [ 75%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.018392 seconds (Warm-up)
Chain 2:                0.016901 seconds (Sampling)
Chain 2:                0.035293 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'lm' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 5e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  500 / 2000 [ 25%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1500 / 2000 [ 75%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.019053 seconds (Warm-up)
Chain 3:                0.019565 seconds (Sampling)
Chain 3:                0.038618 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'lm' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 8e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  500 / 2000 [ 25%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1500 / 2000 [ 75%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.019692 seconds (Warm-up)
Chain 4:                0.017254 seconds (Sampling)
Chain 4:                0.036946 seconds (Total)
Chain 4: 
print(fit)
Inference for Stan model: lm.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

           mean se_mean   sd  2.5%   25%   50%  75% 97.5% n_eff Rhat
intercept  0.35    0.01 0.26 -0.15  0.21  0.35 0.50  0.90  2029    1
beta      -0.06    0.00 0.18 -0.39 -0.16 -0.06 0.05  0.30  2112    1
sigma      0.70    0.01 0.23  0.41  0.55  0.65 0.80  1.29  1633    1
lp__      -0.95    0.04 1.51 -4.99 -1.55 -0.55 0.12  0.72  1209    1

Samples were drawn using NUTS(diag_e) at Mon Jul  6 08:41:28 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Advanced options

Details can be found in the documentation for rstan_create_package() so we only mention some of these briefly here: