Gamma2.f<-function(x) {

#
# Descrição e detalhes:
# Função retorna o número de pares concordantes (C), discordantes (D), a estatística gamma, o seu respectivo erro
# padrão assintótico (obtido pelo método Delta) e o nível descritivo do teste H0: gamma=0.
# A estatística deve ser utilizada em tabelas de contingência IJ quando ambos fatores são ordinais e espera-se um
# padrão monótono (crescente ou decrescente) entre os fatores. Nesse caso, a estatística é particularmente interessante
# pois funciona como um teste de independência bem mais poderoso que o qui-quadrado, pois a estatística qui-quadrado
# testaria se as (I-1)x(J-1) interações são iguais a zero, enquanto que o gamma testa apenas 1 parâmetro.
# Ver Agresti (2002, págs.58-59). Sobre a implementação computacional, veja Thompson (2003, pág.39).
#
# Argumentos obrigatórios:
# x: matriz de contagens, sugere-se utilizar o crosstabs no S-Plus ou o xtabs no R.
#
# Autor: Frederico Zanqueta Poleto <fred@poleto.com>, arquivo disponível em http://www.poleto.com
#
# Referências:
# AGRESTI, A. (2002). Categorical Data Analysis. 2ª ed. John Wiley & Sons, New Jersey.
# THOMPSON, L. A. (2003). S-PLUS (and R) Manual to Accompany Agresti’s Categorical Data Analysis (2002). 2ª ed.
#    [Não publicado, disponível em http://math.cl.uh.edu/~thompsonla/Splusdiscrete2.pdf]
#
# Exemplos:
# Gamma2.f( crosstabs(count~idade+pressao,data=df1) ) #S-Plus
# Gamma2.f( xtabs(count~idade+pressao,data=df1) ) #R
#

# Check for using S-PLUS and output is from crosstabs
if(is.null(version$language) && inherits(x, "crosstabs")) { oldClass(x)<-NULL;
attr(x, "marginals")<-NULL}

n <- nrow(x)
m <- ncol(x)
pi.c<-pi.d<-matrix(0,nr=n,nc=m)

row.x<-row(x)
col.x<-col(x)

for(i in 1:(n)){
	for(j in 1:(m)){
		pi.c[i, j]<-sum(x[row.x<i & col.x<j]) + sum(x[row.x>i & col.x>j])
		pi.d[i, j]<-sum(x[row.x<i & col.x>j]) + sum(x[row.x>i & col.x<j])
	}
}

C <- sum(pi.c*x)/2
D <- sum(pi.d*x)/2

psi<-2*(D*pi.c-C*pi.d)/(C+D)^2
sigma<-sqrt(sum(x*psi^2)-sum(x*psi)^2) # Agresti, p. 76

gamma <- (C - D)/(C + D)
pv<-2*pnorm(-abs(gamma/sigma))

list(gamma = gamma, C = C, D = D, sigma=sigma, pvalue= pv)

}