rethinking (version 2.13)

quap: Compute quadratic approximate posterior distribution


Find mode of posterior distribution for arbitrary fixed effect models and then produce an approximation of the full posterior using the quadratic curvature at the mode.


quap( flist , data , start , method="BFGS" , hessian=TRUE , debug=FALSE , verbose=FALSE , ... )
map( flist , data , start , method="BFGS" , hessian=TRUE , debug=FALSE , verbose=FALSE , ... )



A formula or alist of formulas that define the likelihood and priors. See details.


A data frame or list containing the data


Optional named list specifying parameters and their initial values


Search method for optim. Defaults to BFGS.


If TRUE, attempts to compute the Hessian


If TRUE, prints various internal steps to help with debugging


Additional arguments to pass to optim


Returns an object of class map with the following slots.


The function call


The estimated posterior modes


Variance-covariance matrix


List returned by optim


The data


Formula list


Minus log-likelihood function that accepts a vector of parameter values


This command provides a convenient interface for finding quadratic approximations of posterior distributions for models defined by a formula or by a list of formulas. This procedure is equivalent to penalized maximum likelihood estimation and the use of a Hessian for profiling, and therefore can be used to define many common regularization procedures. The point estimates returned correspond to a maximum a posterior, or MAP, estimate. However the intention is that users will use extract.samples and other functions to work with the full posterior.

flist should be a either a single formula that defines the likelihood function or rather a list of formulas that define the likelihood and priors for parameters. See examples below.

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.

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 or mu <- invlink(a + b*x). The names "link" and "invlink" must be recognized by map. It currently recognizes log/exp and logit/logistic. Any other link function can be coded directly into the likelihood formula. For example y ~ dfoo(invlink(mu)).

The start list is optional. For each parameter with a defined prior, an initial value will be sampled from the prior. Explicit named initial values can also be provided in this list.

Methods are defined for coef, summary, logLik, vcov, nobs, deviance, link, DIC, and show.

See Also



Run this code
flist <- alist(
    dist ~ dnorm( mu , sigma ) ,
    mu <- a+b*speed ,
    c(a,b) ~ dnorm(0,1) , 
    sigma ~ dexp(1)
fit <- quap( flist , start=list(a=40,b=0.1,sigma=20) , data=cars )

# regularized logistic regression example
y <- c( rep(0,10) , rep(1,10) )
x <- c( rep(-1,9) , rep(1,11) )

flist0 <- alist(
    y ~ dbinom( 1 , p ) ,
    logit(p) <- a + b*x

flist1 <- alist(
    y ~ dbinom( 1 , p ),
    logit(p) <- a + b*x ,
    c(a,b) ~ dnorm(0,10)

# without priors, same problem as:
# glm3a <- glm( y ~ x , family=binomial , data=list(y=y,x=x) )
fit3a <- quap( flist0 , data=list(y=y,x=x) , start=list(a=0,b=0) )

# now with regularization
fit3b <- quap( flist1 , data=list(y=y,x=x) , start=list(a=0,b=0) )

# vector parameters
fit4 <- quap(
        admit ~ dbinom( apps , p ),
        logit(p) <- a[dept_id] + b*male,
        a[dept_id] ~ dnorm(0,4),
        b ~ dnorm(0,1)
        admit = UCBadmit$admit,
        apps = UCBadmit$applications,
        male = ifelse( UCBadmit$applicant.gender=="male" , 1 , 0 ),
        dept_id = coerce_index(UCBadmit$dept)

# }

Run the code above in your browser using DataCamp Workspace