[R] the chi-square test for trend

[Eric Lecoutre]

Hi,

Please find here a function to compute Mantel-Haenszel Chi squared as well
as Cochran-Armitage test for trend (suitable when one of the dimensions is 2).
Both functions take as argument a contingency table.
By the way, the MH Chi² is very easy to compute as it is defined as n*rho.

Eric Lecoutre

------------------------------------------------------------------
tablepearson=function(x,scores.type="table")
{

# Statistic
	sR=scores(x,1,scores.type)
	sC=scores(x,2,scores.type)
	n=sum(data)
	Rbar=sum(apply(x,1,sum)*sR)/n
	Cbar=sum(apply(x,2,sum)*sC)/n
	ssr=sum(x*(sR-Rbar)^2)
	ssc=sum(t(x)* (sC-Cbar)^2)
	tmpij=outer(sR,sC,FUN=function(a,b) return((a-Rbar)*(b-Cbar)))
	ssrc= sum(x*tmpij)
	v=ssrc
	w=sqrt(ssr*ssc)
	r=v/w
# ASE
	bij=outer(sR,sC, FUN=function(a,b)return((a-Rbar)^2*ssc + (b-Cbar)^2*ssr))	
	tmp1=1/w^2
	tmp2=x*(w*tmpij - (bij*v)/(2*w))^2
	tmp3=sum(tmp2)
	ASE=tmp1*sqrt(tmp3)
# Test
	var0= (sum(x*tmpij) - (ssrc^2/n))/ (ssr*ssc)
	tb=r/sqrt(var0)
	p.value=2*(1-pnorm(tb))
# Output
	out=list(estimate=r,ASE=ASE,statistic=tb,p.value=p.value,name="Pearson 
Correlation",bornes=c(-1,1))
	class(out)="ordtest"
	return(out)
}

tableChisqMH=function(x)
{
	n=sum(x)
	G2=n*(tablepearson(x)^2)
	dll=1
	p.value=1-pchisq(G2,dll)
	out=list(estimate=G2,dll=dll,p.value=p.value,dim=dim(x),name="Mantel-Haenszel 
Chi-square")
	return(out)

}



tabletrend=function(x,transpose=FALSE)
{
	if (any(dim(x)==2))
	{
	if (transpose==TRUE) {
	x=t(x)
	}
	
	if (dim(x)[2]!=2){stop("Cochran-Armitage test for trend must be used with 
a (R,2) table. Use transpose argument",call.=FALSE) }
	
	nidot=apply(x,1,sum)
	n=sum(nidot)

	Ri=scores(x,1,"table")
	Rbar=sum(nidot*Ri)/n
	
	s2=sum(nidot*(Ri-Rbar)^2)
	pdot1=sum(x[,1])/n
	T=sum(x[,1]*(Ri-Rbar))/sqrt(pdot1*(1-pdot1)*s2)
	p.value.uni=1-pnorm(abs(T))
	p.value.bi=2*p.value.uni
	out=list(estimate=T,dim=dim(x),p.value.uni=p.value.uni,p.value.bi=p.value.bi,name="Cochran-Armitage 
test for trend")
	return(out)
	
	}
	else {stop("Cochran-Armitage test for trend must be used with a (2,C) or a 
(R,2) table",call.=FALSE) }
}

-----------------------------------------------------------------------

__________________________________________________

Eric Lecoutre           Informaticien/Statisticien
Institut de Statistique                        UCL

                               (+32) (0)10 47 30 50
                            lecoutre at stat.ucl.ac.be
     http://www.stat.ucl.ac.be/ISpersonnel/lecoutre
__________________________________________________
Le vrai danger, ce n’est pas quand les ordinateurs
penseront comme des hommes, c’est quand les hommes
penseront comme des ordinateurs.     Sydney Harris