hwe.ibf.plot: Plot of an Hardy-Weinberg Testing Analysis
In sergioventurini/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

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)

sergioventurini/HWEintrinsic documentation built on Oct. 19, 2023, 12:40 a.m.