[R] the chi-square test for trend

Mitsuo Igarashi mitsu5 at ruby.famille.ne.jp
Sat Mar 29 01:42:19 CET 2003


Hello Eric. Thank you so much.

I have just started to learn R. I am a bit bewildered to look 
at these programs. 
For example,I can not find "scores() and scores.type" in 
help("scores") and so on. Then I can not run these functions.

Will you or anybody please to teach me how to run the functions?
--------========----------
Mitsuo Igarashi
mitsu5 at ruby.famille.ne.jp


Eric Lecoutre <lecoutre at stat.ucl.ac.be> wrote:

> 
> 
> 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弾st pas quand les ordinateurs
> penseront comme des hommes, c弾st quand les hommes
> penseront comme des ordinateurs.     Sydney Harris



More information about the R-help mailing list