Plot of an Hardy-Weinberg Testing Analysis

Share:

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

1
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@unibocconi.it

References

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

See Also

hwe.ibf, hwe.ibf.mc.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
# 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)