############################################################ # C_alpha test: using both asymptotic Normal and permutation # distributions to calculate p-values. # Wei Pan, weip@biostat.umn.edu # Oct 24, 2010 ############################################################# ### Reference: Basu and Pan (2011). Comparison of Statistical Tests for Disease Association with Rare Variants. Genetic Epi. Downloadable from http://www.sph.umn.edu/biostatistics/research.asp ###### Copyright: Wei Pan, Division of Biostatistics, University of Minnesota ###### weip@biostat.umn.edu 10/24/10 ############################################################# # Input: # X: nobs by nSNP genotype matrix; each element=0, 1 or 2; # Y: binary phenotype of length nobs; each element=0 or 1; # B: #permutations used to calculate p-value. # Output: # a vector of (TestStat, AsymptoticPvalue, PermutationPvalue) ############################################################# sum1<-function(x) {sum(x>0)} CalphaTi<-function(yi, ni, p0=1/2){ (yi-ni*p0)^2 - ni*p0*(1-p0) } ClaphaVi<-function(ni, p0=1/2){ Vi=0 for(u in 0:ni) Vi=Vi + (((u-ni*p0)^2 - ni*p0*(1-p0))^2)*dbinom(u, ni, p0) } Calpha<-function(Y, X, B=500, p0=0.5){ nSNP=ncol(X) ni<-apply(X, 2, sum1) yi<-apply(X[Y>0,], 2, sum1) CalphaT<-V<-0 for(i in 1:length(ni)){ CalphaT = CalphaT + (yi[i]-ni[i]*p0)^2 - ni[i]*p0*(1-p0) for(u in 0:ni[i]) V=V + (((u-ni[i]*p0)^2 - ni[i]*p0*(1-p0))^2)*dbinom(u, ni[i], p0) } ####permutation: Ts<-rep(0, B) for(b in 1:B){ Yb<-sample(Y, length(Y)) yi<-apply(X[Yb>0,], 2, sum1) for(i in 1:length(ni)){ Ts[b] = Ts[b] + (yi[i]-ni[i]*p0)^2 - ni[i]*p0*(1-p0) } } if (V==0) Apv=1 else Apv=1-pchisq(CalphaT^2/V, df=1) c(CalphaT, V, Apv, sum(Ts^2>CalphaT^2)/B ) }