|
I am amazed by the prompt responses to my mail requesting the program routine for calculating kappa and weighted kappa in S-PLUS. A BIG Thank You to all of you who responded.
Neema
I have listed the responses below:
1. Chuck Cleland -
I have attached a file containing a kappa function I picked up long ago. I can't remember where it came from.
<<skappa.q>>
***************************************************************************
2.James L. Pratt, SONUS Pharmaceuticals, Inc.
Here is a program written to cpmuted Weighted kappas for any size matrix. I've sent it to MathSoft, hoping they'll add to future versions of S-Plus. It doesn't compute 95% confidence intervals, but does compute/return asymptotic standard errors from which confidence intervals may be obatined (I've not looked into it).
Below is the basic function I've used to compute kappa statistics in my
situation. It runs okay, but is slowed down by my call to crosstabs within
kappa.counts function call to implicitly compute counts variable. I think
the code coould be made to run more efficiently, but have not taken the time
to do so. Also, when one reviews Fleiss' book, will notice there are other
weighting schemes. The one used here is for linear weights. Perhaps the
creation of the weight matrix w can be accomplished without for loops.
kappa.counts<-function(counts,wt=T){
##counts is a k-by-k contingency table
n.tot<-sum(counts)
k<-dim(counts)[1]
if(wt){
w<-matrix(1,nrow=k,ncol=k)
for (i in 1:(k-1)){
for (j in (i+1):k){
w[i,j]<-1-abs(i-j)/(k-1)
w[j,i]<-w[i,j]
}
}
}
else {
w<-matrix(0,nrow=k,ncol=k)
for (i in 1:k) w[i,i]<-1
}
print(w) #w is weight matrix, identity matrix if wt=F
P<-counts/n.tot #k-by-k matrix of proportions Pij
Pi<-rowSums(P) #k-by-1 vector of marginal proportions Pi.
Pj<-colSums(P) #k-by-1 vector of marginal proportions P.j
#observed proportion of agreement (weighted if needed)
po<-sum(P*w)
# * is element-wise multiplication, %*% is linear algebra multiplication
# t() is transpose function
# expected proportion of agreement (weighted if needed)
pe<-sum((Pj%*%t(Pi))*w)
kappa<-(po-pe)/(1-pe)
#outer(x,y,"+") creates outer-sum of two matrices x and y, similar to a
#cross-product, adding elemtents rather than multiplying them
# part 1 of numerator for variance estimate
var1<-sum(P*(w-matrix(outer(Pj%*%t(w),Pi%*%w,
FUN="+"),ncol=k,byrow=F)*(1-kappa))^2)
# part 2 of numerator for variance estimate
var2<-(kappa-pe*(1-kappa))^2
#denominator for variance estiamte
den<-(n.tot*(1-pe)^2)
#asymptotic se
ase<-sqrt((var1-var2)/den)
# part 1 of numerator for variance estimate under null hyp
var1.0<-sum((Pj%*%t(Pi))*(w-matrix(outer(Pj%*%t(w),Pi%*%w,
FUN="+"),ncol=k,byrow=T))^2)
#asymptotic se under null hyp
ase.0<-sqrt((var1.0-pe^2)/den)
print(counts)
cbind(po,pe,kappa,var1,var2,den,ase,z.value=kappa/ase,
p.value=2*(1-pnorm(abs(kappa/ase))),ase.0,z.value.0=kappa/ase.0,
p.value.0=2*(1-pnorm(abs(kappa/ase.0))))
}
kappa.counts(crosstabs(~factor(NEA,levels=0:3)+factor(NEB,levels=0:3),
data="">
drop.unused.levels=F,margin=list("N/Total" = integer(0))))
It is important that the counts argument to kappa.counts be a square matrix.
This is why have chosen the drop.unused.levels=F in the call.
***************************************************
3. recently found a kappa function in the archives. It's quite nice, as it
calculates both weighted and unweighted kappas. You can find the function
at:
http://www.biostat.wustl.edu/hyperlists/s-news/199810/msg00089.html
<http://www.biostat.wustl.edu/hyperlists/s-news/199810/msg00089.html>
sally
*******************************************************
4. Not sure if this is what you want--it does the overall kappa and individual
stats but you would have to add weight assignments and CI. The input
file is not a matrix, but a two column file with the raw data for
categories (it creates the matrix from the input file). In my case I was
looking at agreement between cover type recorded on the ground (column
1), and cover type recorded in an image classification (column 2).
Credit to Robin Reich at Colorado State, who wrote the code.
Sandra L. Haire
USGS-Biological Resources Division
"kap"<-
function(data)
{
data.tbl<-table(data[,1],data[,2])
rt <- matrix(apply(data.tbl, 1, sum), ncol = 1)
ct <- matrix(apply(data.tbl, 2, sum), nrow = 1)
gt <- sum(data.tbl)
dgt <- sum(diag(data.tbl))
pm <- rt %*% ct
kappa <- 1 - (gt * (gt - dgt))/(sum(pm) - sum(diag(pm)))
nr<-length(data.tbl[,1])
kap1<-matrix(0, ncol=1, nrow=nr)
overall<-dgt/gt
for(i in 1:nr) {
kap1[i]<-(gt*data.tbl[i,i]-rt[i]*ct[i])/(gt*rt[i]-rt[i]*ct[i])
}
cat("error matrix \n")
print(data.tbl)
cat("\n")
cat("\n Overall accuracy ",round(overall,4), "\n")
cat("\n Overall Kappa is ", round(kappa, 4), "\n")
cat("Individual Kappa Statistics \n")
kap1
}
The information or documents contained in this electronic mail message is intended to be privileged and confidential information only for the use of the individual or entity named herein or above. If you are not the intended recipient, you are hereby notified that any dissemination, distribution or copy of this message or document is strictly prohibited. If you have received this electronic mail transmission in error, notify the sender by reply e-mail and delete all copies from your system. Erlanger Health System is not responsible for errors in this electronic mail message. Any personal comments made do not necessarily reflect the views of Erlanger Health System.
skappa.q
Description: Binary data
|