Statistical methods for analyzing results of preference study I've chosen to compare the results of three different machine-learning models outputting the mouth animation of a talking 3D model for my thesis using a user preference study.
Users will see two videos at a time, each from one of the models, and indicate which they think look the most natural, or if they are neutral. The videos are synchronized. For each user, one comparison will be made for each sentence (out of a set of sentences) and for each pair of models to be compared.
Now, I could arrange the data as simply A vs. B, B vs. C, and A vs. C. Then, let's say A was preferred a majority of the time over B (a majority out of a total of sentences*subjects comparisons between them). What test should I use to see if this is statistically significant? Could I even just compare them all at the same time by treating the "wins" as "points" to get three sums or means representing "scores" and then use a repeated measures ANOVA followed by e.g. Tukey HSD?
A lot of methods (Tukey's HSD, Student's t-test...) have the assumption that the score in A needs to be independent to the score in B, but since this is a preference study, they aren't really independent, right? E.g. if there's a total of 20 un-neutral comparisons, if A wins 15, then B wins 5.
I would greatly appreciate any insight!
 A: Right, I'd treat this as a multilevel logistic regression. Formally:
$y_{ijk} = \beta_0 + \beta_1V_{ijk} + \alpha_i^p + \alpha_i^{p*v}V_{ijk} + \alpha_j^s + \alpha_j^{s*v}V_{ijk} + \alpha_{ij}^{p*s*v}V_{ijk} + \epsilon_{ijk}$
where you're modeling the response (y) to the $i^{th}$ sentence from the $j^{th}$ subject to the $k^{th}$ ml model.
$\beta_0$ is the overall intercept (i.e. mean response)
$\beta_1$ reflect the change expected overall moving from one model to another (this depends on how you encode the models - here I'm assuming dummy coding)
$V_{ijk}$ is the ml model $k$ for subject $j$, sentence $i$
The $\alpha_i^P$'s and $\alpha_j^s$'s represent "random intercepts", or how much each individual's ($p$) or sentence's ($i$) mean differs from the overall mean.
The $\alpha_i^{p*v}$ $\alpha_j^{s*v}$ represent "random slopes", or how much each individual or sentence shows a different pattern of preferences among the ML models.
Finally, the $\alpha_{ij}^{p*s*v}$ represents an interaction between the possible subject-sentence pairings.
To estimate this kind of model in R, you would use the following code:
library(lme4)
m <- glmer(y ∼ v + (v|subject) + (v|sentence) + (v|subject:sentence), 
           data=df, family='binomial')

Here's a some other good references on these models:
https://www.annualreviews.org/doi/abs/10.1146/annurev-psych-122414-033702
http://jakewestfall.org/publications/JWK.pdf
https://www.sciencedirect.com/science/article/pii/S0749596X08000193
Edit
Based on your comments below, and points I initially overlooked, I think the model should be reparameterized slightly. Basically, I don't think you have the data to estimate the final $\alpha_{ij}^{p*s*v}$ from above. Additionally, I overlooked the point that there's no clearly defined 'baseline' in your case, since the levels of your $V_{ijk}$ should be m1_vs_m2, m1_vs_m3, and m2_vs_m3, and it isn't clear how to interpret an estimated difference between m1_vs_m2 and m2_vs_m3.
Relatedly, I was a little bit vague about the meaning of $y_{ijk}$. This should be coded as 0/1, where 1 means the right-hand model was selected. For completeness, this means that for the category m1_vs_m2, a 0 would indicate that m1 was selected, and a 1 would indicated m2 was selected; for the category m1_vs_m3, a 0 would indicate that m1 was selected, and a 1 would indicated m3 was selected; for the category m2_vs_m3, a 0 would indicate that m2 was selected, and a 1 would indicated m3 was selected.
Your data would be set up to look something like this:
  y           v subject sentence
1 1       m1vm2       1        8
2 1       m1vm2       1        1
3 1       m1vm3       1        4
4 0       m2vm3       1        6
5 1       m1vm2       1        1
6 0       m1vm2       1       12

The model using lme4 syntax would look like this:
library(lme4)
m <- glmer(y ∼ v + (v|subject) + (v|sentence), 
           data=df, family='binomial')

The (v|subject) and (v|sentence) terms are the parts of the model that take care of your concerns of non-independence. You're partitioning the variance into variance attributable to characteristic differences between subjects and sentences. Basically, you're allowing the intercept and the effect of v to vary by each of those factors. 
For interpretation, you're going to want to examine the intercept of this model and the coefficients for the levels of v. Assuming m1_vs_m2 is your baseline, the intercept would reflect the log of the odds that m2 is preferred. In a null-hypothesis testing framework, if the intercept is statistically significant, then that corresponds to a preference for m2 over m1. The coefficients for v will represent deflections from this intercept for your other two comparison categories.
I'm not sure what software you were planning to use to analyze your data, but that's a bit beyond the scope of the original question. If you think this is a reasonable approach, I'd then post on stackoverflow to figure out how to program this in your language of choice.
