Dear all: This is a newbie''s question. Will anyone who knows the package to do the program: " Mantel-Haenszel extension of the chi-square test for trend (Schlesselman 1982, pp. 203-206) " please tell the location? I appreciate in advance. --------========---------- Mitsuo Igarashi mitsu5 at ruby.famille.ne.jp
Can''t help you there but please see Mark Myatt''s Home page http://www.myatt.demon.co.uk/ for a lot of useful information and utilities for epidemiological work. ----- Original Message ----- From: "Mitsuo Igarashi" <mitsu5 at ruby.famille.ne.jp> To: <r-help at stat.math.ethz.ch> Sent: Friday, March 28, 2003 10:04 AM Subject: [R] the chi-square test for trend> Dear all: > > This is a newbie''s question. > > Will anyone who knows the package to do the program: > " Mantel-Haenszel extension of the chi-square test for trend > (Schlesselman 1982, pp. 203-206) " > > please tell the location? > > I appreciate in advance. > > > --------========---------- > Mitsuo Igarashi > mitsu5 at ruby.famille.ne.jp > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://www.stat.math.ethz.ch/mailman/listinfo/r-help >
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
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