Beyond Item Response Theory

My previous four posts have focused on basic item response models¹. These models are commonly found in an educational testing context but are less often seen in clinical research settings where questionnaires are used for the measuring and assessment of individuals’ conditions. Due to the higher level of complexity in the study design of clinical studies, IRT models are often discarded for the analyses of questionnaires. Simple sum scores of the questionnaires, instead of the more robust estimates derived from IRT models, are taken directly to represent the quantity of the latent constructs. This has undesirable effects of introducing bias and overconfidence and ignoring uncertainty in measuring these latent constructs.

There is no need to trade measurement inefficiency for study design complexity. Complex models arising from complicated study designs can be naturally extended to incorporate IRT models that deal with latent construct measurement, thus enhancing the quality of the study. In this post, we will walk through such an example by building up a model for analyzing individuals’ changes—in terms of latent variables—over time. A Bayesian framework will be adopted, and we will be working with Stan code directly. As always, we begin with simulations. But since the model we are building is quite complicated—as any model trying to mirror reality is—I’ll first provide some context.

Treatment of Alcohol Dependence

Suppose some researchers were carrying out a study to examine the effectiveness of treatments on alcohol dependence. Individuals were recruited and randomly assigned to one of the three treatment conditions. All treatments lasted for three months, during which the participants were assessed for the effectiveness of the treatments to track their clinical progress. The assessments included collecting data such as measures of treatment outcomes (e.g., days of heavy drinking in the past 14 days) and mediating factors on treatment effectiveness suggested by previous studies (e.g., self-efficacy). It was hypothesized that treatments for alcohol dependence work partially through raising individuals’ self-efficacy in controlling alcohol use. Hence, during the assessments, questionnaires on self-efficacy in alcohol control were also administered to measure this potential mediating factor.

The DAG above explicates the causal assumptions of this fictitious study². Here’s the description of the variables in the DAG.

• $E$: Participants’ self-efficacy on alcohol use control

  Since self-efficacy $E$ is not directly observed, it is represented as a circled node in the DAG.

• $R$: Item responses collected through self-efficacy questionnaires

  To measure the unobserved self-efficacy $E$, tools like questionnaires are required to measure such a latent construct. $R$ stands for the responses collected through the questionnaire. These responses would allow the estimation of the variable $E$ for each participant. Note that the item parameter $I$ is left out for simplicity. If present, it would point to $R$ as item estimates also affect the responses $R$.

• $A$: Participants’ age

• $T$: Treatment condition received by a participant

• $D$: Latent treatment outcome

  $D$ is the latent quantity that underlies the (observed) treatment outcome.

• $D^{\ast}$: Treatment outcome. Here, it is the days of heavy drinking in the past 14 days.

The arrows among the nodes in the DAG indicate the directions of influence. Therefore, the graph is basically saying that the treatments affect the outcomes through two pathways. One direct, and the other indirect, through self-efficacy. Age also has direct influences on self-efficacy and the treatment outcome. The labels on the edges mark the regression coefficients, which are the parameters of interest for our later simulations and model testing.

The DAG presented above ignores the time dimension for simplicity. The second DAG below includes it. To avoid cluttering the graph, only three, instead of four, time points are shown. The subscripts on the variables mark the time points. $t=0$ indicates the baseline (i.e., the first) assessment. A cautionary note here is that age only directly influences self-efficacy and the latent treatment outcome at baseline ($A \rightarrow E_0, D_0$). At subsequent time points, they are influenced by age only indirectly through $E_0$ and $D_0$. This slight complication will become clearer in the following description of the simulation.

Simulating Treatments

With the background provided, let’s now simulate the above scenario. The code for the simulation, as well as later model fitting, will be inevitably long. It’s a natural consequence of realistic stats modeling. This post is thus less of a pedagogical step-by-step guide and more of a conceptual demonstration of practical real-world data analyses. That said, intricate details of code implementation won’t be hidden and are laid out as is. If anything seems mysterious, the key to demystifying it is to run some code. Experimentation is an ally that you should always trust.

Efficacy and Treatment Outcome

Let’s first sketch out the overall scenario by simulating 30 participants for each treatment condition. The age, gender, and the assigned treatment condition of each participant are saved in the variable A, G, and Tx, respectively. I also simulated a “baseline” for each of the participants. These baselines are the quantification of individual differences and can be more or less thought of as the amount of efficacy possessed by the individuals before receiving any treatment.

 1# remotes::install_github("liao961120/stom")
 2library(stom)
 3set.seed(1977)
 4
 5Ntx = 3        # number of treatments
 6Ns = Ntx * 30  # number of subjects
 7Nt = 4         # number of time points
 8
 9Tx = rep(1:Ntx, each = Ns / Ntx)  # Treatment condition for each subj
10G = rep(1:2, times = Ns / 2)      # Gender of each subj
11A = rtnorm( Ns, m = 36, lower = 18, upper = 80, s = 20 )  # Age
12tau = 0.8                 # std for Subject baseline Efficacy
13subj = rnorm(Ns, 0, tau)  # Subject baseline Efficacy (subject intercept)
14
15# Transform Age to a reasonable scale
16minA = min(A)
17A = (A - minA) / 10  # 1 unit = 10 years in original age

The simulation of the participants’ age deserves some elaboration. These ages are drawn from a truncated normal distribution, which is essentially a normal distribution with an upper and a lower bound applied. The boundaries are set here such that the participants come from a normal distribution with a mean age of 36 and a standard deviation of 20, with ages below 18 or over 80 left out from the sample³. After drawing samples for the ages, the ages are scaled such that (1) the youngest participant has a scaled age of zero, and (2) a unit of difference in the scaled age corresponds to a 10-year difference in the raw age. This is done to align the scale of the age to the scales of the other variables. Otherwise, the original large scale of the age would make the effects difficult to interpret.

Now, let’s lay out the parameters for generating participants’ self-efficacy ($E$) and treatment outcomes ($D$):

 1B_AE = .1             # Effect of Age on Efficacy
 2B_TE = matrix(c(      # Effect of Treatment on Efficacy
 3      .3, .7,  1.3,
 4      .3,  1,   .7
 5    ), byrow=T, nrow=2
 6)
 7B_AD = .2             # Effect of Age on Outcome
 8B_ED = 1              # Effect of Efficacy on Outcome
 9B_TD = c(0, 0, 0)     # Effect of Treatment on Outcome
10
11delta = -1.8   # Global Intercept for "E model"
12alpha = 0      # Global Intercept for "D model"

As mentioned previously, it is assumed that the treatments differentially affect efficacy according to gender (i.e., the treatment effect on efficacy interacts with gender). This is specified in B_TE ($\beta_{TE}$) as a 2-by-3 matrix. The first row of B_TE corresponds to male and the second to female. The three columns correspond to Treatment 1 (Tx == 1), 2 (Tx == 2), and 3 (Tx == 3) respectively. For the direct treatment effects on the outcomes ($\beta_{TD}$), I set them all to zero and assume no interaction with the gender, to keep things simple.

With all these parameters prepared, we can build up the generative process of efficacy and the outcomes.

As shown in the code chunk below, we first set up the time points t. t starts with zero because it is assumed that the first assessment took place right before the treatments started. Hence, any treatment effects at this point of measure should be zero. t = 0 does us the favor as it naturally cancels out the treatment effects $\beta_{TD}$ and $\beta_{TE}$.

Let’s now simulate self-efficacy $E$. Before receiving the treatments, the efficacy is assumed to be slightly influenced by age, through B_AE. In addition, baseline individual differences are modeled through the random subject intercept subj. After receiving the treatments, treatment effects get added on through the term B_TE. Efficacy is then modeled as the sum of the aforementioned effects.

After generating participants’ efficacy, we can again follow the arrows of the DAGs to write down the generative process of $D^{\ast}$, coded as D_latent below.

Finally, the treatment outcomes here, days of heavy drinking in the past 14 days, are assumed to follow a binomial distribution $D \sim \text{Binomial}(14, p)$. The latent score D_latent is reversed with the negative sign since we want higher latent scores to result in fewer heavy drinking days. Next, we simply map the latent scores to probabilities with the logistic function.

 1t = 0:(Nt - 1)  # time points of measure
 2
 3# E model (causes of E)
 4E = sapply(t, function(time) {
 5    b_TE = sapply( 1:Ns, function(i) B_TE[ G[i],Tx[i] ] )
 6    delta + subj + B_AE * A + b_TE * time
 7})
 8
 9# D model (causes of D)
10D_latent = sapply(t, function(time) {
11   alpha + B_TD[Tx]*time + B_AD * A + B_ED * E[, time + 1]
12})
13D = sapply(t, function(time) {
14   mu = -D_latent[, time+1]
15   rbinom( Ns, size=14, prob=logistic(mu) )
16})

Now, let’s collect the values we have simulated into a long-form data frame. Each row of the data frame corresponds to a response ($D$) of a participant collected at a specific time point. I will name this data frame the (treatment-)outcome-level data frame. Later, we will see another data frame that records item-level responses.

 1# Outcome-level responses (subject-time)
 2dO = expand.grid( Sid=1:Ns, time=t, KEEP.OUT.ATTRS=F )
 3dO$A = with( dO, A[Sid] )
 4dO$Tx = with( dO, Tx[Sid] )
 5dO$G = with( dO, G[Sid] )
 6dO$D = NA
 7dO$D_latent = NA
 8dO$E = NA
 9for ( i in 1:nrow(dO) ) {
10   s = dO$Sid[i]
11   t_ = dO$time[i] + 1
12   dO$E[i] = E[s, t_]
13   dO$D_latent[i] = D_latent[s, t_]
14   dO$D[i] = D[s, t_]
15}
16str(dO)

'data.frame':   360 obs. of  8 variables:
 $ Sid     : int  1 2 3 4 5 6 7 8 9 10 ...
 $ time    : int  0 0 0 0 0 0 0 0 0 0 ...
 $ A       : num  1.57 2.58 1.21 0.92 0 ...
 $ Tx      : int  1 1 1 1 1 1 1 1 1 1 ...
 $ G       : int  1 2 1 2 1 2 1 2 1 2 ...
 $ D       : int  13 9 6 11 12 13 7 12 13 12 ...
 $ D_latent: num  -2.534 -1.132 0.286 -0.942 -2.003 ...
 $ E       : num  -2.8466 -1.6485 0.0451 -1.1261 -2.0029 ...

Item Responses

Item responses are generated similarly to $D$ since they both depend on efficacy $E$. Let’s assume the questionnaire for measuring participants’ efficacy contains 20 items (Ni = 20), and that the easiness of the items is equally spaced and ranges from -.6.3 to 6.3.

With the item easiness and a baseline choice preference kappa set, we now can generate the item responses $R$ from the participants’ efficacy.

 1Ni = 20  # number of items
 2I = seq(-4.5, 4.5, length = Ni)  # item easiness (sums to zero)
 3kappa = logit(cumsum(simplex(c(1, 2, 4, 4, 2, 1))))
 4kappa = kappa[-length(kappa)]
 5
 6# Item-level responses (subject-item-time)
 7dI = expand.grid( Sid=1:Ns, Iid=1:Ni, time=t, KEEP.OUT.ATTRS=F )
 8for (i in 1:nrow(dI)) {
 9   dI$R[i] = with(dI, {
10      rordlogit(E[Sid[i], time[i] + 1] + I[Iid[i]], kappa = kappa)
11   })
12}
13str(dI)

'data.frame':   7200 obs. of  4 variables:
 $ Sid : int  1 2 3 4 5 6 7 8 9 10 ...
 $ Iid : int  1 1 1 1 1 1 1 1 1 1 ...
 $ time: int  0 0 0 0 0 0 0 0 0 0 ...
 $ R   : int  1 1 1 1 1 1 1 1 1 1 ...

Wrapping up

When a simulation gets massive, it is always better to pack up all the code into a function. There will certainly be times when we have to modify the structure of, or the parameter values in, the simulation. A function helps a lot in these situations.

The simulation code we have gone through so far is condensed into the sim_data() function in simulation.R. You can just load it and run the simulation:

1source("simulation.R")
2d = sim_data()  # Change param vals by overwriting default args
3str(d)

List of 3
 $ dat   :List of 18
  ..$ Ns      : num 90
  ..$ Ntx     : num 3
  ..$ Nt      : num 4
  ..$ Nk      : num 6
  ..$ Ni      : num 20
  ..$ NI      : num 7200
  ..$ Sid_I   : int [1:7200] 1 2 3 4 5 6 7 8 9 10 ...
  ..$ Iid_I   : int [1:7200] 1 1 1 1 1 1 1 1 1 1 ...
  ..$ time_I  : int [1:7200] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ R       : int [1:7200] 1 1 1 1 1 1 1 1 1 1 ...
  ..$ NO      : num 360
  ..$ Sid_O   : int [1:360] 1 2 3 4 5 6 7 8 9 10 ...
  ..$ time_O  : int [1:360] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ G       : int [1:360] 1 2 1 2 1 2 1 2 1 2 ...
  ..$ A       : num [1:360] 1.947 5.137 0.634 2.086 5.84 ...
  ..$ Tx      : int [1:360] 1 1 1 1 1 1 1 1 1 1 ...
  ..$ D       : int [1:360] 5 9 13 6 8 14 6 14 9 11 ...
  ..$ D_latent: num [1:360] -0.534 -2.038 -1.235 -0.354 -0.725 ...
 $ params:List of 13
  ..$ alpha  : num 0
  ..$ delta  : num -1.8
  ..$ B_AE   : num 0.1
  ..$ B_TE   : num [1:2, 1:3] 0.3 0.3 0.7 1 1.3 0.7
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "male" "female"
  .. .. ..$ : NULL
  ..$ B_AD   : num 0.2
  ..$ B_ED   : num 1
  ..$ B_TD   : num [1:3] 0 0 0
  ..$ E      : num [1:90, 1:4] -0.923 -3.065 -1.362 -0.772 -1.893 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:90] "male" "female" "male" "female" ...
  .. .. ..$ : NULL
  ..$ I      : num [1:20] -4.5 -4.03 -3.55 -3.08 -2.61 ...
  ..$ sigma_I: num 2.8
  ..$ kappa  : num [1:5] -2.56 -1.3 0 1.3 2.56
  ..$ tau    : num 0.8
  ..$ subj   : num [1:90] 0.682 -1.779 0.374 0.82 -0.677 ...
 $ others:List of 3
  ..$ minA    : num 18.5
  ..$ D       : int [1:90, 1:4] 5 9 13 6 8 14 6 14 9 11 ...
  ..$ D_latent: num [1:90, 1:4] -0.534 -2.038 -1.235 -0.354 -0.725 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:90] "male" "female" "male" "female" ...
  .. .. ..$ : NULL

Notice the structure of the returned value by sim_data(). It is a nested named list with three elements at the top level: $dat, $params, and $others.

• $dat holds the data for model fitting. It has this specific structure to match stan’s data block. It makes more sense later when we see stan’s counterpart of the code.

• $params are the true parameter values used in the simulation. They are returned so we can later compare the true parameter values to those estimated by the model.

• $others holds other information that does not belong to the former two but is nevertheless needed. For instance, the minimum age minA is saved in $others so that raw ages in new data can be scaled accordingly if one would like to make predictions with new datasets.

Model Formulation

Now we have the simulated data, let’s begin constructing the model. In a Bayesian framework, model construction feels pretty much like a simulation (assume you are using Stan, not some wrappers around it) since both are essentially describing the same data-generating process. Bayes, as complicated as it may seem, really just comes down to four simple components:

1. Data-generating process
2. Prior distribution
3. Data
4. Posterior distribution

The data-generating process describes the relations between the parameters, that is, how the parameters arise from others. The priors give preliminary information on the parameters (the priors could be thought of as part of the data-generating process as well). With these prepared, the Bayesian machine incorporates the information from the data to update the priors—according to the data-generating process—and arrives at the posterior distribution.

The posterior is nothing else but a joint probability distribution of all parameters, after combining the information from the priors, the data-generating process, and the data. It is what we believe the parameters are, given the data and our assumptions. Everything we need for inference derives from the posterior.

So we’re now ready to formulate our model (data-generating process) in Bayesian terms. This massive model can be chopped into three parts. Let’s divide and conquer.

Causes of Efficacy

The first part of the model deals with the generation of Efficacy. The formulation is shown below.

The first two lines of the formulation describe how efficacy $E$ is generated from the effects of participant baseline $\gamma_{Sid}$, age $\beta_{AE}$, and treatment $\beta_{TE}$. The last three lines specify the prior distributions for the relevant parameters. The “+” sign on the prior distribution of $\tau$ is an abbreviation for the truncated normal distribution with negative values removed. Since $\tau$ is a standard deviation, it must not be negative. The truncation at zero imposes this constraint on $\tau$.

Causes of Treatment Outcome

With $E$ generated, we can then describe the process generating the treatment outcome $D^{\ast}$. The outcomes in the current example are the number of heavy-drinking days in the past 14 days. Hence, it is a count variable with an upper bound of 14, which makes the binomial an ideal distribution for modeling the outcome. The binomial and the logit link are used here to map the count outcomes to the underlying latent variable $D$. As seen in the simulation, a negative sign is needed here to reverse the relation between $D^{\ast}$ and $D$ such that a higher value of $D$ corresponds to a fewer count of heavy-drinking days. The third formula lays out the generative process of $D$, linking the influences of the treatment, age, and efficacy to $D$.

A special kind of prior distribution—exemplified by the distribution of $\beta_{TD}$ here—deserves some attention. The distribution is assumed to be a normal distribution with an unknown mean $\mu_{\beta_{TD}}$ and unknown standard deviation $\sigma_{\beta_{TD}}$. Their priors are then given by the next two lines. This nested structure of prior distributions is known as hierarchical priors. By specifying such a hierarchical structure in the priors, one can partial-pool information across the three treatment effects, thus reducing variation among the $\beta_{TD}$ parameters.

Partial pooling is used here for two reasons. First, it makes sense theoretically to partial-pool the direct treatment effects on the outcome since different treatments are still similar in some sense that information across all treatments is useful and should be shared. Second, partial pooling has the effect of reducing variation among the pooled parameters. Without such a soft constraint on the $\beta_{TD}$ parameters, the model here will be unidentified (i.e., multiple sets of solutions exist).

Causes of Item Response

The final piece of the model describes the generation of participants’ responses to the items measuring efficacy. The submodel here has a structure nearly identical to the rating scale model described in the previous post. The only difference here is that a sum-to-zero constraint is imposed on the item parameters, which allows the item easiness to be anchored and thus comparable to the true parameter values from the simulation. The item parameters are also partially pooled here, achieved through the hierarchical prior on $\sigma_I$.

Stan: A Brief Introduction

We have completed our formal model-building. Now, we have to code it in Stan, a language for implementing Bayesian models.

Stan is hard for anyone first meeting it. Sadly, I think there is really no “introductory” material on Stan because Bayesian statistics is already an advanced topic. How, then, could anyone get started with Bayes if everything is advanced and hard?

We need scaffolds. The scaffolds for learning statistics are simulations and programming skills that license them. Building up simulations allows one to experiment with stat models and therefore provides opportunities for understanding. So to learn Stan, one should probably begin with examples. The Stan User’s Guide provides tons of example models. Start with simple models or any one that you’re familiar with. Simulate data based on the model’s assumption, construct the model by modifying the example code, and fit the model to see if it correctly recovers the parameters. If succeeded, extend the simulation and the model for the next round of testing. If failed, simplify the model (and the simulation, if needed) to pinpoint the potential problems.

Interestingly, the process of learning Stan is no different from using Stan. In both situations, one is unable to proceed without scaffolds. Stan gives you full flexibility, and with this modeling power, one immediately finds out there are infinite ways to go wrong. Simulation is the compass that keeps us oriented while navigating through the mess of modeling. That said, some acquaintance with the structure of Stan files does help when getting started.

Stan File Structure

The template below sketches the most common structure of a Stan program. A Stan program consists of several “blocks”. Three of the most basic blocks are shown here: the data, the parameters, and the model blocks. The data block simply contains variable-declaration code for the data. The parameters block declares variables that hold the parameters. The model block is where the data-generating process and the priors get specified.

1data {
2   // Declare varaiables for data (passed in from R)
3}
4parameters {
5   // Declare model parameters
6}
7model {
8   // Describe data-generating process (including priors)
9}

With the blocks written up, one then compiles the stan program and feeds it the data. Stan would then construct and sample from the posterior. When all of these are done, posterior samples along with diagnostic information of sample quality are returned.

That’s it. Don’t read too much while learning Stan. Getting the hands dirty is a much better way. So now let’s proceed to our massive model.

Stan Model

I’ll present our model coded in Stan one block at a time so that we don’t get overwhelmed. The stan file gets executed from the top to the bottom, and the blocks have to be defined in the given order as shown in the above template. Through this structure, the parameters and the model blocks have access to the data variables defined in the data block, and the model block has access to the parameter and data variables defined in previous blocks. Let’s first look at the data block of our model.

The data block

 1data {
 2    int Ns;  // num of subjects
 3    int Ntx; // num of treatments
 4    int Nt;  // num of time points
 5    int Nk;  // num of Likert scale choices
 6    int Ni;  // num of items in the self-efficacy scale
 7
 8    // Item-level responses (NI=Ns*Ni*Nt)
 9    int NI;
10    array[NI] int<lower=1,upper=Ns> Sid_I;     // Subject ID
11    array[NI] int<lower=1,upper=Ni> Iid_I;     // Item ID
12    array[NI] int<lower=0,upper=Nt-1> time_I;  // time point of obs.
13    array[NI] int<lower=1,upper=Nk> R;         // Responses on Efficacy scale
14
15    // Outcome-level responses (NO=Ns*Nt)
16    int NO;
17    array[NO] int<lower=1,upper=Ns> Sid_O;     // Subject ID
18    array[NO] int<lower=0,upper=Nt-1> time_O;  // time point of obs.
19    array[NO] real<lower=0,upper=20> A;        // Age scaled: (A-min(A))/10 
20    array[NO] int<lower=1> Tx;                 // Treatment received
21    array[NO] int<lower=1,upper=2> G;          // Gender
22    array[NO] int<lower=0,upper=14> D;         // Binomial outcome (heavy-drinking in past 14 days)
23}

In our data block, we simply define all variables that hold information about the data we pass in through R. Thus, the variables here correspond exactly to the $dat element of the simulated data list.

Stan is a statically typed language, meaning that the type must be specified for any variables defined. In addition, boundaries could also be imposed on the typed variables. This is done through the angle brackets <> after the type definition, as shown in some of the variables above. For the data block, the constraint on the boundaries guards against unexpected input data format. For the parameters block, these constraints are often used to truncate the prior distributions later assigned to the parameters.

parameters and transformed parameters

 1parameters {
 2    // IRT model params
 3    real<lower=0> sigma_I;
 4    vector[Ni-1] zI_raw;
 5    ordered[Nk-1] kappa;
 6
 7    // E model params
 8    matrix[2,Ntx] B_TE;  // Treatment on Efficacy (indirect effect)
 9    real B_AE;           // Age on Efficacy
10    real delta;          // global intercept (E linear model)
11    real<lower=0> tau;   // subj baseline std (E linear model)
12    vector[Ns] zSubj;    // subject baseline
13    
14    // D model params
15    vector[Ntx] zB_TD;       // Treatment on Outcome (direct effect)
16    real mu_TD;              // Common mean among direct treatment effects 
17    real<lower=0> sigma_TD;  // Common std among direct treatment effects
18    real B_AD;               // Age on outcome
19    real B_ED;               // Efficacy on outcome
20    real alpha;              // global intercept (D linear model)
21}
22transformed parameters {
23    // IRT item params (sum-to-zero contrained)
24    vector[Ni] I = sigma_I * append_row( -sum(zI_raw), zI_raw );
25
26    // subject baseline efficacy
27    vector[Ns] subj;
28    subj = tau * zSubj;
29
30    matrix[Ns,Nt] E;
31    // Transformed E
32    for ( i in 1:NO ) {
33        int sid = Sid_O[i];
34        int time = time_O[i];
35        E[sid,time+1] = delta + subj[sid] + B_AE*A[i] + B_TE[G[i],Tx[i]]*time;
36    }
37
38    // Direct treatment effect
39    vector[Ntx] B_TD = fma(zB_TD, sigma_TD, mu_TD);
40}

The parameters block defines variables for the parameters in a stan model. This should be enough, theoretically. However, Bayesian models are often very complex (such as our massive model). This can lead to complications in posterior geometries such that Stan’s algorithm would have difficulty exploring the distribution. To deal with this issue, it is often required to apply mathematical tricks known as reparameterization to modify the posterior such that it becomes easier for Stan to sample from.

This is where the transformed parameters block comes into play. In general, when reparameterizing our model, we define those reparameterized parameters in the parameters block. These parameters, however, are often in forms that are unintuitive for us. We thus utilize the transformed parameters block to transform those parameters back into the forms we are familiar with. The transformed parameters get returned from Stan as if they are parameters after sampling. But they are not actual parameters per se but only functions of parameters. That said, once we have the transformed parameters block set up and the sampler working, it’s safe to simply treat them all as parameters. Understanding the details of the transform block helps when we have trouble coding the model.

Here, we utilize the transform block to convert the non-centered parameters (zB_TD and zE) back to the familiar centered forms (subj, B_TD and E). In addition to non-centered reparameterizations, another transformation is performed here to reconstruct the full vector of the sum-to-zero constrained item easiness parameters.

model

 1model {
 2    // Priors for IRT parameters
 3    zI_raw ~ std_normal();
 4    sigma_I ~ exponential(1);
 5    kappa ~ std_normal();
 6
 7    // Priors for causes of E (T -> E)
 8    to_vector(B_TE) ~ normal(0, 1.5);
 9    B_AE ~ std_normal();
10    delta ~ normal(0, 1.5);
11    tau ~ std_normal();
12    zSubj ~ std_normal();
13
14    // Priors for causes of D (T -> D <- E)
15    zB_TD ~ std_normal();
16    mu_TD ~ std_normal();
17    sigma_TD ~ std_normal();
18    B_AD ~ std_normal();
19    alpha ~ normal(0, 1.5);
20    B_ED ~ std_normal();
21
22    // Causes of E (see transformed parameters block)
23    
24    // Causes of D
25    vector[NO] mu;
26    for ( i in 1:NO ) {
27        int sid, time;
28        sid = Sid_O[i];
29        time = time_O[i];
30        mu[i] = alpha + B_TD[Tx[i]]*time + B_AD*A[i] + B_ED*E[sid,time+1];
31    }
32    D ~ binomial_logit( 14, -mu );
33
34    // Measurement model (IRT)
35    vector[NI] phi;
36    for ( i in 1:NI )
37        phi[i] = E[Sid_I[i],time_I[i]+1] + I[Iid_I[i]];
38    R ~ ordered_logistic( phi, kappa );
39}

We are finally at the model block, often the most terrifying part of a stan program. It becomes less terrifying if one knows that the code in the model block can be classified into two functional groups. The first is the sampling statements. These statements tell Stan to update the posterior distribution, with information from the priors, the data, or the combination of both. The posterior gets updated by adding values to itself. This means that the order of the sampling statements does not influence the output (ignoring differences resulting from numerical instability). In general, although it is legitimate to switch the order of the sampling statements, we usually organize them in ways that make our modeling logic easier to grasp.

The function of the second type of code in the model block is computation. Often, before reaching a sampling statement, some intermediate computation is required. For instance, D ~ binomial_logit( 14, -mu ) depends on the variable mu. So mu has to be computed before one can invoke the sampling statement.

For our massive model here, I first define the prior distributions in the first three blocks, with those related grouped into the same block of code. I then proceed to specify the relations between the parameters and the data. For efficacy $E$, the relevant relationships are already defined in the transformed parameters block, so only a line of comment is left as a placeholder. Below $E$ is the block for the “causes of D”, which defines the upstream variables influencing D (through mu). Finally, the last block defines the generative process for the item responses.

R Interface

There are interfaces written in many interpreted languages built to interact with Stan. In R, the most commonly used are rstan and cmdstanr. I will be using cmdstanr.

Here’s the code for calling Stan to fit our model. The model takes about 5 minutes to run on an M1 Macbook Air and 16 minutes on a much older Windows laptop. So it would be wise to save the fitted model once Stan finishes sampling.

 1source("simulation.R")
 2set.seed(1977)
 3d = sim_data()
 4
 5# Compile stan file
 6m = cmdstanr::cmdstan_model("m1.stan")
 7
 8# Fit model
 9fit = m$sample(
10    data=d$dat,
11    chains=3, 
12    parallel_chains=3
13)
14
15# Save fitted model for later use
16fit$save_object("m1.RDS")

I’ve put the fitted model m1.RDS on Google Drive so you don’t have to refit the model. All the scripts in this post are also available on GitHub.

Model Checking

With our model fitted, we can now check if it correctly recovers the parameter. But before doing so, we need to first assess the quality of the posterior samples.

Quality of Posterior Samples

When fitting Bayesian models, at least two chains (i.e., two independent sets of samples) are required to evaluate the quality of the posterior samples. These chains could be visualized in trace plots, which record the parameter values in the order they are sampled. Thus, a curve in a trace plot is essentially the trace of a sampler exploring the posterior (in one of the dimensions). For high-quality sampling, traces of different samplers should be independent and fluctuate at similar ranges of values. This results in chains that mix well, and the trace plots would then look like caterpillars. Conversely, if the posterior geometries are difficult to explore, the samplers may be stuck in local regions. Or, they may explore different regions of the posterior (e.g., when the posterior is multi-modal), resulting in a complete separation of the chains in the trace plots. When these happen, the chains may disagree with each other, leading to the non-convergence of samples from different chains.

Non-convergence indicates that inferences based on the posterior samples are problematic and shouldn’t be trusted. Upon such circumstances, the modeler should try to locate the sources of non-convergence and fix the model. Misspecification of the generative process, loose priors, unidentifiability, and unexpected data are all potential causes of non-convergence.

To evaluate how well the chains mix, we look at a statistic known as $\hat{R}$ (rhat). Typically, an $\hat{R}$ value less than 1.05 for a particular parameter indicates a good enough mixing of chains. Let’s load our fitted model and look at the sampling quality. I will use some handy wrappers around the functions from cmdstanr through my package stom.

1source("simulation.R")
2set.seed(1977)
3d = sim_data()
4
5m = readRDS("m1.RDS")   # Load fitted model
6s = stom::precis(m, 5)  # Compute posterior summary
7head(s)

# A tibble: 6 x 8
  variable    mean    sd    q5    q95  rhat ess_bulk ess_tail
  <chr>      <dbl> <dbl> <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 sigma_I    2.56  0.379  2.02  3.26   1.01     330.     543.
2 zI_raw[1] -1.68  0.248 -2.11 -1.29   1.01     349.     567.
3 zI_raw[2] -1.43  0.208 -1.79 -1.10   1.01     360.     565.
4 zI_raw[3] -1.17  0.175 -1.48 -0.891  1.01     349.     525.
5 zI_raw[4] -1.05  0.156 -1.31 -0.802  1.01     355.     635.
6 zI_raw[5] -0.826 0.126 -1.04 -0.626  1.01     354.     698.

The function stom::precis computes parameter summaries from the posterior samples, including the $\hat{R}$ statistic. Let’s plot the $\hat{R}$ values for all of the parameters to see if any of them are unexpectedly large.

1plot(s$rhat)

The $\hat{R}$ values look fine. The largest one is 1.0106079, which is much lower than 1.05. Just to be careful, let’s inspect the trace plots of those parameters with the largest $\hat{R}$ values and the lowest effective sample size (ESS).

1library(bayesplot)
2color_scheme_set("viridis")
3pars = c(
4   s$variable[s$rhat > 1.006],
5   s$variable[s$ess_bulk < 400],
6   s$variable[s$ess_tail < 400] 
7)
8mcmc_trace(m$draws(), pars=pars)

The trace plots seem fine too. Now we can proceed to check the recovery of the parameters.

Parameter Recovery

Let’s first look at how well the IRT submodel works. I plot the estimated (y-axis) against the true (x-axis) parameter values. The grey dashed line shows the identity. Points lying on this line indicate perfect recovery through the posterior mean. The vertical pink bars cover 90% of posterior density around the mean.

 1par( mfrow = c(1, 2) )
 2for ( p in c("I", "kappa", "E") ) {
 3    if ( p == "E" ) par(mfrow = c(1, 1))    
 4    mt = d$params[[p]]        # True value
 5    mm = get_pars(s, p)$mean  # Posterior mean
 6    u = get_pars(s, p)$q95    # Posterior .05 quantile 
 7    l = get_pars(s, p)$q5     # Posterior .95 quantile 
 8    plot( mt, mm , main=p, ylim=c(-4.6,4.6),
 9          xlab="True", ylab="Estimated")
10    for ( i in seq_along(mm) )
11        lines( c(mt[i],mt[i]), c(u[i],l[i]), lwd=3, col=col.alpha(2) )
12    abline(0,1, lty="dashed", col="grey")
13}

Let’s next look at the parameters for modeling the causes of $E$ and $D$. This time I plot the true parameter values as solid black dots. The pink open circles are the posterior mean estimates, and the bars mark the 90% density around the mean.

 1beta = c( "B_TE", "B_AE", "B_AD", "B_ED", "B_TD", "delta", "alpha", "tau", "sigma_I" )
 2b_true = lapply( beta, function(p) d$params[[p]] ) |> unlist()
 3b_est = lapply( beta, function(p) get_pars(s, p)$mean ) |> unlist()
 4b_est_upp = lapply( beta, function(p) get_pars(s, p)$q5 ) |> unlist()
 5b_est_low = lapply( beta, function(p) get_pars(s, p)$q95 ) |> unlist()
 6
 7plot( b_true, pch=19, ylim=c(-2.6, 3.2), ylab="Parameter value", xlab="" )
 8abline(h = 0, lty="dashed", col="grey")
 9points( b_est, col=2 )
10for ( i in seq_along(b_est) )
11    lines( c(i,i), c(b_est_upp[i],b_est_low[i]), lwd=4, col=col.alpha(2) )
12for ( v in c(6:9,12) )
13    abline( v=v+.5, lty="dashed", col="grey" )
14for ( v in 7:9 )
15    mtext( beta[v-5], at=v )
16mtext( beta[1], at = length(get_pars(d,"B_TE")) * .5 )
17mtext( beta[5], at = 11 )
18for ( i in 6:9 )
19    mtext( beta[i], at = i+7 )

With all these checks, we can see that the model does pretty well in recovering the parameters! So now, we’re ready to utilize the posterior samples for inference.

Posterior Predictions

Bayesian models are generative in the sense that they can generate samples, much like nature generates stochastic observations of various kinds. This can be done because, instead of providing point estimates, a Bayesian model returns a full distribution of parameter estimates (i.e., the posterior distribution). A distribution contains information of uncertainty. Thus, the samples computed from the posterior also retain such uncertainty.

Given the generative nature of Bayesian models, we can utilize the posterior distribution to compute various predictions. These predictions carry forward the uncertainty of the posterior and are thus known as posterior predictive distributions. Based on the purpose of the inference, we can design how the predictions are computed. Here, I provide two examples. The first computes the predicted trajectory of the clinical outcomes for all six treatment-gender combinations. The second compares the effect of the treatments in terms of the expected reduction in days of heavy drinking relative to the baseline treatment. Since everything derived from the posterior is a distribution, we have to plot, rather than tabulate, our predictions.

Before starting, let me first demonstrate the computation of a single prediction from the posterior. The rest should be more straightforward since posterior predictive distributions are simply the generalization from a single prediction to multiple predictions using all posterior samples.

A Single Posterior Prediction

1post = stom::extract2(m)
2str(post$delta())

 num [1:3000] -2.51 -2.12 -1.76 -2.16 -1.86 ...

1post$delta(idx = c(2, 4, 6))

[1] -2.11809 -2.15574 -1.69995

stom::extract2() extracts the posterior and exposes its parameters as R6 methods. To access the values of a parameter, we use the syntax ${param}(). In our example here, doing so would return 3000 samples since 1000 samples were drawn for each of the 3 chains from our model. We can control the samples drawn by passing an optional argument idx into the parameter method. For instance, post$delta(2) would return the second set of posterior samples.

What I would like to do here is to compute the latent scores $D$ underlying the treatment outcome for a single participant at time 1, 2, 3, and 4. I will use subject 31 (Sid == 31) as an example here. So let’s first prepare the data.

1Sid = 31  # subject ID
2idx_subj = d$dat$Sid_O == Sid
3A = d$dat$A[idx_subj][1]     # Age of subject 31
4G = d$dat$G[idx_subj][1]     # Gender of subject 31
5Tid = d$dat$Tx[idx_subj][1]  # Treatment received by subject 31
6c(A, G, Tid)

[1] 0.987944 1.000000 2.000000

Now, let’s extract 1 of the 3000 sets of the posterior samples.

1idx = 2999  # posterior sample drawn
2alpha = post$alpha(idx)
3delta = post$delta(idx)
4b_AD  = post$B_AD(idx)
5b_ED  = post$B_ED(idx)
6b_AE  = post$B_AE(idx)
7b_TD  = post$B_TD(idx)
8b_TE  = post$B_TE(idx)
9subj  = post$subj(idx)

With the parameters and the data prepared, we can compute our prediction.

1time = 1:4
2D = sapply( time, function(t) {
3    E = delta + subj[Sid] + b_AE*A + b_TE[G,Tid]*(t-1)
4    alpha + b_TD[Tid]*(t-1) + b_AD*A + b_ED*E
5})
6D

[1] -0.8849468 -0.1329965  0.6189538  1.3709041

Okay, now we have a prediction from a single set of posterior samples. In general, we would like to have predictions based on all the posterior samples. To do this, we simply wrap the code above in a function and repeat it for all of the posterior samples.

 1predict_obs = function(Tid, Sid, A, G, idx) {
 2    alpha = post$alpha(idx)
 3    delta = post$delta(idx)
 4    b_AD  = post$B_AD(idx)
 5    b_ED  = post$B_ED(idx)
 6    b_AE  = post$B_AE(idx)
 7    b_TD  = post$B_TD(idx)
 8    b_TE  = post$B_TE(idx)
 9    subj  = post$subj(idx)
10    
11    time = 1:4
12    D = sapply( time, function(t) {
13        E = delta + subj[Sid] + b_AE*A + b_TE[G,Tid]*(t-1)
14        alpha + b_TD[Tid]*(t-1) + b_AD*A + b_ED*E
15    })
16    D
17}
18
19# A single prediction based on the 2999th set of posterior samples
20( predict_obs(Tid, Sid, A, G, idx=2999) )

[1] -0.8849468 -0.1329965  0.6189538  1.3709041

1# Compute over the full posterior
2D = sapply( 1:3000, function(i) {
3    predict_obs(Tid, Sid, A, G, idx=i)
4})
5str(D)

 num [1:4, 1:3000] -0.319 0.397 1.114 1.83 -0.891 ...

Counterfactual Predictions

In the demonstration above, we compute the predictions based on one of the observations in our data, that is, the data from subject 31. We used his gender, age, and treatment condition as the data, and his estimated baseline efficacy subj[Sid], to compute the prediction. But we can do something more interesting by going beyond our current data to compute counterfactual predictions. For instance, we can arbitrarily determine the treatment received and then compute the prediction. Doing so effectively asks the question: what would the outcome be if subject 31 had been treated with Treatment 3 instead of 2? Likewise, the age and gender of the participant can be changed to ask similar questions. We can go even further to simulate a completely new population, with a different set of baseline subject efficacy parameters, if we have some idea of what such a population is like (e.g., the degree of variation in baseline conditions).

Another thing worth noting is the time variable. In our data and the above prediction, time is treated as discrete time points (1-4). But we can treat time as continuous when we compute predictions, simply by using a smaller time step, such as time = seq(1,4,.01). Doing so results in a smoother trajectory—as we’ll see later—for the predictions across time.

Visualizing Treatment Effects

Now, let’s use the full posterior to compute the posterior predictive distributions. The posterior predictive distributions computed here are counterfactuals: all 90 subjects are used repeatedly for computing each of the six treatment-gender combinations. The treatment received and gender of the subjects are set to match the group’s condition when computing predictions for that combination.

The functions for plotting and calculating the relevant quantities are available in the script func_post_predict.R. The manipulation of the posterior samples in the script may seem sophisticated. But they are all based on the core function predict_obs(), which is nearly the same as its simplified version presented previously.

Trajectories By Treatment and Gender

The six charts below plot the trajectories of treatment outcomes for each of the treatment-gender conditions. The charts on the left (blue) are the trajectories for male participants, and those on the right (red) are the trajectories for female participants. The three rows counting from top to bottom correspond to the three treatment conditions, respectively.

Note that these trajectories are expected days of heavy drinking. Since we are modeling the outcomes as counts from the binomial distribution, the “expectation of the counts” is the product of $N$ and $P$, i.e., 14 * logistic(-D), where D is the latent quantity computed from predict_obs(). The means of these expected counts computed from the full posterior are shown as the thick curves in the charts below. The shaded regions cover 90% density of the expected trajectory distributions. The grey lines plot the empirical observation under that condition⁴.

 1source("func_post_predict.R")
 2post = stom::extract2(m)
 3
 4mar = par("mar")
 5mar[1] = 2.5
 6mar[3] = 1
 7mar[4] = 1.5
 8par( mfrow = c(3, 2), mar = mar )
 9plot_model_prediction(Tid=1, G=1, col_main=4)  # Treatment 3
10plot_model_prediction(Tid=1, G=2, col_main=2)  # Treatment 3
11plot_model_prediction(Tid=2, G=1, col_main=4)  # Treatment 2
12plot_model_prediction(Tid=2, G=2, col_main=2)  # Treatment 2
13plot_model_prediction(Tid=3, G=1, col_main=4)  # Treatment 1
14plot_model_prediction(Tid=3, G=2, col_main=2)  # Treatment 1

From these charts, we can see that Treatment 3 works best for male participants. For females, however, it is Treatment 2 that is the best. These differences might not be apparent enough in the form of trajectories. So next, let’s directly compare the treatments by inspecting their contrasts.

Treatment Contrasts

The code below computes the relevant contrasts for later uses. Recall that all 90 participants are used repeatedly for each treatment-gender condition. Since 3000 is indivisible by 90, I discarded 30 posterior samples and used the first 2970 samples just to keep things simpler. A more rigorous but quite computationally expensive alternative is to use all 3000 instead of 33 (2970/90) posterior samples for computing predictions for a single participant under each treatment-gender condition.

 1Sids = rep( 1:d$dat$Ns, length=2970 ) |> sample()
 2times = seq(1, 4, .01)
 3# table(Sids)
 4
 5cf = list(
 6    T1 = list( G1=c(), G2=c(), all=c() ),
 7    T2 = list( G1=c(), G2=c(), all=c() ),
 8    T3 = list( G1=c(), G2=c(), all=c() )
 9)
10for ( tx in 1:3 ) {
11    for ( g in 1:2 )
12        cf[[tx]][[g]] = sapply( 1:2970, function(i) {
13            D = predict_obs(Tid=tx, Sid=Sids[i], A=NULL, G=g, idx=i) 
14            14 * logistic(-D)
15        })
16    cf[[tx]][[3]] = sapply( 1:2970, function(i) {
17        D = predict_obs(Tid=tx, Sid=Sids[i], A=NULL, G=NULL, idx=i) 
18        14 * logistic(-D)
19    })
20}

With the data prepared, we can now directly compare the treatments. For illustrative purposes, let’s first ignore gender and see how it might lead us to suboptimal decisions.

 1T3T1_noG = treatment_contrast(3,1,3)
 2T2T1_noG = treatment_contrast(2,1,3)
 3
 4p1 = plot_treatment_contrast_noG(T3T1_noG, "Treatment 3 - Treatment 1")
 5p2 = plot_treatment_contrast_noG(T2T1_noG, "Treatment 2 - Treatment 1")
 6plot_grid(p1, p2,
 7          align = 'h',
 8          axis = "lb",
 9          ncol=1
10)

The figure above compares Treatment 3 to Treatment 2, with Treatment 1 acting as the common baseline. The horizontal axis shows the difference in the treatment outcomes between the target and baseline treatments. A negative difference indicates that the target treatment results in fewer expected days of heavy drinking and is hence better compared to the baseline.

Judging from the figure alone, one might conclude that Treatment 3 is better than Treatment 2 and supports prescribing Treatment 3 for all patients. However, since we simulated the data, we know this is not true because Treatment 3 and 2 work differentially for males and females. So let’s replot the figure but consider gender this time.

 1T3T1 = treatment_contrast2(3,1)
 2T2T1 = treatment_contrast2(2,1)
 3p3 = plot_treatment_contrast(T3T1, "Treatment 3 - Treatment 1")
 4p4 = plot_treatment_contrast(T2T1, "Treatment 2 - Treatment 1")
 5
 6plot_grid(p3, p4,
 7          align = 'h',
 8          axis = "lb",
 9          ncol=1
10)

Now, it becomes apparent how the treatments work differentially for males and females. Treatment 3 seemingly gives better results overall because it works particularly well for males. However, it is much less effective for females. Female patients will benefit much more if they receive Treatment 2.

Can We Do Better?

Previously, we picture a rather “static” mechanism that underlies the working of the treatments (see the DAG below). In particular, we assume that self-efficacy mediates between the treatment and the clinical outcomes. In other words, it is assumed that the treatments unidirectionally affect self-efficacy and self-efficacy unidirectionally affects the treatment outcomes. However, this assumption becomes a bit oversimplified once we think harder about how people’s lives unfold during treatments. The treatments may have raised participants’ self-efficacy in alcohol control, leading to improved clinical outcomes. Yet the improvement in clinical outcomes very likely also raises subsequent self-efficacy. So instead of a unidirectional influence from treatment to efficacy to clinical outcomes, a bi-directional interaction seems to be more plausible, with previous clinical outcomes and self-efficacy acting as feedback to later statuses.

Such a scenario might resemble something sketched in the DAG below. Here, each state of the efficacy ($E_t$) and the clinical outcomes ($D_t$) at time $t$ influences the states at $t+1$ directly. $E_0$ and $D_0$ represent the initial states before receiving treatments. These initial states may be influenced by several variables (e.g., age, gender, ethnicity, etc.), as illustrated by the greyed nodes on the left-most part of the graph. $T$ represents the treatment and is assumed to influence only $E$ here for the sake of simplicity.

In addition to the variables we have modeled in our unidirectional causal model, the DAG here introduces additional variables, $S_i$. These variables indicate extraneous influences on the clinical outcomes at time $t$. For instance, they might be the stress level of each participant at time $t$. Stressful life events could worsen clinical outcomes and often come and go without signs. When these random events occur and lead to decreases in participants’ clinical outcomes, the negative effect is likely to carry over to future states, such as lowering participants’ self-efficacy subsequently. The arrows $D_t \rightarrow E_{t+1}$ and $E_t \rightarrow D_{t+1}$ would allow one to model such effects.

The assumption of this dynamic feedback interaction leads to two questions. First, such a fine-grained interaction across time requires denser longitudinal data for the model to make sense. When the dataset is sparse in terms of time, the link between the current and the next state might be too weak as multiple events may have occurred in between and obscure their relationship. Thus, experience sampling methods or ecological momentary assessments, where data are collected daily, are preferable for this model. Second, the model is only conceptual at best so far. To prove that this dynamic feedback model can work, we would need to run simulations, construct the model, and test it, as always. This would probably take weeks to months for me. Once I’m done, I might document the process in another post someday.