ac1_agreement_statistics.R

gistfile1.r

# AC1 statistic for 2 raters special case
# table = k x k table which represents table(rater1,rater2), must have equal number of rows and columns
# N = population size which will be stick in standard error correction, N=Inf is no correction.
# conflev = Confidence Level associated with the confidence interval (0.95 is the default value)

AC1 <- function(table,conflev=0.95,N=Inf,print=TRUE){
if(dim(table)[1] != dim(table)[2]){
stop('The table should have the same number of rows and columns!')
}
n <- sum(table)
f <- n/N
pa <- sum(diag(table))/n # formula 18
q <- ncol(table) # number of categories
pkk <- diag(table)/n
pak <- sapply(1:q,function(i)sum(table[i,]))/n
pbk <- sapply(1:q,function(i)sum(table[,i]))/n
pik <- (pak + pbk)/2
pegama <- (sum(pik*(1-pik)))/(q-1)
gama <- (pa - pegama)/(1 - pegama) # AC1 statistics
# 2 raters special case variance
pkl <- table/n
soma <- 0;
for(k in 1:q){
for(l in 1:q){
soma <- soma + (pkl[k,l]*((1-(pik[k]+pik[l])/2)^2))
}
}
vgama <- ((1-f)/(n*(1-pegama)^2)) * (pa*(1-pa) - 4*(1-gama)*((1/(q-1))*sum(pkk*(1-pik)) - pa*pegama) + 4*((1-gama)^2) * ((1/((q-1)^2))*soma - pegama^2))
epgama <- sqrt(vgama)# AC1 standard error
lcb <- max(0,gama - epgama*qnorm(1-(1-conflev)/2,0,1)) # lower confidence bound
ucb <- min(1,gama + epgama*qnorm(1-(1-conflev)/2,0,1)) # upper confidence bound
if(print==TRUE){cat('Raw agreement:',pa,'Chance-independent agreement:',pegama,'\n')
cat('Agreement coeficient (AC1):',gama,'AC1 standard error:',epgama,'\n') 
cat(conflev*100,'% Confidence Interval (AC1): (',lcb,',',ucb,')\n')
}
invisible(c(pa,pegama,gama,epgama,lcb,ucb))
}

 table3 <- matrix(c(118,2,5,0),nrow=2,ncol=2)
# AC1(table3)
# x <- AC1(table3,print=F)
# print(x)