| hwe.ibf.plot {HWEintrinsic} | 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 |
|
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 |
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 |
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
See Also
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)