############################################################################
# PROGRAM FOR COMPUTING KAPPA ESTIMATE AND ITS ASYMPTOTIC VARIANCE         #
# USING CLASSIC APPROACH AND U-STATISTICS BASED METHOD UNDER COMPLETE DATA #
# Developed by Yan Ma                                                      #
# 10/16/2006                                                               #
############################################################################

###############################################################################
# All you need to have before this computation are : m, n, g, and pair        #
# m:data matrix with n rows and (2*g) cols                                    #
# n: number of subjects                                                       #
# g: number of categories                                                     #
# pair:number of pairs of kappas                                              #
# x12, y12: responses from the first pair of judges X and Y                   #
# x34, y34: responses from the second pair of judges X and Y                  #
# x56, y56: responses from the third pair of judges X and Y                   #
# w: weight matrix                                                            #
# p12:vector of cell and marginal proportions from the first pair of judges   #
# p34:vector of cell and marginal proportions from the second pair of judges  #
# p56:vector of cell and marginal proportions from the third pair of judges   #
# k:The classic kappa estimates                                               #
# uk:The u-statistics based kappa estimates                                   #
# stdk:The asymptotic standard errors of classic kappa estimates              #
# stduk:The asymptotic standard errors of u-statistics based kappa estimates  #      
###############################################################################


m=read.table("e:/yan ma/kappa/HIV data/5dian/Cgrp2/cgrp2.txt")
g=6
n=24
pair=3
###########################################################################################
# This is to show part of the data set, in which the first 6 cols are obs from one judge 
# and the last 6 cols are from the other judge.
# >m
# 0	0	0	0	0	1	0	0	0	0	0	1
# 0	0	0	0	0	1	0	0	0	0	0	1
# 0	0	0	0	0	1	0	0	0	0	0	1
# 1	0	0	0	0	0	1	0	0	0	0	0
# 0	0	0	0	0	1	0	0	0	0	0	1
# 1	0	0	0	0	0	1	0	0	0	0	0
# 0	1	0	0	0	0	0	1	0	0	0	0
# 1	0	0	0	0	0	1	0	0	0	0	0
# 0	0	0	0	1	0	0	0	0	0	0	1
# 0	1	0	0	0	0	0	1	0	0	0	0
# ......................
# ......................
############################################################################################


#####k12
x12=m[1:n, 1:g]
x12=as.matrix(x12)


y12=m[1:n, (g+1):(2*g)]
y12=as.matrix(y12)

xy12=t(x12)%*%y12
xy12=as.vector(t(xy12))

##p12={marginal proportions and cell proportions}
p12=matrix(nrow=(2*g+g*g), ncol=1)
##p12[1:g, ]=marginal proportions of X: p121., p122., ...,p12g.
for (i in 1:g)
{
p12[i,]=mean(x12[,i])
}


##p12[(g+1):(2*g), ]=marginal proportions of y:p12.1,p12.2,...,p12.g
for (i in (g+1):(2*g))
{
p12[i,]=mean(y12[,(i-g)])
}	

###product terms: p1211,p1212,...,p121g, p1221,...,p12gg
for (i in 1:(g*g))
{
p12[(i+2*g),]=xy12[i]/n
}	



#weight functions  
#w: identity matrix
w=matrix(0,ncol=g, nrow=g)
diag(w)=1



pi12=0
for (i in 1:g)
{
for (j in 1:g)
{
pi12=pi12+w[i,j]*p12[((i-1)*g+2*g+j),]
}
}

pie12=0
for (i in 1:g)
{
for (j in 1:g)
{
pie12=pie12+w[i,j]*p12[i,]*p12[(j+g),]
}
}
k12=(pi12-pie12)/(1-pie12)
k12



###k34

x34=m[(n+1):(2*n), 1:g]
y34=m[(n+1):(2*n), (g+1):(2*g)]
x34=as.matrix(x34)
y34=as.matrix(y34)

xy34=t(x34)%*%y34
xy34=as.vector(t(xy34))

p34=matrix(nrow=(2*g+g*g), ncol=1)
for (i in 1:g)
{
p34[i,]=mean(x34[,i])
}

for (i in (g+1):(2*g))
{
p34[i,]=mean(y34[,(i-g)])
}	

for (i in 1:(g*g))
{
p34[(i+2*g),]=xy34[i]/n
}	

pi34=0
for (i in 1:g)
{
for (j in 1:g)
{
pi34=pi34+w[i,j]*p34[((i-1)*g+2*g+j),]
}
}

pie34=0
for (i in 1:g)
{
for (j in 1:g)
{
pie34=pie34+w[i,j]*p34[i,]*p34[(j+g),]
}
}
k34=(pi34-pie34)/(1-pie34)
k34




###k56

x56=m[(2*n+1):(3*n), 1:g]
y56=m[(2*n+1):(3*n), (g+1):(2*g)]
x56=as.matrix(x56)
y56=as.matrix(y56)

xy56=t(x56)%*%y56
xy56=as.vector(t(xy56))


p56=matrix(nrow=(2*g+g*g), ncol=1)
for (i in 1:g)
{
p56[i,]=mean(x56[,i])
}

for (i in (g+1):(2*g))
{
p56[i,]=mean(y56[,(i-g)])
}	

for (i in 1:(g*g))
{
p56[(i+2*g),]=xy56[i]/n
}	

pi56=0
for (i in 1:g)
{
for (j in 1:g)
{
pi56=pi56+w[i,j]*p56[((i-1)*g+2*g+j),]
}
}

pie56=0
for (i in 1:g)
{
for (j in 1:g)
{
pie56=pie56+w[i,j]*p56[i,]*p56[(j+g),]
}
}
k56=(pi56-pie56)/(1-pie56)
k56






dev=matrix(0,nrow=pair, ncol=pair*(2*g+g*g))
for (i in 1:g)
{
dev[1,i]=-sum(w[i,]*p12[(g+1):(2*g),])*(1-pi12)/(1-pie12)^2
dev[1,(i+g)]=-sum(w[,i]*p12[1:g,])*(1-pi12)/(1-pie12)^2
dev[2,(i+2*g+g*g)]=-sum(w[i,]*p34[(g+1):(2*g),])*(1-pi34)/(1-pie34)^2
dev[2,(i+3*g+g*g)]=-sum(w[,i]*p34[1:g,])*(1-pi34)/(1-pie34)^2
dev[3,(i+4*g+2*g*g)]=-sum(w[i,]*p56[(g+1):(2*g),])*(1-pi56)/(1-pie56)^2
dev[3,(i+5*g+2*g*g)]=-sum(w[,i]*p56[1:g,])*(1-pi56)/(1-pie56)^2
}
for (i in 1:(g*g))
{
dev[1,(i+2*g)]=as.vector(t(w))[i]/(1-pie12)
dev[2,(i+4*g+g*g)]=as.vector(t(w))[i]/(1-pie34)
dev[3,(i+6*g+2*g*g)]=as.vector(t(w))[i]/(1-pie56)
}
dev



##v12=(x12, y12, x12*y12)
v12=matrix(nrow=n,ncol=2*g+g*g)
for (i in 1:g)
{
v12[,i]=x12[,i]
v12[, (i+g)]=y12[,i]
for (j in 1:g)
{
v12[,(j+(i-1)*g+2*g)]=x12[,i]*y12[,j]
}
}

##vp12=v12-p12
vp12=matrix(nrow=n, ncol=2*g+g*g)
for (i in 1:n)
{
	vp12[i,]=v12[i,]-p12
}



####
v34=matrix(nrow=n,ncol=2*g+g*g)
for (i in 1:g)
{
v34[,i]=x34[,i]
v34[, (i+g)]=y34[,i]
for (j in 1:g)
{
v34[,(j+(i-1)*g+2*g)]=x34[,i]*y34[,j]
}
}
vp34=matrix(nrow=n, ncol=2*g+g*g)
for (i in 1:n)
{
	vp34[i,]=v34[i,]-p34
}

############
v56=matrix(nrow=n,ncol=2*g+g*g)
for (i in 1:g)
{
v56[,i]=x56[,i]
v56[, (i+g)]=y56[,i]
for (j in 1:g)
{
v56[,(j+(i-1)*g+2*g)]=x56[,i]*y56[,j]
}
}

vp56=matrix(nrow=n, ncol=2*g+g*g)
for (i in 1:n)
{
	vp56[i,]=v56[i,]-p56
}


vp=matrix(nrow=3*(2*g+g*g)*n, ncol=1)
for (i in 1:n)
{
	vp[((i-1)*3*(2*g+g*g)+1):(i*3*(2*g+g*g)),1]=c(vp12[i,],vp34[i,],vp56[i,])
}
vp

sig=matrix(nrow=3*(2*g+g*g), ncol=3*(2*g+g*g)*n)
for (i in 1:n)
{
	sig[,((i-1)*3*(2*g+g*g)+1):(i*3*(2*g+g*g))]=vp[((i-1)*3*(2*g+g*g)+1):(i*3*(2*g+g*g)),]%*%t(vp[((i-1)*3*(2*g+g*g)+1):(i*3*(2*g+g*g)),])
}
sig

sigm=matrix(nrow=3*(2*g+g*g), ncol=3*(2*g+g*g))
sigm=sig[,1:(3*(2*g+g*g))]
for (i in 1:(n-1))
{

sigm=sigm+sig[,(i*3*(2*g+g*g)+1):((i+1)*3*(2*g+g*g))]

}
sigm
sigma=sigm/n
sigma
sigk=dev%*%sigma%*%t(dev)
sigk
##standard error
sqrt(diag(sigk)/n)

 

#############alternative U-stat
eu121=rep(0,n)
for (i in 1:g)
{
for (j in 1:g)
{
eu121=eu121+w[i,j]*(x12[,i]*y12[,j]-x12[,i]*p12[(j+g),]-y12[,j]*p12[i,]+p12[((i-1)*g+2*g+j),])
}
}
eu121=eu121/2


#E(U121ij|Z12i)
eu122=rep(0,n)
for (i in 1:g)
{
for (j in 1:g)
{
eu122=eu122+w[i,j]*(x12[,i]*p12[(j+g),]+y12[,j]*p12[i,])
}
}
eu122=1-1/2*eu122


eu341=rep(0,n)
for (i in 1:g)
{
for (j in 1:g)
{
eu341=eu341+w[i,j]*(x34[,i]*y34[,j]-x34[,i]*p34[(j+g),]-y34[,j]*p34[i,]+p34[((i-1)*g+2*g+j),])
}
}
eu341=eu341/2

eu342=rep(0,n)
for (i in 1:g)
{
for (j in 1:g)
{
eu342=eu342+w[i,j]*(x34[,i]*p34[(j+g),]+y34[,j]*p34[i,])
}
}
eu342=1-1/2*eu342




eu561=rep(0,n)
for (i in 1:g)
{
for (j in 1:g)
{
eu561=eu561+w[i,j]*(x56[,i]*y56[,j]-x56[,i]*p56[(j+g),]-y56[,j]*p56[i,]+p56[((i-1)*g+2*g+j),])
}
}
eu561=eu561/2

eu562=rep(0,n)
for (i in 1:g)
{
for (j in 1:g)
{
eu562=eu562+w[i,j]*(x56[,i]*p56[(j+g),]+y56[,j]*p56[i,])
}
}
eu562=1-1/2*eu562
eu=cbind(eu121,eu122,eu341,eu342,eu561,eu562)
eu
eu=t(eu)
eu
dim(eu)




#####
u121=matrix(0,ncol=(n-1),nrow=(n-1))
for (i in 1:(n-1))
{
	for (j in (i+1):n)
	{
       for (k in 1:g)
        {
         for (q in 1:g)
       {
	u121[i,(j-1)]=u121[i,(j-1)]+w[k,q]*(x12[i,k]-x12[j,k])*(y12[i,q]-y12[j,q])  
}
}
}
}

u121=sum(u121)/2
u121

u122=matrix(0,ncol=(n-1),nrow=(n-1))
for (i in 1:(n-1))
{
	for (j in (i+1):n)
	{
for (k in 1:g)
{
for (q in 1:g)
{
	u122[i,(j-1)]=u122[i,(j-1)]+w[k,q]*(x12[i,k]*y12[j,q]+x12[j,k]*y12[i,q])
}	
}
}
}
u122=n*(n-1)/2-sum(u122)/2
u122

###
u341=matrix(0,ncol=(n-1),nrow=(n-1))
for (i in 1:(n-1))
{
	for (j in (i+1):n)
	{
for (k in 1:g)

{
 for (q in 1:g)
{
	u341[i,(j-1)]=u341[i,(j-1)]+w[k,q]*(x34[i,k]-x34[j,k])*(y34[i,q]-y34[j,q])    
}
}
}
}

u341=sum(u341)/2
u341

u342=matrix(0,ncol=(n-1),nrow=(n-1))
for (i in 1:(n-1))
{
	for (j in (i+1):n)
	{
for (k in 1:g)
{
for (q in 1:g)
{

	u342[i,(j-1)]=u342[i,(j-1)]+w[k,q]*(x34[i,k]*y34[j,q]+x34[j,k]*y34[i,q])
	
}
}
}
}

u342=n*(n-1)/2-sum(u342)/2
u342


###
u561=matrix(0,ncol=(n-1),nrow=(n-1))
for (i in 1:(n-1))
{
	for (j in (i+1):n)
	{
for (k in 1:g)
{
for (q in 1:g)
{
	u561[i,(j-1)]=u561[i,(j-1)]+w[k,q]*(x56[i,k]-x56[j,k])*(y56[i,q]-y56[j,q])    
}
}
}
}
u561=sum(u561)/2
u561

u562=matrix(0,ncol=(n-1),nrow=(n-1))
for (i in 1:(n-1))
{
	for (j in (i+1):n)
	{
for (k in 1:g)
{
for (q in 1:g)
{
	u562[i,(j-1)]=u562[i,(j-1)]+w[k,q]*(x56[i,k]*y56[j,q]+x56[j,k]*y56[i,q])
	
}
}
}
}
u562=n*(n-1)/2-sum(u562)/2
u562



thta=2/(n*(n-1))*matrix(c(u121,u122,u341,u342,u561, u562), nrow=6,ncol=1)


#u-stat

usig=matrix(nrow=2*pair, ncol=n*2*pair)
for (i in 1:n)
{
	usig[,((i-1)*(2*pair)+1):(i*(2*pair))]=(eu[,i]-thta)%*%t(eu[,i]-thta)
}

usig

usigm=matrix(nrow=2*pair, ncol=2*pair)
usigm=usig[,1:(2*pair)]
for (i in 1:(n-1))
{

usigm=usigm+usig[,(i*(2*pair)+1):((i+1)*(2*pair))]

}
usigm
usigm=4/n*usigm
usigm


udev=matrix(nrow=pair,ncol=2*pair)
udev[1,]=c(1/thta[2,1], -thta[1,1]/(thta[2,1]^2),0,0,0,0)
udev[2,]=c(0,0,1/thta[4,1], -thta[3,1]/(thta[4,1]^2),0,0)
udev[3,]=c(0,0,0,0,1/thta[6,1], -thta[5,1]/(thta[6,1]^2))
udev
usigk=udev%*%usigm%*%t(udev)



uk12=thta[1,]/thta[2,]
uk34=thta[3,]/thta[4,]
uk56=thta[5,]/thta[6,]

k=c(k12, k34, k56)
uk=c(uk12, uk34, uk56)
#The classic kappa estimates:
k
#The U-statistics based kappa estimates:
uk
#The asymptotic standard error of classic kappa estimates:
stdk=sqrt(diag(sigk)/n)
stdk
#The asymptotic standard error of u-statistics based kappa estimates:
stduk=sqrt(diag(usigk)/n)
stduk