MCMC for binary regression

MCMC for binary regression

mcmc_bin(data, formula, nsim = 1000, burnin = round(0.1 * nsim),
  lag = 10, type = "logit", sample_c = TRUE, sample_d = TRUE,
  sigma_beta = 3, mean_c = FLAMES:::mean_sd_beta(a = 1, b = 2)$mean,
  sd_c = FLAMES:::mean_sd_beta(a = 1, b = 2)$sd,
  mean_d = FLAMES:::mean_sd_beta(a = 2, b = 1)$mean,
  sd_d = FLAMES:::mean_sd_beta(a = 2, b = 1)$sd, a_lambda = 0.01,
  b_lambda = 0.99, var_df = 0.02, var_c = ifelse(sample_d, 0.005,
  0.02), var_d = ifelse(sample_c, 0.005, 0.02), var_lambda = 0.05,
  method = "ARMS", const = 1, const_beta = 1, const_c = 1,
  const_d = 1, const_df = 1, const_lambda = 1, fitm = FALSE)

Arguments

data	Dataset to be used
formula	Simple formula
nsim	Sample size required for MCMC
burnin	Burn in for MCMC
lag	Lag for MCMC
type	"logit", "probit", "cauchit", "robit", "cloglog" or "loglog"
sample_c	Should c be sampled?
sample_d	Should d be sampled?
sigma_beta	Variance of beta prior
mean_c	Mean for c a priori
sd_c	Standard deviation for c a priori
mean_d	Mean for d a priori
sd_d	Standard deviation for d a priori
a_lambda	Inferior limit for lambda
b_lambda	Superior limit for lambda
var_df	Variance to sample 1-exp(-df/const)
var_c	Variance to sample from c (if sample_d = TRUE) otherwise log(c/(1-c))
var_d	Variance to sample d (if sample_c = TRUE) otherwise log(c/(1-c))
var_lambda	Variance to sample lambda
method	"metropolis" or "ARMS"
const	A constant to help on sampling degrees of freedom \(\tilde{df} = df/c\)
const_beta	A constant to tunning the acceptance rate (default = 1)
const_c	A constant to tunning the acceptance rate (default = 1)
const_d	A constant to tunning the acceptance rate (default = 1)
const_df	A constant to tunning the acceptance rate (default = 1)
const_lambda	A constant to tunning the acceptance rate (default = 1)
fitm	Should return fit measures? "full" for measures based on full dataset or "bootstrap" to use a bootstrap technique

Value

Chains of all parameters

Examples

if (FALSE) {
##-- Seed ----
set.seed(123456)

##-- Data ----
n <- 1000
n_cov <- 2

##-- Covariates
X <- matrix(rnorm(n*n_cov), ncol = n_cov)

##-- Coefficients
betas <- c(0, -1, 0.5)
XBeta <- cbind(1, X)%*%betas

##-- c parameter
c1 <- 0.20
d1 <- 0.95

type_data = "cloglog"

##-- p and y
p <- FLAMES:::inv_link(x = XBeta, type = type_data, df = df)*(d1-c1) + c1
y <- rbinom(n = n, size = 1, prob = p)

bd <- data.frame(y = y, X)

##-- MCMC
nsim <- 2000
burnin <- 10000
lag <- 10

f <- y ~ X1 + X2
type <- type_data

##-- ARMS ~ 4 minutes (soon in c++)
out_arms <- mcmc_bin(data = bd, formula = f,
                     nsim = nsim, burnin = burnin, lag = lag,
                     type = type, sample_c = TRUE, sample_d = TRUE,
                     method = "ARMS")

##-- ARMS ~ 55 seconds (soon in c++)
out_met <- mcmc_bin(data = bd, formula = f,
                    nsim = nsim, burnin = burnin, lag = lag,
                    type = type, sample_c = TRUE, sample_d = TRUE,
                    method = "metropolis")

##-- GLM
out_glm <- glm(formula = f, data = bd, family = "binomial")

coef(out_glm)
coef(out_arms)
coef(out_met)
betas

par(mfrow = c(2, 5))
plot(out_arms, ask = F)
plot(out_met, ask = F)

summary(out_arms)
summary(out_met)
}