hwx.test: 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.

Usage

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)

Arguments

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

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

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 enumeratio

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:
$ PvaluesThe four computed P values corresponding to the test statistics: LLR, Prob, U and Chisq in that order.
$ observedThe four observed statistics in the same order as above
$ ntrialsThe number of tables examined during the calculation if done by Monte Carlo
$ tableCountThe total number of tables if done by full enumeration
$ genotypesThe input matrix of genotype counts
$ allelesThe allele counts $m$ corresponding to the input genotype counts
$ statNameWhich statistic to use for the histogram and in the p.value item
$ methodWhich method was used, exact or monte
$ detailAn integer indicating how much detail to print. Use 0 for no printing
$ SEvector with the standard error for each stat. Only applicable with Monte Carlo tests

References

The methods are described by http://dx.doi.org/10.1534/genetics.109.108977{Engels, 2009. Genetics 183:1431}.

Examples

Run this code

# 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