Plate/Experiment model

To generate the Stan code for modeling experiment and plate effects, I used the following model.

Y ~ (1|experiment/plate)

I made this model under the assumption that the relationship between plates and experiments replicates the relationsip between “cask” and “batch” in the Pastes dataset in lme4.

See the following “Pastes” data from lme4. This was the only example I found in lme4 that had only nested categorical variables like experiment and plate, no continuous regressors.

Strength of a chemical paste product; its quality depending on the delivery batch, and the cask within the delivery. 3 casks selected at random from each delivery were sampled and the samples were kept for reference. Ten of the delivery batches were sampled at random and two analytical tests carried out on each of the 30 samples”.

• strength is the response
• batch is the delivery batch from which the sample was sample. A factor with 10 levels: ‘A’ to ‘J’.
• cask is a cask within the delivery batch from which the sample was chosen. A factor with 3 levels: ‘a’ to ‘c’.
• sample is the sample of paste whose strength was assayed, two assays per sample. A factor with 30 levels: ‘A:a’ to ‘J:c’.

I used this example to reason through how we model experiments and plates. batches correspond to experiment and casks correspond to plates.

Note that casks are encoded as factors, with levels ‘a’, ‘b’, and ‘c’ nested completely in batches.

levels(.data$cask)

[1] "a" "b" "c"

In contrast, the samples column gives a unique id for each cask, in this case created by concatenating the batch and cask columns.

levels(.data$sample)

 [1] "A:a" "A:b" "A:c" "B:a" "B:b" "B:c" "C:a" "C:b" "C:c" "D:a" "D:b" "D:c" "E:a" "E:b" "E:c" "F:a" "F:b" "F:c" "G:a"
[20] "G:b" "G:c" "H:a" "H:b" "H:c" "I:a" "I:b" "I:c" "J:a" "J:b" "J:c"

In other words, both ways of encoding the nested variable are present, with levels nested within the parent variable, and with unique labels for each instance of the nested variable.

The latter ‘unique’ case matches the labeling of our plates with unique ZIDs. The former case of ‘within-parent-levels’ matches labeling of row and column indexes as integers between 1-12.

This example in lme4 comes with 4 equivilent formulas.

You’ll note that brms generates the same code for each of them.

code1 == code2; code2 == code3; code3 == code4

[1] TRUE
[1] TRUE
[1] TRUE

The code generated – I model experiments and plates using this code

code2

// generated with brms 1.9.0
functions { 
} 
data { 
  int<lower=1> N;  // total number of observations 
  vector[N] Y;  // response variable 
  // data for group-level effects of ID 1 
  int<lower=1> J_1[N]; 
  int<lower=1> N_1; 
  int<lower=1> M_1; 
  vector[N] Z_1_1; 
  // data for group-level effects of ID 2 
  int<lower=1> J_2[N]; 
  int<lower=1> N_2; 
  int<lower=1> M_2; 
  vector[N] Z_2_1; 
  int prior_only;  // should the likelihood be ignored? 
} 
transformed data { 
} 
parameters { 
  real temp_Intercept;  // temporary intercept 
  real<lower=0> sigma;  // residual SD 
  vector<lower=0>[M_1] sd_1;  // group-level standard deviations 
  vector[N_1] z_1[M_1];  // unscaled group-level effects 
  vector<lower=0>[M_2] sd_2;  // group-level standard deviations 
  vector[N_2] z_2[M_2];  // unscaled group-level effects 
} 
transformed parameters { 
  // group-level effects 
  vector[N_1] r_1_1 = sd_1[1] * (z_1[1]); 
  // group-level effects 
  vector[N_2] r_2_1 = sd_2[1] * (z_2[1]); 
} 
model { 
  vector[N] mu = rep_vector(0, N) + temp_Intercept; 
  for (n in 1:N) { 
    mu[n] = mu[n] + (r_1_1[J_1[n]]) * Z_1_1[n] + (r_2_1[J_2[n]]) * Z_2_1[n]; 
  } 
  // priors including all constants 
  target += student_t_lpdf(sigma | 3, 0, 10)
    - 1 * student_t_lccdf(0 | 3, 0, 10); 
  target += student_t_lpdf(sd_1 | 3, 0, 10)
    - 1 * student_t_lccdf(0 | 3, 0, 10); 
  target += normal_lpdf(z_1[1] | 0, 1); 
  target += student_t_lpdf(sd_2 | 3, 0, 10)
    - 1 * student_t_lccdf(0 | 3, 0, 10); 
  target += normal_lpdf(z_2[1] | 0, 1); 
  // likelihood including all constants 
  if (!prior_only) { 
    target += normal_lpdf(Y | mu, sigma); 
  } 
} 
generated quantities { 
  // actual population-level intercept 
  real b_Intercept = temp_Intercept; 
} 

J_1 is just integer indices for batch.

And J_2 is just integer indices for unique cask id