rethinking (version 1.59)

map2stan: Build RStan models from formulas

Description

Compiles lists of formulas, like those used in map, into Stan model code. Allows for arbitary fixed effect and mixed effect regressions. Computes DIC and WAIC. Allows for simple imputation of missing values.

Usage

map2stan( flist , data , start , pars , constraints=list() , types=list() , 
  sample=TRUE , iter=2000 , warmup=floor(iter/2) , chains=1 , debug=FALSE , 
  verbose=FALSE , WAIC=TRUE , cores=1 , rng_seed , ... )

Arguments

flist

A formula or list of formulas that define the likelihood and priors. Can also pass in a map model fit. See details.

data

A data frame or list containing the data

start

Optional named list specifying parameters and their initial values

pars

Optional: character vector of parameters to return samples for

constraints

Optional: named list of custom parameter constraints, using Stan notation

types

Optional: named list of custom parameter types, using Stan notation

sample

If FALSE, builds Stan code without sampling

iter

Number of iterations of sampling. By default, half of these iterations are warmup.

warmup

Number of warmup iterations. By default, half of iter.

chains

Number of independent chains to sample from

debug

If TRUE, prints various internal steps to help with debugging

verbose

If TRUE, prints extra progress messages.

WAIC

When TRUE, computes WAIC after sampling, storing the result

cores

Number of processor cores to distribute chains over, using parallel.

...

Additional arguments to pass to stan

Value

Returns an object of class map2stan with the following slots.

call

The function call

model

Stan model code

stanfit

stanfit object returned by stan

coef

The posterior means

vcov

Minimal variance-covariance matrix, just holding diagonal variances

data

The data

start

List of starting values that were used in sampling

pars

Parameter names monitored in samples

formula

Formula list from call

formula_parsed

List of parsed formula information. Useful mainly for debugging.

Details

This command provides a convenient interface for building arbitary fixed effect and mixed effect generalized linear models, as defined by a list of formulas. Syntax is similar to map, but also allowing multivariate priors corresponding to varying (aka random) effects, as well as simple imputation schemes.

flist should be either (1) a single formula that defines the likelihood function or rather a list of formulas that define the likelihood and linear models and priors for parameters (see examples below) or (2) a previously fit map model.

Likelihood formulas take the form y ~ dfoo(bar), where y is the outcome variable, dfoo is a density function such as dnorm, and bar is a parameter of the density.

Prior formulas take the same form, but the outcome should be a parameter name. Identical priors can be defined for multiple parameters by using c(par1,par2,...) on the left hand side of the formula. See example below.

A special case of prior formula is for varying effects. For single varying effects, such as varying intercepts alone, all that is needed is to define a prior and mark it as conditional on a grouping variable in the data. For example: aj[id] ~ dnorm(0,sigma_id) specifies a vector of varying effects aj, one for each unique value in id. For correlated varying effects, such as both varying intercepts and slopes, a parameter vector is specified and a multivariate prior is used instead. For example: c(aj,bj)[id] ~ dmvnorm(0,Sigma_id) specifices varying intercepts aj and varying slopes bj.

Linear models can be specified as formulas of the form mu <- a + b*x for a direct link. To use a link function, use the form link(mu) <- a + b*x. The name "link" must be recognized by map2stan. It currently recognizes log and logit.

Imputation of missing values is available by specifying distributions for predictor variables that contain NA values. map2stan will split the variable into observed values and a vector of parameters used to estimate the missing values, assigning the same distribution to each. See the example.

The start list is optional. When missing from the list, for each parameter with a defined prior, an initial value will be sampled from the prior. Sampled initial values will try to respect parameter constraints. For varying effect parameter vectors, initial values will always be set to zero. Specific initial values can be specified in the start list. See examples below.

The Stan model code includes a generated quantities block that computes the deviance for each iteration of parameter samples. When sampling completes, map2stan computes DIC, the deviance information criterion, from the samples. DIC information is available from show and DIC, as well as being attributes of the returned object.

WAIC can be computed with WAIC, or by setting WAIC=TRUE when calling map2stan. This is currently the default. WAIC is calculated entirely after Stan completes sampling.

Methods are defined for extract.samples, link, sim, ensemble, compare, coef, summary, logLik, vcov, nobs, deviance, plot, pairs, and show.

See Also

resample, map, stan, link, sim, glimmer

Examples

# NOT RUN {
library(rethinking)
data(chimpanzees)

# don't want any variables with NAs
d <- list( 
    pulled_left = chimpanzees$pulled_left ,
    prosoc_left = chimpanzees$prosoc_left ,
    condition = chimpanzees$condition ,
    actor = as.integer( chimpanzees$actor ) ,
    blockid = as.integer( chimpanzees$block )
)

# RStan fit
m2 <- map2stan(
    alist(
        pulled_left ~ dbinom(1,theta),
        logit(theta) <- a + bp*prosoc_left + bpc*condition*prosoc_left ,
        a ~ dnorm(0,10),
        bp ~ dnorm(0,10),
        bpc ~ dnorm(0,10)
    ) ,
    data=d, chains=2, cores=1 )

precis(m2)
summary(m2)
plot(m2)
pairs(m2)

# now RStan fit of model with varying intercepts on actor
m3 <- map2stan(
    alist(
        pulled_left ~ dbinom(1,theta),
        logit(theta) <- a + aj[actor] + bp*prosoc_left + bpc*condition*prosoc_left,
        aj[actor] ~ dnorm( 0 , sigma_actor ),
        a ~ dnorm(0,10),
        bp ~ dnorm(0,10),
        bpc ~ dnorm(0,10),
        sigma_actor ~ dcauchy(0,1)
    ) ,
    data=d,
    iter=5000 , warmup=1000 , chains=2 , cores=1 )

precis(m3)
plot(m3)
pairs(m3)

# varying intercepts on actor and experimental block
m4 <- map2stan(
    alist(
        pulled_left ~ dbinom(1,theta),
        logit(theta) <- a + aj + ak + bp*prosoc_left + bpc*condition*prosoc_left,
        aj[actor] ~ dnorm( 0 , sigma_actor ),
        ak[blockid] ~ dnorm( 0 , sigma_block ),
        a ~ dnorm(0,10),
        bp ~ dnorm(0,10),
        bpc ~ dnorm(0,10),
        sigma_actor ~ dcauchy(0,1),
        sigma_block ~ dcauchy(0,1)
    ) ,
    data=d,
    iter=5000 , warmup=1000 , chains=2 , cores=1 )

precis(m4)
summary(m4)
plot(m4)

# compare posterior means
coeftab(m2,m3,m4)
plot(coeftab(m2,m3,m4))

# show WAIC for m2,m3,m4
compare(m2,m3,m4)
plot(compare(m2,m3,m4))

###########
# varying slopes models

# varying slopes on actor
# also demonstrates use of multiple linear models
# see Chapter 13 for discussion
m5 <- map2stan(
    alist(
        # likeliood
        pulled_left ~ dbinom(1,p),

        # linear models
        logit(p) <- A + (BP + BPC*condition)*prosoc_left,
        A <- a + a_actor[actor],
        BP <- bp + bp_actor[actor],
        BPC <- bpc + bpc_actor[actor],

        # adaptive prior
        c(a_actor,bp_actor,bpc_actor)[actor] ~
                                dmvnorm2(0,sigma_actor,Rho_actor),

        # fixed priors
        c(a,bp,bpc) ~ dnorm(0,1),
        sigma_actor ~ dcauchy(0,2),
        Rho_actor ~ dlkjcorr(4)
    ) , data=d , iter=5000 , warmup=1000 , chains=3 , cores=3 )

# same model but with non-centered parameterization
# see Chapter 13 for explanation and more elaborate example

m6 <- map2stan(
    alist(
        # likeliood
        pulled_left ~ dbinom(1,p),

        # linear models
        logit(p) <- A + (BP + BPC*condition)*prosoc_left,
        A <- a + a_actor[actor],
        BP <- bp + bp_actor[actor],
        BPC <- bpc + bpc_actor[actor],

        # adaptive prior - non-centered
        c(a_actor,bp_actor,bpc_actor)[actor] ~
                                dmvnormNC(sigma_actor,Rho_actor),

        # fixed priors
        c(a,bp,bpc) ~ dnorm(0,1),
        sigma_actor ~ dcauchy(0,2),
        Rho_actor ~ dlkjcorr(4)
    ) , data=d , iter=5000 , warmup=1000 , chains=3 , cores=3 )

###########
# Imputation example

# simulate data:
#  linear regression with two predictors
#  both predictors have valules missing at random
N <- 100
N_miss <- 10
x1 <- rnorm( N )
x2 <- rnorm( N )
y <- rnorm( N , 2*x1 - 0.5*x2 , 1 )
x1[ sample(1:N,size=N_miss) ] <- NA
x2[ sample(1:N,size=N_miss) ] <- NA

# formula with distributions assigned to both predictors
f <- alist(
    y ~ dnorm( mu , sigma ),
    mu <- a + b1*x1 + b2*x2,
    x1 ~ dnorm( mu_x1, sigma_x1 ),
    x2 ~ dnorm( mu_x2, sigma_x2 ),
    a ~ dnorm( 0 , 100 ),
    c(b1,b2) ~ dnorm( 0  , 10 ),
    c(mu_x1,mu_x2) ~ dnorm( 0 , 100 ),
    c(sigma_x1,sigma_x2) ~ dcauchy(0,2),
    sigma ~ dcauchy(0,2)
)

m <- map2stan( f , data=list(y=y,x1=x1,x2=x2) , sample=TRUE )

# show observed outcomes against retrodicted outcomes
# cases with missing values shown with red posterior intervals
v <- link(m)
mu <- apply( v , 2 , mean )
ci <- apply( v , 2 , PI )
plot( y ~ mu )
cicols <- ifelse( is.na(x1) | is.na(x2) , "red" , "gray" )
for( i in 1:N ) lines( ci[,i] , rep(y[i],2) , col=cicols[i] )
# }