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.
hwx.test(c, method = "auto", cutoff = 1e+07, B = 1e+05,
statName = "LLR", histobins = 0, histobounds = c(0, 0), showCurve = T,
safeSecs = 100, detail = 2)
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.
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
.
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.
The number of trials to perform if Monte Carlo method is used
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.
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.
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.
whether to show a blue curve indicating the asymptotic (chi squared) distribution. This only works for LLR
and Chisq
After this many seconds the calculation will be aborted. This is a safety valve to prevent attempts to compute impossibly large sets of tables.
Determines how much detail is printed. If it is set to 0, nothing is printed (useful if you use hwx.test
programmatically.)
Returns a list of class hwtest
which includes the following items:
The four computed P values corresponding to the test statistics: LLR
, Prob
, U
and Chisq
in that order.
The four observed statistics in the same order as above
The number of tables examined during the calculation if done by Monte Carlo
The total number of tables if done by full enumeration
The input matrix of genotype counts
The allele counts \(m\) corresponding to the input genotype counts
Which statistic to use for the histogram and in the p.value
item
Which method was used, “exact” or “monte”
An integer indicating how much detail to print. Use 0 for no printing
vector with the standard error for each stat. Only applicable with Monte Carlo tests
The methods are described by Engels, 2009. Genetics 183:1431.
# NOT RUN {
# 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)
# }
Run the code above in your browser using DataLab