R/ustat.R
hweustat.RdEstimates double reduction and tests for equilibrium while accounting for double reduction. It does this using an approach called "U-process minimization", where we minimize a function of a U-statistic that should be 0 at equilibrium given the true double reduction rate.
hweustat(nvec, thresh = NULL, effdf = TRUE)A vector containing the observed genotype counts,
where nvec[[i]] is the number of individuals with genotype
i-1. This should be of length ploidy+1.
The threshold for ignoring the genotype. We keep
genotypes such that nvec >= thresh.
Setting this to 0 uses all genotypes. Setting this to
NULL uses a heuristic that works well in practice.
A logical. Should we use the ad-hoc
"effective degrees of freedom" (TRUE) or not (FALSE)?
A list with some or all of the following elements:
alphaThe estimated double reduction parameter(s).
In diploids, this value is NULL.
chisq_hweThe chi-square test statistic for testing against the null of equilibrium.
df_hweThe degrees of freedom associated with
chisq_hwe.
p_hweThe p-value against the null of equilibrium.
This is a two-step estimator, where we first obtain a consistent estimate of the double reduction parameter, use this to estimate the covariance of estimators, then use this to obtain our final estimate of the double reduction parameter.
set.seed(1)
ploidy <- 6
size <- 1000
r <- 0.1
alpha <- 0.1
qvec <- hwefreq(r = r, alpha = alpha, ploidy = ploidy)
nvec <- c(rmultinom(n = 1, size = size, prob = qvec))
hweustat(nvec = nvec)
#> $alpha
#> [1] 0.09900657
#>
#> $chisq_hwe
#> [1] 2.936928
#>
#> $df_hwe
#> [1] 2
#>
#> $p_hwe
#> [1] 0.2302789
#>