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

Modeling
1

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

1 Like
2

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
3

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)

4

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

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