Hello,
I ran a linear mixed model with repeated measures using PROC MIXED. I used the RCORR statement to obtain the intraclass correlation coefficient (ICC) for the non-independent variable 'score' across readers. Is there a function within PROC MIXED that outputs the 95% CI and p-value for the ICC?
PROC MIXED DATA=final COVTEST;
CLASS subject_id reader (REF='1');
MODEL score = reader / S CL;
REPEATED / SUBJECT=subject_id TYPE=cs R RCORR;
RUN;
Thank you for your help!
OK. The attachment is the original paper you mentioned .
The following marco is I wrote before. Hope could give you a hand.
%macro intracc(data=,target=,rater=, depvar= ,df=,fas_pps=,label=);
proc glm data=&data outstat=_stats_ noprint ;
class &target &rater;
model &depvar = &target &rater ;
run;
proc sort data=_stats_;
by _name_ _SOURCE_;
run;
data want;
retain msw msb wms ems edf bms bdf jms jdf k;
set _stats_;
by _name_;
if upcase(_type_)='SS1' then delete;
if upcase(_source_)='ERROR' then do;
ems=ss/df;
edf=df;
end;
if upcase(_source_)="%upcase(&target)" then do;
bms=ss/df;
msb=bms;
bdf=df;
end;
if upcase(_source_)="%upcase(&rater)" then do;
jms=ss/df;
jdf=df;
k=df+1;
end;
if last._name_ then do;
msw=((ems*edf)+(jms*jdf))/(edf+jdf);
wms=msw;
n=bdf+1;
sffixed=(bms-ems)/(bms+((k-1)*ems)); /* ICC(3,1);*/
/*95% Confidence Interval for sfsingle - ICC(3,1)
Suppose that each of N subjects (targets) yields K observations.*/
F0=bms/ems;
V1=(n-1)*(k-1);
V2=n-1;
FL=F0/quantile('F',.975,V2,V1);
FU=F0*quantile('F',.975,V1,V2);
sffixed_lower=(FL-1)/(FL+k-1); /* ICC(3,1) lower CI*/
sffixed_upper=(FU-1)/(FU+k-1); /* ICC(3,1) upper CI*/
output;
end;
label sffixed="Shrout-Fleiss reliability: fixed set";
run;
data report;
set want;
length fas_pps label proj value $ 80;
fas_pps="%upcase(&fas_pps.)";proj="&label."; label proj='指标';
label='值';value=put(sffixed,10.3 -l);output;
label='95%置信区间';value=cats('(',put(sffixed_lower,10.3),',',put(sffixed_upper,10.3),')');output;
keep proj label value fas_pps;
run;
title;
proc append base=output data=report force;run;
%mend intracc;
Hello @Ksharp ,
Thank you very much for sharing these resources!
I was able to run the %INTRACC macro to obtain the 6 ICC measures from Shrout and Fleiss. The macro to obtain the 95% CI does not align with the ICC (3,1) measure; are you aware of any existing code to calculate this? I tried calculating it with the code below (to align with Shrout and Fleiss 1979 seminal paper), however with no success.
DATA new;
f=finv(.99975, 135, 270);
fo=(0.1670116195/0.0242011619);
fl=(fo/f)*(135);
fu=((fo)*(f))*((135)*(2));
lowerci=(fl-1)/(fl+(2));
upperci=(fu-1)/(fu+(2));
PUT f= fo= fl= fu= lowci= upperci=;
RUN;
What does "with no success" mean?
Your PUTs a variable named LOWCI. You create with a calculation LOWERCI.
Change the PUT to use the correct variable and that will appear in the output requested. Since you have to look in the LOG to see the output did you not see the note
NOTE: Variable lowci is uninitialized.
@Adele3 wrote in a different thread:
Hello @FreelanceReinh,
Your response to this post has been very helpful as I am trying to calculate the 95% CI for the Shrout Fleiss ICC (3,1).
When calculating the relevant F quantiles using the FINV function, it is not clear what to input as the numerator degrees of freedom and denominator degrees of freedom. Are you able to provide clarity?
@Adele3 wrote:
(...) The macro to obtain the 95% CI does not align with the ICC (3,1) measure; ... I tried calculating it with the code below (to align with Shrout and Fleiss 1979 seminal paper), however with no success.
DATA new; f=finv(.99975, 135, 270); fo=(0.1670116195/0.0242011619); fl=(fo/f)*(135); fu=((fo)*(f))*((135)*(2)); lowerci=(fl-1)/(fl+(2)); upperci=(fu-1)/(fu+(2)); PUT f= fo= fl= fu= lowci= upperci=; RUN;
The formulas for lowerci and upperci match those from the Shrout and Fleiss 1979 paper, with k = 3 (i.e., three raters). Your definition of f suggests that n = 136 (i.e., 136 targets). It also suggests that your confidence level is 1−a=0.9995, which seems unusual to me. So you don't want the most common a=0.05, i.e., a 95% confidence interval?
Your definition of fo suggests that BMS=0.1670116195 and EMS=0.0242011619 (in the paper's notation).
The formulas for FL and FU in the paper (equations (8) and (9) on p. 424), however, do not match your definitions: They are using two different F quantiles, whereas you use the same and then you multiply with the some degrees of freedom (why?).
So I would have rather expected to see something like this:
n=136; k=3; alpha=0.05; fl=fo/finv(1-alpha/2, n-1, (n-1)*(k-1)); fu=fo*finv(1-alpha/2, (n-1)*(k-1), n-1);
Hello @FreelanceReinh,
Thank you for taking the time to help and answer my question. I am able to obtain the CI's with the help of your code! Many thanks again!
OK. The attachment is the original paper you mentioned .
The following marco is I wrote before. Hope could give you a hand.
%macro intracc(data=,target=,rater=, depvar= ,df=,fas_pps=,label=);
proc glm data=&data outstat=_stats_ noprint ;
class &target &rater;
model &depvar = &target &rater ;
run;
proc sort data=_stats_;
by _name_ _SOURCE_;
run;
data want;
retain msw msb wms ems edf bms bdf jms jdf k;
set _stats_;
by _name_;
if upcase(_type_)='SS1' then delete;
if upcase(_source_)='ERROR' then do;
ems=ss/df;
edf=df;
end;
if upcase(_source_)="%upcase(&target)" then do;
bms=ss/df;
msb=bms;
bdf=df;
end;
if upcase(_source_)="%upcase(&rater)" then do;
jms=ss/df;
jdf=df;
k=df+1;
end;
if last._name_ then do;
msw=((ems*edf)+(jms*jdf))/(edf+jdf);
wms=msw;
n=bdf+1;
sffixed=(bms-ems)/(bms+((k-1)*ems)); /* ICC(3,1);*/
/*95% Confidence Interval for sfsingle - ICC(3,1)
Suppose that each of N subjects (targets) yields K observations.*/
F0=bms/ems;
V1=(n-1)*(k-1);
V2=n-1;
FL=F0/quantile('F',.975,V2,V1);
FU=F0*quantile('F',.975,V1,V2);
sffixed_lower=(FL-1)/(FL+k-1); /* ICC(3,1) lower CI*/
sffixed_upper=(FU-1)/(FU+k-1); /* ICC(3,1) upper CI*/
output;
end;
label sffixed="Shrout-Fleiss reliability: fixed set";
run;
data report;
set want;
length fas_pps label proj value $ 80;
fas_pps="%upcase(&fas_pps.)";proj="&label."; label proj='指标';
label='值';value=put(sffixed,10.3 -l);output;
label='95%置信区间';value=cats('(',put(sffixed_lower,10.3),',',put(sffixed_upper,10.3),')');output;
keep proj label value fas_pps;
run;
title;
proc append base=output data=report force;run;
%mend intracc;
Hello @Ksharp,
Thank you very much for taking the time to provide your help and sharing code; I am able to obtain the CI's! Your code helped me better understand the methodology behind the formula by Shrout and Fleiss.
Many thanks again!
SAS Innovate 2025 is scheduled for May 6-9 in Orlando, FL. Sign up to be first to learn about the agenda and registration!
Learn how use the CAT functions in SAS to join values from multiple variables into a single value.
Find more tutorials on the SAS Users YouTube channel.
Ready to level-up your skills? Choose your own adventure.