bdw.reg {BDgraph}

Bayesian estimation of (zero-inflated) Discrete Weibull regression

Description

Bayesian estimation of the parameters for Discrete Weibull (DW) regression. The conditional distribution of the response given the predictors is assumed to be DW with parameters q and beta, dependent on the predictors, and, with an additional parameter pi under zero inflation.

Usage

bdw.reg( data, formula = NA, iter = 5000, burnin = NULL, 
         dist.q = dnorm, dist.beta = dnorm, 
         par.q = c( 0, 1 ), par.beta = c( 0, 1 ), par.pi = c( 1, 1 ), 
         initial.q = NULL, initial.beta = NULL, initial.pi = NULL, 
         ZI = FALSE, scale.proposal = NULL, adapt = TRUE, print = TRUE )

Arguments

data

data.frame or matrix corresponding to the data, containing the variables in the model.

formula

object of class formula as a symbolic description of the model to be fitted. For the case of data.frame, it is taken as the model frame (see model.frame).

iter

number of iterations for the sampling algorithm.

burnin

number of burn-in iterations for the sampling algorithm.

dist.q

Prior density for the regression coefficients associated to the parameter q. The default is a Normal distribution (dnorm). Any density function which has two parameters and can support the log = TRUE flag can be used, e.g. dnorm, dlnorm, dunif etc.

dist.beta

Prior density for the regression coefficients associated to the parameter beta. The default is a Normal distribution (dnorm). Any density function which has two parameters and can support the log = TRUE flag can be used, e.g. dnorm, dlnorm, dunif etc.

par.q

vector of length two corresponding to the parameters of dist.q.

par.beta

vector of length two corresponding to the parameters of dist.beta.

par.pi

vector of length two corresponding to the parameters of the beta prior density on pi.

initial.q, initial.beta, initial.pi

vector of initial values for the regression coefficients and for pi (if ZI = TRUE).

ZI

logical: if FALSE (default), the conditional distribution of the response given the predictors is assumed to be DW with parameters q and beta. If TRUE, a zero-inflated DW distribution will be applied.

scale.proposal

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

adapt

logical: if TRUE (default), the proposals will be adapted. If FALSE, no adapting will be applied.

print

logical: if TRUE (default), tracing information is printed.

Details

The regression model uses a logit link function on q and a log link function on beta, the two parameters of a DW distribution, with probability mass function given by

DW(y) = q^{y^\beta} - q^{(y+1)^\beta}, y = 0, 1, 2, \ldots

For the case of zero inflation (ZI = TRUE), a zero-inflated DW is considered:

f(y) = ( 1 - pi) I(y = 0 ) + pi DW(y)

where 0 \leq pi \leq 1 and I(y = 0 ) is an indicator for the point mass at zero for the response y.

Value

sample

MCMC samples

q.est

posterior estimates of q

beta.est

posterior estimates of beta

pi.est

posterior estimates of pi

accept.rate

acceptance rate of the MCMC algorithm

Author(s)

Veronica Vinciotti, Reza Mohammadi [e-mail], and Pariya Behrouzi

References

Vinciotti, V., Behrouzi, P., and Mohammadi, R. (2022) Bayesian structural learning of microbiota systems from count metagenomic data, arXiv preprint, doi:10.48550/arXiv.2203.10118

Peluso, A., Vinciotti, V., and Yu, K. (2018) Discrete Weibull generalized additive model: an application to count fertility, Journal of the Royal Statistical Society: Series C, 68(3):565-583, doi:10.1111/rssc.12311

Haselimashhadi, H., Vinciotti, V. and Yu, K. (2018) A novel Bayesian regression model for counts with an application to health data, Journal of Applied Statistics, 45(6):1085-1105, doi:10.1080/02664763.2017.1342782

See Also

bdgraph.dw, bdgraph, ddweibull, bdgraph.sim

Examples

## Not run: 
# - - Example 1

q    = 0.6
beta = 1.1
n    = 500

y = BDgraph::rdweibull( n = n, q = q, beta = beta )

output = bdw.reg( data = y, y ~ ., iter = 5000 )

output $ q.est
output $ beta.est

traceplot( output $ sample[ , 1 ], acf = T, pacf = T )
traceplot( output $ sample[ , 2 ], acf = T, pacf = T )

# - - Example 2

q    = 0.6
beta = 1.1
pii  = 0.8
n    = 500

y_dw = BDgraph::rdweibull( n = n, q = q, beta = beta )
z = rbinom( n = n, size = 1, prob = pii ) 
y = z * y_dw

output = bdw.reg( data = y, iter = 5000, ZI = TRUE )

output $ q.est
output $ beta.est
output $ pi.est

traceplot( output $ sample[ , 1 ], acf = T, pacf = T )
traceplot( output $ sample[ , 2 ], acf = T, pacf = T )
traceplot( output $ sample[ , 3 ], acf = T, pacf = T )

# - - Example 3

theta.q    = c( 0.1, -0.1, 0.34 )  # true parameter
theta.beta = c( 0.1, -.15, 0.5  )  # true parameter

n  = 500

x1 = runif( n = n, min = 0, max = 1.5 )
x2 = runif( n = n, min = 0, max = 1.5 )

reg_q = theta.q[ 1 ] + x1 * theta.q[ 2 ] + x2 * theta.q[ 3 ]
q     = 1 / ( 1 + exp( - reg_q ) )

reg_beta = theta.beta[ 1 ] + x1 * theta.beta[ 2 ] + x2 * theta.beta[ 3 ]
beta     = exp( reg_beta )

y = BDgraph::rdweibull( n = n, q = q, beta = beta )

data = data.frame( x1, x2, y ) 

output = bdw.reg( data, y ~. , iter = 5000 )

# - - Example 4

theta.q    = c( 1, -1, 0.8 )  # true parameter
theta.beta = c( 1, -1, 0.3 )  # true parameter
pii = 0.8

n  = 500

x1 = runif( n = n, min = 0, max = 1.5 )
x2 = runif( n = n, min = 0, max = 1.5 )

reg_q = theta.q[ 1 ] + x1 * theta.q[ 2 ] + x2 * theta.q[ 3 ]
q     = 1 / ( 1 + exp( - reg_q ) )

reg_beta = theta.beta[ 1 ] + x1 * theta.beta[ 2 ] + x2 * theta.beta[ 3 ]
beta     = exp( reg_beta )

y_dw = BDgraph::rdweibull( n = n, q = q, beta = beta )
z    = rbinom( n = n, size = 1, prob = pii ) 
y    = z * y_dw

data = data.frame( x1, x2, y ) 

output = bdw.reg( data, y ~. , iter = 5000 )

## End(Not run)