Skip to content

Commit

Permalink
version 0.82
Browse files Browse the repository at this point in the history
  • Loading branch information
@lee7801 @gaborcsardi
lee7801 authored and gaborcsardi committed Feb 26, 2011
1 parent eb5f7a7 commit a34552b
Show file tree
Hide file tree
Showing 12 changed files with 383 additions and 259 deletions.
9 changes: 5 additions & 4 deletions DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -1,13 +1,14 @@
Package: SKAT
Type: Package
Title: SNP-set (Sequence) Kernel Association Test
Version: 0.81
Date: 2011-12-18
Version: 0.82
Date: 2011-02-26
Author: Seunggeun Lee, Larisa Miropolsky and Micheal Wu
Maintainer: Seunggeun (Shawn) Lee <phila78@gmail.com>
Description: Kernel based SNP set test
License: GPL (>= 2)
Depends: R (>= 2.13.0)
Packaged: 2013-01-15 03:46:01 UTC; sglee
Packaged: 2013-03-02 14:33:42 UTC; sglee
NeedsCompilation: yes
Repository: CRAN
Date/Publication: 2013-01-15 07:26:17
Date/Publication: 2013-03-02 17:05:32
22 changes: 11 additions & 11 deletions MD5
View file
Original file line number Diff line number Diff line change
@@ -1,20 +1,20 @@
f930ee455f9d3230b885b53c1b2a3ea1 *DESCRIPTION
8e398189e90bb9a837d92cc679f848b8 *DESCRIPTION
eeb7248f5e33d7458779501c9a9bbb73 *NAMESPACE
05b68595c51d4529303e7d32ccccfe47 *R/Function.R
411be20546043ef1b62084eebb547ced *R/Function.R
1dfb97d6d0ed6d4e023825ed0b22627a *R/Function_Power_Resampling.R
1fb2369f9bea6eb49ee1d0e8e1c34678 *R/KMTest_Linear.R
b00fbd5d834c1a0e42b9aab70657636e *R/KMTest_Linear.R
aeb01912fa7696ed7dae5b8012f1249a *R/KMTest_Logistic.R
73d981729a5570a4d4bcf7484c2e90f8 *R/KMTest_Logistic_VarMatching.R
3c1e8814623ff4d6177505c67b420f5a *R/KMTest_Optimal.R
e1155fbf7a2b4ae1fcf8bb07c496934d *R/KMTest_Logistic_VarMatching.R
19c391154790744bd551e3918d461b7f *R/KMTest_Optimal.R
5ff61f1e2188ea06c83a4b5f9b3345c8 *R/KMTest_Optimal_VarMatching.R
a270d1f8802d88edb99d9e9ed343d164 *R/Kernel.R
8bd5979c851ff3711c2e0b2a923a573c *R/Main.R
dedd98339da2922e136233176d41fdfa *R/Main_SSD.R
66ce4e65dc9e3840c24b919a3c6d2b5e *R/Null_Model.R
5704ca2797f81bf19b6e6ed7dfd96e6e *R/Power.R
05548a552bae2a0b51dbff39386e3fc7 *R/Power_Corr.R
b32e19496b2aa6411253eb0b6ec77aeb *R/Power.R
d2e084f5648696600b31e629ffe01786 *R/Power_Corr.R
192e270bf689a25409a2b16e5aaaab5e *R/SKAT_CompQuadForm.R
50a74b01e9676420fe143e0e7da8d215 *R/SSD.R
b8f1362618c38bdc46c2fc36ed2a4bf2 *R/SSD.R
0cefb76544272ec9f29b9ecee879af42 *data/SKAT.example.rda
3781bcbcae5bc581403b3a91960cebd5 *data/SKAT.haplotypes.rda
56ee07d0e42edb2859c948aed2354370 *data/datalist
Expand All @@ -23,7 +23,7 @@ ba60f2917612cca53407a1933d988f39 *inst/doc/Example1.SetID
e3b99ed0a36686c879c840e293e7687d *inst/doc/Example1.bim
9a8d91e7b3c8d88b3d465b2c5ee5e614 *inst/doc/Example1.fam
3fbd65beed55ea38efbe61d12b85d148 *inst/doc/SKAT.Rnw
9bc8443e1bc25652ca0c270a1eec4022 *inst/doc/SKAT.pdf
85c171b89352e4dddc26b4bdd063fc4c *inst/doc/SKAT.pdf
c15e095daf297308d467f36902f60e62 *man/Close_SSD.rd
d49220922e600dfd17ba774ebcab9aae *man/Generate_SSD_SetID.rd
f9800e7e2bc262cf021d87dd53df0876 *man/Get_Genotypes_SSD.rd
Expand All @@ -38,7 +38,7 @@ d92c7f5e671eb8b752fb6417c8866a6a *man/Power_Continuous.rd
3aba9ec25d378071cedfcb7fc8bb7554 *man/SKAT.SSD.All.rd
8d815f91cd850341802b9ecc334a1569 *man/SKAT.example.rd
32db14db92e122e11566931363c6ded8 *man/SKAT.haplotypes.rd
f9c8c73df507bcac10c74d19d6153502 *man/SKAT.rd
1002219e19e470781e07086b7eb2c633 *man/SKAT.rd
f1c70002116e5b2b7c888e38df99950e *man/SKAT_Null_Model.rd
e28440439fe71f4e6500ad5b642b2d96 *man/SKAT_Null_Model_MomentAdjust.rd
6abae0d1663d5a07b8e6cf01898f810c *man/SSD_FILE_OPEN.rd
Expand All @@ -47,7 +47,7 @@ c026767413f7c44c9eb102f2b06adc6f *src/DArray.h
562431ac9e3186630ba12aa47abf580f *src/NPsort.h
738c83d5696bc74457128b8e219e89fa *src/bed_reader.cpp
947d99c46efedcb27f4fb034ca117f0e *src/bed_reader.h
81d80f21da444ea9a62e61789d781fa0 *src/error_messages.h
40e3d557b770197e0e877d6e8076f651 *src/error_messages.h
0196768b51ff2bb278448a989fd33171 *src/interface.cpp
f3cc776422beec8e9bdf508cd990334f *src/interface_to_R.cpp
739723792bc105af0e258f3f265f7522 *src/interface_to_R.h
Expand Down
41 changes: 35 additions & 6 deletions R/Function.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,7 +122,7 @@ Get_Liu_PVal<-function(Q, W, Q.resampling = NULL){

Q.Norm<-(Q.all - param$muQ)/param$sigmaQ
Q.Norm1<-Q.Norm * param$sigmaX + param$muX
p.value<- 1-pchisq(Q.Norm1, df = param$l,ncp=param$d)
p.value<- pchisq(Q.Norm1, df = param$l,ncp=param$d, lower.tail=FALSE)

p.value.resampling = NULL
if(length(Q.resampling) > 0){
Expand All @@ -149,7 +149,7 @@ Get_Liu_PVal.MOD<-function(Q, W, Q.resampling = NULL){

Q.Norm<-(Q.all - param$muQ)/param$sigmaQ
Q.Norm1<-Q.Norm * param$sigmaX + param$muX
p.value<- 1-pchisq(Q.Norm1, df = param$l,ncp=param$d)
p.value<- pchisq(Q.Norm1, df = param$l,ncp=param$d, lower.tail=FALSE)

p.value.resampling = NULL
if(length(Q.resampling) > 0){
Expand All @@ -167,13 +167,32 @@ Get_Liu_PVal.MOD.Lambda<-function(Q.all, lambda){

Q.Norm<-(Q.all - param$muQ)/param$sigmaQ
Q.Norm1<-Q.Norm * param$sigmaX + param$muX
p.value<- 1-pchisq(Q.Norm1, df = param$l,ncp=param$d)
p.value<- pchisq(Q.Norm1, df = param$l,ncp=param$d, lower.tail=FALSE)

return(p.value)

}


Get_Liu_PVal.MOD.Lambda.Zero<-function(Q, muQ, muX, sigmaQ, sigmaX, l, d){


Q.Norm<-(Q - muQ)/sigmaQ
Q.Norm1<-Q.Norm * sigmaX + muX

temp<-c(0.05,10^-10, 10^-20,10^-30,10^-40,10^-50, 10^-60, 10^-70, 10^-80, 10^-90, 10^-100)
#qchisq(temp, df=1000000000,lower.tail=FALSE)
out<-qchisq(temp,df = l,ncp=d, lower.tail=FALSE)
#cat(c(Q.Norm1,l,d, out))
#cat("\n")
IDX<-max(which(out < Q.Norm1))

pval.msg<-sprintf("Pvalue < %e", temp[IDX])
return(pval.msg)

}


Get_Davies_PVal<-function(Q, W, Q.resampling = NULL){

K<-W/2
Expand All @@ -195,7 +214,8 @@ Get_Davies_PVal<-function(Q, W, Q.resampling = NULL){
}


re<-list(p.value = re$p.value[1], param=param,p.value.resampling = p.value.resampling )
re<-list(p.value = re$p.value[1], param=param,p.value.resampling = p.value.resampling
, pval.zero.msg=re$pval.zero.msg )
return(re)
}

Expand All @@ -219,7 +239,7 @@ Get_Lambda<-function(K){
stop("No Eigenvalue is bigger than 0!!")
}
lambda<-lambda1[IDX2]
lambda
return(lambda)

}

Expand Down Expand Up @@ -302,8 +322,17 @@ Get_PValue.Lambda<-function(lambda,Q){
p.val[i]<-p.val.liu[i]
}
}

p.val.msg = NULL
#cat(p.val[1])
if(p.val[1] == 0){

param<-Get_Liu_Params_Mod_Lambda(lambda)
p.val.msg<-Get_Liu_PVal.MOD.Lambda.Zero(Q[1], param$muQ, param$muX, param$sigmaQ, param$sigmaX, param$l, param$d)

}

return(list(p.value=p.val, p.val.liu=p.val.liu, is_converge=is_converge))
return(list(p.value=p.val, p.val.liu=p.val.liu, is_converge=is_converge, pval.zero.msg=p.val.msg))

}

Expand Down
16 changes: 12 additions & 4 deletions R/KMTest_Linear.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,21 +61,25 @@ KMTest.linear.Linear = function(res,Z,X1, kernel, weights, s2, method,res.out,n.
if( method == "liu" ){

out<-Get_Liu_PVal(Q, W.1, Q.res)
pval.zero.msg=NULL

} else if( method == "liu.mod" ){

out<-Get_Liu_PVal.MOD(Q, W.1, Q.res)
pval.zero.msg = NULL

} else if( method == "davies" ){

out<-Get_Davies_PVal(Q, W.1, Q.res)
pval.zero.msg = out$pval.zero.msg

} else {
stop("Invalid Method!")
}


re<-list(p.value = out$p.value, p.value.resampling = out$p.value.resampling, Test.Type = method, Q = Q, param=out$param )
re<-list(p.value = out$p.value, p.value.resampling = out$p.value.resampling
, Test.Type = method, Q = Q, param=out$param ,pval.zero.msg=pval.zero.msg )

return(re)

Expand Down Expand Up @@ -123,22 +127,26 @@ SKAT.linear.Other = function(res,Z,X1, kernel, weights = NULL, s2, method,res.ou
if( method == "liu" ){

out<-Get_Liu_PVal(Q, W,Q.res)
pval.zero.msg=NULL

} else if( method == "liu.mod" ){

out<-Get_Liu_PVal.MOD(Q, W,Q.res)
out<-Get_Liu_PVal.MOD(Q, W,Q.res)
pval.zero.msg = NULL

} else if( method == "davies" ){

out<-Get_Davies_PVal(Q, W1,Q.res)
out<-Get_Davies_PVal(Q, W1,Q.res)
pval.zero.msg = out$pval.zero.msg

} else {
stop("Invalid Method!")
}



re<-list(p.value = out$p.value, p.value.resampling = out$p.value.resampling, Test.Type = method, Q = Q, param=out$param )
re<-list(p.value = out$p.value, p.value.resampling = out$p.value.resampling
, Test.Type = method, Q = Q, param=out$param, pval.zero.msg=pval.zero.msg )

return(re)

Expand Down
50 changes: 33 additions & 17 deletions R/KMTest_Logistic_VarMatching.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,15 +63,19 @@ KMTest.logistic.Linear.VarMatching = function(res, Z, X1, kernel, weights = NULL
p.value.resampling<-re$p.value[-1]
}

p.value.noadj= re$p.value.noadj[1]
if(length(Q.all) > 1){
p.value.noadj=NULL
if(!is.null(re$p.value.noadj)){
p.value.noadj= re$p.value.noadj[1]
if(length(Q.all) > 1){

p.value.noadj.resampling<-re$p.value.noadj[-1]
p.value.noadj.resampling<-re$p.value.noadj[-1]
}
}

re<-list(p.value = p.value, p.value.resampling = p.value.resampling
, p.value.noadj = p.value.noadj, p.value.noadj.resampling = p.value.noadj.resampling
, Test.Type = method, Q = Q, Q.resampling = Q.res, param=NULL)
, Test.Type = method, Q = Q, Q.resampling = Q.res, param=NULL
, pval.zero.msg=re$pval.msg)

return(re)
}
Expand All @@ -87,26 +91,35 @@ SKAT_PValue_Logistic_VarMatching<-function(Q, Z1, p_all, Q.sim, type="Other"){
# aSKAT pvalue
Q.Norm<-(Q - param$muQ)/sqrt(param$varQ)
Q.Norm1<-Q.Norm * sqrt(2*param$df) + param$df
p.value<- 1-pchisq(Q.Norm1, df = param$df, ncp=0)
p.value<- pchisq(Q.Norm1, df = param$df, ncp=0, lower.tail=FALSE)

df1=param$df
if(is.null(Q.sim) && param$n.lambda==1){
re<-Get_Satterthwaite(param$muQ, param$varQ)
Q.Norm1<-Q / re$a
p.value<- 1-pchisq(Q.Norm1, df = re$df, ncp=0)

p.value<- pchisq(Q.Norm1, df = re$df, ncp=0, lower.tail=FALSE)
df1=re$df
#cat("df1:", re$df, "\n")

}
pval.msg = NULL
if(p.value[1] == 0){
pval.msg<-Get_Liu_PVal.MOD.Lambda.Zero(Q.Norm1, muQ=0, muX=0, sigmaQ=1, sigmaX=1, l=df1, d=0)

}


# SKAT pvalue
param.noadj<-param$param.noadj
Q.Norm<-(Q - param.noadj$muQ)/param.noadj$sigmaQ
Q.Norm1<-Q.Norm * param.noadj$sigmaX + param.noadj$muX
p.value.noadj<- 1-pchisq(Q.Norm1, df = param.noadj$l,ncp=0)

p.value.noadj = NULL
if(!is.null(param$param.noadj)){
param.noadj<-param$param.noadj
Q.Norm<-(Q - param.noadj$muQ)/param.noadj$sigmaQ
Q.Norm1<-Q.Norm * param.noadj$sigmaX + param.noadj$muX
p.value.noadj<- pchisq(Q.Norm1, df = param.noadj$l,ncp=0, lower.tail=FALSE)
}
#cat("df2:", param.noadj$l, "\n")

out<-list(p.value=p.value, p.value.noadj=p.value.noadj, param=param)
out<-list(p.value=p.value, p.value.noadj=p.value.noadj, param=param, pval.msg=pval.msg)

return(out)

Expand Down Expand Up @@ -181,13 +194,16 @@ SKAT_Logistic_VarMatching_GetParam1<-function(Z1, p_all, Q.sim, type="Other"){

}


################################
# ver 0.82 changed
SKAT_Logistic_VarMatching_GetParam1_OnlySim<-function(Z1, p_all, Q.sim){


out.svd = Get_Lambda_U_From_Z(Z1)
lambda<-out.svd$lambda
#out.svd = Get_Lambda_U_From_Z(Z1)
#lambda<-out.svd$lambda

lambda<-Get_Lambda(t(Z1) %*% Z1)


muQ = sum(lambda)
varQ.sim<-var(Q.sim)

Expand Down
30 changes: 23 additions & 7 deletions R/KMTest_Optimal.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,7 @@ SKAT_Optiaml_Each_Q<-function(param.m, Q.all, r.all, lambda.all, method=NULL){

# get pvalue
Q.Norm<-(Q - muQ)/sqrt(varQ) * sqrt(2*df) + df
pval[,i]<- 1-pchisq(Q.Norm, df = df)
pval[,i]<- pchisq(Q.Norm, df = df, lower.tail=FALSE)
# will be changed later

if(!is.null(method)){
Expand Down Expand Up @@ -142,18 +142,25 @@ SKAT_Optimal_Integrate_Func_Davies<-function(x,pmin.q,param.m,r.all){

}

SKAT_Optimal_PValue_Davies<-function(pmin.q,param.m,r.all){
# add pmin on 02-13-2013
SKAT_Optimal_PValue_Davies<-function(pmin.q,param.m,r.all, pmin=NULL){

#re<-try(integrate(SKAT_Optimal_Integrate_Func_Davies, lower=0, upper=30, subdivisions=500, pmin.q=pmin.q,param.m=param.m,r.all=r.all,abs.tol = 10^-15), silent = TRUE)

re<-try(integrate(SKAT_Optimal_Integrate_Func_Davies, lower=0, upper=40, subdivisions=1000, pmin.q=pmin.q,param.m=param.m,r.all=r.all,abs.tol = 10^-25), silent = TRUE)

if(class(re) == "try-error"){
re<-SKAT_Optimal_PValue_Liu(pmin.q,param.m,r.all)
re<-SKAT_Optimal_PValue_Liu(pmin.q,param.m,r.all, pmin)
return(re)
}

pvalue<-1-re[[1]]
if(!is.null(pmin)){
if(pmin *length(r.all) < pvalue){
pvalue = pmin *length(r.all)
}
}


return(pvalue)

Expand Down Expand Up @@ -183,13 +190,20 @@ SKAT_Optimal_Integrate_Func_Liu<-function(x,pmin.q,param.m,r.all){

}


SKAT_Optimal_PValue_Liu<-function(pmin.q,param.m,r.all){
# add pmin on 02-13-2013
SKAT_Optimal_PValue_Liu<-function(pmin.q,param.m,r.all, pmin=NULL){

re<-integrate(SKAT_Optimal_Integrate_Func_Liu, lower=0, upper=40, subdivisions=2000
,pmin.q=pmin.q,param.m=param.m,r.all=r.all,abs.tol = 10^-25)

pvalue<-1-re[[1]]

if(!is.null(pmin)){
if(pmin *length(r.all) < pvalue){
pvalue = pmin *length(r.all)
}
}

return(pvalue)

}
Expand Down Expand Up @@ -268,19 +282,21 @@ SKAT_Optimal_Get_Pvalue<-function(Q.all, Z1, r.all, method){
param.m<-SKAT_Optimal_Param(Z1,r.all)
Each_Info<-SKAT_Optiaml_Each_Q(param.m, Q.all, r.all, lambda.all, method=method)
pmin.q<-Each_Info$pmin.q
pmin<-Each_Info$pmin
pval<-rep(0,n.q)

if(method == "davies" || method=="optimal" || method=="optimal.mod" || method=="optimal.adj"){

for(i in 1:n.q){
pval[i]<-SKAT_Optimal_PValue_Davies(pmin.q[i,],param.m,r.all)
pval[i]<-SKAT_Optimal_PValue_Davies(pmin.q[i,],param.m,r.all, pmin[i])

}


} else if(method =="liu" || method =="liu.mod" || method=="optimal.moment" || method=="optimal.moment.adj" ){

for(i in 1:n.q){
pval[i]<-SKAT_Optimal_PValue_Liu(pmin.q[i,],param.m,r.all)
pval[i]<-SKAT_Optimal_PValue_Liu(pmin.q[i,],param.m,r.all, pmin[i])
}

} else {
Expand Down

0 comments on commit a34552b

Please sign in to comment.