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
[EMAIL PROTECTED]
Eric Lecoutre <[EMAIL PROTECTED]> 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
> [EMAIL PROTECTED]
> 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
______________________________________________
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
