Description Usage Arguments Details Value Author(s) See Also Examples
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.
1 2 3 |
flist |
A formula or list of formulas that define the likelihood and priors. Can also pass in a |
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 |
iter |
Number of iterations of sampling. By default, half of these iterations are warmup. |
warmup |
Number of warmup iterations. By default, half of |
chains |
Number of independent chains to sample from |
debug |
If |
verbose |
If |
WAIC |
When |
cores |
Number of processor cores to distribute chains over, using |
... |
Additional arguments to pass to |
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
.
Returns an object of class map2stan
with the following slots.
call |
The function call |
model |
Stan model code |
stanfit |
|
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. |
Richard McElreath
resample
, map
, stan
, link
, sim
, glimmer
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 | ## 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] )
## End(Not run)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.