Where are "all" the variances in mixed? Manually calculating Intraclass Correlation Coefficients

[Stefan Kreisel]

SORRY - THIS IS WORK IN PROGRESS ...


Intraclass Correlation Coefficients (ICC) are hard to understand for us plain commoners, especially if the focus is not primarily on "classical" reliability,

\[ \text{ICC(3,1) by Shrout or ICC(C,1) by McCraw = } {\sigma^2_\text{Subject} \over \sigma^2_\text{Subject} + \sigma^2_\text{Error}} \]

but rather on the "expected trial-to-trial noise in the data" as in a test-re-test setting = seeing how much the scores agree with one another when repeating trials.
In this case it is suggested that one also includes systematic error in the denominator due to trials, so that

\[ \text{ICC(2,1) by Shrout or ICC(A,1) by McCraw = } {\sigma^2_\text{Subject} \over \sigma^2_\text{Subject} + \sigma^2_\text{Trials} + \sigma^2_\text{Error}} \]

Moreover, one quickly gets pointed towards a concept called "agreement" or "absolute reliability" (vs. "relative reliability"); in any case, more often than not, this concept is termed "Standard Error of Measurement" (SEM; not to be confused with the Standard Error of the Mean).

We can calculate the SEM by different methods, but most involve the ICC or components thereof - such as variance components.

1. \[ {\text{Pooled standard deviation} * \sqrt{\text{1-ICC }}} \]
2. \[ {\sqrt{\text{Error Mean Square} \text{ = MS}_\text{Error}}} \] or \[ {\sqrt{\sigma^2_\text{Error}}} \]

We can use the inbuilt Stata commands such as icc to get all the Intraclass Correlation Coefficients we want. If we want to

\[ \text{ICC(2,1)} = {\text{MS}_\text{Subject} - \text{MS}_\text{Error} \over \text{MS}_\text{Subject} + \text{(Number of trials - 1)} * \text{MS}_\text{Error}} \]

* Example generated by -dataex-. To install: ssc install dataex
clear
input int outcome byte(person_row _person_id trial_column)
166 1 1 1
168 2 2 1
160 3 3 1
150 4 4 1
147 5 5 1
146 6 6 1
156 7 7 1
155 8 8 1
160 1 1 2
172 2 2 2
142 3 3 2
159 4 4 2
135 5 5 2
143 6 6 2
147 7 7 2
168 8 8 2
end

[Erik Ruzek]

A couple of quick questions/thoughts.

1. Do you have more than two trials? Mixed models, estimated by maximum likelihood, work best when the sample size for a grouping variable is larger. They can be estimated w/ 5-10 groups, but the variance estimate will be less precise (huge confidence intervals). Also, you will want to use restricted maximum likelihood (use the mixed option reml). Bayesian estimation might work with 5 or fewer, but will be fairly dependent on a prior for the variance.
2. Check out this post on cross-validated. It addresses some of the issues you raise. It also suggests a potential solution if you have 5 or fewer trials - treat trial as a fixed factor. This has the effect of giving you Shrout's ICC(3,1).
3. Can you clarify your question?

[Stefan Kreisel]

Sorry - kids wanted dinner
I'll start from and repeat the beginning, makes more sense ...

Intraclass Correlation Coefficients (ICC) are hard to understand for us plain commoners, especially if the focus is not primarily on "classical" reliability,

\[ \text{Formula 1 - } \text{ICC(3,1) by Shrout or ICC(C,1) by McCraw = } {\sigma^2_\text{Subject} \over \sigma^2_\text{Subject} + \sigma^2_\text{Error}} \]

but rather on the "expected trial-to-trial noise in the data" as in a test-re-test setting = seeing how much the scores agree with one another when repeating trials.
In this case it is suggested (see this excellent paper) that one also includes systematic error in the denominator due to trials, so that

\[ \text{Formula 2 - } \text{ICC(2,1) by Shrout or ICC(A,1) by McCraw = } {\sigma^2_\text{Subject} \over \sigma^2_\text{Subject} + \sigma^2_\text{Trials} + \sigma^2_\text{Error}} \]

Moreover, one quickly gets pointed towards a concept called "agreement" or "absolute reliability" (vs. "relative reliability"); in any case, more often than not, this concept is termed "Standard Error of Measurement" (SEM; not to be confused with the Standard Error of the Mean).

We can calculate the SEM by different methods, but most involve the ICC or components thereof - such as variance components.

1. \[ {\text{Pooled standard deviation} * \sqrt{\text{1-ICC }}} \]
2. \[ {\sqrt{\text{Error Mean Square} \text{ = MS}_\text{Error}}} \] or \[ {\sqrt{\sigma^2_\text{Error}}} \]

We can use the inbuilt Stata commands such as icc to get all the Intraclass Correlation Coefficients we want, with confidence intervals and all.
We can even use anova or better yet wsanova (net install sg103.pkg) to churn out mean squares and plug those into formulas like below by hand - just to see how things "really" work; its nice for didactics as well.

\[ {\sigma^2_\text{Subject} \over \sigma^2_\text{Subject} + \sigma^2_\text{Error}} = {\text{MS}_\text{Subject} - \text{MS}_\text{Error} \over \text{MS}_\text{Subject} + \text{(Number of trials - 1)} * \text{MS}_\text{Error}} \]

But at times one might want to use mixed to do the same job. Why? Because more flexible and can handle missings.

A lot of text for following short query. Using the data below I run.

mixed outcome trial_column || person_row : , reml var

and get

[...]

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
person_row: Identity         |
                  var(_cons) |   66.42859    53.0146      13.90072     317.448
-----------------------------+------------------------------------------------
               var(Residual) |   57.10714   30.52504      20.03107    162.8084
------------------------------------------------------------------------------

var(_cons) is exactly σ²_Subject and var(Residual) is exactly σ²_Error; with that information, plug it into Formula 1 above and you get an ICC(3,1) aka ICC(C,1) of 0.53772775.
Of course you could have just pressed postestimation estat icc.

But what if you want to calculate the ICC(2,1) aka ICC(A,1) - if you are more interested in test-re-test "reliability" or the SEM (the 2nd way of calculating SEM provided above is for SEM "consistency"; SEM "agreement" asks for systematic error due to variability within trials and not the random part (the var(Residual)), which requires you to provide information on σ²_Trials (see here), as Formula 2 suggests.

σ²_Trials is easily retrievable using anova and some algebra. But is in in mixed?

But is systematic error due to variability within trials anywhere to be found in mixed output? Is it even possible given how linear mixed-effects models are set up? Does it require some fancy rearrangement of the mixed command?

I am dearly realizing that I never had formal stats training ...

Cheers - Stefan

* Example generated by -dataex-. To install: ssc install dataex
* Data from Table 1, Example data set Trial B - from excellent paper
clear
input int outcome byte(person_row _person_id trial_column)
166 1 1 1
168 2 2 1
160 3 3 1
150 4 4 1
147 5 5 1
146 6 6 1
156 7 7 1
155 8 8 1
160 1 1 2
172 2 2 2
142 3 3 2
159 4 4 2
135 5 5 2
143 6 6 2
147 7 7 2
168 8 8 2
end

[/QUOTE]

[Erik Ruzek]

But is systematic error due to variability within trials anywhere to be found in mixed output? Is it even possible given how linear mixed-effects models are set up? Does it require some fancy rearrangement of the mixed command?

The short answer is yes. It comes from the realization that in your design, person and occasions are cross-classified. Each person has exposure to the same set of trials. Likewise, any trial number is likely to appear in all (or most) of the person ids. How do you estimate such models in Stata?

mixed outcome || _all:R.trial_column || person_row: , reml var

This gives you separate, uncorrelated random variances for trial and person. You can then estimate ICC(3,1) from these variances. As I said in my previous post, in order to get a precise variance estimate for trial, you would want to have many trials (10-20 would be a good start, but 20+ would be best). Presently, with just two trials in your example data, the variance for trials is effectively 0. Note that estat icc does not work after fitting a cross-classified model, so you have to estimate ICC(3,1) by hand.

Yulia Marchenko (StataCorp) wrote an excellent primer on cross-classified modeling in Stata back in 2006. Since that was written, mixed has replaced xtmixed in Stata.

[Stefan Kreisel]

Thanks Erik Ruzek / after re-reading my post, I was under the impression that I got a bit carried away length-wise and nobody would answer :--)

Specifically (acknowledging your comment on the sparsity of my data, but its there for didactical purposes only): isn't

mixed outcome || _all:R.trial_column || person_row: , reml var

doing the same thing as

mixed outcome || person_row: , reml var

?

The Random-effects Parameters minus the fact that there's one for _all (which is, as you point out, quasi zero ...) are (and should be identical).

What I was trying to get at is: can ICC(2,1) be estimated with parameters produced by mixed?
That would necessitate that one can get at the σ²_Trial (or its Mean Square equivalent), as in:

anova outcome trial_column person_row , repeated(person_row)

                           Number of obs =      16     R-squared     =  0.7728
                           Root MSE      = 7.55693     Adj R-squared =  0.5132

                  Source |  Partial SS    df       MS           F     Prob > F
             ------------+----------------------------------------------------
                   Model |        1360     8         170       2.98     0.0842
                         |
             trial_col~n |       30.25     1       30.25       0.53     0.4904
              person_row |     1329.75     7  189.964286       3.33     0.0677
                         |
                Residual |      399.75     7  57.1071429   
             ------------+----------------------------------------------------
                   Total |     1759.75    15  117.316667

or equivalently, but formatted in a more accessible way:

wsanova outcome trial_column , id(person_row) between(_person_id)

That output does suggest that there is within-trial variability -> trial_column MS = 30.25

Thanks for the literature reference!

Cheers

[Erik Ruzek]

What I was trying to get at is: can ICC(2,1) be estimated with parameters produced by mixed?

Perhaps I'm not understanding the intent of your question, but in post #4, I showed that you could get ICC(2,1) from the following model, given that it reports a variance for trial_column:

mixed outcome || _all:R.trial || person_row: , reml  // no variance

However, it does so provided you have 1) non-trivial between-trial variance and 2) more than a handful of trials. At this point, we cannot conclude whether you have 1 because unfortunately you do not have 2, at least in the example data.

I acknowledge that the ANOVA results give you a partial sum of squares for trial_column but it is teeny tiny relative to the SS of the other components (and not significant if you put stock into such metrics). Effectively it is telling you the same thing as the mixed results - not enough information. Note that if you have unbalanced data (some people complete fewer trials than others), then the maximum likelihood (full ml or reml) estimators are more efficient and you should use mixed.

[Stefan Kreisel]

Ah, OK - then your

You can then estimate ICC(3,1) from these variances.

in post #4 would be a typo, I guess.

I can see where you are getting at and acknowledge that given

mixed outcome || _all:R.trial || person_row: , reml var

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
_all: Identity               |
             var(R.trial_~n) |   3.62e-15          .             .           .
-----------------------------+------------------------------------------------
person_row: Identity         |
                  var(_cons) |   68.10715   52.51827      15.02521      308.72
-----------------------------+------------------------------------------------
               var(Residual) |      53.75     26.875      20.17334    143.2119
------------------------------------------------------------------------------

the _all: = 3.62e-15 is really small!

I - nonetheless - would like to understand why doing it per-anova, one does get a non-close-to-zero value with which one can continue to calculate, and get somewhat non-equivalent ICC's for 3,1 vs 2,1.
As this is more of a didactical piece, bare with me

Given that (and using the anova results from post #5)

• σ²_{Subject aka person_row} = abs(MS_{person_row} - MS_Residual) / Number of trials
• σ²_Trial = abs(MS_{trial_column} - MS_Residual) / Number of subjects per trial
• σ²_Error = MS_Residual

we get

• σ²_{Subject aka person_row} = abs(189.96 - 57.11) / 2 = 66.43
• σ²_Trial = abs(30.25 - 57.11) / 8 = 3.36
• σ²_Error = 57.11

Then

• ICC(3,1) would be (see post #3, Formula 1) -> 66.43 / (66.43 + 57.11) = 0.5377
• ICC(2,1) would be (see post #3, Formula 2) -> 66.43 / (66.43 + 3.36 + 57.11) = 0.5234

Clearly that difference is negligible in real life, but ...

So what you are saying is that if a "real" dataset had more trials and more variability, then mixed would do an analogous job toanova / wsanova, etc. and provide variance components including σ²_Trial to calculate ICC / SEM by hand? Correct?

Regards

PS Real datasets, say checking some physical parameter over time to measure reliability / agreement over time, never have more than a couple of repeat trials ...

[Erik Ruzek]

Yes, I mistyped the ICC type in post #4.

So what you are saying is that if a "real" dataset had more trials and more variability, then mixed would do an analogous job to anova / wsanova, etc. and provide variance components including σ²_Trial to calculate ICC / SEM by hand? Correct?

Indeed.

PS Real datasets, say checking some physical parameter over time to measure reliability / agreement over time, never have more than a couple of repeat trials

And hence why in certain fields a ton of time and money is spent on developing and calibrating instruments to be highly accurate. I work in psychology and education where we measure fuzzy things like self-concept, motivation, and achievement. Accordingly, we have come to realize that to measure fuzzy things well, one needs a lot of observations of the phenomenon. This leads to multi-item survey instruments, achievement tests with lots of items, many repeated observations of individuals, etc. And on the statistical side, latent variables and measurement models.

[Stefan Kreisel]

Thanks again Erik Ruzek !

Just to validate for anyone else perhaps stumbling across this post. I've re-run the above using a "real world" dataset including 4 trials (in this case: yearly ratings of walking speed - baseline, 1, 2, 3 years thereafter; ie trial_column) in 225 subjects (ie person_row); there are quite a few missings per person, but also a substantial number of subjects with measurements on all 4 occasions. I can't share the data because it's not public domain, but: here are the results.

Part of the exercise of finding a way to get at variance components (via mixed) and then calculating ICC and SEM by hand is the fact that 1.) the Stata command icc doesn't work with missings 2.) using anova (BTW wsanova doesn't work with larger number of cases and demands a set matsize - which can't be done in newer versions of Stata) is that with unbalanced designs getting variances from mean squares is tricky.

mixed outcome trial_column || person_row : , reml var


Mixed-effects REML regression                   Number of obs      =       746
Group variable: person_row                      Number of groups   =       225

                                                Obs per group: min =         1
                                                               avg =       3.3
                                                               max =         4


                                                Wald chi2(1)       =      0.47
Log restricted-likelihood =    125.054          Prob > chi2        =    0.4947

------------------------------------------------------------------------------
     outcome |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
trial_column |   .0041471   .0060735     0.68   0.495    -.0077567    .0160509
       _cons |   1.368669   .0169906    80.55   0.000     1.335368     1.40197
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
person_row: Identity         |
                  var(_cons) |   .0099891   .0019467      .0068178    .0146356
-----------------------------+------------------------------------------------
               var(Residual) |   .0335563   .0020411      .0297851     .037805
------------------------------------------------------------------------------
LR test vs. linear regression: chibar2(01) =    46.59 Prob >= chibar2 = 0.0000

and

mixed outcome  || _all:R.trial_column || person_row : , var reml

Mixed-effects REML regression                   Number of obs      =       746

-----------------------------------------------------------
                |   No. of       Observations per Group
 Group Variable |   Groups    Minimum    Average    Maximum
----------------+------------------------------------------
           _all |        1        746      746.0        746
     person_row |      225          1        3.3          4
-----------------------------------------------------------

                                                Wald chi2(0)       =         .
Log restricted-likelihood =  131.84613          Prob > chi2        =         .

------------------------------------------------------------------------------
     outcome |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |    1.37975   .0162026    85.16   0.000     1.347994    1.411507
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
_all: Identity               |
             var(R.trial_~n) |   .0006772   .0007156      .0000854    .0053723
-----------------------------+------------------------------------------------
person_row: Identity         |
                  var(_cons) |   .0101357   .0019429      .0069613    .0147577
-----------------------------+------------------------------------------------
               var(Residual) |   .0329062   .0020064      .0291996    .0370834
------------------------------------------------------------------------------
LR test vs. linear regression:       chi2(2) =    52.31   Prob > chi2 = 0.0000

So for

• σ²_{Subject aka person_row} take person_row: Identity var(_cons) from the first mixed= 0.0099891
• σ²_Error take var(Residual) from the first mixed= 0.0335563
• σ²_Trial take _all: Identity var(R.trial_~n) from the second mixed= 0.0006772

the ICCs are

• ICC(3,1) = 0.0099891 / (0.0099891 + 0.0335563) = 0.229
• ICC(2,1) = 0.0099891 / (0.0099891 + 0.0006772 + 0.0335563) = 0.226

So just out of curiosity, I tried to replicate those values with the R package library(psych) and its command ICC. The advantage of going down that road - albeit a lot less fun than trying to figure out how things work in Stata - is: R's ICC 1.) handles missings and 2.) spits out variance components if one knows how. So.

> givemesomeresults <- ICC(x)


                         type      ICC
Single_raters_absolute   ICC(1,1)  0.2288503
Single_random_raters     ICC(2,1)  0.2318366
Single_fixed_raters      ICC(3,1)  0.2354842
Average_raters_absolute  ICC(1,k)  0.5427653
Average_random_raters    ICC(2,k)  0.5469425
Average_fixed_raters     ICC(3,k)  0.5519854

Now

• post-mixed manual ICC(3,1) ≠ R's ICC (3,1) as in 0.229 ≠ 0.2354842, and
• post-mixed manual ICC(2,1) ≠ R's ICC (2,1) as in 0.226 ≠ 0.2318366.

but close enough!

And whats cool is:

> givemesomeresults$lme

             variance  Percent
ID       0.0101356931  0.23183659
Items    0.0006772061  0.01548993
Residual 0.0329062276  0.75267349
Total    0.0437191268  1.00000000

• σ²_{Subject aka person_row} from mixed= 0.0099891; from R = ID = 0.0101356931
• σ²_Error from mixed= 0.0335563; from R = Residual = 0.0329062276
• σ²_Trial from mixed=0.0006772; from R = Items = 0.0006772061

So again not equal, but close enough for me to be using Stata mixed for ICC and SEM from now on (and understanding what's happening, even more important ...).

But - HaHa! If you just run

mixed outcome  || _all:R.trial_column || person_row : , var reml

you'll see: identical output from Stata and R! Just realized that - take it for what it's worth.

Best

[Erik Ruzek]

• post-mixed manual ICC(2,1) ≠ R's ICC (2,1) as in 0.226 ≠ 0.2318366.

I think you are making an incorrect assumption about which variances you use in these calculations. When you run the model,

mixed outcome || _all:R.trial_column || person_row : , var reml

you pull all the variances for the ICC(2,1) calculation from the output of that model. Hence,

di .0101357/(.0006772+.0101357+.0329062)
.23183689

is equivalent up to six decimal places as what is reported by R.

Pulling them from the prior model doesn't make sense because it hasn't accounted for the total variance in the same way. They are different models for the total variation in outcome.

[Stefan Kreisel]

Yep. That was my surprise when

mixed outcome || _all:R.trial_column || person_row : , var reml

put it all together. Mystery solved.

[Stefan Kreisel]

Perhaps one follow-up.

The formulas for ICC derived from the variances as presented in post 2 for ICC(3,1) and ICC(2,1) are the usual ones you find in textbooks.
As I inadvertently did in post #9 was to take σ²_{Subject aka person_row} and σ²_Error from mixed outcome trial_column || person_row : , reml var and mix it with σ²_Trial from mixed outcome || _all:R.trial_column || person_row : , var reml. As Erik Ruzek pointed out: that's a no-go: use the variances from one formula only and don't mix.
Now I can calculate ICC(3,1) (but not ICC(2,1) - because I need a parameter for σ²_Trial) from both mixed outcome trial_column || person_row : , reml var using Formula 1 form post #2 and mixed outcome || _all:R.trial_column || person_row : , var reml using Formula 2 form post #2 - because both provide information for σ²_{Subject aka person_row} and σ²_Error.
But: as I demonstrated by error in post #9 - that gives (albeit only slightly) differing results for the ICC(3,1).

Intuitively I would use mixed outcome trial_column || person_row : , reml var (and not mixed outcome || _all:R.trial_column || person_row : , var reml) if I wanted to calculate ICC(3,1), because all I need is σ²_{Subject aka person_row} and σ²_Error and the parameters here are "complete" for that task.

If there really is only one correct road to Rome, it would also imply that in the textbook formulas for ICC(3,1) and ICC(2,1) - σ²_{Subject aka person_row} and σ²_Error actually are not identical, but slightly different as in

• ICC(3,1) = σ²_{Subject aka person_row}#1 / ( σ²_{Subject aka person_row}#1 + σ²_Error#1)
• ICC(2,1) = σ²_{Subject aka person_row}#2 / ( σ²_{Subject aka person_row}#2 + σ²_Trial + σ²_Error#2)

depending on which estimation they stem from.

[daniel klein]

For what it's worth, kappaetc (SSC) can estimate ICCs, too. Using the example data in #1

* Example generated by -dataex-. To install: ssc install dataex
* Data from Table 1, Example data set Trial B - from excellent paper
clear
input int outcome byte(person_row _person_id trial_column)
166 1 1 1
168 2 2 1
160 3 3 1
150 4 4 1
147 5 5 1
146 6 6 1
156 7 7 1
155 8 8 1
160 1 1 2
172 2 2 2
142 3 3 2
159 4 4 2
135 5 5 2
143 6 6 2
147 7 7 2
168 8 8 2
end

// set-up
keep outcome person_row trial_column
reshape wide outcome , i(person_row) j(trial_column)

kappaetc outcome1 outcome2 , icc(3) var
kappaetc outcome1 outcome2 , icc(2) var

display r(sigma2_s) / ( r(sigma2_s) + r(sigma2_r) + r(sigma2_e) )

yields

. kappaetc outcome1 outcome2 , icc(3) var

Interrater reliability                           Number of subjects =       8
Two-way mixed-effects model                     Ratings per subject =       2
------------------------------------------------------------------------------
               |   Coef.     F     df1     df2      P>F   [95% Conf. Interval]
---------------+--------------------------------------------------------------
      ICC(3,1) |  0.5377   3.33     7.00    7.00   0.068    0.0000     0.8865
---------------+--------------------------------------------------------------
      sigma2_s | 66.4286
      sigma2_e | 57.1071
------------------------------------------------------------------------------
Confidence interval is clipped at the lower limit.

. kappaetc outcome1 outcome2 , icc(2) var

Interrater reliability                           Number of subjects =       8
Two-way random-effects model                    Ratings per subject =       2
------------------------------------------------------------------------------
               |   Coef.     F     df1     df2      P>F   [95% Conf. Interval]
---------------+--------------------------------------------------------------
      ICC(2,1) |  0.5377   3.33     7.00    7.00   0.068    0.0000     0.8912
---------------+--------------------------------------------------------------
      sigma2_s | 66.4286
      sigma2_r |  0.0000 (replaced)
      sigma2_e | 57.1071
------------------------------------------------------------------------------
Confidence interval is clipped at the lower limit.

. 
. display r(sigma2_s) / ( r(sigma2_s) + r(sigma2_r) + r(sigma2_e) )
.55274889

Methods and formulas are given in

Gwet, K. L. (2021). Handbook of Inter-Rater Reliability. Volume II: Analyses of Quantitative Ratings. (5th ed.) Gaithersburg, MD: Advanced Analytics, LLC.
Gwet, K. L. (2014). Handbook of Inter-Rater Reliability. (4th ed.) Gaithersburg, MD: Advanced Analytics, LLC.

[Stefan Kreisel]

Wonderful! Thanks and regards. Stefan