Variance function of cobin

Description

B''(B'^{-1}(\mu))

Usage

Vft(mu)

Arguments

Value

B''(B'^{-1}(\mu))

Cumulant (log partition) function of cobin

Description

B(x) = \log\{(\exp(x)-1)/x)\}

Usage

bft(x)

Arguments

Value

B(x) = \log\{(\exp(x)-1)/x)\}

First derivative of cobin cumulant (log partition) function

Description

B'(x) = 1/(1-\exp(-x))-1/x. When g is canonical link of cobin, this is same as g^{-1}

Usage

bftprime(x)

cobitlinkinv(x)

Arguments

Value

B'(x) = 1/(1-\exp(-x))-1/x.

Inverse of first derivative of cobin cumulant (log partition) function

Description

Calculates (B')^{-1}(y) using numerical inversion (Newton-Raphson), where B'(x) = 1/(1-\exp(-x))-1/x. This is the cobit link function g, the canonical link function of cobin.

Usage

bftprimeinv(y, x0 = 0, tol = 1e-08, maxiter = 100)

cobitlink(y, x0 = 0, tol = 1e-08, maxiter = 100)

Arguments

Value

(B')^{-1}(y)

Second derivative of cobin cumulant (log partition) function

Description

B''(x) = 1/x^2 + 1/(2-2*\cosh(x)) used Taylor series expansion for x near 0 for stability

Usage

bftprimeprime(x)

Arguments

Value

B''(x) = 1/x^2 + 1/(2-2*\cosh(x))

Third derivative of cobin cumulant (log partition) function

Description

B'''(x) = 1/(4*\tanh(x/2)*\sinh(x/2)^2) - 2/x^3 used Taylor series expansion for x near 0 for stability

Usage

bftprimeprimeprime(x)

Arguments

Value

B'''(x) = 1/(4*\tanh(x/2)*\sinh(x/2)^2) - 2/x^3

cobin family class

Description

Specifies the information required to fit a cobin generalized linear model with known lambda parameter, using glm().

Usage

cobinfamily(lambda = stop("'lambda' must be specified"), link = "cobit")

Arguments

Value

An object of class "family", a list of functions and expressions needed by glm() to fit a cobin generalized linear model.

cobin generalized linear (mixed) models

Description

Fit Bayesian cobin regression model under canonical link (cobit link) with Markov chain Monte Carlo (MCMC). It supports both fixed-effect only model

y_i \mid x_i \stackrel{ind}{\sim} cobin(x_i^T\beta, \lambda^{-1}),

for i=1,\dots,n, and random intercept model (v 1.0.x only supports random intercept),

y_{ij} \mid x_{ij}, u_i \stackrel{ind}{\sim} cobin(x_{ij}^T\beta + u_i, \lambda^{-1}), \quad u_i\stackrel{iid}{\sim} N(0, \sigma_u^2)

for i=1,\dots,n (group), and j=1,\dots,n_i (observation within group). See dcobin for details on cobin distribution.

Usage

cobinreg(
  formula,
  data,
  link = "cobit",
  contrasts = NULL,
  priors = list(beta_intercept_scale = 100, beta_scale = 100, beta_df = Inf),
  nburn = 1000,
  nsave = 1000,
  nthin = 1,
  MH = FALSE,
  lambda_fixed = NULL
)

Arguments

Details

The prior setting can be controlled with "priors" argument. Prior for regression coefficients are independent normal or t prior centered at 0. "priors" is a named list of:

• beta_intercept_scale, Default 100, the scale of the intercept prior

• beta_scale, Default 100, the scale of nonintercept fixed-effect coefficients

• beta_df, Default Inf, degree of freedom of t prior. If beta_df=Inf, it corresponds to normal prior

• lambda_grid, Default 1:70, candidate for lambda (integer)

• lambda_logprior, Default p(\lambda)\propto \lambda \Gamma(\lambda+1)/\Gamma(\lambda+5), log-prior of lambda. Default choice arises from beta negative binomial distribution; (\lambda-1)\mid \psi \sim negbin(2,\psi), \psi\sim Beta(2,2).

if random intercept model, u ~ InvGamma(a_u,b_u) with

• a_u, Default 1, first parameter of Inverse Gamma prior of u

• b_u, Default 1, second parameter of Inverse Gamma prior of u

Value

Returns list of

if random effect model, also returns

Examples

requireNamespace("betareg", quietly = TRUE)
library(betareg) # for dataset example
data("GasolineYield", package = "betareg")

# basic model 
out1 = cobinreg(yield ~ temp, data = GasolineYield, 
               nsave = 2000, link = "cobit")
summary(out1$post_save)
plot(out1$post_save)

# random intercept model
out2 = cobinreg(yield ~ temp + (1 | batch), data = GasolineYield, 
               nsave = 2000, link = "cobit")
summary(out2$post_save)
plot(out2$post_save)

Density of Irwin-Hall distribution

Description

Irwin-Hall distribution is defined as a sum of m uniform (0,1) distribution. Its density is given as

f(x;m) = \frac{1}{(m-1)!}\sum_{k=0}^{m} (-1)^k {m \choose k} \max(0,x-k)^{m-1}, 0 < x < m

The density of Bates distribution, defined as an average of m uniform (0,1) distribution, can be obtained from change-of-variable (y = x/m),

h(y;m) = \frac{m}{(m-1)!}\sum_{k=0}^{m} (-1)^k {m \choose k} \max(0,my-k)^{m-1}, 0 < y < 1

Usage

dIH(x, m, log = FALSE)

Arguments

Details

Due to alternating series representation, m > 80 may yield numerical issues

Value

(log) density evaluated at x

Examples

m = 8
xgrid= seq(0, m, length = 500)
hist(colSums(matrix(runif(m*1000), nrow = m, ncol = 1000)), freq = FALSE) 
lines(xgrid, dIH(xgrid, m, log = FALSE))
# Bates distribution
xgrid= seq(0, 1, length = 500)
hist(colMeans(matrix(runif(m*1000), nrow = m, ncol = 1000)), freq = FALSE) 
lines(xgrid, m*dIH(xgrid*m, m, log = FALSE))

Density function of cobin (continuous binomial) distribution

Description

Continuous binomial distribution with natural parameter \theta and dispersion parameter 1/\lambda, in short Y \sim cobin(\theta, \lambda^{-1}), has density

p(y; \theta, \lambda^{-1}) = h(y;\lambda) \exp(\lambda \theta y - \lambda B(\theta)), \quad 0 \le y \le 1

where B(\theta) = \log\{(e^\theta - 1)/\theta\} and h(y;\lambda) = \frac{\lambda}{(\lambda-1)!}\sum_{k=0}^{\lambda} (-1)^k {\lambda \choose k} \max(0,\lambda y-k)^{\lambda-1}. When \lambda = 1, it becomes continuous Bernoulli distribution.

Usage

dcobin(x, theta, lambda, log = FALSE)

Arguments

Details

For the evaluation of h(y;\lambda), see ?cobin::dIH.

Value

density of cobin(\theta,\lambda^{-1})

Examples

xgrid = seq(0, 1, length = 500)
plot(xgrid, dcobin(xgrid, 0, 1), type="l", ylim = c(0,3)) # uniform 
lines(xgrid, dcobin(xgrid, 0, 3))
plot(xgrid, dcobin(xgrid, 2, 3), type="l")
lines(xgrid, dcobin(xgrid, -2, 3))

Density function of micobin (mixture of continuous binomial) distribution

Description

Micobin distribution with natural parameter \theta and dispersion psi, denoted as micobin(\theta, \psi), is defined as a dispersion mixture of cobin:

Y \sim micobin(\theta, \psi) \iff Y | \lambda \sim cobin(\theta, \lambda^{-1}), (\lambda-1) \sim negbin(2, \psi)

so that micobin density is a weighted sum of cobin density with negative binomial weights.

Usage

dmicobin(x, theta, psi, r = 2, log = FALSE, l_max = 70)

Arguments

Value

density of micobin(\theta, \psi)

Examples

hist(rcobin(1000, 2, 3), freq = FALSE)
xgrid = seq(0, 1, length = 500)
lines(xgrid, dcobin(xgrid, 2, 3))

Find the MLE of cobin GLM

Description

Find the maximum likelihood estimate of a cobin generalized linear model with unknown dispersion. This is a modification of stats::glm to include estimation of the additional parameter, lambda, for a cobin generalized linear model, in a similar manner to the MASS::glm.nb. Note that MLE of regression coefficient does not depends on lambda.

Usage

glm.cobin(
  formula,
  data,
  weights,
  subset,
  na.action,
  start = NULL,
  etastart,
  mustart,
  control = glm.control(...),
  method = "glm.fit",
  model = TRUE,
  x = FALSE,
  y = TRUE,
  contrasts = NULL,
  ...,
  lambda_list = 1:70,
  link = "cobit",
  verbose = TRUE
)

Arguments

Details

Since dispersion parameter lambda is discrete, it does not provide standard error of lambda. With cobit link, we strongly encourage Bayesian approaches, using cobin::cobinreg() function.

Value

The object is like the output of glm but contains additional components, the ML estimate of lambda and the log-likelihood values for each lambda in the lambda_list.

Examples

 

requireNamespace("betareg", quietly = TRUE)
library(betareg)# for dataset example
data(GasolineYield, package = "betareg")
# cobin regression, frequentist
out_freq = glm.cobin(yield ~ temp, data = GasolineYield, link = "cobit")
summary(out_freq) 
# cobin regression, Bayesian
out = cobinreg(yield ~ temp, data = GasolineYield, 
               nsave = 10000, link = "cobit")
summary(out$post_save)
plot(out$post_save)

micobin generalized linear (mixed) models

Description

Fit Bayesian micobin regression model under canonical link (cobit link) with Markov chain Monte Carlo (MCMC). It supports both fixed-effect only model

y_i \mid x_i \stackrel{ind}{\sim} micobin(x_i^T\beta, \psi),

for i=1,\dots,n, and random intercept model (v 1.0.x only supports random intercept),

y_{ij} \mid x_{ij}, u_i \stackrel{ind}{\sim} micobin(x_{ij}^T\beta + u_i, \psi), \quad u_i\stackrel{iid}{\sim} N(0, \sigma_u^2)

for i=1,\dots,n (group), and j=1,\dots,n_i (observation within group). See dmicobin for details on micobin distribution.

Usage

micobinreg(
  formula,
  data,
  link = "cobit",
  contrasts = NULL,
  priors = list(beta_intercept_scale = 100, beta_scale = 100, beta_df = Inf),
  nburn = 1000,
  nsave = 1000,
  nthin = 1,
  psi_fixed = NULL
)

Arguments

Details

The prior setting can be controlled with "priors" argument. Prior for regression coefficients are independent normal or t prior centered at 0. "priors" is a named list of:

• beta_intercept_scale, Default 100, the scale of the intercept prior

• beta_scale, Default 100, the scale of nonintercept fixed-effect coefficients

• beta_df, Default Inf, degree of freedom of t prior. If beta_df=Inf, it corresponds to normal prior

• lambda_max, Default 70, upper bound for lambda (integer)

• psi_ab, Default c(2,2), beta shape parameters for \psi (length 2 vector).

if random intercept model, u ~ InvGamma(a_u,b_u) with

• a_u, Default 1, first parameter of Inverse Gamma prior of u

• b_u, Default 1, second parameter of Inverse Gamma prior of u

Value

Returns list of

if random effect model, also returns

Examples

requireNamespace("betareg", quietly = TRUE)
library(betareg)# for dataset example
data("GasolineYield", package = "betareg")

# basic model 
out1 = micobinreg(yield ~ temp, data = GasolineYield, 
               nsave = 2000, link = "cobit")
summary(out1$post_save)
plot(out1$post_save)

# random intercept model
out2 = micobinreg(yield ~ temp + (1 | batch), data = GasolineYield, 
               nsave = 2000, link = "cobit")
summary(out2$post_save)
plot(out2$post_save)

Cumulative distribution function of cobin (continuous binomial) distribution

Description

Continuous binomial distribution with natural parameter \theta and dispersion parameter 1/\lambda, in short Y \sim cobin(\theta, \lambda^{-1}), has density

p(y; \theta, \lambda^{-1}) = h(y;\lambda) \exp(\lambda \theta y - \lambda B(\theta)), \quad 0 \le y \le 1

where B(\theta) = \log\{(e^\theta - 1)/\theta\} and h(y;\lambda) = \frac{\lambda}{(\lambda-1)!}\sum_{k=0}^{\lambda} (-1)^k {\lambda \choose k} \max(0,\lambda y-k)^{\lambda-1}. When \lambda = 1, it becomes continuous Bernoulli distribution.

Usage

pcobin(q, theta, lambda)

Arguments

Value

c.d.f. of cobin(\theta,\lambda^{-1})

Examples

xgrid = seq(0, 1, length = 500)
out = pcobin(xgrid, 1, 2)
plot(ecdf(rcobin(10000, 1, 2)))
lines(xgrid, out, col = 2)

Cumulative distribution function of micobin (mixture of continuous binomial) distribution

Description

Micobin distribution with natural parameter \theta and dispersion psi, denoted as micobin(\theta, \psi), is defined as a dispersion mixture of cobin:

Y \sim micobin(\theta, \psi) \iff Y | \lambda \sim cobin(\theta, \lambda^{-1}), (\lambda-1) \sim negbin(2, \psi)

so that micobin cdf is a weighted sum of cobin cdf with negative binomial weights.

Usage

pmicobin(q, theta, psi, r = 2, l_max = 70)

Arguments

Value

c.d.f. of micobin(\theta, \psi)

Examples

xgrid = seq(0, 1, length = 500)
out = pmicobin(xgrid, 1, 1/2)
plot(ecdf(rmicobin(10000, 1, 1/2)))
lines(xgrid, out, col = 2)

Random variate generation for cobin (continuous binomial) distribution

Description

Continuous binomial distribution with natural parameter \theta and dispersion parameter 1/\lambda, in short Y \sim cobin(\theta, \lambda^{-1}), has density

p(y; \theta, \lambda^{-1}) = h(y;\lambda) \exp(\lambda \theta y - \lambda B(\theta)), \quad 0 \le y \le 1

where B(\theta) = \log\{(e^\theta - 1)/\theta\} and h(y;\lambda) = \frac{\lambda}{(\lambda-1)!}\sum_{k=0}^{\lambda} (-1)^k {\lambda \choose k} \max(0,\lambda y-k)^{\lambda-1}. When \lambda = 1, it becomes continuous Bernoulli distribution.

Usage

rcobin(n, theta, lambda)

Arguments

Details

The random variate generation is based on the fact that cobin(\theta, \lambda^{-1}) is equal in distribution to the sum of \lambda cobin(\theta, 1) random variables, scaled by \lambda^{-1}. Random variate generation for continuous Bernoulli is done by inverse cdf transform method.

Value

random samples from cobin(\theta,\lambda^{-1}).

Examples

hist(rcobin(1000, 2, 3), freq = FALSE)
xgrid = seq(0, 1, length = 500)
lines(xgrid, dcobin(xgrid, 2, 3))

Sample Kolmogorov-Gamma random variables

Description

A random variable X follows Kolmogorov-Gamma(b,c) distribution, in short KG(b,c), if

X \stackrel{d}{=} \dfrac{1}{2\pi^2}\sum_{k=1}^\infty \dfrac{\epsilon_k}{k^2 + c^2/(4\pi^2)}, \quad \epsilon_k\stackrel{iid}{\sim} Gamma(b,1)

where \stackrel{d}{=} denotes equality in distribution. The random variate generation is based on alternating series method, a fast and exact method (without infinite sum truncation) implemented in cpp. This function only supports integer b, which is sufficient for cobin and micobin regression models.

Usage

rkgcpp(n, b, c)

Arguments

Value

It returns n independent Kolmogorov-Gamma(b[i],c[i]) samples. If input b or c is scalar, it is assumed to be length n vector with same entries.

Examples

rkgcpp(100, 1, 2)
rkgcpp(100, 1, rnorm(100))
rkgcpp(100, rep(c(1,2),50), rnorm(100))

Random variate generation for micobin (mixture of continuous binomial) distribution

Description

Micobin distribution with natural parameter \theta and dispersion psi, denoted as micobin(\theta, \psi), is defined as a dispersion mixture of cobin:

Y \sim micobin(\theta, \psi) \iff Y | \lambda \sim cobin(\theta, \lambda^{-1}), (\lambda-1) \sim negbin(2, \psi)

Usage

rmicobin(n, theta, psi, r = 2)

Arguments

Value

random samples from micobin(\theta,\psi).

Examples

hist(rmicobin(1000, 2, 1/3), freq = FALSE)
xgrid = seq(0, 1, length = 500)
lines(xgrid, dmicobin(xgrid, 2, 1/3))

spatial cobin regression model

Description

Fit Bayesian spatial cobin regression model under canonical link (cobit link) with Markov chain Monte Carlo (MCMC).

y(s_{i}) \mid x(s_{i}), u(s_i) \stackrel{ind}{\sim} cobin(x(s_{i})^T\beta + u(s_i), \lambda^{-1}), \quad u(\cdot)\sim GP

for i=1,\dots,n. See dcobin for details on cobin distribution. It currently only supports mean zero GP with exponential covariance

cov(u(s_i), u(s_j)) = \sigma_u^2\exp(-\phi_u d(s_i,s_j))

where \phi_u corresponds to inverse range parameter.

Usage

spcobinreg(
  formula,
  data,
  link = "cobit",
  coords,
  NNGP = FALSE,
  contrasts = NULL,
  priors = list(beta_intercept_scale = 10, beta_scale = 2.5, beta_df = Inf),
  nngp.control = list(n.neighbors = 15, ord = order(coords[, 1])),
  nburn = 1000,
  nsave = 1000,
  nthin = 1
)

Arguments

Details

The prior setting can be controlled with "priors" argument. Prior for regression coefficients are independent normal or t prior centered at 0. "priors" is a named list of:

• beta_intercept_scale, Default 100, the scale of the intercept prior

• beta_scale, Default 100, the scale of nonintercept fixed-effect coefficients

• beta_df, Default Inf, degree of freedom of t prior. If beta_df=Inf, it corresponds to normal prior

• lambda_grid, Default 1:70, candidate for lambda (integer)

• lambda_logprior, Default p(\lambda)\propto \lambda \Gamma(\lambda+1)/\Gamma(\lambda+5), log-prior of lambda. Default choice arises from beta negative binomial distribution; (\lambda-1)\mid \psi \sim negbin(2,\psi), \psi\sim Beta(2,2).

• logprior_sigma.sq, Default half-Cauchy on the sd(u) =\sigma_u, log prior of var(u)=\sigma_u^2

• phi_lb, lower bound of uniform prior of \phi_u (inverse range parameter of spatial random effect). Can be same as phi_ub

• phi_ub, lower bound of uniform prior of \phi_u (inverse range parameter of spatial random effect). Can be same as phi_lb

Value

Returns list of

if NNGP = TRUE, also returns

Examples

 # Please see https://github.com/changwoo-lee/cobin-reproduce

spatial micobin regression model

Description

Fit Bayesian spatial micobin regression model under canonical link (cobit link) with Markov chain Monte Carlo (MCMC).

y(s_{i}) \mid x(s_{i}), u(s_i) \stackrel{ind}{\sim} micobin(x(s_{i})^T\beta + u(s_i), \psi), \quad u(\cdot)\sim GP

for i=1,\dots,n. See dmicobin for details on micobin distribution. It currently only supports mean zero GP with exponential covariance

cov(u(s_i), u(s_j)) = \sigma_u^2\exp(-\phi_u d(s_i,s_j))

where \phi_u corresponds to inverse range parameter.

Usage

spmicobinreg(
  formula,
  data,
  link = "cobit",
  coords,
  NNGP = FALSE,
  contrasts = NULL,
  priors = list(beta_intercept_scale = 10, beta_scale = 2.5, beta_df = Inf),
  nngp.control = list(n.neighbors = 15, ord = order(coords[, 1])),
  nburn = 1000,
  nsave = 1000,
  nthin = 1
)

Arguments

Details

The prior setting can be controlled with "priors" argument. Prior for regression coefficients are independent normal or t prior centered at 0. "priors" is a named list of:

• beta_intercept_scale, Default 100, the scale of the intercept prior

• beta_scale, Default 100, the scale of nonintercept fixed-effect coefficients

• beta_df, Default Inf, degree of freedom of t prior. If beta_df=Inf, it corresponds to normal prior

• lambda_max, Default 70, upper bound for lambda (integer)

• psi_ab, Default c(2,2), beta shape parameters for \psi (length 2 vector).

• logprior_sigma.sq, Default half-Cauchy on the sd(u) =\sigma_u, log prior of var(u)=\sigma_u^2

• phi_lb, lower bound of uniform prior of \phi_u (inverse range parameter of spatial random effect). Can be same as phi_ub

• phi_ub, lower bound of uniform prior of \phi_u (inverse range parameter of spatial random effect). Can be same as phi_lb

Value

Returns list of

if NNGP = TRUE, also returns

Examples

 # Please see https://github.com/changwoo-lee/cobin-reproduce