bdw {BDWreg}

Bayesian parameter estimation for discrete Weibull regression

Description

Bayesian estimation of the parameters for discrete Weibull regression. The conditional distribution of the response given the predictors is assumed to be DW with parameters q and beta dependent on the predictors.

Usage

bdw(data, formula = NA, reg.q = FALSE, reg.b = FALSE,
  logit = TRUE, initial = c(.5,.5), iteration = 25000,
  v.scale = 0.1,RJ = FALSE,  dist.q = imp.d, dist.b = imp.d,
  q.par = c(0, 0), b.par = c(0, 0), penalized = FALSE,
  dist.l = imp.d, l.par = c(0, 0), bi.period = 0.25, cov = 1,
  sampling = c("bin"), est = Mode, fixed.l = -1,jeffrey = FALSE)

Arguments

data

A data frame containing the variables in the model. If not found in data, the variables are taken from the environment (formula).

formula

An object of class "formula". A symbolic description of the model to be fitted.

reg.q

Logical flag. If TRUE, the model includes a dependency of q on the predictors x, as explained in 'logit' argument.

reg.b

Logical flag. If TRUE, the model includes a dependency of beta on the predictors x, given by:

log(\beta)= \gamma_0+\gamma_1X_1+\ldots+\gamma_pX_p.

If FALSE, beta(x) is constant.

logit

Logical flag. If TRUE, the model includes a dependency of q on predictors via a logit transformation as explained below.

The conditional distribution of Y (response) given x (predictors) is assumed a DW(q(x),beta(x)).

If logit=TRUE & (reg.q=TRUE)

log(q/(1-q))=\theta_0+\theta_1 X_1+\ldots+\theta_pX_p.

If logit=FALSE & (reg.q=TRUE)

log(-log(q))=\theta_0+\theta_1 X_1+\ldots+\theta_pX_p.

initial

Vector of initial values for parameters. In all cases DW(q,B), DW(regQ,B), DW(q,rebB) and DW(regQ,regB) the first parameters correspond to q or corresponding regression coefficients and next is beta or corresponding regression coefficients. If penalized=TRUE then an extra value must be added for tuning parameter.

iteration

Number of MCMC iterations.

v.scale

The scale of the proposal function. Setting to lower values results in an increase in the acceptance rate of the sampler.

RJ

Logical flag. If TRUE, Reversible-Jump sampling is used to draw samples from the posterior. Otherwise, a Metropolis-Hastings sampler applies.

dist.q

Density function. Prior density for theta. If not set, an improper prior applies. Any two parameter density function is allowed, e.g. dnorm, dlaplace, dunif etc. Any customized density function must support log=TRUE flag.

dist.b

Density function. Prior density for gamma. If not set, an improper prior applies. Any two parameter density function is allowed, e.g. dnorm, dlaplace, dunif etc. Any customized density function must support log=TRUE flag.

q.par

A vector of length two corresponding to the parameters of 'dist.q'. The default is set to c(0,0).

b.par

A vector of length two corresponding to the parameters of 'dist.b'. The default is set to c(0,0).

penalized

logical flag. If TRUE, an hyper-parameter inducing shrinkage is considered. In this case, prior must be set to an informative prior, e.g. Gaussian, Laplace. See also 'l.par' and 'dist.l' below.

dist.l

Density function. Hyper prior for penalty term. If not set, an improper prior is used. Any non-negative two parameter density function is allowed, e.g. dgamma, dbeta,... Any customized density function must support log=TRUE flag.

l.par

A vector of length two corresponding to the parameters of the hyper-prior 'dist.l'. The default is set to c(0,0).

bi.period

A numeric value in (0,1) indicating the burn-in period of the MCMC chain. The default is set to 0.25 meaning that 25% of values remove from the beginning of the output chain. See 'sampling'.

cov

A value in {1,2,3,4}. If set to 1 then the adaptive-MH is performed; 2: an independent uniform proposal; 3: an independent Laplace proposal and 4: an independent Gaussian proposal. The default is 1. If cov=1 then the scale of the proposal is adapted from the data.

sampling

Choose between independent (indp), systematic (syst) and burn-in (bin). If set to indp then the chain is not ordered! The default is 'bin'. Sampling interval for Systematic sampling is calculated from iteration(1-bi.period). Similarly for indp the number of samples is computed from iteration(1-bi.period).

est

Statistic. The statistic that is used to estimate the parameters from marginal densities. Possible values are: mode, mean, median or any customized univariate measure of location. The default is 'mode'.

fixed.l

A positive number. Set to a positive value corresponding to a parameter that does not need estimation. Set to any negative value to disable this option.

jeffrey

A logical flag. Set to a TRUE to use Jeffrey prior. Notice that MCMC based on Jeffrey can take considerably long time and the results in the most cases are worse than using an improper flat (=1) prior.

Value

res

Estimation of the parameters from marginal densities using the statistic specified in 'est'.

chain

A list of values including sampler configurations.

chain

Including estimation of marginal densities

acceptance.rate

Acceptance rate of sampler

RejAccChain

A zero-one vector reporting rejection (0) and acceptance (1) of the samples

error

Total number of errors during the sampling procedure

minf

Minimum loglikelihood among all iterations

minState

Coefficients corresponding to minimum loglikelihood

lb

The number of gamma parameters

lq

The number of theta parameters

model.chain

If RJ=TRUE then this chain contains acceptance (+values) and rejection (-values) of models

duration

Time elapsed until completion of the sampling procedure

Author(s)

Hamed Haselimashhadi <[e-mail]>

References

Haselimashhadi, Vinciotti and Yu (2015), A new Bayesian regression model for counts in medicine.

See Also

plot.bdw, summary.bdw, bdw.mc

Examples

set.seed(123)
#========== example 1 - estimating DW parameters under logit transformation ==========
q = .41    # <<< true parameters
b = 1.1    # <<< true parameters
y = BDWreg:::rdw(n = 200,q = q,beta = b) #<<< generating data
result  = bdw(data = y, v.scale = .10,initial = c(.5,.5),iteration = 8000 )
plot(result)
summary(result)


## Not run: 
  #==== example 2 - estimating logit-DW(regQ,beta) parameters using RJ ======
  set.seed(1234)
  n = 500
  x1 = runif(n = n, min = 0, max = 1.5)
  x2 = runif(n = n, min = 0, max = 1.5)

  theta0 = .6   #<<< true parameter
  theta1 =  0   #<<< true parameter
  theta2 = .34  #<<< true parameter

  lq = theta0 +  x1*theta1 + x2*theta2

  q = exp(lq - log(1+exp(lq)) )
  beta = 1.5

  y = c()
  for(i in 1:n){
    y[i] = BDWreg:::rdw(1,q = q[i],beta = beta)
  }

  data = data.frame(x1,x2,y) # <<<- data
  result2 = bdw(data = data              ,
               formula = y~.             ,
               RJ = TRUE                 ,
               initial = rep(.5,4)       ,
               iteration = 25000         ,
               reg.b = FALSE,reg.q = TRUE,
               v.scale = .1              ,
               q.par = c(0,1)            ,
               b.par = c(0,1)            ,
               dist.q = dnorm            ,
               dist.b = dnorm
  )
  plot(result2)
  summary(result2)

## End(Not run)