BridgeChangeReg: Bridge Regression with Change-point

In jongheepark/BridgeChange: Bayesian Change-point models with Bridge prior for analyzing time-series and panel data.

Bridge Regression with Change-point

Description

Univariate response linear change-point model with Bridge prior.

Usage

BridgeChangeReg(
  formula,
  data = parent.frame(),
  n.break = 0,
  standardize = TRUE,
  intercept = TRUE,
  mcmc = 100,
  burn = 100,
  verbose = 100,
  thin = 1,
  reduce.mcmc = NULL,
  a = NULL,
  b = NULL,
  c0 = 0.1,
  d0 = 0.1,
  nu.shape = 2,
  nu.rate = 2,
  known.alpha = FALSE,
  alpha.start = 1,
  alpha.limit = FALSE,
  alpha.MH = TRUE,
  beta.start = NULL,
  Waic = FALSE,
  marginal = FALSE
)

Arguments

data	Data.frame object. Design matrix. Columns correspond to variables.
n.break	The number of change-point(s). If n.break = 0, the model corresponds to the usual regression. Default is n.break = 0.
standardize	If TRUE, y and X are both scaled to have zero mean and unit variance. Default is TRUE. We recommend standardize = TRUE unless the original data are already scaled.
intercept	A boolean. If TRUE, estimate intercept by MLE. The estimated intercept is used to detect breaks. This option does not affect the result when n.break = 0. We recommend intercept = TRUE when the number of break is not zero.
mcmc	The number of MCMC iterations after burn-in.
verbose	A switch which determines whether or not the progress of the sampler is printed to the screen. If verbose is greater than 0, the iteration number and the posterior density samples are printed to the screen every verboseth iteration.
thin	The thinning interval used in the simulation. The number of MCMC iterations must be divisible by this value.
reduce.mcmc	The number of reduced MCMC iterations for marginal likelihood computations. If reduce.mcmc = NULL, mcmc/thin is used.
a	a is the shape1 beta prior for transition probabilities. By default, the expected duration is computed and corresponding a and b values are assigned. The expected duration is the sample period divided by the number of states.
b	b is the shape2 beta prior for transition probabilities. By default, the expected duration is computed and corresponding a and b values are assigned. The expected duration is the sample period divided by the number of states.
c0	Scale parameter for Gamma distribution. Used for the prior of σ^2. Default is 0.1.
d0	Shape parameter for Gamma distribution. Used for the prior of σ^2. Default is 0.1.
nu.shape	Shape parameter for Gamma distribution. Used for the prior for τ. Default is 2.0.
nu.rate	Rate parameter for Gamma distribution. Used for the prior for τ. Default is 2.0.
known.alpha	If TRUE, a user must specify a numeric value [0, 2] in alpha.start. Default is FALSE and therefore α will be estimated.
alpha.start	Starting value for alpha. When known.alpha = TRUE, alpha is fixed to the value of this argument. Default is 1.
alpha.limit	If TRUE, alpha is sampled from (0,1), otherwise alpha is sampled between (0,2). Default is FALSE.
alpha.MH	If TRUE, alpha is updated by the Metropolis–Hastings algorithm. If FALSE the Griddy gibbs sampler is used instead. Default is TRUE.
beta.start	Starting values of beta. If NULL, randomly choose beta.start from OLS or standard normal distribution. Default is NULL.
Waic	If TRUE, WAIC is computed after the parameter estimation. Default is FALSE.
marginal	If TRUE, the marginal likelihood is computed based on Chib's method. Default is FALSE.
fomula	Inherited from lm. For example, Y ~ X + Z.
burnin	The number of burn-in iterations for the sampler.

Value

An mcmc object that contains the posterior sample of coefficients and variances as columns. Rows correspond to mcmc samples. This object can be summarized by functions provided by the coda package. The object contains an attribute "intercept" that stores mcmc samples for intercept and an attribute state storage matrix that contains posterior samples of hidden states.

Examples

####################################
## no change
####################################
set.seed(1999)

## simulate data
K <- 100
n <- 100
X <- matrix(rnorm(n*K), n, K)
sig2 <- 4
beta.true <- matrix(rnorm(K), K, 1)*5
beta.true[sample(1:K, K/2, replace=FALSE)] <- 0
Y <- X%*%beta.true + rnorm(n, 0, sqrt(sig2))

#########################
## HMBB estimation
#########################
## Be aware that 100 runs are too short for analysis. 
mcmc = 100; burn = 100; verbose = 100; thin = 1;
formula <- Y ~ X
reg.cp0 <- BridgeChangeReg(formula=formula, n.break = 0, Waic = TRUE)
reg.cp1 <- BridgeChangeReg(formula=formula, n.break = 1, Waic = TRUE)
reg.cp2 <- BridgeChangeReg(formula=formula, n.break = 2, Waic = TRUE)

## model selection by WAIC
waic <- WaicCompare(list(reg.cp0, reg.cp1, reg.cp2), print = TRUE)
plotWaic(waic)

## obtain beta.hat
beta.est <- apply(reg.cp0[, grep("beta", colnames(reg.cp0))], 2, mean)

## plot
plot(beta.true, beta.est,  xlab="TRUE", ylab="EST", type = "n",
       xlim = range(beta.est), ylim = range(beta.true), asp = 1)
abline(a=0, b=1, col="red", lty = 3, lwd = 1.5)
points(beta.est, beta.true, col="darkblue")

####################################
## change-point case
####################################
set.seed(11199)
## simulate data
K <- 100
n <- 100
X <- matrix(rnorm(n*K), n, K)
sig2 <- 4
beta.true <- matrix(NA, 2, K)

beta.true[1,] <- matrix(rnorm(K, 1, 1), K, 1)*2
beta.true[2,] <- matrix(rnorm(K, -1, 1), K, 1)
# beta.true[1, sample(1:K, K/2, replace=FALSE)] <- 0
# beta.true[2, sample(1:K, K/2, replace=FALSE)] <- 0
mu1 <- X[1:(n/2), ]%*%beta.true[1,]
mu2 <- X[((n/2)+1):n, ]%*%beta.true[2,] 
Y  <- c(rnorm(n/2, mu1, sqrt(sig2)), rnorm(n/2, mu2, sqrt(sig2)))

## fit he model
formula <- Y ~ X
fit.cp0 <- BridgeChangeReg(formula=formula, mcmc=1000, burn=1000, n.break = 0, Waic = TRUE)
fit.cp1 <- BridgeChangeReg(formula=formula, mcmc=1000, burn=1000, n.break = 1, Waic = TRUE)
fit.cp2 <- BridgeChangeReg(formula=formula, mcmc=1000, burn=1000, n.break = 2, Waic = TRUE)

## model selection by WAIC
## WAIC prefers the no break model but the ground truth is the one break model
(waic <- WaicCompare(list(fit.cp0, fit.cp1, fit.cp2), print = TRUE))

## plot generated data
par(mfrow=c(1,3))
plot(1:length(beta.true[1,]), beta.true[1,],
    xlab="predictor", ylab="coefficients",
    ylim=range(beta.true), type='n', main="Parameter changes")
points(beta.true[1,], col="red", pch="1", cex=1)
points(beta.true[2,], col="blue", pch="2", cex=1)
plot(Y)
plotWaic(waic)

## true break at t = 50
par(mfrow=c(1, 2))
MCMCpack::plotState(fit.cp1, legend.control =c(60, 0.85), main="One break")
MCMCpack::plotState(fit.cp2, legend.control =c(60, 0.85), main="Two breaks")

## obtain beta.hat
beta.est <-  matrix(apply(fit.cp1[, grep("beta", colnames(fit.cp1))], 2, mean), 2, , byrow=TRUE)

## plot
par(mfrow=c(1, 2))
## regime 1
plot(beta.true[1,], beta.est[1,],  xlab="TRUE", ylab="EST", type = "n",main="Regime 1",
       xlim = range(beta.est), ylim = range(beta.true), asp = 1)
abline(a=0, b=1, col="red", lty = 3, lwd = 1.5)
points(beta.true[1,], beta.est[1,], col="darkblue")

## regime 2
plot(beta.true[2,], beta.est[2,],  xlab="TRUE", ylab="EST", type = "n", main="Regime 2",
       xlim = range(beta.est), ylim = range(beta.true), asp = 1)
abline(a=0, b=1, col="red", lty = 3, lwd = 1.5)
points(beta.true[2,], beta.est[2,], col="darkblue")

####################################
## dotplot over time
## time-varying movements of selected covariates
####################################
## all covariates
dotplotRegime(fit.cp1, hybrid=FALSE, location.bar=12, x.location="default",
              text.cex=0.8, main="Time-varying Movements of All Covariates")

jongheepark/BridgeChange documentation built on Jan. 12, 2023, 4:48 p.m.