# hwe.ibf.plot: Plot of an Hardy-Weinberg Testing Analysis In HWEintrinsic: Objective Bayesian Testing for the Hardy-Weinberg Equilibrium Problem

 hwe.ibf.plot R Documentation

## Plot of an Hardy-Weinberg Testing Analysis

### Description

Plot of the null posterior probability of a Hardy-Weinberg testing problem based on intrinsic priors as described in Consonni et al. (2011).

### Usage

``````hwe.ibf.plot(y, t.vec, M = 1e+05, bf = FALSE)
``````

### Arguments

 `y` an object of class "HWEdata". `t.vec` vector of training sample sizes. `M` number of Monte Carlo iterations. `bf` logical: if TRUE the plot reports the Bayes factor based on intrinsic priors.

### Details

This function allows to create a plot of the null posterior probability versus a given set of training sample sizes. It simply performs a repeated analysis using `hwe.ibf.mc` on each of the training sample sizes contained in `t.vec`.

### Value

`hwe.ibf.plot` returns as the output an invisible list with the following components:

 `mc` matrix containing the Monte Carlo estimates of the Bayes factor and the null posterior probability for each training sample size in `t.vec`. `std` vector containing the standard Bayes factor and the null posterior probability for the data in hand. `exact` matrix containing the exact values of the Bayes factor and the null posterior probability for each training sample size in `t.vec`; available only for the two alleles case.

### Note

The Bayes factor computed here is for the unrestricted model (`M_1`) against the Hardy-Weinberg case (`M_0`).

### Author(s)

Sergio Venturini sergio.venturini@unicatt.it

### References

Consonni, G., Moreno, E., and Venturini, S. (2011). "Testing Hardy-Weinberg equilibrium: an objective Bayesian analysis". Statistics in Medicine, 30, 62–74. https://onlinelibrary.wiley.com/doi/10.1002/sim.4084/abstract

`hwe.ibf`, `hwe.ibf.mc`.

### Examples

``````# The following code reproduces Figure 4 in Consonni et al. (2011) #
## Not run:
# ATTENTION: it may take a long time to run! #

data(simdata)
n <- sum(dataset1@data.vec, na.rm = TRUE)
f <- c(.1,.5,1)
t <- round(f*n)
p11 <- p21 <- seq(0,1,length.out=100)
ip <- array(NA,c(length(f),length(p11),length(p21)))
for (k in 1:length(f)) {
ip[k,,] <- outer(X = p11, Y = p21, FUN = Vectorize(ip.tmp), t[k])
print(paste(k," / ",length(f),sep=""), quote = FALSE)
}

r <- 2
R <- r*(r + 1)/2
l <- 4
tables <- matrix(NA, nrow = R, ncol = l)
tables[, 1] <- dataset1@data.vec
tables[, 2] <- dataset2@data.vec
tables[, 3] <- dataset3@data.vec
tables[, 4] <- dataset4@data.vec
lik <- array(NA, c(l, length(p11), length(p21)))
M <- 300000
par(mfrow = c(4, 4))
for (k in 1:l) {
y <- new("HWEdata", data = tables[, k])
lik[k,,] <- lik.multin(y, p11, p21)

nlev <- 10
for (q in 1:length(f)) {
contour(p11, p21, ip[q,,], xlab = expression(p[11]),
ylab = expression(p[21]), nlevels = nlev, col = gray(0),
main = "", cex.axis = 1.75, cex.lab = 1.75, labcex = 1.4)
lines(p11^2, 2*p11*(1 - p11), lty = "longdash", col = gray(0), lwd = 2)
contour(p11, p21, lik[k,,], nlevels = nlev, add = TRUE,
col = gray(.7), labcex = 1.2)
abline(a = 1, b = -1, lty = 3, col = gray(.8))
}
hwe.ibf.plot(y = y, t.vec = seq(1,n,1), M = M)
}

## End(Not run)
``````

HWEintrinsic documentation built on Sept. 8, 2023, 5:56 p.m.