[S] Kappa

Dr. L. Y. Hin (z044106@mailserv.cuhk.edu.hk)
Fri, 16 Oct 1998 02:13:26 +0800


Dear S-plusers,

Thanks to Melania who kindly offered the following function:

function(x,y,n,weight=T,file="kappa.list",t1="x and y")
{
############################################################################
#
## Reference: J.L.Fleiss , Statistical Method for rates and proportions
Wiley
## Chapter 13
## used weights : Fleiss and Cohen 1973
############################################################################
#
## x - categories through the first method
## y - categories through the second method
## n - number of categories
## when the categories are ordered (ex. stage)
## the weight should be = T (ref: A.Agresti:Categorical Data Analysis
## page 367, 10.5.2 Generalized Kappa measures)
#################################################################
by_table(x,y)
#n_max(ncol(by),nrow(by))

### makes the table a square matrix
square_function(a,M)
{
#################
# adds extra row and columns such that
# the matrix is MxM
###################
nM_ncol(a)
mM_nrow(a)
if (dimnames(a)[[2]][1]!=1)
{a_cbind(0,a)
nM_nM+1
dimnames(a)[[2]][1]_1
}
if (dimnames(a)[[1]][1]!=1)
{a_rbind(0,a)
mM_mM+1
dimnames(a)[[1]][1]_1
}
if (dimnames(a)[[2]][nM]!=M)
{a_cbind(a,0)
nM_nM+1
dimnames(a)[[2]][nM]_M
}
if (dimnames(a)[[1]][mM]!=M)
{a_rbind(a,0)
mM_mM+1
dimnames(a)[[1]][mM]_M
}
if (nM!=M)
{for (i in 1:M)
{ if (dimnames(a)[[2]][i]!=i)
{ dummy_vector(mode="numeric",length=mM)
l_ncol(a)
t_i-1
u_i+1
v_l+1
nam1_dimnames(a)[[2]][1:t]
nam2_dimnames(a)[[2]][i:l]
a_cbind(a[,1:t],dummy,a[,i:l])
dimnames(a)[[2]][i]_i
dimnames(a)[[2]][1:t]_nam1
dimnames(a)[[2]][u:v]_nam2
}
}
}
if (mM!=M)
{for (i in 1:M)
{if (dimnames(a)[[1]][i]!=i)
{dummy_vector(mode="numeric",length=M)
l_nrow(a)
t_i-1
u_i+1
v_l+1
nam1_dimnames(a)[[1]][1:t]
nam2_dimnames(a)[[1]][i:l]
a_rbind(a[1:t,],dummy,a[i:l,])
dimnames(a)[[1]][i]_i
dimnames(a)[[1]][1:t]_nam1
dimnames(a)[[1]][u:v]_nam2
}
}
}
return(a)
}

a_square(by,n)

p_a/sum(a)
w_matrix(0,ncol=n,nrow=n)
for (i in 1:n)
{w[i,i]_1
}
if (weight)
{for (i in 1:n)
{for (j in 1:n)
{w[i,j]_1-((i-j)**2/(n-1)**2)
}
}
}
po_sum(w*p)
hori_apply(p,1,sum)
vert_apply(p,2,sum)
pe_matrix(0,ncol=1,nrow=n)
w1_w%*%vert
w2_hori%*%w
pv_matrix(0,ncol=n,nrow=1)
for (i in 1:n)
{pe_cbind(pe,vert[i]*hori)
pv_rbind(pv,w1[i]+w2)
}
pe_pe[,-1]
pv_pv[-1,]
pv_w-pv
pv_sum(pe*pv*pv)
pe_sum(pe*w)
kw_(po-pe)/(1-pe)
se_(sqrt(pv-pe**2))/((1-pe)*sqrt(sum(a)))
ci1_kw-1.96*se
ci2_kw+1.96*se
storage.mode(kw)_"single"
storage.mode(ci1)_"single"
storage.mode(ci2)_"single"
storage.mode(se)_"single"
sink(file,append=T)
cat(t1,"\n")
cat("\n","\n")
cat("A Kappa of 1 shows perfect concordance","\n")
cat("Below 0.40 is a poor concordance","\n")
cat("Negative values show that agreement is weaker than expected by
chance","\n")
cat("The table used :","\n")
for (i in 1:n)
{cat(" ",a[i,],"\n")
}
cat("Kappa=",kw,"\n")
cat("Weights=","\n")
for (i in 1:n)
{cat(" ",w[i,],"\n")
}
cat("Standard error=",se,"\n")
cat("95% confidence intervals :",ci1," ; ",ci2,"\n")
cat("\n","\n","\n","\n")
sink()
return(kw,se)
}

Gratitude goes to Lynd for the following function:

"cohen.kappa"<-
function(xtab)
{
#computes Cohen's kappa and variance estimates, 95% CI
#from a symmetric agreement table xtab
#return value is a list with elements kap and vark
#See Bishop et al. Disc. Mult. Analysis pp. 395-397
#DISCLAIMER: The contents of the materials made available here are
#not provided #with a warranty of any kind. LYND BACON & ASSOCIATES,
#LTD DOES NOT MAKE AND #HEREBY DISCLAIMS ANY EXPRESS OR IMPLIED
#WARRANTIES INCLUDING, BUT NOT LIMITED #TO, THE WARRANTIES OF
#NON-INFRINGEMENT, MERCHANTABILITY OR FITNESS FOR A #PARTICULAR
#PURPOSE, OR ARISING FROM A COURSE OF DEALING, USAGE OR TRADE
#PRACTICE.

ci <- vector(length = 2)
kap <- NULL
vark <- NULL
totn <- NULL
if(nrow(xtab) != ncol(xtab)) {
cat("\n freq. table not symmetric.\n")
return(kap, vark, totn)
}
# compute obs props:
totn <- sum(xtab)
obs <- xtab/totn # calc marginals, expected diag props.
pi <- apply(xtab, 1, sum)/totn
pj <- apply(xtab, 2, sum)/totn
exp <- outer(pi, pj)
theta.2 <- sum(diag(exp))
theta.1 <- sum(diag(obs))
kap <- (theta.1 - theta.2)/(1 - theta.2)
theta.3 <- diag(obs) %*% (pi + pj)
theta.4 <- sum(obs * outer(pi, pj, FUN = "+")^2)
# cat("\n thetas 1,2,3,4 \n", theta.1, theta.2, theta.3, theta.4, "\n")
# calc the var est
vark <- (theta.1 * (1 - theta.1))/(1 - theta.2)^2
vark <- vark + (2 * (1 - theta.1) * (2 * theta.1 * theta.2 - theta.3))/(1
- theta.2)^3
vark <- vark + ((1 - theta.1)^2 * (theta.4 - 4 *theta.2^2))/(1 -theta.2)^4
vark <- as.vector(vark/totn)
ci[1] <- kap - 1.96 * sqrt(vark)
ci[2] <- kap + 1.96 * sqrt(vark)
return(kap, vark, totn, ci)
}

Hope these functions are useful reference to others.

Sincerely,

Lin

--------------------------------------------------------------------
Dr. L. Y. Hin
Department of Obstetrics & Gynaecology
Prince of Wales Hospital
Hong Kong SAR
Tel. (852) 2632 2807
Fax. (852) 2636 0008
--------------------------------------------------------------------

-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news