12 Measurement

12.1 Repeatability

Repeatability or test/retest reliabity is the agreement of measurements across technical replications.

12.1.1 Mean/difference - Bland/Altman plots

The Bland/Altman (Bland and Altman 1999) or Tukey mean/difference plot (Tukey et al. 1977) is a plot of agreement between two measured quantities. Here I use the mricloudpy package to read in data and convert it to a dataframe. Here, we have two measures of the same data. B/A plots typically add a 1.96 sd bar to detect outlying differences.

import statsmodels.api as sm
import numpy as np

import sys 
import os
sys.path.append("/home/bcaffo/sandboxes/MRICloudPy/mricloudpy/")

from mricloudpy.mricloudpy import Data
example = "/home/bcaffo/sandboxes/MRICloudPy/mricloudpy/sample_data/"
d = Data(example)

import_data: Data files found 
['/home/bcaffo/sandboxes/MRICloudPy/mricloudpy/sample_data/kirby127a_3_1_ax_283Labels_M2_corrected_stats.txt', '/home/bcaffo/sandboxes/MRICloudPy/mricloudpy/sample_data/kirby142a_3_1_ax_283Labels_M2_corrected_stats.txt', '/home/bcaffo/sandboxes/MRICloudPy/mricloudpy/sample_data/kirby239a_3_1_ax_283Labels_M2_corrected_stats.txt', '/home/bcaffo/sandboxes/MRICloudPy/mricloudpy/sample_data/kirby346a_3_1_ax_283Labels_M2_corrected_stats.txt', '/home/bcaffo/sandboxes/MRICloudPy/mricloudpy/sample_data/kirby422a_3_1_ax_283Labels_M2_corrected_stats.txt', '/home/bcaffo/sandboxes/MRICloudPy/mricloudpy/sample_data/kirby492a_3_1_ax_283Labels_M2_corrected_stats.txt', '/home/bcaffo/sandboxes/MRICloudPy/mricloudpy/sample_data/kirby501a_3_1_ax_283Labels_M2_corrected_stats.txt', '/home/bcaffo/sandboxes/MRICloudPy/mricloudpy/sample_data/kirby505a_3_1_ax_283Labels_M2_corrected_stats.txt', '/home/bcaffo/sandboxes/MRICloudPy/mricloudpy/sample_data/kirby656a_3_1_ax_283Labels_M2_corrected_stats.txt', '/home/bcaffo/sandboxes/MRICloudPy/mricloudpy/sample_data/kirby679a_3_1_ax_283Labels_M2_corrected_stats.txt', '/home/bcaffo/sandboxes/MRICloudPy/mricloudpy/sample_data/kirby742a_3_1_ax_283Labels_M2_corrected_stats.txt', '/home/bcaffo/sandboxes/MRICloudPy/mricloudpy/sample_data/kirby800a_3_1_ax_283Labels_M2_corrected_stats.txt', '/home/bcaffo/sandboxes/MRICloudPy/mricloudpy/sample_data/kirby814a_3_1_ax_283Labels_M2_corrected_stats.txt']
import_data: Importing...

visit1 = Data("mricloud/visit1/").df
visit2 = Data("mricloud/visit2/").df

import_data: Data files found 
['mricloud/visit1/142_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit1/239_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit1/346_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit1/422_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit1/492_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit1/501_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit1/505_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit1/656_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit1/679_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit1/742_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit1/800_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit1/814_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit1/815_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit1/849_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit1/906_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit1/913_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit1/916_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit1/934_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit1/959_kirby_3_1_ax_283Labels_M2_corrected_stats.txt']
import_data: Importing...

import_data: Data files found 
['mricloud/visit2/142_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit2/239_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit2/346_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit2/422_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit2/492_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit2/501_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit2/505_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit2/656_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit2/679_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit2/742_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit2/800_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit2/814_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit2/815_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit2/849_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit2/906_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit2/913_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit2/916_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit2/934_kirby_3_1_ax_283Labels_M2_corrected_stats.txt', 'mricloud/visit2/959_kirby_3_1_ax_283Labels_M2_corrected_stats.txt']
import_data: Importing...

import statsmodels.api as sm
import numpy as np
import matplotlib.pyplot as plt

csf1 = visit1[(visit1['Type'] == 1) & (visit1['Level'] == 1) & (visit1['Object'] == 'CSF')]['Volume']
csf2 = visit2[(visit2['Type'] == 1) & (visit2['Level'] == 1) & (visit2['Object'] == 'CSF')]['Volume']
plt.scatter(csf1, csf2)

<matplotlib.collections.PathCollection at 0x78d037e1bca0>

sm.graphics.mean_diff_plot(csf1, csf2)

Often the B/A plot is done on the log scale. This plots the log ratio on the y axis versus the geometric mean on the x axis.

sm.graphics.mean_diff_plot(np.log(csf1), np.log(csf2))

12.2 ICC

The intra-class correlation coefficient is a measure of agreement. It measures the ratio of the inter-subject variation to the total variation (intra and inter). I like to think of ICC as a random effect model. If $Y_{i j}$ is measurement $j$ on subject $i$ then consider the random effect model

$Y_{i j} = μ + U_{i} + ϵ_{i j} U_{i} \sim N (0, σ_{u}^{2}) ϵ_{i j} \sim N (0, σ^{2})$

Then, the ICC is defined as

$\frac{σ_{U}^{2}}{σ_{U}^{2} + σ^{2}} .$

Note, this ICC model applies even if there’s more than 2 measurements per subject.

Consider two subjects, however. There’s an easy moment estimator in that $Y_{i 2} - Y_{i 1} = ϵ_{i 2} - ϵ_{i 1}$ Thus, the variance of the differences is an estimator of $2^2. Similarly, $(Y_{i 1} + Y_{i 2}) / 2 = U_{i} + (ϵ_{i 2} + ϵ_{i 1}) / 2.$ Thus, the variance of the average is an estimate of $σ_{u}^{2} + σ^{2} / 4$ . Thus, we have two equations and two uknowns. This solution has the benefit that it doesn’t depend on the normality of the random effects. However, it can produce negative estimates. Another approach simply uses maximum likelihood.

import pingouin as pg
import pandas as pd
import statsmodels.formula.api as smf


csfdf1 = visit1[(visit1['Type'] == 1) & (visit1['Level'] == 1) & (visit1['Object'] == 'CSF')]
csfdf2 = visit2[(visit2['Type'] == 1) & (visit2['Level'] == 1) & (visit2['Object'] == 'CSF')]

csfdf = pd.concat( [csfdf1, csfdf2] )
csfdf['logvolume'] = np.log(csfdf['Volume'])

md = smf.mixedlm("logvolume ~ 1", csfdf, groups=csfdf["ID"]).fit()
md.summary()

sigmasq = md.scale