All arguments are the same as in and passed intact to the ccfast.
See help for this function.
y
character name of the vector of case-control status. Cases are denoted as 1 and controls as 0.
data
An object of gwaa.data-class
snpsubset
Index, character or logical vector with subset of SNPs to run analysis on.
If missing, all SNPs from data are used for analysis.
idsubset
Index, character or logical vector with subset of IDs to run analysis on.
If missing, all people from data are used for analysis.
times
If more then one, the number of replicas to be used in derivation of
empirical genome-wide significance. See emp.qtscore, which
calls qtscore with times>1 for details
quiet
do not print warning messages
bcast
If the argument times > 1, progress is reported once in bcast replicas
Details
In the analysis of empirical significance, first time the function
ccfast is called and result object is
saved. Later, the function ccfast is called
times times with replace=FALSE in order to generate
the distribution under the null. Each call, minimal P-value is extracted
and compared with original P-values. For a particular SNP, empirical
P-value is obtained as a proportion of times minimal Ps from resampled data
was less then the original P.
The list elements effB, effAB and effBB are the ones obtained from the
analysis of the original (not permuted) data set
Value
Object of class scan.gwaa-class
Author(s)
Yurii Aulchenko
See Also
ccfast,
emp.qtscore,
scan.gwaa-class
Examples
require(GenABEL.data)
data(srdta)
a<-ccfast("bt",data=srdta,snps=c(500:800))
plot(a)
# this does not make sense, as the whole experiment must be analysed, not a small region!
# also, times = 20 is way too small (should be at least 200)
b<-emp.ccfast("bt",data=srdta,snps=c(500:800),bcast=10, times = 20)
plot(b)
# compare qvalues and empirical P
qv<-qvaluebh95(a[,"P1df"])$qval
qv
b[,"P1df"]
plot(qv,b[,"P1df"],xlim=c(0,1),ylim=c(0,1))
abline(a=0,b=1)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(GenABEL)
Loading required package: MASS
Loading required package: GenABEL.data
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/GenABEL/emp.ccfast.Rd_%03d_medium.png", width=480, height=480)
> ### Name: emp.ccfast
> ### Title: Genome-wide significance for a case-control GWA scan
> ### Aliases: emp.ccfast
> ### Keywords: htest
>
> ### ** Examples
>
> require(GenABEL.data)
> data(srdta)
> a<-ccfast("bt",data=srdta,snps=c(500:800))
Warning in ccfast("bt", data = srdta, snps = c(500:800)) :
11 people (out of 2500 ) excluded as not having cc status
> plot(a)
> # this does not make sense, as the whole experiment must be analysed, not a small region!
> # also, times = 20 is way too small (should be at least 200)
> b<-emp.ccfast("bt",data=srdta,snps=c(500:800),bcast=10, times = 20)
Warning in ccfast(y = y, data = data, snpsubset = snpsubset, idsubset = idsubset, :
11 people (out of 2500 ) excluded as not having cc status
50%100%> plot(b)
> # compare qvalues and empirical P
> qv<-qvaluebh95(a[,"P1df"])$qval
> qv
[1] 0.9826648 0.9826648 0.9826648 0.9826648 0.9826648 0.9837415 0.9826648
[8] 0.9935462 0.9826648 0.9826648 0.9826648 0.9826648 0.9826648 0.9826648
[15] 0.9826648 0.9826648 0.9826648 0.9826648 0.9602959 0.9826648 0.9602959
[22] 0.9826648 0.9973161 0.9826648 0.9826648 0.9826648 0.9935462 0.9826648
[29] 0.9826648 0.9826648 0.9973161 0.9826648 0.9973161 0.9826648 0.9826648
[36] 0.9602959 0.9826648 0.9826648 0.9602959 0.9826648 0.9602959 0.9602959
[43] 0.9826648 0.9973161 0.9826648 0.9826648 0.9935462 0.9826648 0.9826648
[50] 0.9826648 0.9826648 0.9826648 0.9973161 0.9826648 0.9602959 0.9602959
[57] 0.9602959 0.9837415 0.9826648 0.9973161 0.9602959 0.9602959 0.9837415
[64] 0.9826648 0.9826648 0.9826648 0.9826648 0.9602959 0.9602959 0.9602959
[71] 0.9826648 0.9826648 0.9973161 0.9602959 0.9602959 0.9826648 0.9826648
[78] 0.9973161 0.9826648 0.9826648 0.9826648 0.9826648 0.9826648 0.9826648
[85] 0.9602959 0.9826648 0.9899597 0.9973161 0.9850907 0.9973161 0.9826648
[92] 0.9720685 0.9826648 0.9826648 0.9837415 0.9826648 0.9826648 0.9826648
[99] 0.9973161 0.9826648 0.9826648 0.9826648 0.9973161 0.9826648 0.9826648
[106] 0.9899597 0.9826648 0.9826648 0.9826648 0.9826648 0.9973161 0.9973161
[113] 0.9973161 0.9850907 0.9985882 0.9826648 0.9973161 0.9834113 0.9973161
[120] 0.9826648 0.9826648 0.9973161 0.9602959 0.9602959 0.9602959 0.9826648
[127] 0.9973161 0.9973161 0.9602959 0.9826648 0.9826648 0.9826648 0.9826648
[134] 0.9826648 0.9602959 0.9602959 0.9826648 0.9826648 0.9826648 0.9602959
[141] 0.9602959 0.9602959 0.9826648 0.9826648 0.9826648 0.9826648 0.9935462
[148] 0.9955128 0.9935462 0.9973161 0.9826648 0.9935462 0.9826648 0.9826648
[155] 0.9826648 0.9602959 0.9826648 0.9602959 0.9602959 0.9826648 0.9826648
[162] 0.9826648 0.9720685 0.9826648 0.9826648 0.9826648 0.9826648 0.9826648
[169] 0.9826648 0.9602959 0.9826648 0.9826648 0.9826648 0.9973161 0.9973161
[176] 0.9602959 0.9826648 0.9602959 0.9602959 0.9826648 0.9720685 0.9826648
[183] 0.9826648 0.9826648 0.9973161 0.9602959 0.9826648 0.9826648 0.9973161
[190] 0.9826648 0.9826648 0.9826648 0.9973161 0.9826648 0.9899597 0.9712952
[197] 0.9602959 0.9602959 0.9602959 0.9602959 0.9826648 0.9602959 0.9826648
[204] 0.9826648 0.9826648 0.9826648 0.9837415 0.9826648 0.9826648 0.9826648
[211] 0.9826648 0.9826648 0.9973161 0.9826648 0.9826648 0.9973161 0.9826648
[218] 0.9826648 0.9973161 0.9602959 0.9826648 0.9602959 0.9985882 0.9826648
[225] 0.9826648 0.9602959 0.9602959 0.9602959 0.9826648 0.9826648 0.9602959
[232] 0.9602959 0.9826648 0.9826648 0.9826648 0.9826648 0.9826648 0.9973161
[239] 0.9973161 0.9602959 0.9826648 0.9826648 0.9826648 0.9826648 0.9826648
[246] 0.9826648 0.9837415 0.9973161 0.9602959 0.9973161 0.9826648 0.9826648
[253] 0.9602959 0.9602959 0.9826648 0.9826648 0.9602959 0.9826648 0.9602959
[260] 0.9973161 0.9826648 0.9826648 0.9826648 0.9826648 0.9826648 0.9826648
[267] 0.9826648 0.9602959 0.9602959 0.9973161 0.9826648 0.9602959 0.9837415
[274] 0.9826648 0.9826648 0.9826648 0.9602959 0.9826648 0.9826648 0.9602959
[281] 0.9826648 0.9826648 0.9602959 0.9602959 0.9826648 0.9826648 0.9973161
[288] 0.9602959 0.9826648 0.9602959 0.9826648 0.9837415 0.9826648 0.9901093
[295] 0.9602959 0.9602959 0.9826648 0.9602959 0.9973161 0.9602959 0.9826648
> b[,"P1df"]
[1] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[16] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[31] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[46] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[61] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[76] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[91] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[106] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[121] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[136] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[151] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[166] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[181] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[196] 1.00 1.00 1.00 1.00 1.00 1.00 0.85 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[211] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[226] 1.00 1.00 0.85 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[241] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[256] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[271] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[286] 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
[301] 1.00
> plot(qv,b[,"P1df"],xlim=c(0,1),ylim=c(0,1))
> abline(a=0,b=1)
>
>
>
>
>
> dev.off()
null device
1
>