# How to describe a split-plot ordinal logistic regression

I am having difficulties reporting a split-plot ordinal logistic regression. I get the concept of ordinal logistic regression involving a categorical representation of a latent continuous variable, but am having difficulty with the notion of the link function and whether the categorical variable is expressed as a probability.

### Description of the Data

I recently ran an experiment on placebo caffeine withdrawal reduction. We found that when people drinking decaf believed they were receiving caffeine they would report their caffeine withdrawal symptoms reducing more than when they they (correctly) believed they were receiving decaf, and both groups reported their symptoms reducing more than a group given tap water. Part of this experiment involved asking all participants at the start of the study how much they would expect water, decaf, and caffeine to reduce their caffeine withdrawal symptoms on a scale from 0-10.

Although people were not assigned to their experimental groups at the time I wanted to test whether there were any differences in average 0-10 ratings of the withdrawal reduction potential of the three beverages in the people who would be allocated to the three experimental groups.

### The Model

To do this I ran a hierarchical Bayesian split-plot analysis, with one three-level between-subject factor based on the experimental group people were allocated to (Told Caffeine vs Told Decaf vs Told Water), one three-level within-subjects factor beverage type (coffee vs decaf vs water) which represented the type of beverage whose withdrawal-reduction potential people were rating.

At first I ran this as a standard Gaussian analysis, with an identity inverse-link function, however when I went to generate posterior predictive plots I found that the model would make predictions below and above the actual range of the data. On the advice of some of the clever folks working on the Stan project I switched to an ordinal logistic regression, and, although it required learning new techniques (see here) it worked well and generated sensible estimates as well as posterior predictions constrained to the 0-10 scale.

My problem is that I am not really familiar with the language of ordinal logistic regression, so actually reporting the analysis I undertook, in a way that is comprehensive and correct (like the answer to one of my previous questions here) has been very challenging. I was hoping to get some help.

### The Model Description

I started by saying that participants were part of a hypothetical population of 24-hour abstinent moderate to heavy coffee drinkers characterised by two things,

1. their group membership (Told Caffeine vs Told Decaf vs Told Water)
2. the expectancies they hold about the withdrawal reduction potential of coffee, decaf and water.

I then say that we make the assumption that these withdrawal reduction expectancies can be expressed on a continuous 0-10 scale which we will model as a latent variable. The more the hypothetical person expects their beverage to reduce their caffeine withdrawal the higher their score on the latent variable.

If we designate P to be an unordered vector identifying each individual in the population, G to be an unordered vector of the group each participant belongs to (Told Caffeine vs Told Decaf vs Told Water), B to be and unordered vector designating the beverage type (coffee vs decaf vs water) about which the each individual P holds the expectancies of caffeine withdrawal, then scores on the latent variable can be modeled as an additive combination of the linear predictors plus error, expressed via the equation

$$\begin{equation*} z_{GxBxP|G[i,j,k|i]} = \beta_0 + \beta_{G[i]} + \beta_{B[j]} + \beta_{GxB[i,j]} + \beta_{P|G[k|i]} + \varepsilon_k \end{equation*}$$

where $$z_{GxBxP|G[i,j,k|i]}$$ is the amount that participant k in group i expects beverage j to reduce their caffeine withdrawal on the latent 0-10 scale; $$\beta_0$$ is the mean amount all participants across all groups expect all beverages to reduce their caffeine withdrawal; $$\beta_{G[i]}$$ is the average deflection from the grand mean for the experimental group i that participant k belongs to; $$\beta_{B[j]}$$ is the deflection from the grand mean for the beverage j that participant k is rating; $$\beta_{GxB[i,j]}$$ is the average deflection for participants in group i rating beverage j; $$\beta_{P|G[k│i]}$$ is the unique deflection for participant k in group i; and $$\varepsilon_k$$ is the amount that participant k’s score deviates from the score predicted by the linear equation.

So as to allow our statistical model of the latent variable $$z_{GxBxP|G[i,j,k|i]}$$ to generate predictions on the 0-10 scale of the data we will use an ordered logistic regression model. In this model each person in the hypothetical population’s latent expectancies of caffeine withdrawal reduction $$z_{GxBxP|G[i,j,k|i]}$$ were split into 101 discrete ordered categories, each category spanning an interval of 0.1 on the 0-10 scale of the latent variable. This model of the latent scale can be represented as

$$\begin{equation*} y_{GxBxP|G[i,j,k|i]} = \left\{ \begin{array}{l} 1\ \text{if} \ z_{GxBxP|G[i,j,k|i]} \in (0.0, 0.1)\\ 2\ \text{if} \ z_{GxBxP|G[i,j,k|i]} \in (0.1, 0.2)\\ 3\ \text{if} \ z_{GxBxP|G[i,j,k|i]} \in (0.2, 0.3)\\ \vdots\\ 99\ \text{if} \ z_{GxBxP|G[i,j,k|i]} \in (9.8, 9.9)\\ 100\ \text{if} \ z_{GxBxP|G[i,j,k|i]} \in (9.9, 10.0)\\ \end{array} \right. \end{equation*}$$

where the probability that the latent variable $$z_{GxBxS|G[i,j,k|i]}$$ falls within the boundaries of either of the 101 categories can be expressed as

$$\begin{equation*} y_{GxBxP|G[i,j,k|i]} = \text{logistic}(z_{GxBxS|G[i,j,k|i]}, \sigma^2) \end{equation*}$$

where $$\sigma^2$$ is the variance of the errors $$\varepsilon_k$$.

Now I know that I am doing several things wrong in this last part (and probably elsewhere) but I don't know where.

For example:

1. is it true that $$\text{logistic}(z_{GxBxS|G[i,j,k|i]}, \sigma^2)$$ is a probability? And, if so, shouldn't the expression be wrapped in a $$Pr(...=...)$$ type statement

$$\begin{equation*} \text{Pr}(y_{GxBxP|G[i,j,k|i]} = \text{category}\ x) = \text{logistic}(z_{GxBxS|G[i,j,k|i]}, \sigma^2) \end{equation*}$$

or something like that?

1. And is there a better way to express the ordinal logistic function around $$y_{GxBxP|G[i,j,k|i]}$$? I have seen it written

$$\begin{equation*} y_{GxBxP|G[i,j,k|i]} = \text{logit}^{-1}(z_{GxBxS|G[i,j,k|i]}, \sigma^2) \end{equation*}$$

but in this sort of formulation the $$Pr(...=...)$$ statement is missing.

1. I also know that there needs to be some discussion of the area between cutpoints under the probability density function, but I might be getting confused with ordinal probit regression now.

As you can probably tell I am getting tied up in knots with all this stuff. I would like to try to explain the model correctly and fully, using the correct mathematical notation.