Hi all! I am trying to use the generalized skewed t distribution priors in BASIE Framework for some randomized controlled trial data that I have collected.
I am wondering if there is a way to directly specify skewed t distribution priors in brms. If not, is there a way to pass the Stan code that BASIE provides to brms so that I can specify the desired prior (with mean, sd, skewness, and degrees of freedom)? Sorry I am fairly new to directly writing Stan codes, so any help would be great, thanks!
BASIE prior (the prior is the mixture of two skewed t-distributions)
Here is an example of my model (using brm_multiple for multiple imputed datasets):
fit <- brm_multiple(
self_integrity_esteem_e ~ self_integrity_esteem_b + treated +
class_size + grade + gender + age_b + student_type + father_edu_f + mother_edu_f +
father_occ_salary + mother_occ_salary + adult_members + siblings + duration_int +
(treated | class_id),
data = dat_mi,
prior = prior("normal(0, 1)", class = "b"),
chains = 4,
seed = 1234
)
Below is the Stan code provided by BASIE:
functions {
// taken from https://discourse.mc-stan.org/t/skewed-distribution/1383/19
// This function already incorporates the v/m adjustments, so sd = sigma and mean=mu, approximately. Doesn't work so great at low q
real sgt_log(vector x, real mu, real s, real l, real p, real q) {
// Skewed generalised t
int N;
real lz1;
real lz2;
real v;
real m;
real r;
real out;
N = dims(x)[1];
lz1 = lbeta(1.0/p,q);
lz2 = lbeta(2.0/p,q-1.0/p);
v = q^(-1.0/p)*((3*l^2+1)*exp(lbeta(3.0/p,q-2.0/p)-lz1)-4*l^2*exp(lz2-lz1)^2)^(-0.5);
m = 2*v*s*l*q^(1.0/p)*exp(lz2-lz1);
out = 0;
for (n in 1:N) {
r = x[n]-mu+m;
if (r<0)
out = out+log(p)-log(2*v*s*q^(1.0/p)*exp(lz1)*(fabs(r)^p /(q*(v*s)^p*(l*(-1)+1)^p)+1)^(1.0/p+q));
else
out = out+log(p)-log(2*v*s*q^(1.0/p)*exp(lz1)*(fabs(r)^p /(q*(v*s)^p*(l*(1)+1)^p)+1)^(1.0/p+q));
}
return out;
}
// Skew t is just sgt with p fixed at 2
real skew_t_log(vector x, real mu, real s, real l, real q) {
// Skewed generalised t
int N;
real lz1;
real lz2;
real v;
real m;
real r;
real out;
real l2;
real simp1;
real simp2;
N = dims(x)[1];
lz1 = lbeta(0.5,q);
lz2 = lbeta(1.0,q-0.5);
v = q^-.5*((3*l^2+1)*exp(lbeta(1.5,q-1)-lz1)-4*l^2*exp(lz2-lz1)^2)^(-0.5);
m = 2*v*s*l*sqrt(q)*exp(lz2-lz1);
out = 0;
// Some of the pieces in the loop can be pulled back out
l2 = log(2);
simp1 = 2*v*s*sqrt(q)*exp(lz1);
simp2 = q*(v*s)^2;
for (n in 1:N) {
r = x[n]-mu+m;
// With p=2, abs(r)^p == r^p so no need for abs
if (r<0)
out = out+l2-log(simp1*(r^2 /(simp2*(l*(-1)+1)^2)+1)^(0.5+q));
else
out = out+l2-log(simp1*(r^2 /(simp2*(l*(1)+1)^2)+1)^(0.5+q));
}
return out;
}
}
data {
real impact_est; // estimated treatment effects
real<lower=0> se_est; // s.e. of effect estimates
real mu; //user-specified prior mean of true effect distribution
real sigma; //user-specified prior std dev
real lambda; //user-specified prior skew
real q; //user-specified prior taily
real p; //user-specified prior plateau
real df;
}
parameters {
vector[1] impact_true;
}
model {
if (p==2) impact_true ~ skew_t(mu, sigma, lambda, q);
else impact_true ~ sgt(mu, sigma, lambda,p, q);
impact_est - impact_true ~ student_t(df,0,se_est);
}