Dear users,
my question is about Dirichlet regression but it can hold for basic regression as well.
I want my predictors’ coefficients (Age and Year) to have a sum-zero constraint.
Specifically, Age and Year are both vectors of coefficients, for each vector, the sum should be zero.
\sum Age = 0 and \sum Year = 0
In my understanding, I should work on the design matrix which I can get in the following code using brsm.
Any suggestion?
Many thanks for the outstanding work of this community.
Get data
#install.packages("brms")
library(brms)
library(rstan)
library(dplyr)
bind <- function(...) cbind(...)
#Data simulation
N <- length(2000:2009)*length(1:10)
df <- data.frame(
y1 = rbinom(N, 10, 0.5), y2 = rbinom(N, 10, 0.7),
y3 = rbinom(N, 10, 0.9), Age = as.factor(1:10),Year = as.factor(rep(2000:2009,each=length(1:10)))
) %>%
mutate(
size = y1 + y2 + y3,
y1 = y1 / size,
y2 = y2 / size,
y3 = y3 / size
)
df$y <- with(df, cbind(y1, y2, y3))
str(df)
Get model
stan_mod <- make_stancode(bind(y1, y2, y3) ~ Age+Year, df, dirichlet(),save_model = "model1.stan")
brms_standata <- make_standata(bind(y1, y2, y3) ~ Age+Year, df, dirichlet(),save_model = "data1.stan")
# design matrix
brms_standata$X_muy2
brms_standata$X_muy3
stan_mod_a <- stan_model(file="model1.stan")
rstan_model = stan("model1.stan",model_code=stan_mod_a,
data=brms_standata)
Stan model - script
// generated with brms 2.18.0
functions {
/* dirichlet-logit log-PDF
* Args:
* y: vector of real response values
* mu: vector of category logit probabilities
* phi: precision parameter
* Returns:
* a scalar to be added to the log posterior
*/
real dirichlet_logit_lpdf(vector y, vector mu, real phi) {
return dirichlet_lpdf(y | softmax(mu) * phi);
}
}
data {
int<lower=1> N; // total number of observations
int<lower=2> ncat; // number of categories
vector[ncat] Y[N]; // response array
int<lower=1> K_muy2; // number of population-level effects
matrix[N, K_muy2] X_muy2; // population-level design matrix
int<lower=1> K_muy3; // number of population-level effects
matrix[N, K_muy3] X_muy3; // population-level design matrix
int prior_only; // should the likelihood be ignored?
}
transformed data {
int Kc_muy2 = K_muy2 - 1;
matrix[N, Kc_muy2] Xc_muy2; // centered version of X_muy2 without an intercept
vector[Kc_muy2] means_X_muy2; // column means of X_muy2 before centering
int Kc_muy3 = K_muy3 - 1;
matrix[N, Kc_muy3] Xc_muy3; // centered version of X_muy3 without an intercept
vector[Kc_muy3] means_X_muy3; // column means of X_muy3 before centering
for (i in 2:K_muy2) {
means_X_muy2[i - 1] = mean(X_muy2[, i]);
Xc_muy2[, i - 1] = X_muy2[, i] - means_X_muy2[i - 1];
}
for (i in 2:K_muy3) {
means_X_muy3[i - 1] = mean(X_muy3[, i]);
Xc_muy3[, i - 1] = X_muy3[, i] - means_X_muy3[i - 1];
}
}
parameters {
vector[Kc_muy2] b_muy2; // population-level effects
real Intercept_muy2; // temporary intercept for centered predictors
vector[Kc_muy3] b_muy3; // population-level effects
real Intercept_muy3; // temporary intercept for centered predictors
real<lower=0> phi; // precision parameter
}
transformed parameters {
real lprior = 0; // prior contributions to the log posterior
lprior += student_t_lpdf(Intercept_muy2 | 3, 0, 2.5);
lprior += student_t_lpdf(Intercept_muy3 | 3, 0, 2.5);
lprior += gamma_lpdf(phi | 0.01, 0.01);
}
model {
// likelihood including constants
if (!prior_only) {
// initialize linear predictor term
vector[N] muy2 = rep_vector(0.0, N);
// initialize linear predictor term
vector[N] muy3 = rep_vector(0.0, N);
// linear predictor matrix
vector[ncat] mu[N];
muy2 += Intercept_muy2 + Xc_muy2 * b_muy2;
muy3 += Intercept_muy3 + Xc_muy3 * b_muy3;
for (n in 1:N) {
mu[n] = transpose([0, muy2[n], muy3[n]]);
}
for (n in 1:N) {
target += dirichlet_logit_lpdf(Y[n] | mu[n], phi);
}
}
// priors including constants
target += lprior;
}
generated quantities {
// actual population-level intercept
real b_muy2_Intercept = Intercept_muy2 - dot_product(means_X_muy2, b_muy2);
// actual population-level intercept
real b_muy3_Intercept = Intercept_muy3 - dot_product(means_X_muy3, b_muy3);
}