R: Test for HW by either full enumeration or Monte Carlo.
hwx.test
R Documentation
Test for HW by either full enumeration or Monte Carlo.
Description
The hwx.test() function is the main function of the HWxtest package. This function produces a valid test for Hardy-Weinberg frequencies for virtually any set of genotype counts. It will use either a full-enumeration method in which all possible tables with the same allele numbers are examined, or a Monte Carlo test where a large number of random tables is examined. To decide which to use, it calls xcountCutoff to determine whether the number of tables to examine is greater than cutoff. If it is, then mtest is used. Otherwise xtest is used. The result is a robust test which will always provide a meaningful and accurate P value. Each table examined is compared with the observed counts according to each of four measures of fit: “LLR”, “Prob”, “U”, or “Chisq” corresponding to the log-likelihood ratio, the null-hypothesis probability, the U-score or the Pearson X^2 value. It can also plot a histogram showing the distribution of any of these statistics.
The genotype counts. You must provide the number of each genotype. So if there are k alleles, you need to include the number of each of the k(k+1)/2 genotypes. The format of x is somewhat flexible: It can be a square matrix, but only the lower-left half is used. It can be a vector of the observations in the order a_11, a_21, a_22, a_31, ..., a_kk. For compatability with the packages genetics and adegenet, it can also be an object of class genind, genotype, or a data.frame. If c contains multiple samples, the parallel package will be used in an attempt to employ multi-cores.
method
Can be “auto”, “exact” or “monte” to indicate the method to use. If “auto”, the hwx.test will first check to see whether the total number of tables exceeds a cutoff specified by the parameter cutoff.
cutoff
If method is set to “auto”, then cutoff is used to decide whether to perform the test via the full enumeration or Monte Carlo method. If the number of tables is less than cutoff, then a full enumeration is performed. Otherwise the method will be Monte Carlo with B random trials.
B
The number of trials to perform if Monte Carlo method is used
statName
can be “LLR”, “Prob”, “U”, or “Chisq” depending on which one is to be ploted. Note that P values for all four are computed regardless of which one is specified with this parameter.
histobins
If 0, no histogram is plotted. If 1 or TRUE a histogram with 500 bins is plotted. If histobins is set to a number greater than 1, a histogram with histobins bins is plotted.
histobounds
A vector containing the left and right boundaries for the histogram's x axis. If you leave this as the default, c(0,0), then hwx.test will compute reasonable bounds to include most of the distribution.
showCurve
whether to show a blue curve indicating the asymptotic (chi squared) distribution. This only works for LLR and Chisq
safeSecs
After this many seconds the calculation will be aborted. This is a safety valve to prevent attempts to compute impossibly large sets of tables.
detail
Determines how much detail is printed. If it is set to 0, nothing is printed (useful if you use hwx.test programmatically.)
Value
Returns a list of class hwtest which includes the following items:
$ Pvalues
The four computed P values corresponding to the test statistics: LLR, Prob, U and Chisq in that order.
$ observed
The four observed statistics in the same order as above
$ ntrials
The number of tables examined during the calculation if done by Monte Carlo
$ tableCount
The total number of tables if done by full enumeration
$ genotypes
The input matrix of genotype counts
$ alleles
The allele counts m corresponding to the input genotype counts
$ statName
Which statistic to use for the histogram and in the p.value item
$ method
Which method was used, “exact” or “monte”
$ detail
An integer indicating how much detail to print. Use 0 for no printing
$ SE
vector with the standard error for each stat. Only applicable with Monte Carlo tests
# Data from Louis and Dempster 1987 Table 2 and Guo and Thompson 1992 Figure 2:
c <- c(0,3,1,5,18,1,3,7,5,2)
hwx.test(c)
# To see a histogram of the LLR statistic:
hwx.test(c, histobins=TRUE)
# For a histogram of the U statistic and other details of the result:
hwx.test(c, statName="U", histobins=TRUE, detail=3)
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(HWxtest)
------------------------------------------------
HWxtest
version 1.1.7
Please see the package vignette for instructions
------------------------------------------------
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/HWxtest/hwx.test.Rd_%03d_medium.png", width=480, height=480)
> ### Name: hwx.test
> ### Title: Test for HW by either full enumeration or Monte Carlo.
> ### Aliases: hwx.test
>
> ### ** Examples
>
> # Data from Louis and Dempster 1987 Table 2 and Guo and Thompson 1992 Figure 2:
> c <- c(0,3,1,5,18,1,3,7,5,2)
> hwx.test(c)
***** Sample of 45 diploids with 4 alleles
Full enumeration of 162365 tables to test for HW
P value (LLR) = 0.012945
P value (Prob) = 0.017442
P value (U) = 0.003343 (test for heterozygote excess)
P value (Chisq) = 0.020170
> # To see a histogram of the LLR statistic:
> hwx.test(c, histobins=TRUE)
***** Sample of 45 diploids with 4 alleles
Full enumeration of 162365 tables to test for HW
P value (LLR) = 0.012945
P value (Prob) = 0.017442
P value (U) = 0.003343 (test for heterozygote excess)
P value (Chisq) = 0.020170
Loading required namespace: ggplot2
> # For a histogram of the U statistic and other details of the result:
> hwx.test(c, statName="U", histobins=TRUE, detail=3)
***** Sample of 45 diploids with 4 alleles
Full enumeration of 162365 tables to test for HW
P value (LLR) = 0.012945
P value (Prob) = 0.017442
P value (U) = 0.003343 (test for heterozygote excess)
P value (Chisq) = 0.020170
Observed Test Statistics:
LLR : -8.591
Prob : 1.808e-06
U : -29.53
Chisq : 14.63
>
>
>
>
>
>
> dev.off()
null device
1
>
providing 
.
