[R] Interrater and intrarater reliability

[Tore Wentzel-Larsen]

Dear R users

The following function is R code for the main compuations in the article:
M. Eliasziw, S Lorraine Young, M Gail Woodbury and Karen Fryday-Field (1994):
Statistical Methodology for the Concurrent Assessment of Intrarater and 
Intrarater Reliability: Using Goniometric Measurements as an Example.
Physical Therapy 74 (8); 777-788

The function gives the estimated inter- and intrarater reliabilities
(rhohat) for fixed and random rater effects, partial intrarater
reliabilities for each rater, F test statistics with corresponding 
p-values, lower bounds for one-sided confiodence intervals and
standard errors of measurement.

The defaults are set up for use in a current project.
Script for testing in the example in the article is included below the function.

Tore Wentzel-Larsen
statistician
Centre for Clinical Research
Haukeland University Hospital
N-5021 Bergen, Norway
email: tore.wentzel-larsen at helse-bergen.no

the function:
---------------------------------------------------------------------------------
relInterIntra<-function(Results,nsubj=40,nrater=3,nmeas=2,raterLabels=c('a','b','c'),rho0inter=0.6,rho0intra=.8,conf.level=.95)
{
# gives the reliability coefficients in the article by Eliasziw et. al. 1994; Phys. Therapy 74.8; 777-788.
# all references in this function are to this article.
# input: Results, data frame representating a data structure as in Table 1 (p. 779), 
#	with consecutive measurements for each rater in adjacent columns.
# rho0inter: null hypothesis value of the interrater reliability coefficient
# rho0intra: null hypothesis value of the intrarater reliability coefficient
# conf.level: confidence level of the one-sided confidence intervals reported for the reliability coefficients
Frame1<-data.frame(cbind(
	rep(1:nsubj,nrater*nmeas),rep(1:nrater,rep(nsubj*nmeas,nrater)),
	rep(rep(1:nmeas,rep(nsubj,nmeas)),nrater),
	matrix(as.matrix(Results),ncol=1)))
names(Frame1)<-c('Subject','Rater','Repetation','Result')
Frame1$Subject<-factor(Frame1$Subject)
Frame1$Rater<-factor(Frame1$Rater,labels=raterLabels)
Frame1$Repetation<-factor(Frame1$Repetation)
nn<-nsubj # this and following two commands: aliases for compatibility with Eliasziw et. al. notation
tt<-nrater
mm<-nmeas
aovFull<-aov(Result~Subject*Rater,data=Frame1)
meanSquares<-summary(aovFull)[[1]][,3]
for (raterAct in 1:tt)
{
raterActCat<-raterLabels[raterAct]
aovAct<-aov(Result~Subject,data=Frame1[Frame1$Rater==raterActCat,])
meanSquares<-c(meanSquares,summary(aovAct)[[1]][2,3])
}
names(meanSquares)<-c('MSS','MSR','MSSR','MSE',paste('MSE',levels(Frame1$Rater),sep=''))

MSS<-meanSquares[1]
MSR<-meanSquares[2]
MSSR<-meanSquares[3]
MSE<-meanSquares[4]
MSEpart<-meanSquares[-(1:4)] # the same for random and fixed, see table 2 (p. 780) and 3 (p. 281)
sighat2Srandom<-(MSS-MSSR)/(mm*tt)
sighat2Rrandom<-(MSR-MSSR)/(mm*nn)
sighat2SRrandom<-(MSSR-MSE)/mm
sighat2e<-MSE # the same for random and fixed, see table 2 (p. 780) and 3 (p. 281)
sighat2Sfixed<-(MSS-MSE)/(mm*tt)
sighat2Rfixed<-(MSR-MSSR)/(mm*nn)
sighat2SRfixed<-(MSSR-MSE)/mm
sighat2e.part<-MSEpart # the same for random and fixed, see table 2 (p. 780) and 3 (p. 281)
rhohat.inter.random<-sighat2Srandom/
	(sighat2Srandom+sighat2Rrandom+sighat2SRrandom+sighat2e)
rhohat.inter.fixed<-(sighat2Sfixed-sighat2SRfixed/tt)/
	(sighat2Sfixed+(tt-1)*sighat2SRfixed/tt+sighat2e)
rhohat.intra.random<-(sighat2Srandom+sighat2Rrandom+sighat2SRrandom)/
	(sighat2Srandom+sighat2Rrandom+sighat2SRrandom+sighat2e)
rhohat.intra.fixed<-(sighat2Sfixed+(tt-1)*sighat2SRfixed/tt)/
	(sighat2Sfixed+(tt-1)*sighat2SRfixed/tt+sighat2e)
rhohat.intra.random.part<-(sighat2Srandom+sighat2Rrandom+sighat2SRrandom)/
	(sighat2Srandom+sighat2Rrandom+sighat2SRrandom+sighat2e.part)
rhohat.intra.fixed.part<-(sighat2Sfixed+(tt-1)*sighat2SRfixed/tt)/
	(sighat2Sfixed+(tt-1)*sighat2SRfixed/tt+sighat2e.part)
Finter<-(1-rho0inter)*MSS/((1+(tt-1)*rho0inter)*MSSR)
Finter.p<-1-pf(Finter,df1=nn-1,df2=(nn-1)*(tt-1))
alpha<-1-conf.level
nu1<-
(nn-1)*(tt-1)*
(
tt*rhohat.inter.random*(MSR-MSSR)+
nn*(1+(tt-1)*rhohat.inter.random)*MSSR+
nn*tt*(mm-1)*rhohat.inter.random*MSE
)^2/
(
(nn-1)*(tt*rhohat.inter.random)^2*MSR^2+
(nn*(1+(tt-1)*rhohat.inter.random)-tt*rhohat.inter.random)^2*MSSR^2+
(nn-1)*(tt-1)*(nn*tt*(mm-1))*rhohat.inter.random^2*MSE^2
)
nu2<-
(nn-1)*(tt-1)*
(
nn*(1+(tt-1)*rhohat.inter.fixed)*MSSR+
nn*tt*(mm-1)*rhohat.inter.fixed*MSE
)^2/
(
(nn*(1+(tt-1)*rhohat.inter.fixed))^2*MSSR^2+
(nn-1)*(tt-1)*(nn*tt*(mm-1))*rhohat.inter.fixed^2*MSE^2
)
F1<-qf(1-alpha,df1=nn-1,df2=nu1)
F2<-qf(1-alpha,df1=nn-1,df2=nu2)
lowinter.random<-nn*(MSS-F1*MSSR)/
(
nn*MSS+F1*(tt*(MSR-MSSR)+
nn*(tt-1)*MSSR+nn*tt*(mm-1)*MSE)
)
lowinter.random<-min(c(lowinter.random,1))
lowinter.fixed<-nn*(MSS-F2*MSSR)/
(
nn*MSS+F2*
(nn*(tt-1)*MSSR+nn*tt*(mm-1)*MSE)
)
lowinter.fixed<-min(c(lowinter.fixed,1))
Fintra<-(1-rho0intra)*MSS/((1+(mm-1)*rho0intra)*MSE*tt)
Fintra.p<-1-pf(Fintra,df1=nn-1,df2=nn*(mm-1))
Fintra.part<-(1-rho0intra)*MSS/((1+(mm-1)*rho0intra)*MSEpart*tt)
Fintra.part.p<-1-pf(Fintra.part,df1=nn-1,df2=nn*(mm-1))
F3<-qf(1-alpha,df1=nn-1,df2=nn*(mm-1))
lowintra<-(MSS/tt-F3*MSE)/(MSS/tt+F3*(mm-1)*MSE)
lowintra<-min(c(lowintra,1))
F4<-qf(1-alpha,df1=nn-1,df2=nn*(mm-1))
lowintra.part<-(MSS/tt-F4*MSEpart)/(MSS/tt+F4*(mm-1)*MSEpart)
for (raterAct in 1:tt) lowintra.part[raterAct]<-min(lowintra.part[raterAct],1)
SEMintra<-sqrt(MSE)
SEMintra.part<-sqrt(MSEpart)
SEMinter.random<-sqrt(sighat2Rrandom+sighat2SRrandom+sighat2e)
SEMinter.fixed<-sqrt(sighat2SRfixed+sighat2e)
rels<-c(rhohat.inter.random,rhohat.intra.random,rhohat.inter.fixed,rhohat.intra.fixed,
rhohat.intra.random.part,rhohat.intra.fixed.part,
Finter,Finter.p,Fintra,Fintra.p,Fintra.part,Fintra.part.p,
lowinter.random,lowinter.fixed,lowintra,lowintra.part,
SEMintra,SEMintra.part,SEMinter.random,SEMinter.fixed)
names(rels)<-c('rhohat.inter.random','rhohat.intra.random','rhohat.inter.fixed','rhohat.intra.fixed',
paste('rhohat.intra.random.part',raterLabels,sep='.'),paste('rhohat.intra.fixed.part',raterLabels,sep='.'),
'Finter','pvalue.Finter','Fintra','pvalue.Fintra',paste('Fintra',raterLabels,sep='.'),paste('pvalue.Fintra',raterLabels,sep='.'),
'lowinter.random','lowinter.fixed','lowintra',paste('lowintra',raterLabels,sep='.'),
'SEMintra',paste('SEMintra.part',raterLabels,sep='.'),'SEMinter.random','SEMinter.fixed')
rels
} # end of function relInterIntra
--------------------------------------------------------------------------------------------------------


testing code for the Goniometer data from the article:
--------------------------------------------------------------------------------------------------------
table4<-matrix(c(
-2,16,5,11,7,-7,18,4,0,0,-3,3,7,-6,1,-13,2,4,-10,8,7,-3,-5,5,0,7,-8,1,-3,
0,16,6,10,8,-8,19,5,-3,0,-2,-1,9,-7,1,-14,1,4,-9,9,6,-2,-5,5,-1,6,-8,1,-3,
1,15,6,10,6,-8,19,5,-2,-2,-2,1,9,-6,0,-14,0,3,-10,8,7,-4,-7,5,-1,6,-8,2,-3,
2,12,4,9,5,-9,17,5,-7,1,-4,-1,4,-8,-2,-12,-1,7,-10,2,8,-5,-6,3,-4,4,-10,1,-5,
1,14,4,7,6,-10,17,5,-6,2,-3,-2,4,-10,-2,-12,0,6,-11,8,7,-5,-8,4,-3,4,-11,-1,-4,
1,13,4,8,6,-9,17,5,-5,1,-3,1,2,-9,-3,-12,0,4,-10,8,7,-5,-7,4,-4,4,-10,0,-5
),ncol=6)
relIIgon<-relInterIntra(Results=table4,nsubj=29,nrater=2,nmeas=3,raterLabels=c('universal','Lamoreux'))
relIIgon
--------------------------------------------------------------------------------------------------------