####kappa.u: A function for calculating kappa statistics based on U-statistics#### ####reference: Ma Y, Tang W, Feng C, and Tu X.M. "Inference for Kappas for Longitudinal Study Data: Applications to Sexual Health Research" Biometrics 2008 64: 781-789. #####reuqired parameters####### ### n: number of subjects ### g: number of categories ### t: number of repeated measurements (t=2 or 3 for the currently version) ### x: a n*t matrix of ratings of the first rater ### y: a n*t matrix of ratings of the second rater ### weight: weighted function, 1 for unweighted kappa; 2 for Cicchetti-Allison weighted kappa; 3 for Fleiss-Cohen weighted kappa #####results##### ### u.kappa: U-statistic based kappa estimate ### standard error: standard error of the kappa estimate ### test statistic: test statistic for testing equal kappas over time ### p-value: p-value for the test of equal kappas kappa.u=function(n, g, t, x, y, weight) { if (weight==1){ w=matrix(0,ncol=g, nrow=g) diag(w)=1 } if (weight==2){ w=matrix(0,ncol=g, nrow=g) for (i in 1:g) { for (j in 1:g) { w[i,j]=1-abs(i-j)/(g-1) } } } if (weight==3){ w=matrix(0,ncol=g, nrow=g) for (i in 1:g) { for (j in 1:g) { w[i,j]=1-(i-j)^2/(g-1)^2 } } } if (t==2){ m=matrix(0,nrow=(n*t), ncol=(2*g)) for (i in (1:n)) { m[i,x[i,1]]=1 m[(n+i),x[i,2]]=1 m[i,(g+y[i,1])]=1 m[(n+i),(g+y[i,2])]=1 } 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=matrix(nrow=(2*g+g*g), ncol=1) for (i in 1:g) { p12[i,]=mean(x12[,i]) } for (i in (g+1):(2*g)) { p12[i,]=mean(y12[,(i-g)]) } for (i in 1:(g*g)) { p12[(i+2*g),]=xy12[i]/n } 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 } 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 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 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 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 thta=2/(n*(n-1))*matrix(c(u121,u122,u341,u342), nrow=4,ncol=1) 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 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 eu=cbind(eu121,eu122,eu341,eu342) eu=t(eu) usig=matrix(nrow=2*t, ncol=n*2*t) for (i in 1:n) { usig[,((i-1)*(2*t)+1):(i*(2*t))]=(eu[,i]-thta)%*%t(eu[,i]-thta) } usigm=matrix(nrow=2*t, ncol=2*t) usigm=usig[,1:(2*t)] for (i in 1:(n-1)) { usigm=usigm+usig[,(i*(2*t)+1):((i+1)*(2*t))] } usigm=4/n*usigm udev=matrix(nrow=t,ncol=2*t) udev[1,]=c(1/thta[2,1], -thta[1,1]/(thta[2,1]^2),0,0) udev[2,]=c(0,0,1/thta[4,1], -thta[3,1]/(thta[4,1]^2)) usigk=udev%*%usigm%*%t(udev) uk12=thta[1,]/thta[2,] uk34=thta[3,]/thta[4,] u.kappa=c(uk12, uk34) se.kappa=sqrt(diag(usigk)/n) #standard error ###Ho: k12=k34 a=matrix(c(1,-1), ncol=2,nrow=1, byrow=T) uk=matrix(nrow=1,ncol=1) uk[1,]=c(uk12-uk34) chi.stat=uk^2/(a%*%(usigk/n)%*%t(a)) pvalue=1-pchisq(chi.stat,1) ukappa=list(c("u.kappa"=u.kappa, "standard error"=se.kappa, "test statistic"=chi.stat, "p-value"=pvalue)) return(ukappa) } if (t==3){ m=matrix(0,nrow=(n*t), ncol=(2*g)) for (i in (1:n)) { m[i,x[i,1]]=1 m[(n+i),x[i,2]]=1 m[(2*n+i),x[i,3]]=1 m[i,(g+y[i,1])]=1 m[(n+i),(g+y[i,2])]=1 m[(2*n+i),(g+y[i,3])]=1 } 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=matrix(nrow=(2*g+g*g), ncol=1) for (i in 1:g) { p12[i,]=mean(x12[,i]) } for (i in (g+1):(2*g)) { p12[i,]=mean(y12[,(i-g)]) } for (i in 1:(g*g)) { p12[(i+2*g),]=xy12[i]/n } 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 } 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 } 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 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 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 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 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 thta=2/(n*(n-1))*matrix(c(u121,u122,u341,u342,u561, u562), nrow=6,ncol=1) 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 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=t(eu) usig=matrix(nrow=2*t, ncol=n*2*t) for (i in 1:n) { usig[,((i-1)*(2*t)+1):(i*(2*t))]=(eu[,i]-thta)%*%t(eu[,i]-thta) } usigm=matrix(nrow=2*t, ncol=2*t) usigm=usig[,1:(2*t)] for (i in 1:(n-1)) { usigm=usigm+usig[,(i*(2*t)+1):((i+1)*(2*t))] } usigm=4/n*usigm udev=matrix(nrow=t,ncol=2*t) 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)) usigk=udev%*%usigm%*%t(udev) uk12=thta[1,]/thta[2,] uk34=thta[3,]/thta[4,] uk56=thta[5,]/thta[6,] u.kappa=c(uk12, uk34, uk56) se.kappa=sqrt(diag(usigk)/n) #standard error ###Ho: k12=k34=k56 a=matrix(c(1,-1,0,1,0,-1), ncol=3,nrow=2, byrow=T) uk=matrix(nrow=1,ncol=2) uk[1,]=c(uk12-uk34, uk12-uk56) chi.stat=uk%*%solve(a%*%(usigk/n)%*%t(a))%*%t(uk) pvalue=1-pchisq(chi.stat,2) ukappa=list(c("kappa at t=1 (k1)"=u.kappa[1], "kappa at t=2 (k2)"=u.kappa[2], "kappa at t=3 (k3)"=u.kappa[3], "standard error of k1"=se.kappa[1], "standard error of k2"=se.kappa[2], "standard error of k3"=se.kappa[3], "test statistic for testing H0:k1=k2=k3"=chi.stat, "p-value for testing H0:k1=k2=k3"=pvalue)) return(ukappa) } }