组内相关系数(The Intraclass Correlation Coefficient)是测量实验可重复性的指标。例如某项实验有n个测量目标值,重复测量了k次。
举个例子有一个nk的矩阵。BMS表示测量目标(row)之间的组间均方差(mean square)。WMS表示测量目标的组内均方差。JMS表示重复测量(column)之间的组间均方差。EMS表示残差的均方差。
以下是常用的三种ICC的计算方式:
假设有一个64的矩阵,
根据公式计算可得,
因此,有:
下面贴出ICC的计算matlab脚本。
function [ICC] = ICC(X, model, rating)
% ---------------------------------------------------------------
% Intraclass Correlation Coefficient
%
% [ICC] = ICC(X, model, rating)
% 返回矩阵 X的ICC值
%
% The input X is an n x k matrix, where n is the number of tar-
% gets and k the number of judges. ICCs are divided into three
% models for two types of ratings.
%
% Model Types:
% '1' = random effects model => a different group of k judges,
% who have been selected from a larger population of jud-
% ges, rate each target
% '2' = mixed effects model => the same group of k judges, who
% have been selected from a larger population of judges,
% rate each target
% '3' = 固定效应模型 => the same group of k judges, who
% are the only judges of interest, rate each targets.
%
% Rating Types:
% '1' = individual rating
% 'k' = mean of several ratings
% ---------------------------------------------------------------
% Script:
[n,k] = size(X);
table = ANOVA_for_ICC(X, 'off');
BMS = table{6,4}; WMS = table{7,4};
JMS = table{3,4}; EMS = table{8,4};
if isequal(rating,'1')
if isequal(model,'1')
ICC = (BMS - WMS) / ...
(BMS + (k-1)*WMS);
elseif isequal(model,'2')
ICC = (BMS - EMS) / ...
(BMS + (k-1)*EMS + k*(JMS-EMS)/n);
elseif isequal(model,'3')
ICC = (BMS - EMS) / ...
(BMS + (k-1)*EMS);
else
error('Model type does not exist.');
end
elseif isequal(rating,'k')
if isequal(model,'1')
ICC = (BMS - WMS) / ...
BMS;
elseif isequal(model,'2')
ICC = (BMS - EMS) / ...
(BMS + (JMS-EMS)/n);
elseif isequal(model,'3')
ICC = (BMS - EMS) / ...
BMS;
else
error('Model type does not exist.');
end
else
error('Rating type does not exist.');
end
% ---------------------------------------------------------------
调用的函数为:
function [table] = ANOVA_for_ICC(X, displayopt)
% ---------------------------------------------------------------
% Analysis of Variance for Ratings (for Intraclass Correlations)
%
% [table] = ANOVA_for_ICC(X, displayopt)
% returns the ANOVA summary table as a cell array
%
% The input X is an n x k matrix, where n is the number of tar-
% gets and k is he number of judges. The option of displaying
% an ANOVA summary table is handled by displayopt, which can be
% set to 'on' (default) or 'off.'
% ---------------------------------------------------------------
% Script:
[m,n] = size(X);
total = m*n;
A = sum(sum(X.^2));
Bc = sum(sum(X,1).^2)/m;
Br = sum(sum(X,2).^2)/n;
D = sum(sum(X))^2/total;
SSbc = Bc-D;
SSwc = A-Bc;
SSbr = Br-D;
SSwr = A-Br;
SSe = A-Br-Bc+D;
DFbc = n-1;
DFwc = n*(m-1);
DFbr = m-1;
DFwr = m*(n-1);
DFe = DFbc*DFbr;
MSbc = SSbc/DFbc;
MSwc = SSwc/DFwc;
MSbr = SSbr/DFbr;
MSwr = SSwr/DFwr;
MSe = SSe/DFe;
Fbc = MSbc/MSe; Pbc = 1-fcdf(Fbc, DFbc, DFe);
Fbr = MSbr/MSe; Pbr = 1-fcdf(Fbr, DFbr, DFe);
table = { 'Source of Variance' 'SS' 'df' 'MS' 'F' 'Prob>F'
'Columns' [] [] [] [] []
' betweeen' SSbc DFbc MSbc Fbc Pbc
' within' SSwc DFwc MSwc [] []
'Rows' [] [] [] [] []
' between' SSbr DFbr MSbr Fbr Pbr
' within' SSwr DFwr MSwr [] []
'Residual' SSe DFe MSe [] [] };
if isequal(nargin, 1) || isequal(displayopt, 'on')
statdisptable(table, 'ANOVA for ICC', ...
'ANOVA Summary Table', '');
end
% ---------------------------------------------------------------
运行结果,如图所示:
FMRI工具包afni中有命令3dICC和3dlme以及3dlmer可计算ICC的值。