I have three judges (A,B,C) that estimate a value for N objects. Each judge estimates the same object three times, so that my data frame looks like

 object    A1    A2    A3    B1    B2    B3    C1    C2    C3
      1 19.93 20.05 20.15 20.01 20.66 19.72 20.53 21.02 20.41
      2 19.78 19.16 19.47 19.90 20.50 19.60 20.41 20.76 20.42
      3 19.22 19.33 19.41 19.82 20.39 19.49 20.33 20.69 20.32
      4 19.21 19.67 19.41 19.76 20.34 19.42 20.30 20.67 20.29

Now I would like to test whether the variation beetween judges is greater than within each judge. The problem is, that each object is different, which means that the rows connot be considered measurements of the same random varaible, but the true value differs from object to object.

Does someone know of a simple test for this problem?

Here is what came to my mind:

  1. Compare the average spans (max - min) per object within each judge with the span between all judges. From a practical point of view, the absolute value of thsi difference is is more important than a significance test, but it would nevertheless be interesting to have a test for statistical significance.

  2. Do a paired t-test for all combinations (AB, AC, BC) with a Bonferroni-correction. This raises the question, how the three observations of each judge are taken into account. By adding differences of all combinations?

  3. Use ICC (intra-class correlation) for the average values of each judge. This approach, however, throws away two third of the data, because it reduces the three values per judge per object to a single one.

  • 1
    $\begingroup$ You could try a random effects only model, so a model with only random effects. $\endgroup$ – user2974951 Aug 29 at 10:00

Let's assume the variance of the measurements can be decomposed the following way:

$$ \sigma^{2}_{Y}=\sigma^{2}_{\mathrm{Judge}} + \sigma^{2}_{\mathrm{Object}} + \sigma^{2}_{\mathrm{Repeats:Judge}} +\sigma^{2}_{\mathrm{Residual}} $$

Where $\sigma^{2}_{\mathrm{Judge}}$ is the variance between judges, $\sigma^{2}_{\mathrm{Object}}$ is the variance between objects, $\sigma^{2}_{\mathrm{Repeats:Judge}}$ is the variance within judges and $\sigma^{2}_{\mathrm{Residual}}$ is the residual variance. These variance components can be estimated using a linear mixed model. But first, let's inspect the data:


We can see that there seems to be a lot of variance between the judges. The estimates of Judge C (blue boxes) are consistently higher than those of judge B (green) and judge A (red) who is lowest. The variance within each judge seems to be lower compared to the variance between judges. But there is heterogeneity there: Judge B's measurements (the green boxplots) seem more variable than those of judge A or C. Finally, the variance between objects also doesn't seem to be that big. Let's fit the model using your data and the R package lme4 now:


# The model with no intercept

mod <- lmer(y ~ -1 + (1|object) + (1|judge/rep), data = dat)


Random effects:
 Groups    Name        Variance  Std.Dev.
 rep:judge (Intercept)   0.08427  0.2903 
 object    (Intercept)   0.02960  0.1720 
 judge     (Intercept) 399.71254 19.9928 
 Residual                0.02768  0.1664 

# Confidence intervals for random effects based on profile likelihood

confint(mod, method = "profile")

             2.5 %     97.5 %
.sig01  0.16717461  0.5998100
.sig02  0.07714821  0.5310402
.sig03 10.62671323 59.1203551
.sigma  0.12845541  0.2273330

The model (fit by REML) results in the following estimates:

  • $\widehat{\sigma^{2}}_{\mathrm{Judge}} = 399.7125$
  • $\widehat{\sigma^{2}}_{\mathrm{Object}} = 0.0296$
  • $\widehat{\sigma^{2}}_{\mathrm{Repeats:Judge}} = 0.0843$
  • $\widehat{\sigma^{2}}_{\mathrm{Residual}} = 0.0277$

We can compare the variances (or standard deviations) by looking at their 95% confidence intervals (see output above). The 95% confidence interval for the standard deviation between judges is $(10.627; 59.120)$ whereas the corresponding 95% confidence interval for the standard deviation within judges is $(0.167; 0.600)$. These confidence intervals do not overlap from which we can conclude that they are not compatible. Hence, these data provide evidence that the within and between judge variances are different.

  • $\begingroup$ Thanks for the detailed explanation and R code. The data are actually approximately periodic in the objects, so a polynomial fit is not appropriate for the raw data. I can apply lmer, however, on the centered data by subtracting the mean over all judgements per object. Unfortunately this results in a "boundary (singular) fit" error. I will post this in another answer to this question and would greatly appreciate a comment form your side. $\endgroup$ – cdalitz Aug 30 at 8:22

Applying the suggestion of @COOLSerdash to the entrie data is problematic, because the values are not linear (nor polynomial), but approximately periodic: expert judgments per object I have tried to remove this effect by subtracting the mean estimate per object ID:

mean.y <- (x$A1 + x$A2 + x$A3 + x$B1 + x$B2 + x$B3 + x$C1 + x$C2 + x$C3) / 9.0
N <- nrow(x)
tmp <- rbind(data.frame(judge="A",
                        rep=c(rep(1,N), rep(2,N), rep(3,N)),
                        delta=c(x$A1-mean.y, x$A2-mean.y, x$A3-mean.y)),
                        rep=c(rep(1,N), rep(2,N), rep(3,N)),
                        delta=c(x$B1-mean.y, x$B2-mean.y, x$B3-mean.y)),
                        rep=c(rep(1,N), rep(2,N), rep(3,N)),
                        delta=c(x$C1-mean.y, x$C2-mean.y, x$C3-mean.y)))

enter image description here

Then the measurement could be approximately modeled as ($N$ is a normal random variable): $$\delta_{judge} = bias_{judge} + N(0,\sigma_{judge}) + N(0, \sigma_{total})$$

Following @COOLSerdash's appraach, this corresponds to the following mixed model in R's lmer syntax:

> mod <- lmer(delta ~ 0 + judge + (1|judge:rep), data = tmp)
> summary(mod)

Random effects:
 Groups    Name        Variance Std.Dev.
 judge:rep (Intercept) 0.09699  0.3114  
 Residual              1.04120  1.0204
Number of obs: 3825, groups:  judge:rep, 9

> confint(mod, method = "profile")

            2.5 %     97.5 %
.sig01  0.1652415 0.44148447
.sigma  0.9979202 1.04371959

As these confidence intervals do not overlap, I conclude that $\sigma_{judge}$ is significantly smaller than $\sigma_{total}$ and the variation within judges is thus smaller than between judges.

Hint: This answer was edited on Sep 26, because in my earlier post from Aug 30, I had used the wrong model specification in R. I could find the error due to the helpful comments of @COOLSerdash.

  • $\begingroup$ I don't understand what you mean by periodic? Could you elaborate, please? $\endgroup$ – COOLSerdash Aug 31 at 6:26
  • 1
    $\begingroup$ @COOLSerdash: I have edited my answer to include a scatterplot of the entire data. The object id actually is a time stamp (frame number in a video). $\endgroup$ – cdalitz Aug 31 at 19:07
  • $\begingroup$ Thanks for adding more information. I don't fully grasp at the moment why the trend is a problem when you are concerned about variances. But if you want to remove the trend, I'd probably try fitting a LOESS model or a spline model for each object and judge. The residuals from these models would be void of the trend. $\endgroup$ – COOLSerdash Sep 1 at 7:46
  • $\begingroup$ @COOLSerdash: As LOESS is a weighted regression over a kNN neighborhood, my method of simply subtracting the mean at the same time (object id) is formally identical to LOESS with a neighborhood of 𝑘=9. The resulting values then should only represent variance in judgements. My question is whether "$\verb#delta ~ -1 + (1|judge:rep)#$" is the correct formula for the statistical model $\delta = bias_{judge} + N(0,\sigma_{judge}) + N(0,\sigma_{residual})$? $\endgroup$ – cdalitz Sep 1 at 9:05
  • $\begingroup$ I don't think the model you fit corresponds to the formula. The formula includes a fixed effects for judge ($bias_{judge}$) and a random effect for judge. The model, however, fits no fixed effects and estimates two variances: the within-judge variance and the residual variance. $\endgroup$ – COOLSerdash Sep 1 at 10:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.