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.
