####################################################################
# PROGRAM FOR COMPUTING KAPPA ESTIMATE AND ITS ASYMPTOTIC VARIANCE #
# USING U-STATISTICS BASED METHOD UNDER MISSING COMPLETE AT RANDOM #
# 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+1) 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                   #
# r1, r2: mssing indicators for the first pair of judges                      #
# r3, r4: mssing indicators for the second pair of judges                     #
# r5, r6: mssing indicators for the third pair of judges                      #
# n12: total number of observed from the first pair of judges                 #
# n34: total number of observed from the second pair of judges                #
# n56: total number of observed from the third pair of judges                 #
# 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   #
# uk:The kappa estimates                                                      #
# stduk:The asymptotic standard errors of kappa estimates                     #      
###############################################################################

m=read.table("e:/yan ma/kappa/HIV data/5dian/missgrp2/merge5grp2.txt")
n=30
g=6
pair=3
#################################################################################################
# This is to show part of the data set, in which the first 6 cols are obs from one judge 
# and the next 6 cols are from the other judge and the last col is missing indicator
# >m
# 0	0	0	0	0	1	0	0	0	0	0	1	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	1	1
# 0	0	0	1	0	0	0	0	0	0	0	1	1
# 0	0	0	0	0	1	0	0	0	0	0	1	1
# .........................
# .........................
##################################################################################################

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


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


r1=m[1:n,(2*g+1)]
r2=r1
n12=sum(r1)

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

p12=matrix(nrow=(2*g+g*g), ncol=1)
for (i in 1:g)
{
p12[i,]=sum(r1*x12[,i])/n12
}
for (i in (g+1):(2*g))
{
p12[i,]=sum(r1*y12[,(i-g)])/n12
}	
for (i in 1:(g*g))
{
p12[(i+2*g),]=xy12[i]/n12
}	




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)


r3=m[(n+1):(2*n),(2*g+1)]
r4=r3
n34=sum(r3)

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

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

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

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


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))

r5=m[(2*n+1):(3*n), (2*g+1)]
r6=r5
n56=sum(r5)

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

p56=matrix(nrow=(2*g+g*g), ncol=1)
for (i in 1:g)
{
p56[i,]=sum(r5*x56[,i])/n56
}
for (i in (g+1):(2*g))
{
p56[i,]=sum(r5*y56[,(i-g)])/n56
}	


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


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


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


#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=r1*r2*(1-1/2*eu122)
eu122


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

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=r3*r4*(1-1/2*eu342)
eu342





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


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=r5*r6*(1-1/2*eu562)
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)]+r1[i]*r2[i]*r1[j]*r2[j]*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)]+r1[i]*r2[i]*r1[j]*r2[j]*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)]+r3[i]*r4[i]*r3[j]*r4[j]*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)]+r3[i]*r4[i]*r3[j]*r4[j]*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)]+r5[i]*r6[i]*r5[j]*r6[j]*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)]+r5[i]*r6[i]*r5[j]*r6[j]*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/((mean(r1))^2),u122/((mean(r1))^2),u341/((mean(r3))^2),u342/((mean(r3))^2),u561/((mean(r5))^2), u562/((mean(r5))^2)), nrow=6,ncol=1)
thta



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

usig12


usigm12=matrix(nrow=2, ncol=2)
usigm12=usig12[,1:2]
for (i in 1:(n-1))
{

usigm12=usigm12+usig12[,(i*2+1):((i+1)*2)]

}
usigm12
usigm12=4*usigm12/sum(r1*r2)
usigm12

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

usig34


usigm34=matrix(nrow=2, ncol=2)
usigm34=usig34[,1:2]
for (i in 1:(n-1))
{

usigm34=usigm34+usig34[,(i*2+1):((i+1)*2)]

}
usigm34
usigm34=4*usigm34/sum(r3*r4)
usigm34


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

usig56


usigm56=matrix(nrow=2, ncol=2)
usigm56=usig56[,1:2]
for (i in 1:(n-1))
{

usigm56=usigm56+usig56[,(i*2+1):((i+1)*2)]

}
usigm56
usigm56=4*usigm56/sum(r5*r6)
usigm56

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

usig1234


usigm1234=matrix(nrow=2, ncol=2)
usigm1234=usig1234[,1:2]
for (i in 1:(n-1))
{

usigm1234=usigm1234+usig1234[,(i*2+1):((i+1)*2)]

}
usigm1234
usigm1234=4*usigm1234/(sum(r1*r2)*sum(r3*r4))*sum(r1*r2*r3*r4)
usigm1234

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

usig1256


usigm1256=matrix(nrow=2, ncol=2)
usigm1256=usig1256[,1:2]
for (i in 1:(n-1))
{

usigm1256=usigm1256+usig1256[,(i*2+1):((i+1)*2)]

}
usigm1256
usigm1256=4*usigm1256/(sum(r1*r2)*sum(r5*r6))*sum(r1*r2*r5*r6)
usigm1256

###
usig3456=matrix(nrow=2, ncol=n*2)
for (i in 1:n)
{
	usig3456[,((i-1)*2+1):(i*2)]=(eu[3:4,i]-thta[3:4,])%*%t(eu[5:6,i]-thta[5:6,])
}

usig3456


usigm3456=matrix(nrow=2, ncol=2)
usigm3456=usig3456[,1:2]
for (i in 1:(n-1))
{

usigm3456=usigm3456+usig3456[,(i*2+1):((i+1)*2)]

}
usigm3456
usigm3456=4*usigm3456/(sum(r3*r4)*sum(r5*r6))*sum(r3*r4*r5*r6)
usigm3456



usigm=matrix(0, nrow=2*pair, ncol=2*pair)
usigm[1:2, 1:2]=usigm12
usigm[3:4, 3:4]=usigm34
usigm[5:6, 5:6]=usigm56

usigm[1:2, 3:4]=usigm1234
usigm[1:2, 5:6]=usigm1256
usigm[3:4, 5:6]=usigm3456
usigm[3:4, 1:2]=t(usigm1234)
usigm[5:6, 1:2]=t(usigm1256)
usigm[5:6, 3:4]=t(usigm3456)


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,]
uk=c(uk12, uk34, uk56)
stduk=sqrt(diag(usigk)/n)
#The kappa estimates:
uk
#The asymptotic standard errors of kappa estimates:
stduk