B Developing a Model
Our final modeling strategy is an evolution from other attempts. The development proceeded through multiple iterations described in chapter 4, but doesn’t tell the full story. We learn more from others when they share what didn’t work along with the final path that did work. There is knowledge to be gained in failed experiments, because then there is one more way to not do something, just like a failing outcome reduces the variance of the Beta distribution.
In the first attempt at modeling, we used a classical GLM to get a baseline understanding of the data, but the fact that some estimates for certain subjects failed due to complete separation reinforced the our adoption of non-classical techniques. Our first Bayesian model was derived from Lee and Wagenmakers (2014) which used nested loops to iterate over subjects and SOA values. The data were required to be stored in a complicated way that made it difficult to comprehend and extend.
We moved on to using arm::bayesglm to remove convergence issues, but we were met with other limitations such as linear parameterization and lack of hierarchical modeling. The book Statistical Rethinking (McElreath 2020) offers a great first introduction to Bayesian multilevel modeling. McElreath’s rethinking package accompanies the book, and offers a compact yet expressive syntax for models that get translated into a Stan model. A model with age group and block can be written using rethinking::ulam as
rethinking::ulam(alist(
k ~ binomial_logit(n, p),
p = exp(b + bG[G] + bT[trt]) * (x - (a + aG[G] + aT[trt])),
a ~ normal(0, 0.06),
aG[G] ~ normal(0, sd_aG),
aT[trt] ~ normal(0, sd_aT),
b ~ normal(3, 1),
bG[G] ~ normal(0, sd_bG),
bT[trt] ~ normal(0, sd_bT),
c(sd_aG, sd_aT, sd_bG, sd_bT) ~ half_cauchy(0, 5)
), data = df, chains = 4, cores = 4, log_lik = TRUE)While learning about multilevel models, we tried writing a package that generates a Stan program based on R formula syntax. At the time the concepts of no-pooling, complete pooling, and partial pooling were vaguely understood, and the package was plagued by the same lack of flexibility that rstanarm and brms have. Then it was discovered that brms and rstanarm already did what we were trying to do, but programming experience was invaluable.
We also tried using lme4, rstanarm, and brms, and learned more about the concepts of fixed and random effects. We noticed that parameterization can have a significant affect on the efficiency of a model and the inferential power of the estimated parameters. When fitting a classical model, there is little difference in estimating a + bx vs. d(x - c) since the latter can just be expanded as -cd + dx which is essentially the same as the first parameterization, but there is a practical difference in the interpretation of the parameters. The second parameterization implies that there is a dependence among the parameters that can be factored out. In the context of psychometric functions, there is a stronger connection between PSS and c and the JND and d. This parameterization made it easier to specify priors and also increased the model efficiency. Of the modeling tools mentioned, only rethinking and Stan allow for arbitrary parameterization.
We finally arrived at a model that worked well, but learned that using a binary indicator variable for the treatment comes with the assumption of higher uncertainty for one of the conditions. The linear model that we arrived at is displayed in equation (B.1).
\[\begin{equation} \theta = \exp(\beta + \beta_G +(\beta_T + \beta_{TG})\times trt) \left[x - (\alpha + \alpha_G + (\alpha_T + \alpha_{TG})\times trt)\right] \tag{B.1} \end{equation}\]
Using an indicator variable in this fashion also introduced an interaction effect into the model that we almost did not account for after switching to using a factor variable. Interaction effects between factors is handled by creating a new factor that is essentially the cross-product of other factor variables. E.g. for factor variables \(x\) and \(y\)
\[ x = \begin{bmatrix} a \\ b \\ c \end{bmatrix}, y = \begin{bmatrix} i \\ j \end{bmatrix}\Longrightarrow x\times y = \begin{bmatrix} ai & aj \\ bi & bj \\ ci & cj \end{bmatrix} \]
The final round of reparameterization came in the form of adopting non-centered parameterization for more efficient models. To us, \(Z \sim N(0, 1^2);\quad X = 3 + 2Z\) is the same as \(X \sim N(3, 2^2)\), but to a computer the process of sampling from \(X\) can be more difficult than sampling from \(Z\) (discussed in chapter 2).