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

**URL:** https://discourse.mc-stan.org/t/survival-analysis-with-brms-whats-the-design-matrix/21925
**Category:** Modeling
**Tags:** brms
**Created:** [April 16, 2021, 3:51pm UTC](https://discourse.mc-stan.org/t/survival-analysis-with-brms-whats-the-design-matrix/21925 "2021-04-16T15:51:13Z")
**Posts on this page:** 4
**Page:** 1

<div class="post-metadata">

### Author: ![fsdias](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/fsdias/32/4198_2.png) [@fsdias](https://discourse.mc-stan.org/u/fsdias)
#### Post date: [April 16, 2021, 3:51pm UTC](https://discourse.mc-stan.org/t/survival-analysis-with-brms-whats-the-design-matrix/21925/1 "2021-04-16T15:51:13Z")

</div>

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 %>% 
  mutate(adopted = ifelse(out_event == "Adoption", 1, 0),
         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

---

<div class="post-metadata">

### Author: ![andrjohns](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andrjohns/32/15297_2.png) [@andrjohns](https://discourse.mc-stan.org/u/andrjohns)
#### Post date: [April 17, 2021, 2:43am UTC](https://discourse.mc-stan.org/t/survival-analysis-with-brms-whats-the-design-matrix/21925/2 "2021-04-17T02:43:10Z")

</div>

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

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

rstan_model = stan(model_code=brms_stancode,
                   data=brms_standata)

```

---

<div class="post-metadata">

### Author: ![fsdias](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/fsdias/32/4198_2.png) [@fsdias](https://discourse.mc-stan.org/u/fsdias)
#### Post date: [April 18, 2021, 12:23pm UTC](https://discourse.mc-stan.org/t/survival-analysis-with-brms-whats-the-design-matrix/21925/3 "2021-04-18T12:23:50Z")

</div>

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)

---

<div class="post-metadata">

### Author: ![andrjohns](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andrjohns/32/15297_2.png) [@andrjohns](https://discourse.mc-stan.org/u/andrjohns)
#### Post date: [April 18, 2021, 1:40pm UTC](https://discourse.mc-stan.org/t/survival-analysis-with-brms-whats-the-design-matrix/21925/4 "2021-04-18T13:40:05Z")

</div>

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

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

```