lik.multin | R Documentation |
This function provides the value of the likelihood function of the full (unrestricted) model, as described in Consonni et al. (2011).
lik.multin(y, p11, p21)
y |
an object of class "HWEdata". |
p11 |
gentotype proportion for the alleles pair |
p21 |
gentotype proportion for the alleles pair |
This function has been used to generate the likelihood contours that appear in Figure 4 of Consonni et al. (2011) (see the example below).
This function returns the numerical value of the likelihood in correspondence of the argument values.
This function provides the likelihood function value only for the two alleles case.
Sergio Venturini sergio.venturini@unicatt.it
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
,
hwe.ibf.plot
.
# 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.