|
Version 1.0 (March 2013)
aickin.alpha
[R] Aickin's Alpha agreement coefficient
The function presented below computes the agreement coefficient α defined inMaximum Likelihood Estimation of Agreement in the Constant Predictive Probability Model, and Its Relation to Cohen's Kappa Mikel Aickin Biometrics 46, 293-302, June 1990.
Menu
Top Syntax
aickin.alpha(n, d=diag(1, nrow=nrow(n), ncol=ncol(n)), epsilon=1e-7, level=0.95)
Top aickin.alpha arguments list
Top Output
Top Example
We illustrate the use of aickin.alpha through the Alcohol consumption classification example presented in Aickin's paper (cited above), Table 4. The R code below first loads the macro definition, then defines both the data and the agreement matrix, and finally calls aickin.alpha function. Note that in that example, the variables registry and interview are considered as agreeing if registry=interview OR (interview=3 AND registry ∈ {2,4}).
> source('c:/users/pbelisle/My Documents/Home/bin/R/fcts/AickinAlpha.R') # load file where aickin.alpha R function is defined
> # Table 3 data > n3 <- c(60, 2, 7, 0, 2, 0, 1, 5, 67, 2, 0, 1, 2, 0, 1, 0, 67, 0, 0, 2, 0, 0, 1, 1, 27, 0, 1, 0, 3, 1, 3, 1, 64, 5, 0, 1, 2, 6, 1, 7, 82, 2, 0, 0, 0, 0, 0, 0, 13) > n3 <- matrix(n3, ncol=sqrt(length(n3))) > d3 <- diag(1, nrow=nrow(n3), ncol=ncol(n3)) > d3[6,5] <- 1 > d3[5,6] <- 1 > aickin.alpha(n3, d3)
Top Code
aickin.alpha <- function(n, d=diag(1, nrow=nrow(n), ncol=ncol(n)), epsilon=1e-7, level=0.95)
{ # Version 1.0 (March 2013) # # This function computes the alpha coefficient described in: # Maximum Likelihood Estimation of Agreement in the Constant Predictive Probability Model, and Its Relation to Cohen-s Kappa # Mikel Aickin # Biometrics 46, 293-302, June 1990. # # # Function developed by # Lawrence Joseph and Patrick Bélisle # Division of Clinical Epidemiology # Montreal General Hospital # Montreal, Qc, Can # # patrick.belisle@rimuhc.ca # http://www.medicine.mcgill.ca/epidemiology/Joseph/PBelisle/Aickin-Alpha-Agreement-R.html # # Please refer to our webpage for details on each argument. if (any(dim(n) != dim(d))) stop("n and d must be of equal dimensions.") if (diff(dim(n)) != 0) stop("n and d must be square matrices.") m <- nrow(n) n <- n + 1/(m^2) J <- rep(1, m) ssize <- sum(n) A <- sum(n*d) p0 <- A/ssize rows.tot <- as.vector(n%*%matrix(J, ncol=1)) cols.tot <- as.vector(matrix(J, nrow=1)%*%n) pr <- rows.tot/ssize pc <- cols.tot/ssize s <- sum(matrix(pc, nrow=m, ncol=m, byrow=T)*pr*d) alpha <- (p0 - s) / (1 - s) continue <- T while (continue) { previous.alpha <- alpha pr.denominator <- ssize * (1 - alpha + alpha * as.vector(d%*%matrix(pc, ncol=1)) / s) pr <- rows.tot/pr.denominator pr[1] <- 1- sum(pr[-1]) pc.denominator <- ssize * (1 - alpha + alpha * as.vector(matrix(pr, nrow=1)%*%d) / s) pc <- cols.tot/pc.denominator pc[1] <- 1- sum(pc[-1]) s <- sum(matrix(pc, nrow=m, ncol=m, byrow=T)*pr*d) alpha <- (p0 - s) / (1 - s) continue <- abs(alpha-previous.alpha) > epsilon } prdiff <- pr[-1] - pr[1] # m-1 x 1 pcdiff <- pc[-1] - pc[1] # m-1 x 1 d2L.da2 <- - ssize * (1-s)/(1-alpha)/((1-alpha)*s+alpha) R <- alpha/s/((1-alpha)*s+alpha) U <- 1/alpha - 1 d2L.dadpri <- - A * pcdiff * ((ssize/A)^2) # m-1 x 1 d2L.dadpcj <- - A * prdiff * ((ssize/A)^2) # m-1 x 1 d2L.dpridprj <- -sum(n[1,])/(pr[1]^2) + A*(2*s*U+1)*R*R*matrix(pcdiff, ncol=1) %*% matrix(pcdiff, nrow=1) # m-1 x m-1 d2L.dpri2 <- -rows.tot[-1]/(pr[-1]^2) - rows.tot[1]/(pr[1]^2) + A*(2*s*U + 1)*R*R*(pcdiff^2) diag(d2L.dpridprj) <- d2L.dpri2 d2L.dpridpcj <- -A*R + A*(2*s*U+1)*R*R* matrix(pcdiff, ncol=1) %*% matrix(prdiff, nrow=1) # m-1 x m-1 d2L.dpridpci <- -2*A*R + A*(2*s*U+1)*R*R*prdiff*pcdiff diag(d2L.dpridpcj) <- d2L.dpridpci d2L.dpcidpcj <- -sum(n[,1])/(pc[1]^2) + A*(2*s*U+1)*R*R*matrix(prdiff, ncol=1) %*% matrix(prdiff, nrow=1) # m-1 x m-1 d2L.dpci2 <- -cols.tot[-1]/(pc[-1]^2) - cols.tot[1]/(pc[1]^2) + A*(2*s*U + 1)*R*R*(prdiff^2) diag(d2L.dpcidpcj) <- d2L.dpci2 # The matrix of second derivatives is composed of vectors and matrices # computed above. The variance of alpha-s estimator is the top-left component of its inverse. # The matrix is symetric and constructed as follows: # # alpha Pr1 Pr2 ... Prq Pc1 Pc2 ... Pcq # # ------------ ----------------------- ------------------- # alpha - d2L.da2 - - t(d2L.dadpri) - - t(d2L.dadpcj) - # ------------ ----------------------- ------------------- # # Pr1 ---------------------- ------------------ # Pr2 - d2L.dpridprj - - d2L.dpridpcj - # ... - - - - # Prq ---------------------- ------------------ # # Pc1 ------------------- # Pc2 - d2L.dpcidpcj - # ... - - # Pcq ------------------- matrix.second.derivatives.top <- matrix(c(d2L.da2, d2L.dadpri, d2L.dadpcj), nrow=1) matrix.second.derivatives.middle <- cbind(matrix(d2L.dadpri, ncol=1), d2L.dpridprj, d2L.dpridpcj) matrix.second.derivatives.bottom <- cbind(matrix(d2L.dadpcj, ncol=1), t(d2L.dpridpcj), d2L.dpcidpcj) matrix.second.derivatives <- rbind(matrix.second.derivatives.top, matrix.second.derivatives.middle, matrix.second.derivatives.bottom) parms.cov.matrix <- solve(-matrix.second.derivatives) alpha.sd <- sqrt(parms.cov.matrix[1]) z <- qnorm((1+level)/2) lcl <- alpha - z*alpha.sd ucl <- alpha + z*alpha.sd list(alpha=alpha, lcl=lcl, ucl=ucl) } Top Download
aickin.alpha is a free R function. Download version 1.0 now. Top SAS version
aickin.alpha is also available as a SAS macro (see %AickinAlpha page). |