#### wSSU-P: #### weighted SSU test with p-value calculated by permutations; #### the weights of the SNPs are inversely proportional to their var #### (stablized by "adding-2" rule and) #### calculated from the control samples only. #### NOTE: wSSU-P is NOT SSUw or SSUw-P; the latter two use the weights #### as the inverse of the variance of the score function from #### both control and case samples. #### ### Reference: Basu and Pan (2011). Comparison of Statistical Tests for Disease ### Association with Rare Variants. Genetic Epi. 35: 606-619. ### Downloadable from http://www.sph.umn.edu/biostatistics/research.asp ###### Copyright: Wei Pan, Division of Biostatistics, University of Minnesota ###### weip@biostat.umn.edu 3/15/11 ############################################################# # Purpose: as Madsen and Browning's weighted sum test, this test assumes that # the likelihood of a variant's being causal is inversely # proportional to its MAF, hence it will be more powerful than # SSU and SSUw in the presence of non-functional variants with # higher MAFs (e.g. non-functional CVs). For more details see the # first paragraph in the right column on page 608 (Basu and Pan 2011). # # Input: # X: nobs by nSNP genotype matrix; each element=0, 1 or 2; # Y: a binary phenotype of length nobs; each element=0 or 1; # B: #permutations used to calculate p-value. # Output: # P-value ############################################################# wSSUP<-function(Y, X, B=500){ n<-length(Y) Xg <- X #calculate weights based on ONLY controls: Xg0<-Xg[Y==0,] freq<-apply(Xg0, 2, sum) p0s<-(freq+1)/(2*nrow(Xg0)+2) wts<-1/sqrt(n*p0s*(1-p0s)) for(i in 1:n) Xg[i,]<-Xg[i,]*wts #########SumSqUs: U<-t(Xg) %*% (Y-mean(Y)) T1<- sum(U^2) #############Permuttaions: Ts<-rep(0, B) for(b in 1:B){ Yb<-sample(Y, n) Xg <- X #calculate weights based on ONLY controls: Xg0<-Xg[Yb==0,] freq<-apply(Xg0, 2, sum) p0s<-(freq+1)/(2*nrow(Xg0)+2) wts<-1/sqrt(n*p0s*(1-p0s)) for(i in 1:n) Xg[i,]<-Xg[i,]*wts #########SumSqUs: Ub<-t(Xg) %*% (Yb-mean(Yb)) Ts[b]<- sum(Ub^2) } sum(Ts>T1)/B }