# Mantel-Haenszel test for linear trend # table: IxJ table of counts for two ordinal variables # output: Pearson correlation (pcor), test statistic (M2), and p-value # for H0: no linear trend # MH.test uses user-input scores # MH.test.mid uses midrank scores MH.test = function(table, rscore, cscore) { # rscore: vector of I scores for row variable # cscore: vector of J scores for column variable dim=dim(table) rbar=sum(margin.table(table,1)*rscore)/sum(table) rdif=rscore-rbar cbar=sum(margin.table(table,2)*cscore)/sum(table) cdif=cscore-cbar ssr=sum(margin.table(table,1)*(rdif^2)) ssc=sum(margin.table(table,2)*(cdif^2)) ssrc=sum(t(table*rdif)*cdif) pcor=ssrc/(sqrt(ssr*ssc)) M2=(sum(table)-1)*pcor^2 pval=1-pchisq(M2,1) return(list(pcor=pcor,M2=M2,pval=pval,rscore=rscore,cscore=cscore)) } MH.test.mid = function(table) { table=as.table(table) rscore=array(data=NA, dim=dim(table)[1]) cscore=array(data=NA, dim=dim(table)[2]) for( i in 1:dim(table)[1]){ if (i==1) rscore[i]=(margin.table(table,1)[i]+1)/2 else rscore[i]=sum(sum(margin.table(table,1)[1:(i-1)])+(margin.table(table,1)[i]+1)/2) } rscore ri=rscore-sum(table)/2 ri=as.vector(ri) for (j in 1:dim(table)[2]){ if (j==1) cscore[j]=(margin.table(table,2)[j]+1)/2 else cscore[j]=sum(sum(margin.table(table,2)[1:(j-1)])+(margin.table(table,2)[j]+1)/2) } cscore ci=cscore-sum(table)/2 ci=as.vector(ci) v=sum(t(table*ri)*ci) rowd=sum(table)^3-sum(margin.table(table,1)^3) cold=sum(table)^3-sum(margin.table(table,2)^3) w=(1/12)*sqrt(rowd*cold) pcor=v/w M2=(sum(table)-1)*pcor^2 pval=1-pchisq(M2,1) return(list(pcor=pcor,M2=M2,pval=pval,rscore=rscore,cscore=cscore)) }