In the following example, we demonstrate how to perform regression
modeling using the JSEP distribution. The code provided below shows how
to set up the regression model, specify the priors, and fit the model
using the bnrm function.
set.seed(123)
x <- runif(50)
e <- rjsep(50, 0, 0.1, 0.8,2) # Simulated errors from JSEP distribution
y <- 0.5 + 0.8*x + e
xy<-cbind(y,x)
df_reg <- data.frame(y, x)
# Basic exploration
plot(x, y, main = "Regression Data", pch = 19, col = "steelblue")
#save pre-computed data
#saveRDS(xy, file.path(precomp_dir,"jsep_data.rds")) # Save pre-computed dataWe specify the regression model using the brms formula
interface. The regression formula is defined as y ~ x,
where y is the response variable and x is the
predictor variable. The priors for the regression parameters are set as
follows: Prior of \(\alpha\) is set to
log-normal(0, 0.5), \(\beta\) is set to
log-normal(1, 0.5), \(\sigma\) is set
to half-normal(0, 1), and the intercept and slope are set to normal(0,
1). The code syntax for defining the model and these priors is shown
below.
# Define regression formula
formula_reg <- brms::bf(y ~ x)
# Priors for regression
prior_reg <- c(
set_prior("lognormal(log(2),0.25)", class = "alpha"),
set_prior("lognormal(log(2),0.25)", class = "beta"),
set_prior("normal(0,1)", class = "sigma"),
set_prior("normal(0,1)", class = "Intercept"),
set_prior("normal(0,1)", class = "b")
)We evaluate our chosen priors via a prior predictive check, which
entails drawing simulated datasets solely from those priors to verify
they can plausibly reproduce the characteristics of the actual data. In
practice, this is done by calling
bnrm(..., sample_prior = "only"), which generates the
prior‐based simulations without fitting the model to the observed
data.
# Fit regression model
fit_ppc_jsep <- bnrm(
formula = formula_reg,
data = df_reg,
family = jsep(),
prior = prior_reg,
chains = 2,
cores = 1,
iter = 2000,
sample_prior = "only", # Prior predictive check
warmup = 1000,
seed = 123,
refresh = 0
)
# prior predictive check fir jesp
ppc_plot_jsep<- pp_check(fit_ppc_jsep, type = "dens_overlay", nsamples = 200)
#saveRDS(ppc_plot_jsep, file.path(precomp_dir,"ppc_plot_jsep.rds")) # Save pre-computed plot
print(ppc_plot_jsep)
# Show prior predictive check plotThe prior predictive check plot overlays the data simulated from our priors onto the real observations, letting us visually judge whether those priors can plausibly reproduce the distribution of the actual data.
We fit the regression model using the bnrm function,
specifying the formula, data, family, and priors defined earlier. The
fitting process uses MCMC sampling to estimate the posterior
distributions of the parameters.
fit_jsep_bnrm <- bnrm(
formula = formula_reg,
data = df_reg,
family = jsep(),
prior = prior_reg,
chains = 2,
cores = 1,
iter = 2000,
warmup = 1000,
seed = 123,
refresh = 0
)
#saveRDS(fit_jsep_bnrm, file.path(precomp_dir, "fit_jsep_bnrm.rds")) # Save pre-computed fitAfter fitting the model, we can summarize the results to see the estimated parameters and their posterior distributions. The summary will include the mean, standard deviation, and credible intervals for each parameter. The convergence diagnostics will also be checked to ensure that the MCMC chains have mixed well and converged to the posterior distribution. Posterior predictive checks will be performed to assess the model fit and the ability of the model to predict new data.
summary(fit_jsep_bnrm)
# Convergence diagnostics
trace_plot <- mcmc_trace(fit_jsep_bnrm, facet_args = list(ncol = 2),pars = c("alpha", "beta", "sigma", "Intercept","b_x"))
print(trace_plot)
# Posterior predictive check
pp_plot <- pp_check(fit_jsep_bnrm)
#saveRDS(pp_plot, file.path(precomp_dir,"pp_plot_jsep.rds")) # Save pre-computed plot
# Model evaluation
loo_result <- loo(fit_jsep_bnrm, moment_match = TRUE)
#saveRDS(loo_result, file.path(precomp_dir,"loo_result_jsep.rds")) # Save pre-computed loo result
print(loo_result)
#> Family: jsep
#> Links: mu = identity; sigma = identity; alpha = identity; beta = identity
#> Formula: y ~ x
#> Data: data (Number of observations: 50)
#> Draws: 2 chains, each with iter = 2000; warmup = 1000; thin = 1;
#> total post-warmup draws = 2000
#>
#> Regression Coefficients:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept 0.47 0.04 0.40 0.54 1.00 743 1096
#> x 0.81 0.05 0.71 0.92 1.00 994 1192
#>
#> Further Distributional Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 0.12 0.03 0.08 0.18 1.00 658 995
#> alpha 0.95 0.12 0.73 1.22 1.00 879 1189
#> beta 2.36 0.52 1.49 3.58 1.00 1188 1433
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).