Nothing
SKAT_META_Optimal_Get_Q<-function(Score, r.all){
n.r<-length(r.all)
Q.r<-rep(0,n.r)
for(i in 1:n.r){
r.corr<-r.all[i]
Q.r[i]<-(1-r.corr) * sum(Score^2) + r.corr * sum(Score)^2
}
Q.r = Q.r /2
re<-list(Q.r=Q.r)
return(re)
}
SKAT_META_Optimal_Get_Q_Res<-function(Score.res, r.all){
n.r<-length(r.all)
p<-dim(Score.res)[1]
Q.r<-matrix(rep(0,n.r*p), ncol=n.r)
for(i in 1:n.r){
r.corr<-r.all[i]
Q.r[,i]<-(1-r.corr) * rowSums(Score.res^2) + r.corr * rowSums(Score.res)^2
}
Q.r = Q.r /2
re<-list(Q.r=Q.r)
return(re)
}
SKAT_META_Optimal_Get_Pvalue<-function(Q.all, Phi, r.all, method, isFast=FALSE){
n.r<-length(r.all)
n.q<-dim(Q.all)[1]
p.m<-dim(Phi)[2]
lambda.all<-list()
for(i in 1:n.r){
r.corr<-r.all[i]
R.M<-diag(rep(1-r.corr,p.m)) + matrix(rep(r.corr,p.m*p.m),ncol=p.m)
L<-chol(R.M,pivot=TRUE)
Phi_rho<- L %*% (Phi %*% t(L))
lambda.all[[i]]<-Get_Lambda(Phi_rho, isFast=isFast)
}
# Get Mixture param
param.m<-SKAT_META_Optimal_Param(Phi,r.all)
Each_Info<-SKAT_Optimal_Each_Q(param.m, Q.all, r.all, lambda.all, method=method)
pmin.q<-Each_Info$pmin.q
pval<-rep(0,n.q)
# added
pmin<-Each_Info$pmin
if(method == "davies" || method=="optimal" || method=="optimal.adj" || method=="optimal.mod"){
for(i in 1:n.q){
pval[i]<-SKAT_Optimal_PValue_Davies(pmin.q[i,],param.m,r.all, pmin[i])
}
} else if(method =="liu" || method =="liu.mod" ){
for(i in 1:n.q){
pval[i]<-SKAT_Optimal_PValue_Liu(pmin.q[i,],param.m,r.all, pmin[i])
}
} else {
stop("Invalid Method:", method)
}
# Check the pval
# Since SKAT-O is between burden and SKAT, SKAT-O p-value should be <= min(p-values) * 2
# To correct conservatively, we use min(p-values) * 3
multi<-3
if(length(r.all) < 3){
multi<-2
}
for(i in 1:n.q){
pval.each<-Each_Info$pval[i,]
IDX<-which(pval.each > 0)
pval1<-min(pval.each) * multi
if(pval[i] <= 0 || length(IDX) < length(r.all)){
pval[i]<-pval1
}
# if pval==0, use nonzero min each.pval as p-value
if(pval[i] == 0){
if(length(IDX) > 0){
pval[i] = min(pval.each[IDX])
}
}
}
return(list(p.value=pval,p.val.each=Each_Info$pval))
}
#
# Function get parameters of optimal test
#
SKAT_META_Optimal_Param<-function(Phi,r.all){
p.m<-dim(Phi)[2]
r.n<-length(r.all)
# ZMZ
Z.item1.1<- Phi %*% rep(1,p.m)
ZZ<-Phi
ZMZ<- Z.item1.1 %*% t(Z.item1.1) / sum(ZZ)
# W3.2 Term : mixture chisq
W3.2.t<-ZZ - ZMZ
lambda<-Get_Lambda(W3.2.t)
# W3.3 Term : variance of remaining ...
W3.3.item<-sum(ZMZ *(ZZ-ZMZ)) * 4
# tau term
z_mean_2<- sum(ZZ)/p.m^2
tau1<- sum(ZZ %*% ZZ) / p.m^2 / z_mean_2
# Mixture Parameters
MuQ<-sum(lambda)
VarQ<-sum(lambda^2) *2 + W3.3.item
KerQ<-sum(lambda^4)/(sum(lambda^2))^2 * 12
Df<-12/KerQ
# W3.1 Term : tau1 * chisq_1
tau<-rep(0,r.n)
for(i in 1:r.n){
r.corr<-r.all[i]
term1<-p.m^2*r.corr * z_mean_2 + tau1 * (1-r.corr)
tau[i]<-term1
}
out<-list(MuQ=MuQ,VarQ=VarQ,KerQ=KerQ,lambda=lambda,VarRemain=W3.3.item,Df=Df,tau=tau,
z_mean_2=z_mean_2, p.m=p.m,
tau.1 = tau1,
tau.2= p.m*z_mean_2 )
#param2<<-out
return(out)
}
#######################################################33
# Linear and Logistic
SKAT_META_Optimal = function(Score, Phi, r.all, method="davies", Score.Resampling, isFast=FALSE){
# if r.all >=0.999 ,then r.all = 0.999
IDX<-which(r.all >= 0.999)
if(length(IDX) > 0){
r.all[IDX]<-0.999
}
p.m<-dim(Phi)[2]
n.r<-length(r.all)
###########################################
# Compute Q.r and Q.r.res
##########################################
out.Q<-SKAT_META_Optimal_Get_Q(Score, r.all)
Q.res=NULL
if(!is.null(Score.Resampling)){
Q.res<-SKAT_META_Optimal_Get_Q_Res(Score.Resampling, r.all)$Q.r
}
Q.all<-rbind(out.Q$Q.r, Q.res)
##################################################
# Compute P-values
#################################################
out<-SKAT_META_Optimal_Get_Pvalue(Q.all, Phi/2, r.all, method, isFast=isFast)
param<-list(p.val.each=NULL,q.val.each=NULL)
param$p.val.each<-out$p.val.each[1,]
param$q.val.each<-Q.all[1,]
param$rho<-r.all
param$minp<-min(param$p.val.each)
id_temp<-which(param$p.val.each == min(param$p.val.each))
id_temp1<-which(param$rho >= 0.999) # treat rho > 0.999 as 1
if(length(id_temp1) > 0){
param$rho[id_temp1] = 1
}
param$rho_est<-param$rho[id_temp]
p.value<-out$p.value[1]
p.value.resampling= NULL
if(!is.null(Q.res)){
p.value.resampling=out$p.value[-1]
}
re<-list(p.value = p.value, param=param, p.value.resampling=p.value.resampling)
return(re)
}
##################################################################
#
# Note: fastOption cannot be used for Optimal test
Met_SKAT_Get_Pvalue<-function(Score, Phi, r.corr, method, Score.Resampling=NULL, isFast=FALSE){
#Score.Resampling1<<-Score.Resampling
p.m<-nrow(Phi)
Q.res = NULL
# if Phi==0
if(sum(abs(Phi)) == 0){
warning("No polymorphic SNPs!",call.=FALSE)
return(list(p.value=1, p.value.resampling= NULL, pval.zero.msg=NULL))
}
if(!is.null(Score.Resampling)){
Score.Resampling<-t(Score.Resampling)
}
if(length(Phi) <=1){
r.corr=0
} else{
if(ncol(Phi) <=10){
if(qr(Phi)$rank <= 1){
r.corr=0
}
}
}
if(length(r.corr) > 1){
re = SKAT_META_Optimal(Score, Phi, r.corr, method=method, Score.Resampling=Score.Resampling)
return(re)
}
if (r.corr == 0){
Q<-sum(Score^2)/2
if(!is.null(Score.Resampling)){
Q.res<-rowSums(Score.Resampling^2)/2
}
} else if (r.corr==1){
Q=SKAT_META_Optimal_Get_Q(Score, r.corr)$Q.r
if(!is.null(Score.Resampling)){
Q.res<-SKAT_META_Optimal_Get_Q_Res(Score.Resampling, r.corr)$Q.r
}
a<- as.matrix(sum(Phi))
re<-Get_Liu_PVal(Q, a, Q.res)
return(re)
} else {
Q=SKAT_META_Optimal_Get_Q(Score, r.corr)$Q.r
if(!is.null(Score.Resampling)){
Q.res<-SKAT_META_Optimal_Get_Q_Res(Score.Resampling, r.corr)$Q.r
}
R.M<-diag(rep(1-r.corr,p.m)) + matrix(rep(r.corr,p.m*p.m),ncol=p.m)
L<-chol(R.M,pivot=TRUE)
Phi<- L %*% (Phi %*% t(L))
}
re<-Get_Davies_PVal(Q, Phi, Q.res, isFast=isFast)
if(length(r.corr)==1){
re$Q = Q
}
return(re)
}
#
# out_type
# C: continuous, D:binary, V: Kinship
#
#
SKAT_RunFrom_MetaSKAT<-function(res,Z, X1, kernel, weights=NULL, s2=NULL, pi_1=NULL, P0=NULL, out_type="C", method, res.out, n.Resampling, r.corr, isFast=FALSE){
if (kernel == "linear.weighted") {
Z = t(t(Z) * (weights))
}
# Get Score
Score = as.vector(t(Z) %*% res)
Score.Resampling=NULL
if(!is.null(Score.Resampling)){
Score.Resampling = t(res.out) %*% Z
}
# Phi
if(out_type=="C"){
Score=Score/ sqrt(s2)
if(!is.null(Score.Resampling)){
Score.Resampling = Score.Resampling / sqrt(s2)
}
Phi = t(Z) %*% Z - (t(Z) %*%X1)%*%solve(t(X1)%*%X1)%*% (t(X1) %*% Z )
} else if(out_type=="D"){
Phi = t(Z) %*% (Z * pi_1) - (t(Z * pi_1) %*%X1)%*%solve(t(X1)%*%(X1 * pi_1))%*% (t(X1) %*% (Z * pi_1))
} else if(out_type=="V"){
Phi = t(Z) %*% (P0 %*% Z) # t(Z) P0 Z
} else {
stop("SKAT_RunFrom_MetaSKAT: no-known out_type!")
}
re = Met_SKAT_Get_Pvalue(Score=Score, Phi=Phi, r.corr=r.corr, method=method, Score.Resampling=Score.Resampling, isFast=isFast)
re$IsMeta=TRUE
return(re)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.