.packageName <- "irr"
"finn" <-
function(ratings, s.levels, model = c("oneway", "twoway")) {
  model <- match.arg(model)
	ratings <- as.matrix(na.omit(ratings))

	ns <- nrow(ratings)
	nr <- ncol(ratings)

  SStotal <- var(as.numeric(ratings))*(ns*nr-1)
  MSr <- var(apply(ratings,1,mean))*nr
	MSw <- sum(apply(ratings,1,var)/ns)
	MSc <- var(apply(ratings,2,mean))*ns
	MSe <- (SStotal-MSr*(ns-1)-MSc*(nr-1))/((ns-1)*(nr-1))

	MSexp <- 1/12*(s.levels^2-1)

	method <- paste("Finn-Coefficient (Model=",model,")",sep="")

	if (model=="oneway") {
		#Asendorpf & Wallbott, S. 245, Fu
		#Finn (1970)
		coeff  <- 1-(MSw/MSexp)
		Fvalue <- MSexp/MSw
	  df1    <- Inf
	  df2    <- ns*(nr-1)
	  p.value <- pf(Fvalue, df1, df2, lower.tail=FALSE)
	}
	else {
		#Asendorpf & Wallbott, S. 246, Fa
		coeff  <- 1-(MSe/MSexp)
		Fvalue <- MSexp/MSe
	  df1    <- Inf
	  df2    <- ns*(nr-1)
	  p.value <- pf(Fvalue, df1, df2, lower.tail=FALSE)
	}

  rval <- structure(list(method = method,
                         subjects = ns, raters = nr,
                         irr.name = "Finn", value = coeff,
                         stat.name = paste("F(Inf,",df2,")",sep=""), statistic = Fvalue, p.value = p.value),
                    class="irrlist")

  return(rval)
}

"icc" <-
function(ratings, model = c("oneway", "twoway"), type = c("consistency", "agreement"), k = 1, r0 = 0, conf.level = .95) {
	ratings <- as.matrix(na.omit(ratings))
  model <- match.arg(model)
  type  <- match.arg(type)
  alpha <- 1-conf.level

  ns <- nrow(ratings)
	nr <- ncol(ratings)

  SStotal <- var(as.numeric(ratings))*(ns*nr-1)
  MSr <- var(apply(ratings,1,mean))*nr
	MSw <- sum(apply(ratings,1,var)/ns)
	MSc <- var(apply(ratings,2,mean))*ns
	MSe <- (SStotal-MSr*(ns-1)-MSc*(nr-1))/((ns-1)*(nr-1))

	#Single Score ICCs
	if (k == 1) {
	    if (model=="oneway") {
			#Asendorpf & Wallbott, S. 245, ICu
			#Bartko (1966), [3]

		    icc.name <- "ICC(1)"
		    coeff  <- (MSr-MSw)/(MSr+(nr-1)*MSw)

		    Fvalue <- MSr/MSw*((1-r0)/(1+(k-1)*r0))
		    df1    <- ns-1
		    df2    <- ns*(nr-1)
		    p.value <- pf(Fvalue, df1, df2, lower.tail=FALSE)*2

		    #confidence interval
		    FL <- (MSr/MSw)/qf(1-alpha/2, ns-1, ns*(nr-1))
		    FU <- (MSr/MSw)*qf(1-alpha/2, ns*(nr-1), ns-1)
		    lbound <- (FL-1)/(FL+(nr-1))
		    ubound <- (FU-1)/(FU+(nr-1))
	    }
	    else if (model=="twoway") {
			if (type == "consistency") {
		    #Asendorpf & Wallbott, S. 245, ICa
				#Bartko (1966), [21]
				#Shrout & Fleiss (1979), ICC(3,1)

		    icc.name <- "ICC(C,1)"
				coeff  <- (MSr-MSe)/(MSr+(nr-1)*MSe)

			  Fvalue <- MSr/MSe*((1-r0)/(1+(k-1)*r0))
			  df1    <- ns-1
			  df2    <- (ns-1)*(nr-1)
			  p.value <- pf(Fvalue, df1, df2, lower.tail=FALSE)*2

			  #confidence interval
			  FL <- (MSr/MSe)/qf(1-alpha/2, ns-1, (ns-1)*(nr-1))
			  FU <- (MSr/MSe)*qf(1-alpha/2, (ns-1)*(nr-1), ns-1)
			  lbound <- (FL-1)/(FL+(nr-1))
			  ubound <- (FU-1)/(FU+(nr-1))
			}
			else if (type == "agreement") {
				#Asendorpf & Wallbott, S. 246, ICa'
				#Bartko (1966), [15]
				#Shrout & Fleiss (1979), ICC(2,1)

				icc.name <- "ICC(A,1)"
				coeff  <- (MSr-MSe)/(MSr+(nr-1)*MSe+(nr/ns)*(MSc-MSe))

				a <- (nr*r0)/(ns*(1-r0))
				b <- 1+(nr*r0*(ns-1))/(ns*(1-r0))
				#v <- (a*MSc+b*MSe)^2/((a*MSc)^2/(nr-1)+(b*MSe)^2/((ns-1)*(nr-1)))

				Fvalue <- MSr/(a*MSc+b*MSe)

				a <- (nr*coeff)/(ns*(1-coeff))
				b <- 1+(nr*coeff*(ns-1))/(ns*(1-coeff))
				v <- (a*MSc+b*MSe)^2/((a*MSc)^2/(nr-1)+(b*MSe)^2/((ns-1)*(nr-1)))

			  df1     <- ns-1
			  df2     <- v
			  p.value <- pf(Fvalue, df1, df2, lower.tail=FALSE)*2

			  #confidence interval (McGraw & Wong, 1996)
				FL <- qf(1-alpha/2, ns-1, v)
			  FU <- qf(1-alpha/2, v, ns-1)
			  lbound <- (ns*(MSr-FL*MSe))/(FL*(nr*MSc+(nr*ns-nr-ns)*MSe)+ns*MSr)
			  ubound <- (ns*(FU*MSr-MSe))/(nr*MSc+(nr*ns-nr-ns)*MSe+ns*FU*MSr)
			}
	  }
	}
	#Average Score ICCs
	else {
	    if (model=="oneway") {
			  #Asendorpf & Wallbott, S. 245, Ru
		    icc.name <- paste("ICC(",k,")",sep="")
		    coeff  <- (MSr-MSw)/MSr

		    Fvalue <- MSr/MSw*(1-r0)
		    df1    <- ns-1
		    df2    <- ns*(nr-1)
		    p.value <- pf(Fvalue, df1, df2, lower.tail=FALSE)*2

		    #confidence interval
		    FL <- (MSr/MSw)/qf(1-alpha/2, ns-1, ns*(nr-1))
		    FU <- (MSr/MSw)*qf(1-alpha/2, ns*(nr-1), ns-1)
		    lbound <- 1-1/FL
		    ubound <- 1-1/FU
	    }
	    else if (model=="twoway") {
			if (type == "consistency") {
				#Asendorpf & Wallbott, S. 246, Ra
		    icc.name <- paste("ICC(C,",k,")",sep="")
				coeff  <- (MSr-MSe)/MSr

			  Fvalue <- MSr/MSe*(1-r0)
			  df1    <- ns-1
			  df2    <- (ns-1)*(nr-1)
			  p.value <- pf(Fvalue, df1, df2, lower.tail=FALSE)*2

			  #confidence interval
			  FL <- (MSr/MSe)/qf(1-alpha/2, ns-1, (ns-1)*(nr-1))
			  FU <- (MSr/MSe)*qf(1-alpha/2, (ns-1)*(nr-1), ns-1)
			  lbound <- 1-1/FL
			  ubound <- 1-1/FU
			}
			else if (type == "agreement") {
		    icc.name <- paste("ICC(A,",k,")",sep="")
				coeff  <- (MSr-MSe)/(MSr+(MSc-MSe)/ns)

				a <- r0/(ns*(1-r0))
				b <- 1+(r0*(ns-1))/(ns*(1-r0))
				#v <- (a*MSc+b*MSe)^2/((a*MSc)^2/(nr-1)+(b*MSe)^2/((ns-1)*(nr-1)))

				Fvalue <- MSr/(a*MSc+b*MSe)

				a <- (nr*coeff)/(ns*(1-coeff))
				b <- 1+(nr*coeff*(ns-1))/(ns*(1-coeff))
				v <- (a*MSc+b*MSe)^2/((a*MSc)^2/(nr-1)+(b*MSe)^2/((ns-1)*(nr-1)))

			  df1    <- ns-1
			  df2    <- v
			  p.value <- pf(Fvalue, df1, df2, lower.tail=FALSE)*2

			  #confidence interval (McGraw & Wong, 1996)
				FL <- qf(1-alpha/2, ns-1, v)
			  FU <- qf(1-alpha/2, v, ns-1)
			  lbound <- (ns*(MSr-FL*MSe))/(FL*(MSc-MSe)+ns*MSr)
			  ubound <- (ns*(FU*MSr-MSe))/(MSc-MSe+ns*FU*MSr)
			}
	  }
	}

  rval <- structure(list(subjects = ns, raters = nr,
                         model = model, type = type, k = k,
                         icc.name = icc.name, value = coeff,
                         r0 = r0, Fvalue = Fvalue, df1 = df1, df2 = df2, p.value = p.value,
                         conf.level = conf.level, lbound = lbound, ubound = ubound),
                    class="icclist")
  return(rval)
}

"kappa2" <-
function(ratings, weight = c("unweighted", "equal", "squared")) {
	ratings <- as.matrix(na.omit(ratings))
	if (is.character(weight))
        weight = match.arg(weight)

	ns <- nrow(ratings)
	nr <- ncol(ratings)

	r1 <- ratings[,1]; r2 <- ratings[,2]

	if (!is.factor(r1)) r1 <- factor(r1)
	if (!is.factor(r2)) r2 <- factor(r2)

	#Find factor levels
	if (length(levels(r1)) >= length(levels(r2))) lev <- c(levels(r1), levels(r2))
	else lev <- c(levels(r2), levels(r1))

	lev <- lev[!duplicated(lev)]
	levels(r1) <- lev; levels(r2) <- lev

  #Compute table
	ttab <- table(r1, r2)

	#Compute weights
	weighttab <- as.matrix(ttab)
	nc <- ncol(weighttab)

	if (is.numeric(weight))
		w <- 1-(weight-min(weight))/(max(weight)-min(weight))
  else if (weight == "equal")
    w <- (nc-1):0/(nc-1)
  else if (weight == "squared")
  	w <- 1 - (0:(nc-1))^2/(nc - 1)^2
  else #unweighted
    w <- c(1, rep(0,nc-1))

	wvec <- c(sort(w, decreasing=FALSE), w[2:length(w)])
	nw <- length(w)
	weighttab <- matrix(0, nrow=nw, ncol=nw)
	for (i in 1:nw) {
		weighttab[i,] <- wvec[(nw-(i-1)):(2*nw-i)]
	}

	agreeP <- sum(ttab*weighttab)/ns

	tm1 <- apply(ttab, 1, sum)
	tm2 <- apply(ttab, 2, sum)

	eij <- outer(tm1, tm2)/ns
	chanceP <- sum(eij*weighttab)/ns

	#Kappa for 2 raters
	value <- (agreeP - chanceP)/(1 - chanceP)

	#Compute statistics
	pe <- sum(tm1*tm2)/ns^2

	p.i <- tm1/ns; p.j <- tm2/ns

	w.i <- apply(rep(tm2/ns,nc)*weighttab,2,sum)
	w.j <- apply(rep(tm1/ns,each=nc)*weighttab,1,sum)

	p.i <- rep(p.i, nc); p.j <- rep(p.j, each=nc)
	w.i <- rep(w.i, nc); w.j <- rep(w.j, each=nc)

	var.matrix <- p.i*p.j*(weighttab-(w.i+w.j))^2

	varkappa <- 1/(ns*(1-pe)^2)*(sum(var.matrix)-pe^2)

	SEkappa <- sqrt(varkappa)
	u <- value/SEkappa
	p.value <- 2 * (1 - pnorm(abs(u)))

  rval <- structure(list(method = paste("Cohen's Kappa for 2 Raters (Weights: ",paste(weight,collapse=","),")",sep=""),
                         subjects = ns, raters = nr,
                         irr.name = "Kappa", value = value,
                         stat.name = "z", statistic = u, p.value = p.value),
                    class="irrlist")
  return(rval)
}

"kappam.fleiss" <-
function(ratings, exact = FALSE, detail = FALSE) {
	ratings <- as.matrix(na.omit(ratings))

	ns <- nrow(ratings)
	nr <- ncol(ratings)

	#Build table
	lev <- levels(as.factor(ratings))

	for (i in 1:ns) {
		frow <- factor(ratings[i,],levels=lev)

		if (i==1)
			ttab <- as.numeric(table(frow))
		else
			ttab <- rbind(ttab, as.numeric(table(frow)))
	}

	ttab <- matrix(ttab, nrow=ns)

	agreeP <- sum((apply(ttab^2,1,sum)-nr)/(nr*(nr-1))/ns)

  if (!exact) {
    method  <- "Fleiss' Kappa for m Raters"
  	chanceP <- sum(apply(ttab,2,sum)^2)/(ns*nr)^2
	} else {
    method  <- "Fleiss' Kappa for m Raters (exact value)"
  	for (i in 1:nr) {
  		rcol <- factor(ratings[,i],levels=lev)

  		if (i==1)
  			rtab <- as.numeric(table(rcol))
  		else
  			rtab <- rbind(rtab, as.numeric(table(rcol)))
  	}

  	rtab <- rtab/ns

  	chanceP <- sum(apply(ttab,2,sum)^2)/(ns*nr)^2 - sum(apply(rtab,2,var)*(nr-1)/nr)/( nr-1)
  }

 	#Kappa for m raters
 	value <- (agreeP - chanceP)/(1 - chanceP)

  if (!exact) {
  	pj2 <- chanceP
  	pj3 <- sum(apply(ttab,2,sum)^3)/(ns*nr)^3

  	varkappa <- (2/(ns*nr*(nr-1)))*((pj2-(2*nr-3)*pj2^2+2*(nr-2)*pj3)/(1-pj2)^2)
  	SEkappa <- sqrt(varkappa)

  	u <- value/SEkappa
  	p.value <- 2 * (1 - pnorm(abs(u)))

    if (detail) {
    	pj  <- apply(ttab,2,sum)/(ns*nr)
    	pjk <- (apply(ttab^2,2,sum)-ns*nr*pj)/(ns*nr*(nr-1)*pj)

    	kappaK    <- (pjk-pj)/(1-pj)
    	varkappaK <- ((1+2*(nr-1)*pj)^2+2*(nr-1)*pj*(1-pj))/(ns*nr*(nr-1)^2*pj*(1-pj))
    	SEkappaK  <- sqrt(varkappaK)

    	uK <- kappaK/SEkappaK
    	p.valueK <- 2 * (1 - pnorm(abs(uK)))

    	tableK <- as.table(round(cbind(kappaK, uK, p.valueK), digits=3))

    	rownames(tableK) <- lev
    	colnames(tableK) <- c("Kappa", "z", "p.value")
    }
  }

  if (!exact) {
    if (!detail) {
      rval <- list(method = method,
                   subjects = ns, raters = nr,
                   irr.name = "Kappa", value = value)
    }
    else {
      rval <- list(method = method,
                   subjects = ns, raters = nr,
                   irr.name = "Kappa", value = value,
                   detail = tableK)
    }
    rval <- c(rval, stat.name = "z", statistic = u, p.value = p.value)
  }
  else {
    rval <- list(method = method,
                 subjects = ns, raters = nr,
                 irr.name = "Kappa", value = value)
  }
  class(rval) <- "irrlist"

  return(rval)
}

"kappam.light" <-
function(ratings) {
	ratings <- as.matrix(na.omit(ratings))

	ns <- nrow(ratings)
	nr <- ncol(ratings)

  for (i in 1:(nr-1))
    for (j in (i+1):nr) {
      if ((i==1) & (j==(i+1))) kappas <- kappa2(ratings[,c(i,j)], weight="u")$value
      else kappas <- c(kappas, kappa2(ratings[,c(i,j)], weight="u")$value)
    }

  value <- mean(kappas)

  #Variance & Computation of p-value
  lev    <- levels(as.factor(ratings))
  levlen <- length(levels(as.factor(ratings)))

  for (nri in 1:(nr-1))
    for (nrj in (nri+1):nr) {
      for (i in 1:levlen)
        for (j in 1:levlen) {
          if (i!=j) {
            r1i <- sum(ratings[,nri]==lev[i])
            r2j <- sum(ratings[,nrj]==lev[j])
            if (!exists("dis")) dis <- r1i*r2j
            else dis <- c(dis,r1i*r2j)
          }
        }
        if (!exists("disrater")) disrater <- sum(dis)
        else disrater <- c(disrater,sum(dis))
        rm(dis)
      }

  B <- length(disrater) * prod(disrater)

  chanceP  <- 1-B/ns^(choose(nr,2)*2)
  varkappa <- chanceP/(ns*(1-chanceP))

	SEkappa <- sqrt(varkappa)
	u <- value/SEkappa
	p.value <- 2 * (1 - pnorm(abs(u)))

  rval <- structure(list(method = "Light's Kappa for m Raters",
                         subjects = ns, raters = nr,
                         irr.name = "Kappa", value = value,
                         stat.name = "z", statistic = u, p.value = p.value),
                    class="irrlist")
  return(rval)
}

"kendall" <-
function(ratings, correct=FALSE) {
	ratings <- as.matrix(na.omit(ratings))

	ns <- nrow(ratings)
	nr <- ncol(ratings)

	#Without correction for ties
	if (!correct) {
		#Test for ties
		TIES = FALSE
		testties <- apply(ratings, 2, unique)
		if (!is.matrix(testties)) TIES=TRUE
		else { if (length(testties) < length(ratings)) TIES=TRUE }

		ratings.rank <- apply(ratings,2,rank)

		coeff.name <- "W"
		coeff <- (12*var(apply(ratings.rank,1,sum))*(ns-1))/(nr^2*(ns^3-ns))
	}
	else { #With correction for ties
		ratings <- as.matrix(na.omit(ratings))

		ns <- nrow(ratings)
		nr <- ncol(ratings)

		ratings.rank <- apply(ratings,2,rank)

		Tj <- 0
		for (i in 1:nr) {
			rater <- table(ratings.rank[,i])
			ties  <- rater[rater>1]
			l 	  <- as.numeric(ties)
			Tj	  <- Tj + sum(l^3-l)
		}

		coeff.name <- "Wt"
		coeff <- (12*var(apply(ratings.rank,1,sum))*(ns-1))/(nr^2*(ns^3-ns)-nr*Tj)
	}

	#test statistics
	Xvalue  <- nr*(ns-1)*coeff
	df1     <- ns-1
	p.value <- pchisq(Xvalue, df1, lower.tail = FALSE)

  rval <- list(method = paste("Kendall's coefficient of concordance",coeff.name),
               subjects = ns, raters = nr,
               irr.name = coeff.name, value = coeff,
               stat.name = paste("Chisq(",df1,")",sep=""), statistic = Xvalue, p.value = p.value)
  if (!correct && TIES) rval <- c(rval, error="Coefficient may be incorrect due to ties")
 	class(rval) <- "irrlist"
  return(rval)
}

"meancor" <-
function(ratings, fisher=TRUE) {
	ratings <- as.matrix(na.omit(ratings))

	ns <- nrow(ratings)
	nr <- ncol(ratings)

	for (i in 1:(nr-1)) for (j in (i+1):nr) {
    if ((i==1) & (j==(i+1))) r <- cor(ratings[,i],ratings[,j])
    else r <- c(r, cor(ratings[,i],ratings[,j]))
	}

  delr <- 0
  
  if (fisher) {
    delr <- length(r) - length(r[(r<1) & (r>-1)])
    #Eliminate perfect correlations (r=1, r=-1)
    r <- r[(r<1) & (r>-1)]
    
 	  rz  <- 1/2*log((1+r)/(1-r))
 	  mrz <- mean(rz)

 	  coeff <- (exp(2*mrz)-1)/(exp(2*mrz)+1)
  	SE    <- sqrt(1/(ns-3))

  	u       <- coeff/SE
  	p.value <- 2 * (1 - pnorm(abs(u)))
  }
  else {
    coeff <- mean(r)
  }

  rval <- list(method = "Mean of bivariate correlations R",
               subjects = ns, raters = nr,
               irr.name = "R", value = coeff)

  if (fisher) rval <- c(rval, stat.name = "z", statistic = u, p.value = p.value)
  if (delr>0) rval <- c(rval, error = paste(delr, ifelse(delr==1, "perfect correlation was", "perfect correlations were"), "dropped before averaging"))
  class(rval) <- "irrlist"

  return(rval)
}

"meanrho" <-
function(ratings, fisher=TRUE) {
	ratings <- as.matrix(na.omit(ratings))

	ns <- nrow(ratings)
	nr <- ncol(ratings)

	#Test for ties
	TIES = FALSE
	testties <- apply(ratings, 2, unique)
	if (!is.matrix(testties)) TIES <- TRUE
	else { if (length(testties) < length(ratings)) TIES <- TRUE }

	ratings.rank <- apply(ratings,2,rank)

	for (i in 1:(nr-1)) for (j in (i+1):nr) {
    if ((i==1) & (j==(i+1))) r <- cor(ratings[,i], ratings[,j], method="spearman")
    else r <- c(r, cor(ratings[,i], ratings[,j], method="spearman"))
	}

  delr <- 0

  if (fisher) {
    delr <- length(r) - length(r[(r<1) & (r>-1)])
    #Eliminate perfect correlations (r=1, r=-1)
    r <-  r[(r<1) & (r>-1)]

 	  rz  <- 1/2*log((1+r)/(1-r))
 	  mrz <- mean(rz)

 	  coeff <- (exp(2*mrz)-1)/(exp(2*mrz)+1)
  	SE    <- sqrt(1/(ns-3))

  	u       <- coeff/SE
  	p.value <- 2 * (1 - pnorm(abs(u)))
  }
  else {
    coeff <- mean(r)
  }

  rval <- list(method = "Mean of bivariate rank correlations Rho",
               subjects = ns, raters = nr,
               irr.name = "Rho", value = coeff)

  if (fisher) rval <- c(rval, stat.name = "z", statistic = u, p.value = p.value)
  if (delr>0) {
    if (TIES) {
      rval <- c(rval, error = paste(delr, ifelse(delr==1, "perfect correlation was", "perfect correlations were"), "dropped before averaging",
                                    "\n Coefficient may be incorrect due to ties"))
    }
    else {
      rval <- c(rval, error = paste(delr, ifelse(delr==1, "perfect correlation was", "perfect correlations were"), "dropped before averaging"))
    }
  }
  
  if ((delr==0) & (TIES)) rval <- c(rval, error = "Coefficient may be incorrect due to ties")
  class(rval) <- "irrlist"

  return(rval)
}

"print.icclist" <-
function(x, ...)
{
  icc.title <- ifelse(x$k==1, "Single Score Intraclass Correlation", "Average Score Intraclass Correlation")
	cat(paste(" ",icc.title,"\n\n",sep=""))
	cat(paste("   Model:", x$model, "\n"))
	cat(paste("   Type :", x$type, "\n\n"))
	cat(paste("   Subjects =", x$subjects, "\n"))
	cat(paste("     Raters =", x$raters, "\n"))
	results <- paste(format.char(x$icc.name, width=11, flag="+"), "=", format(x$value, digits=3))
	cat(results)
  cat("\n\n F-Test, H0: r0 =",x$r0,"\n")
	Ftest <- paste(format.char(paste("F(",x$df1,",",format(x$df2, digits=3),")",sep=""), width=11, flag="+"), "=", format(x$Fvalue, digits=3),
	               ", p =", format(x$p.value, digits=3), "\n\n")
  cat(Ftest)
	cat(" ", round(x$conf.level*100,digits=1), "%-Confidence Interval for ICC Population Values:\n", sep="")
	cat(paste("  ", round(x$lbound, digits=3), " < ICC < ", round(x$ubound, digits=3), "\n", sep=""))
}

"print.irrlist" <-
function(x, ...)
{
  cat(" ", x$method, "\n\n",sep="")
	cat(paste(" Subjects =", x$subjects, "\n"))
	cat(paste("   Raters =", x$raters, "\n"))
	results <- paste(format.char(x$irr.name, width=9, flag="+"), "=", format(x$value, digits=3),"\n")
	cat(results)
	if (!is.null(x$statistic)) {
  	statistic <- paste(format.char(x$stat.name, width=9, flag="+"), "=", format(x$statistic, digits=3), "\n")
  	cat("\n", statistic, sep="")
  	cat(paste("  p-value =", format(x$p.value, digits=3), "\n"))
	}
  if (!is.null(x$detail)) {
    cat("\n")
    print(x$detail)
  }
  if (!is.null(x$error)) cat("\n ", x$error, "\n", sep="")
}

"robinson" <-
function(ratings) {
	ratings <- as.matrix(na.omit(ratings))
	ns <- nrow(ratings)
	nr <- ncol(ratings)

	SStotal <- var(as.numeric(ratings))*(ns*nr-1)
	SSb <- var(apply(ratings,1,mean))*nr*(ns-1)
	SSw <- var(apply(ratings,2,mean))*ns*(nr-1)
	SSr <- SStotal-SSb-SSw

	coeff  <- SSb/(SSb+SSr)

  rval <- structure(list(method = "Robinson's A",
                         subjects = ns, raters = nr,
                         irr.name = "A", value = coeff),
                    class="irrlist")
  return(rval)
}

