# Survival analysis with "brms" - what's the design matrix?

Hi,
I’m trying to run a survival model in “brms” with the attached dataset:

``````AustinCats <- read.csv("AustinCats.csv")
d <- AustinCats
rm(AustinCats)

d <-
d %>%
mutate(black = ifelse(color == "Black", "black", "other"))

d <-
d %>%
censored = ifelse(out_event != "Adoption", 1, 0))

model <-
brm(data = d,
family = exponential,
days_to_event | cens(censored) ~ 0 + black,
prior(normal(0, 1), class = b),
iter = 2000, warmup = 1000, chains = 4, cores = 4
)
``````

The corresponding Stan code is:

``````stancode(model)

// generated with brms 2.15.0
functions {
}
data {
int<lower=1> N;  // total number of observations
vector[N] Y;  // response variable
int<lower=-1,upper=2> cens[N];  // indicates censoring
int<lower=1> K;  // number of population-level effects
matrix[N, K] X;  // population-level design matrix
int prior_only;  // should the likelihood be ignored?
}
transformed data {
}
parameters {
vector[K] b;  // population-level effects
}
transformed parameters {
}
model {
// likelihood including constants
if (!prior_only) {
// initialize linear predictor term
vector[N] mu = X * b;
for (n in 1:N) {
// apply the inverse link function
mu[n] = exp(-(mu[n]));
}
for (n in 1:N) {
// special treatment of censored data
if (cens[n] == 0) {
target += exponential_lpdf(Y[n] | mu[n]);
} else if (cens[n] == 1) {
target += exponential_lccdf(Y[n] | mu[n]);
} else if (cens[n] == -1) {
target += exponential_lcdf(Y[n] | mu[n]);
}
}
}
// priors including constants
target += normal_lpdf(b | 0, 1);
}
generated quantities {
}
``````

My question concerns `matrix[N, K] X; // population-level design matrix`:

1. What is the X matrix for in the model?
2. How do I build it so that I can reproduce the analysis using rstan?

Thanks

1 Like

You can replicate the `brms` analysis in rstan by extracting the Stan code and pre-processed data via:

``````brms_stancode = stancode(model)
brms_standata = standata(model)

rstan_model = stan(model_code=brms_stancode,
data=brms_standata)
``````
2 Likes

Thanks @andrjohns

I was able to figure out that I need to use model.matrix() to create X, but I can’t figure out how .

Failed attempts:

X<-model.matrix(days_to_event, color_id, data=d)
X<-model.matrix(days_to_event, as.factor(color_id), data=d)
X<-model.matrix( as.factor(color_id), days_to_event,data=d)
X<-model.matrix( color_id, days_to_event,data=d)

To recreate `X` using `model.matrix`, you would call:

``````model.matrix(~ 0 + black, data=d)
``````
3 Likes