This function fits Tweedie compound Poisson linear models using Markov Chain Monte Carlo methods.

1 2 3 4 5 6 7 8 9 | ```
bcplm(formula, link = "log", data, inits = NULL,
weights, offset, subset, na.action, contrasts = NULL,
n.chains = 3, n.iter = 2000, n.burnin = floor(n.iter / 2),
n.thin = max(1, floor(n.chains * (n.iter - n.burnin) / n.sims)),
n.sims = 1000, n.report = 2, prior.beta.mean = NULL,
prior.beta.var = NULL, bound.phi = 100, bound.p = c(1.01, 1.99),
tune.iter = 5000, n.tune = floor(tune.iter/100),
basisGenerators = c("tp", "bsp", "sp2d"), doFit = TRUE, ...)
``` |

`formula` |
an object of class |

`link` |
a specification for the model link function. This can be either a literal character string or a numeric number. If it is a character string, it must be one of "log", "identity", "sqrt" or "inverse". If it is numeric, it is the same as the |

`inits` |
a list of initial values to be used for each chain. It must be of length |

`data, subset, weights, na.action, offset, contrasts` |
further model specification arguments as in |

`n.chains` |
an integer indicating the number of Markov chains (default: |

`n.iter` |
the number of total iterations per chain (including burn in; default: |

`n.burnin` |
the length of burn in, i.e. number of iterations to discard at the beginning. Default
is |

`n.thin` |
thinning rate. Must be a positive integer. Set |

`n.sims` |
The approximate number of simulations to keep after thinning (all chains combined). |

`n.report` |
if greater than zero, fitting information will be printed out |

`prior.beta.mean` |
a vector of prior means for the fixed effects. Default is a vector of zeros. |

`prior.beta.var` |
a vector of prior variances for the fixed effects. Default is a vector of |

`bound.phi` |
a numeric value indicating the upper bound of the uniform prior for the dispersion parameter. The default is |

`bound.p` |
a vector of lower and upper bounds for the index parameter |

`tune.iter` |
the number of iterations used for tuning the proposal variances used in the Metropolis-Hastings updates. These iterations will not be included in the final output. Default is |

`n.tune` |
a positive integer (default: |

`basisGenerators` |
a character vector of names of functions that generate spline bases. See |

`doFit` |
if |

`...` |
not used. |

This function provides Markov chain Monte Carlo [MCMC] methods for fitting Tweedie compound Poisson linear models within the Bayesian framework. Both generalized linear models and mixed models can be handled. In computing the posterior distribution, the series evaluation method (see, e.g., `dtweedie`

) is employed to evaluate the compound Poisson density.

In the Bayesian model, prior distributions have to be specified for all parameters in the model. Here, Normal distributions are used for the fixed effects (*β*), a Uniform distribution for the dispersion parameter (*φ*), a Uniform distribution for the index parameter (*p*). If a mixed model is specified, prior distributions must be specified for the variance component. If there is one random effect in a group, the inverse Gamma (scale = `0.001`

, shape = `0.001`

) is specified as the prior. If there is more than one random effects in a group, the inverse Wishart (identity matrix as the scale and the dimension of the covariance matrix as the shape) is specified as the prior.

Prior means and variances of the fixed effects can be supplied using the argument `prior.beta.mean`

and `prior.beta.var`

, respectively. The prior distribution of *φ* is uniform on (0, `bound.phi`

). And the bounds of the Uniform for *p* can be specified in the argument `bound.p`

. See details in section 'Arguments'.

In implementing the MCMC, a Gibbs sampler is constructed in which parameters are updated one at a time given the current values of all the other parameters. Specifically, we use the random-walk Metropolis-Hastings algorithm in updating each parameter except for the variance components, which can be simulated directly due to conjugacy.

Before the MCMC, there is a tuning process where the proposal variances of the (truncated) Normal proposal distributions are updated according to the sample variances computed from the simulations in each tuning loop. The goal is to make the acceptance rate roughly between 40% and 60% for univariate M-H updates. The argument `tune.iter`

determines how many iterations are used for the tuning process, and `n.tune`

determines how many loops these iterations should be divided into. These iterations will not be used in the final output.

The simulated values of all model parameters are stored in the `sims.list`

slot of the returned `bcplm`

object. It is a list of `n.chains`

matrices and each matrix has approximately `n.sims`

rows. The `sims.list`

slot is further coerced to be of class `"mcmc.list"`

so that various methods from the `coda`

package can be directly applied to get Markov chain diagnostics, posterior summary and plots. See `coda`

for available methods.

`bcplm`

returns an object of class `"bcplm"`

. See `bcplm-class`

for details of the return values as well as various methods available for this class.

Yanwei (Wayne) Zhang actuary_zhang@hotmail.com

Zhang, Y (2013). Likelihood-based and Bayesian Methods for Tweedie Compound Poisson Linear Mixed Models, *Statistics and Computing*, 23, 743-757.

The users are recommended to see the documentation for `bcplm-class`

, `cpglm`

, `cpglmm`

, `mcmc`

, and `tweedie`

for related information.

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 | ```
## Not run:
# fit the FineRoot data with Bayesian models
# Bayesian cpglm
set.seed(10)
fit1 <- bcplm(RLD ~ factor(Zone) * factor(Stock),
data = FineRoot, tune.iter = 2000,
n.iter = 6000, n.burnin = 1000, n.thin = 5)
gelman.diag(fit1$sims.list)
# diagnostic plots
acfplot(fit1$sims.list, lag.max = 20)
xyplot(fit1$sims.list)
densityplot(fit1$sims.list)
summary(fit1)
plot(fit1)
# now fit the Bayesian model to an insurance loss triangle
# (see Peters et al. 2009)
fit2 <- bcplm(increLoss ~ factor(year) + factor(lag),
data = ClaimTriangle, n.iter = 12000,
n.burnin = 2000, n.thin = 10, bound.p = c(1.1, 1.95))
gelman.diag(fit2$sims.list)
summary(fit2)
# mixed models
set.seed(10)
fit3 <- bcplm(RLD ~ Stock * Zone + (1|Plant),
data = FineRoot, n.iter = 15000,
n.burnin = 5000, n.thin = 10)
gelman.diag(fit3$sims.list)
summary(fit3)
## End(Not run)
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.