Paired difference test for nonhomogenous variances

Question

Concrete scenario: A random sample of adults in Idaho rate each of 'kombucha', 'apple juice' and 'grape juice' on a scale from 0 to 1. Does the population of Idaho assign significantly different mean ratings to these items?

Characteristics of the data (as I understand them):

1. The samples are paired, not independent (an individual's rating of apple juice and grape juice may be correlated)
2. The variance of the samples is not necessarily equal (kombucha may be more polarizing than apple juice)
3. The distributions aren't necessarily normal (kombucha ratings may have a bimodal distribution, for example)
4. The ratings are continuous

Tests that I've looked at:

• Friedman's test ~ Sounds like it would work, but the ranking step throws away information
• Welch's t-test ~ Seems to assume that the samples are independent
• Wilcoxon signed-rank test ~ Assumes that the sample variances are equal

TL;DR: What's a good paired difference test for a single sample where each entity produces two (or more) differently distributed continuous measures?

Answer 1

I think one important question is whether you want to consider subject as a random effect or a fixed effect.

You raise the possibility of fixed effects (treating subjects as blocks) by considering the Friedman test.

Otherwise - treating them as a random affects - it's essentially repeated measures. I'll focus on this approach for the moment.

The question then comes about how to treat the possible issues of bounded variables, heteroskedasticity and non-normality (in approximate order of importance).

One approach might be to consider generalized linear mixed models.

Usually with continuous proportions I'd lean toward thinking about beta distributions, which don't fit into the usual exponential family framework, but if we think for a moment about quasi-binomial models, we have something that will deal with bounded variables (i.e. we can at least model the mean reasonably) and we have a variance that's proportional to mean*(1-mean) (as is the case with the beta, which has variance = mean $ \times$ (1-mean) $\times$ $ \frac{1}{\alpha+\beta+1}$). This would however, restrict us to having a constant $\alpha+\beta$ for all observations (which doesn't strike me as a great issue for a first attempt at a model).

So that would probably be how I'd start thinking about this - approximating a beta-like model via a quasi-binomial GLMM, and we have a model which has the potential to deal with many of the issues and still leaving us comparing means.

If that wasn't adequate (e.g. remaining unmodelled heteroskedasticity), I'd probably start looking at making a Bayesian model via MCMC to deal with all the quirks in this.