mcmc_bin: MCMC for binary regression

In DouglasMesquita/robit: FLAMES: Flexible Link functions with AsyMptotES

MCMC for binary regression

Description

MCMC for binary regression

Usage

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

## Not run: 
##-- 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)

## End(Not run)

DouglasMesquita/robit documentation built on May 28, 2022, 8:30 a.m.