Measuring intra-observer correlation We're far from what you'd call expert statisticians, so we'd like to ask some help in an ongoing medical research.
We have 60 subjects. A researcher took at least 3 times the measurements of 8 of their values of interest, giving more than 60×8×3=1440 data points.
As a first step we'd like to determine the intra-observer correlation (Edit: mistake, this should have been variability) of these measurements, but we hardly even know how to get started.
Asking Google and Wikipedia gave us a few options:
Pearson: this expects pairs of data points, but we have 3 or more measurements for each subject and each attribute. We could arrange them in pairs but we'd rather find a method that requires less manual work
Intra-Class Correlation: The concept of a "group" is not so evident here. How would we group the measurements in this case?
If the above isn't clear, we'd gladly elaborate more.
Thank you all in advance!
 A: To simplify things, we focus on one of the eight measurements and call it $Y$.
How would we measure intra-observer variability if we had only three values of a single object? Simple answer: Look at the variance (or simpler the standard deviation) of the three values. 
Now, how should we proceed with 60 such objects? A straightforward estimate of the intra observer variability is obtained by averaging all 60 variances obtained as described above. (Again, its square root, the "average" standard deviation" is easier to interpret.)
Let's illustrate this in R using three fake objects as toy example:
Y <- c(1, 1.1, 0.9, 1.5, 1.4, 1.7, 0.5, 0.7, 0.7)
Object <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
sqrt(mean(tapply(Y, Object, var))) # Gives 0.1247

Is this method too simple? Here is an alternative based on a random effects model. Such models are often used to find intraclass correlation coefficients (ICC), so they must be of any value somehow! We model average $Y$ by a global mean $\mu$ plus a random intercept $a_i$ adjusting the global mean to the mean measurement of the $i$th object. Formally, measurement $j$ of object $i$ is decomposed as
$$
   Y_{ij} = \mu + a_i + \varepsilon_{ij}.
$$
The error term $\varepsilon_{ij}$ is assumed to have constant variance $\sigma^2$ across objects and also across measurements by object. We may consider $\sigma$ (or its square) as the intra-observer variability. An estimate $\hat\sigma$ of $\sigma$ is found in R:
library(lme4)
lmer(Y~1 + (1|Object)) 

The output contains the residual standard deviation $\hat \sigma = 0.1247$. Thanks to the balanced design (always three measurements per object), this is exactly the same value as obtained by the very simple method explained above. Important: Computing ICC involves looking at standard deviations of random effects. We don't need this here. We just need to check the standard deviation of the residuals in such model.
