Mozgostroje

By now basically everyone (here, here, here, here and here, and there is likely even more out there) who writes a blog and knows how to do a statistical analysis has analysed data from a recent replication study and from the original study (data repository is here).

The study consists of two experiments. Let's focus on Experiment 1 here. The experiment tests a treatment and a control group. The performance is measured by six likert-scale items. The scale of each item has 9 levels. All responses are averaged together and we obtain a single composite score for each group. We are interested whether the treatment works, which would show up as a positive difference between the composite score of the treatment and the control group. Replication study did the same with more subjects.

Let's perform the original analysis to see why this dataset is so "popular".

Legend: O - original study, R - replication, C - control group, T - treatment group

In the original study the difference between the treatment and control is significantly greater than zero. In the replication, it is not. However the ratings in the replication study are overall higher. The author of the original study therefore raised a concern that no difference was obtained in replication because of the presence of ceiling effects.

How do we show that there are ceiling efects in the replication? The authors and bloggers presented various arguments that support some conclusion (mostly that there are no ceiling effects). Ultimately ceiling effects are a matter of degree and since no one knows how to quantify them the whole discussion of the validity of the conclusions from the replication study is heading into an inferential limbo.

My point here is that if the analysis computed the proper effect size - the causal effect size, we would avoid these kinds of arguments and discussions. Actually, this is fairly easy to see. If we interpret T-C estimate causally then it says that by giving people the treatment we increase their score by 0.83 points on average. So when their base score is 9.5 points we predict that the treatment will increase it to 10.33. This however can't happen because the scale goes only from 1 to 10. A quick fix is to say that the resulting score will be a minimum of the raw score and the upper ceiling value 10. This implies that a group with a base score 9.1 will end up with a score 9.93 which is again implausible. We need some function that maps from $[-\infty,\infty]$ (logit scale) into $[1,10]$ in a continuous fashion. The logistic function is the usual choice. (Probit function is also popular.) This is what it looks like:

We measure and report the effect size on the logit scale. To obtain precise prediction we can always transform into raw score. Then a claim that a treatment increases logit score, say, by 2 units does not produce any counter-intuitive results. This is illustrated above for two groups with base logit score of 1 and 7 respectively. The increase by 2 units on the logit scale leads to big jump in the raw score in the former case while in the latter case the raw score remains the same.

The author of the original study voiced concern that the ceiling effects may not be detectable in the composite score and we should focus on individual items. Let's look at the distribution of the scores of the individual items in the original study and in the replication.

We see that RC shows higher scores than OC on all items. We see also that the response distribution is multimodal (see the resume item). Furthermore the outcome variable is discrete. In contrast, the logistic distribution is unimodal and continuous.

As a consequence we need a better model. We will use ordered logit model which can be seen as an extension of the simple logit model. The probability of a response $k=1,\dots,K$ (in our case $K=10$) is $P(x=k) = \mathrm{logistic}(\beta-c_k)-\mathrm{logistic}(\beta-c_{k-1})$. $c_0=-\infty$ and $c_K=\infty$. The remaining coefficients $c_i$ are estimated from the data. $\beta$ is the item difficulty. $c_k$ can be interpreted as the response thresholds. For instance if $\beta$ equals $c_5$ then a response higher than 5 is given half of the times. The coefficients $\beta, c_i$ give only ordinal information and we may wish to fix one of them to a constant value (say $\beta=0$). If we have two groups we can fix the difficulty for the control group to zero and estimate the difficulty in the treatment group. This quantity then gives the causal effect size.

As with the logistic model above, the ordered logit accounts for ceiling effects in accord with our intuition. We already saw one aspect of how this works - big differences in logit score get squashed into small differences in the raw score. How does this work in the other way? At first blush small differences in the raw score at the ceiling should be transformed into big differences in the logit score. However, while the raw score differences quickly get small at the ceiling the estimation error does not decrease so quickly. As a consequence we get a poor signal-to-noise ratio at the ceiling and what we obtain is an uncertain logit-score estimate. Intuitively, this is what we would expect from a ceiling effect - in its presence it is difficult to say whether a reliable difference between two scores exists and what its magnitude is.

I hope you are appetized now, so let's see what the ordered logit model tells us.

To make it precise we will use the following model.

$y_{n,m}^g \sim \mathrm{orderedLogit}(\beta_m +\delta_g,c)$

$y_{n,m}^g$ is the response of subject $n$ to item $m$, i.e. a value in $[1,10]$. $g$ indexes the four groups which are of main interest to us. These are OC,OT,RC and RT. $\beta_m$ gives the difficulty for each of the six items. $c$ is the threshold vector. $\delta_g$ gives the ability of each group. We fix $\delta_{\mathrm{oc}}$ to zero.

The estimated values are displayed below. Reacall that the zero on the vertical axis shows the difficulty of the dog item and the OC group. We see that the trolley item is most difficult while the kitten item is easiest. The estimates of group ability show identical pattern as the raw score estimates. In particular, there is a control-treatment difference of 0.67 [0.21,1.13] in the original study while the estimate from the replication study -0.02 [-0.22,0.19] excludes the possibility of such a big difference in ability.

These results confirm the original analysis. However, in this case, there can be no discussion of ceiling effects. Unlike with the raw score difference, the ceiling effect does not bias the the value towards zero but rather increase the uncertainty. We can demonstrate this by comparing the raw score difference with the logit score difference in a simulation.

Assume that the true values of $c$ are the ones we just estimated above. Assume that the true effect of the treatment is equal to the one measured in the original study $\delta_T-\delta_C = \delta_{OT}-\delta_{OC}=0.67$. We systematically vary the ability of the control group C. For each ability value, we simulate data from the model for single item and 200 subjects (100 in each group). We then estimate the $\delta_T-\delta_C$. We can then ask how does the magnitude of the base ability ($\delta_C$) influence the ability of our method to detect the effect size. In particular we are interested in cases where $\delta_C$ is very high (or very low).

First, generate data.

Here we use a simplified model - one item, two groups.

Run estimation. I will spare you the output here.

The figure sums up this post. The estimates from the two methods (grey and magenta) share the y-axis, such that the vertical position of ground true treatment efect (red) and zero (black) is the same. The horizontal axis shows the ability of the control group ($\delta_C$). The ceiling effect is located on the right. Raw score difference (magenta) suffers from ceiling effect. To put the ceiling effects into context the blue and green circles show the neg-difficulty of each item in the original and replication study respectively. In this way we can say what kind of bias due to ceiling effect we would expect if we use items with the degree of difficulty observed in the original and replication study. The items in the original study (possibly with the exception of the kitten item) do not suffer from ceiling effects. However, several of the items in the replication study are located in the region with noticeable bias. At the same time it is clear that even if the replication study used the kitten item only, the ceiling effects would cut the effect size into half, but would not set it to zero. So ceiling effect yes, but not of a magnitude that could explain the zero raw score difference obtained in the replication.

However, this is not my main message. The figure shows that we should avoid the raw score difference in the first place. The magnitude of $T-C$ is biased towards zero. Even worse $T-C$ falsely claims high precision in the boundary regions - in the regions where the estimate is biased. In comparison, magnitude of $\delta_T-\delta_C$ is not biased by boundary effect. (The jumping around at the ends of the black line is due to discretization error introduced in the data generation proces.) Rather the boundary effect translates into worse precision of the estimate. Note also that the 95% interval of $\delta_T-\delta_C$ always includes the true effect size while $T-C$ incorrectly excludes it in the boundary regions.

What will you do next time when you analyse data where the values are constrained by boundaries? You can use some linear method such as the raw group difference or Anova and pray that your results are not be biased. This post showed you that you should fit a model from the logit family instead.