MCMC inference for LMLS

A Markov chain Monte Carlo (MCMC) sampler for location-scale regression models from the lmls() function. The sampler uses Gibbs updates for the location coefficients and the Riemann manifold Metropolis-adjusted Langevin algorithm (MMALA) from Girolami and Calderhead (2011) with the Fisher-Rao metric tensor for the scale coefficients. The priors for the regression coefficients are assumed to be flat.

To find the optimal step size for the MMALA updates, the dual averaging algorithm from Nesterov (2009) is used during a warm-up phase.

Usage

mcmc(m, num_samples = 1000, num_warmup = 1000, target_accept = 0.8)

Arguments

m

A location-scale regression model from the lmls() function.

num_samples

The number of MCMC samples after the warm-up. Defaults to 1000.

num_warmup

The number of MCMC samples for the warm-up. Defaults to 1000.

target_accept

The target acceptance rate for the dual averaging algorithm used for the warm-up. Defaults to 0.8.

Value

An lmls S3 object, see lmls(). The entry mcmc with the matrices of MCMC samples is added to the object as a list with the names locationand scale.

References

Girolami, M. and Calderhead, B. (2011), Riemann manifold Langevin and Hamiltonian Monte Carlo methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73: 123-214. doi: 10.1111/j.1467-9868.2010.00765.x

Nesterov, Y. (2009), Primal-dual subgradient methods for convex problems. Mathematical Programming, 120: 221–259. doi: 10.1007/s10107-007-0149-x

Examples

library(lmls)
m <- lmls(y ~ poly(x, 2), ~ x, data = abdom, light = FALSE)
m <- mcmc(m)
summary(m, type = "mcmc")
#> 
#> Call:
#> lmls(location = y ~ poly(x, 2), scale = ~x, data = abdom, light = FALSE)
#> 
#> Deviance residuals:
#>      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
#> -3.095000 -0.654600 -0.008431  0.029320  0.636100  4.016000 
#> 
#> Location coefficients (identity link):
#>                Mean    2.5%     50%   97.5%
#> (Intercept)  226.70  225.66  226.70  227.80
#> poly(x, 2)1 2159.39 2131.80 2159.41 2188.45
#> poly(x, 2)2  -99.58 -122.48  -99.32  -76.61
#> 
#> Scale coefficients (log link):
#>                Mean    2.5%     50% 97.5%
#> (Intercept) 1.37052 1.18277 1.37071 1.576
#> x           0.04192 0.03507 0.04204 0.049
#> 
#> Residual degrees of freedom: 605
#> Log-likelihood: -2398.5
#> AIC: 4807.01
#> BIC: 4829.07
#> 
plot(m$mcmc$scale[, 2], type = "l")