######################Multivariate score, SSU, SSUw, UminP, Sum, ###################### adaptive Sum (aSum), aSum-P tests. ######################modified from the program of Fang HAN (Email: fanghan@biostat.umn.edu) on 1/22/10 ####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 library("mvtnorm") # Introduction to input parameters: # Y: disease lables; =0 for controls, =1 for cases; # X: genotype data; each row for a subject, and each column for an SNP; # the value of each element is the # of the copies for an allele; # alpha0 = 0.1: the cut-off value; if the marginal p-value of an SNP is less # than alpha0 and its marginal regression coefficient is negative, # we flip its coding; # B=1000: the number of permutations; ############################ # Output is the p-values of the seven tests in the order of: # 1. Score # 2. SSU # 3. SSUw # 4. UminP # 5. Sum # 6. aSum-P # 7. aSum ############################################ #######Some back-up programs ############################################ PowerUniv <- function(U,V){ n <- dim(V)[1] x <- as.numeric(max(abs(U))) TER <- as.numeric(1-pmvnorm(lower=c(rep(-x,n)),upper=c(rep(x,n)),mean=c(rep(0,n)),sigma=V)) return(TER) } ##################Univariate score test: ################## Added by wei Pan, 10/22/10: ###To speed up the marginal test, use the score test, rather than the Wald's! UscoTest<-function(Y,X){ U<-sum((X-mean(X))*(Y-mean(Y))) V<- mean(Y)*(1-mean(Y))*sum( (X-mean(X))^2) if (abs(U)<1e-20) Ts<-0 else Ts<-U^2/V if (is.na(Ts) || is.infinite(Ts) || is.nan(Ts)) Ts<-0 c(U, 1-pchisq(Ts, df=1)) } SumTest <- function(Y,X,alpha0){ pv<-NULL Uus<-NULL for(i in 1:ncol(X)){ #fit <- glm(Y~X[,i], family=binomial(logit)) ###"if ()" added by WP, 10/22/10, since some X[,j] can be ALL 0: #if (is.na(fit$coefficients[-1])){ # beta <- c(beta, 0); pv <- c(pv,1); # } else{ #beta <- c(beta,fit$coefficients[-1]) #pv <- c(pv,as.numeric(summary(fit)$coefficients[,4][-1])) #} ####Score test to speed up!!! scoRes<-UscoTest(Y,X[,i]) Uus<-c(Uus, scoRes[1]) pv <- c(pv, scoRes[2]) } Xg <- X Xbar<-mean(Xg) Xgb<-Xg Xgb<-Xgb-Xbar U<-t(Xg) %*% (Y-mean(Y)) CovS<- mean(Y)*(1-mean(Y))*(t(Xgb) %*% Xgb) a<-rep(1, length(U)) a[Uus<0 & pv1e-10, diagCovS, 1e-10) Tg2<- t(U) %*% diag(1/diagCovS) %*% U ##distr of Tg1 is sum of cr Chisq_1: cr<-eigen(CovS %*% diag(1/diagCovS), only.values=TRUE)$values ##approximate the distri by alpha Chisq_d + beta: alpha2<-sum(cr*cr*cr)/sum(cr*cr) beta2<-sum(cr) - (sum(cr*cr)^2)/(sum(cr*cr*cr)) d2<-(sum(cr*cr)^3)/(sum(cr*cr*cr)^2) alpha2<-as.real(alpha2) beta2<-as.real(beta2) d2<-as.real(d2) pTg2<-as.numeric(1-pchisq((Tg2-beta2)/alpha2, d2)) ##########score test: ##gInv of CovS: CovS.edecomp<-eigen(CovS) CovS.rank<-sum(abs(CovS.edecomp$values)> 1e-8) inveigen<-ifelse(abs(CovS.edecomp$values) >1e-8, 1/CovS.edecomp$values, 0) P<-solve(CovS.edecomp$vectors) gInv.CovS<-t(P) %*% diag(inveigen) %*% P Tscore<- t(U) %*% gInv.CovS %*% U pTscore<-as.numeric( 1-pchisq(Tscore, CovS.rank) ) #### V1 %*% t(V1)=CovS: # CovS.ev<-ifelse(abs(CovS.edecomp$values) >1e-8, CovS.edecomp$values, 0) # V1<-CovS.edecomp$vectors %*% diag(sqrt(CovS.ev)) ##univariate/marginal tests: Tus<-as.vector(abs(U)/sqrt(diagCovS)) ####added by WP: 10/24/10: nSNP=length(U) Vs <- matrix(c(rep(0,nSNP^2)),nrow=nSNP) for(i in 1:nSNP){ for(j in 1:nSNP){ if (abs(CovS[i,j])>1e-20) Vs[i,j] <- CovS[i,j]/sqrt(CovS[i,i]*CovS[j,j]) else Vs[i,j] <- 1e-20 } } pTus <- as.numeric(PowerUniv(Tus,Vs)) ##########SumTest######################## a<-rep(1, length(U)) Tsum<- sum(t(a)%*%U)/(sqrt(as.numeric(t(a) %*% CovS %*% (a)))) pTsum <- as.numeric( 1-pchisq(Tsum^2, 1) ) ##########aSumTest fit0 <- SumTest(Y,Xg,alpha0) u <- fit0[1,1] pv <- fit0[1,2] v <- fit0[1,3] fit <- rSumTest(Y,Xg,B,alpha0) u0 <- fit[1,1] v0 <- fit[1,2] a <- fit[1,3] b <- fit[1,4] pv0 <- fit[,5] aSumP <- sum(pv>pv0)/length(pv0) aSum <- as.numeric( 1-pchisq(abs(((u-u0)^2/v0-b)/a), 1) ) return(c(pTscore, pTg1, pTg2, pTus, pTsum, aSumP, aSum)) }