Package 'HWEintrinsic' reference manual

Title:	Objective Bayesian Testing for the Hardy-Weinberg Equilibrium Problem
Description:	General (multi-allelic) Hardy-Weinberg equilibrium problem from an objective Bayesian testing standpoint. This aim is achieved through the identification of a class of priors specifically designed for this testing problem. A class of intrinsic priors under the full model is considered. This class is indexed by a tuning quantity, the training sample size, as discussed in Consonni, Moreno and Venturini (2010). These priors are objective, satisfy Savage's continuity condition and have proved to behave extremely well for many statistical testing problems.
Authors:	Sergio Venturini
Maintainer:	Sergio Venturini <[email protected]>
License:	GPL (>= 2) \| file LICENSE
Version:	1.2.3
Built:	2024-12-01 08:41:13 UTC
Source:	CRAN

Objective Bayesian Testing for the Hardy-Weinberg Equilibrium Problem

Description

General (multi-allelic) Hardy-Weinberg equilibrium problem from an objective Bayesian testing standpoint. This aim is achieved through the identification of a class of priors specifically designed for this testing problem. A class of intrinsic priors under the full model is considered. This class is indexed by a tuning quantity, the training sample size, as discussed in Consonni, Moreno and Venturini (2010). These priors are objective, satisfy Savage's continuity condition and have proved to behave extremely well for many statistical testing problems.

Details

The package is loaded with the usual library(HWEintrinsic) command. The most important functions are hwe.ibf, hwe.ibf.mc, and hwe.ibf.plot. It contains also some data sets which have been extensively used in the literature.

Author(s)

Sergio Venturini

Maintainer: Sergio Venturini <[email protected]>

References

Graph of the Conditional Intrinsic Prior for a Two Alleles Dataset

Description

This function produces the 3D plot for the conditional intrinsic prior based on a sample of two alleles data as described in Consonni et al. (2011).

Usage

  cip.2(t, p, k = 30)
cip.2(t, p, k = 30)

Arguments

`t`	training sample size.
`p`	allele frequency used to condition the prior upon.
`k`	number of grid points for the alleles frequencies at which the prior is evaluated.

Value

No object is returned by this function.

Note

This function provides the plot of the conditional intrinsic prior only for two alleles data.

Author(s)

Sergio Venturini [email protected]

References

Examples

## Not run: 
# ATTENTION: the following code may take a long time to run! #

### The following code reproduces Figure 1 in Consonni et al. (2011) ###
p <- 0.5
t <- 5
cip.2(t, p, 30)
t <- 30
cip.2(t, p, 50)

## End(Not run)
## Not run: 
# ATTENTION: the following code may take a long time to run! #

### The following code reproduces Figure 1 in Consonni et al. (2011) ###
p <- 0.5
t <- 5
cip.2(t, p, 30)
t <- 30
cip.2(t, p, 50)

## End(Not run)

Utility Function

Description

This function provides the calculations needed to plot the conditional intrinsic prior based for a two alleles dataset as described in Consonni et al. (2011).

Usage

  cip.tmp(p11, p21, t, p)
cip.tmp(p11, p21, t, p)

Arguments

`p11`	gentotype proportion for the alleles pair $A_{1}A_{1}$ .
`p21`	gentotype proportion for the alleles pair $A_{2}A_{1}$ .
`t`	training sample size.
`p`	allele frequency used to condition the prior upon.

Value

cip.tmp returns the value of the conditional intrinsic prior evaluated at the arguments values.

Author(s)

Sergio Venturini [email protected]

References

Guo and Thompson (1992) Eight Alleles Simulated Data

Description

A sample of size $n = 30$ of genotype frequencies from a population of $r = 8$ alleles simulated under the Hardy-Weinberg equilibrium when the underlying gene frequencies are (.2, .2, .2, .2, .05, .05, .05, .05); see Guo and Thompson (1992), Example 2.

Usage

data(GuoThompson8)data(GuoThompson8)

Format

An object of class HWEdata.

Source

Guo, S.W. and Thompson, E.A. (1992), "Performing the Exact Test of Hardy-Weinberg Proportion for Multiple Alleles". Biometrics, Vol. 48, No. 2, 361–372.

References

Examples

# Example 1 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #


data(GuoThompson8)
plot(GuoThompson8)
n <- sum([email protected], na.rm = TRUE)
out <- hwe.ibf.mc(GuoThompson8, t = n/2, M = 100000, verbose = TRUE)
summary(out, plot = TRUE)

## End(Not run)

# Example 2 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #

M <- 300000
f <- seq(.1, 1, .05)
n <- sum([email protected], na.rm = TRUE)
out <- hwe.ibf.plot(y = GuoThompson8, t.vec = round(f*n), M = M)

## End(Not run)
# Example 1 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #


data(GuoThompson8)
plot(GuoThompson8)
n <- sum(GuoThompson8@data.vec, na.rm = TRUE)
out <- hwe.ibf.mc(GuoThompson8, t = n/2, M = 100000, verbose = TRUE)
summary(out, plot = TRUE)

## End(Not run)

# Example 2 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #

M <- 300000
f <- seq(.1, 1, .05)
n <- sum(GuoThompson8@data.vec, na.rm = TRUE)
out <- hwe.ibf.plot(y = GuoThompson8, t.vec = round(f*n), M = M)

## End(Not run)

Guo and Thompson (1992) Nine Alleles Data

Description

These data are extracted (Guo and Thompson, 1992) from the Rhesus data in Cavalli-Sforza and Bodmer (1971). They consists of information on 8297 individuals with $r = 9$ groups.

Usage

data(GuoThompson9)data(GuoThompson9)

Format

An object of class HWEdata.

Source

Guo, S.W. and Thompson, E.A. (1992), "Performing the Exact Test of Hardy-Weinberg Proportion for Multiple Alleles". Biometrics, Vol. 49, No. 2, 361–372.

References

Cavalli-Sforza, L. and Bodmer, W. (1971), "The Genetics of Human Populations". W.H. Freeman and Company, San Francisco. 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

# Example 1 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #

data(GuoThompson9)
plot(GuoThompson9)
n <- sum([email protected], na.rm = TRUE)
out <- hwe.ibf.mc(GuoThompson9, t = n/2, M = 100000, verbose = TRUE)
summary(out, plot = TRUE)

## End(Not run)

# Example 2 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #

M <- 300000
f <- seq(.1, 1, .05)
n <- sum([email protected], na.rm = TRUE)
out <- hwe.ibf.plot(y = GuoThompson9, t.vec = round(f*n), M = M)

## End(Not run)
# Example 1 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #

data(GuoThompson9)
plot(GuoThompson9)
n <- sum(GuoThompson9@data.vec, na.rm = TRUE)
out <- hwe.ibf.mc(GuoThompson9, t = n/2, M = 100000, verbose = TRUE)
summary(out, plot = TRUE)

## End(Not run)

# Example 2 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #

M <- 300000
f <- seq(.1, 1, .05)
n <- sum(GuoThompson9@data.vec, na.rm = TRUE)
out <- hwe.ibf.plot(y = GuoThompson9, t.vec = round(f*n), M = M)

## End(Not run)

Standard Bayes Factor for the Hardy-Weinebrg Testing Problem

Description

This function provides the calculations for obtaining the standard Bayes factor for the Hardy-Weinberg testing problem. It implements a common default prior (constant) for both the null and the alternative models.

Usage

  hwe.bf(y)
hwe.bf(y)

Arguments

`y`	an object of class "HWEdata".

Value

hwe.bf returns the standard Bayes Factor value for the Hardy-Weinberg testing problem (see the references for the details).

Note

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

Author(s)

Sergio Venturini [email protected]

References

Examples

# Example 1 #
data(GuoThompson8)
plot(GuoThompson8)
hwe.bf(GuoThompson8)

# Example 2 #
data(LouisDempster)
plot(LouisDempster)
hwe.bf(LouisDempster)
# Example 1 #
data(GuoThompson8)
plot(GuoThompson8)
hwe.bf(GuoThompson8)

# Example 2 #
data(LouisDempster)
plot(LouisDempster)
hwe.bf(LouisDempster)

Testing Hardy-Weinberg Equilibrium Using an Intrinsic Prior Approach

Description

This function implements the exact calculation of the Bayes factor based on intrinsic priors for the Hardy-Weinberg testing problem as described in Consonni et al. (2011).

Usage

hwe.ibf(y, t)
hwe.ibf(y, t)

Arguments

`y`	an object of class "HWEdata".
`t`	training sample size.

Details

This function implements the exact formula for the Bayes factor based on intrinsic priors.

Value

hwe.ibf returns the value of the Bayes factor based on intrinsic priors.

Note

The Bayes factor computed here is for the unrestricted model ( $M_1$ ) against the Hardy-Weinberg case ( $M_0$ ). This function provides the output only for the two alleles case.

Author(s)

Sergio Venturini [email protected]

References

Examples

## Not run: 
# ATTENTION: the following code may take a long time to run! #

data(Lindley)
hwe.ibf.exact <- Vectorize(hwe.ibf, "t")
f <- seq(.05, 1, .05)
n <- sum([email protected], na.rm = TRUE)

# Dataset 1 #
plot(dataL1)
npp.exact <- 1/(1 + hwe.ibf.exact(round(f*n), y = dataL1))
npp.std <- 1/(1 + hwe.bf(dataL1))
plot(f, npp.exact, type="l", lwd = 2, xlab = "f = t/n",
	ylab = "Null posterior probability")
abline(h = npp.std, col = gray(.5), lty = "longdash")

# Dataset 2 #
plot(dataL2)
npp.exact <- 1/(1 + hwe.ibf.exact(round(f*n), y = dataL2))
npp.std <- 1/(1 + hwe.bf(dataL2))
plot(f, npp.exact, type="l", lwd = 2, xlab = "f = t/n",
	ylab = "Null posterior probability")
abline(h = npp.std, col = gray(.5), lty = "longdash")

# Dataset 3 #
plot(dataL3)
npp.exact <- 1/(1 + hwe.ibf.exact(round(f*n), y = dataL3))
npp.std <- 1/(1 + hwe.bf(dataL3))
plot(f, npp.exact, type="l", lwd = 2, xlab = "f = t/n",
	ylab = "Null posterior probability")
abline(h = npp.std, col = gray(.5), lty = "longdash")

# Dataset 4 #
plot(dataL4)
npp.exact <- 1/(1 + hwe.ibf.exact(round(f*n), y = dataL4))
npp.std <- 1/(1 + hwe.bf(dataL4))
plot(f, npp.exact, type="l", lwd = 2, xlab = "f = t/n",
	ylab = "Null posterior probability")
abline(h = npp.std, col = gray(.5), lty = "longdash")

## End(Not run)
## Not run: 
# ATTENTION: the following code may take a long time to run! #

data(Lindley)
hwe.ibf.exact <- Vectorize(hwe.ibf, "t")
f <- seq(.05, 1, .05)
n <- sum(dataL1@data.vec, na.rm = TRUE)

# Dataset 1 #
plot(dataL1)
npp.exact <- 1/(1 + hwe.ibf.exact(round(f*n), y = dataL1))
npp.std <- 1/(1 + hwe.bf(dataL1))
plot(f, npp.exact, type="l", lwd = 2, xlab = "f = t/n",
	ylab = "Null posterior probability")
abline(h = npp.std, col = gray(.5), lty = "longdash")

# Dataset 2 #
plot(dataL2)
npp.exact <- 1/(1 + hwe.ibf.exact(round(f*n), y = dataL2))
npp.std <- 1/(1 + hwe.bf(dataL2))
plot(f, npp.exact, type="l", lwd = 2, xlab = "f = t/n",
	ylab = "Null posterior probability")
abline(h = npp.std, col = gray(.5), lty = "longdash")

# Dataset 3 #
plot(dataL3)
npp.exact <- 1/(1 + hwe.ibf.exact(round(f*n), y = dataL3))
npp.std <- 1/(1 + hwe.bf(dataL3))
plot(f, npp.exact, type="l", lwd = 2, xlab = "f = t/n",
	ylab = "Null posterior probability")
abline(h = npp.std, col = gray(.5), lty = "longdash")

# Dataset 4 #
plot(dataL4)
npp.exact <- 1/(1 + hwe.ibf.exact(round(f*n), y = dataL4))
npp.std <- 1/(1 + hwe.bf(dataL4))
plot(f, npp.exact, type="l", lwd = 2, xlab = "f = t/n",
	ylab = "Null posterior probability")
abline(h = npp.std, col = gray(.5), lty = "longdash")

## End(Not run)

Testing Hardy-Weinberg Equilibrium Using an Intrinsic Prior Approach

Description

This function implements the Monte Carlo estimation of the Bayes factor based on intrinsic priors for the Hardy-Weinberg testing problem as described in Consonni et al. (2011).

Usage

hwe.ibf.mc(y, t, M = 10000, verbose = TRUE)
hwe.ibf.mc(y, t, M = 10000, verbose = TRUE)

Arguments

`y`	an object of class "HWEdata".
`t`	training sample size.
`M`	number of Monte Carlo iterations.
`verbose`	logical; if TRUE the function prints the detailed calculation progress.

Details

This function implements a Monte Carlo approximation using importance sampling of the Bayes factor based on intrinsic priors.

Value

hwe.ibf.mc returns an object of the class "HWEintr".

Note

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

Author(s)

Sergio Venturini [email protected]

References

Examples

# Example 1 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #


data(GuoThompson9)
plot(GuoThompson9)
n <- sum([email protected], na.rm = TRUE)
out <- hwe.ibf.mc(GuoThompson9, t = n/2, M = 100000, verbose = TRUE)
summary(out, plot = TRUE)

## End(Not run)

# Example 2 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #

M <- 300000
f <- seq(.1, 1, .05)
n <- sum([email protected], na.rm = TRUE)
out <- hwe.ibf.plot(y = GuoThompson9, t.vec = round(f*n), M = M)

## End(Not run)
# Example 1 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #


data(GuoThompson9)
plot(GuoThompson9)
n <- sum(GuoThompson9@data.vec, na.rm = TRUE)
out <- hwe.ibf.mc(GuoThompson9, t = n/2, M = 100000, verbose = TRUE)
summary(out, plot = TRUE)

## End(Not run)

# Example 2 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #

M <- 300000
f <- seq(.1, 1, .05)
n <- sum(GuoThompson9@data.vec, na.rm = TRUE)
out <- hwe.ibf.plot(y = GuoThompson9, t.vec = round(f*n), M = M)

## End(Not run)

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)
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 [email protected]

References

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([email protected], 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] <- [email protected]
tables[, 2] <- [email protected]
tables[, 3] <- [email protected]
tables[, 4] <- [email protected]
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)
# 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)

Class "HWEdata". Data specification for the Hardy-Weinberg Testing Problem Using the Intrinsic Prior Approach.

Description

This class encapsulates the data specification for a Bayesian objective analysis via intrinsic priors of the Hardy-Weinberg Testing Problem as described in Consonni et al. (2011).

Objects from the Class

Objects can be created by calls of the form new("HWEdata", data), where data are the data in vector form.

Slots

data.mat:: Object of class "matrix"; data in matrix form.
size:: Object of class "numeric"; number of alleles included in the data.
data.vec:: Object of class "numeric"; data in vector form.

Methods

plot: signature(x = "HWEdata", y = "missing"): Provides a pictorial representation for a sample of genotype counts.
summary: signature(object = "HWEdata"): Extracts the contents of an HWEdata object.

Author(s)

Sergio Venturini [email protected]

References

Examples

data.tmp <- c(3, 9, 8)
dataset <- new("HWEdata", data = data.tmp)
data.tmp <- c(3, 9, 8)
dataset <- new("HWEdata", data = data.tmp)

Class "HWEintr". Result for the Hardy-Weinberg Testing Problem Using the Intrinsic Prior Approach.

Description

This class encapsulates the results of a Bayesian objective analysis via intrinsic priors for the Hardy-Weinberg Testing Problem as described in Consonni et al. (2011).

Objects from the Class

Objects can be created by calls of the form new("HWEintr", bf, npp, draws, data), but most often as the result of a call to hwe.ibf or to hwe.ibf.mc.

Slots

bf:: Object of class "numeric"; Bayes factor based on intrinsic priors.
npp:: Object of class "numeric"; posterior probability of the null Hardy-Weinberg model.
draws:: Object of class "numeric"; individual terms of the Monte Carlo sum using importance sampling.
data.mat:: Object of class "matrix"; original data in matrix form.

Methods

plot: signature(x = "HWEintr", y = "missing"): Provides a graphical representation of the estimates.
summary: signature(object = "HWEintr"): Summarizes the information about the estimates.

Author(s)

Sergio Venturini [email protected]

References

Graph of the (Unconditional) Intrinsic Prior for a Two Alleles Dataset

Description

This function produces the 3D-plot for the (unconditional) intrinsic prior based on a sample of two alleles data as described in Consonni et al. (2011).

Usage

  ip.2(t, k = 30)
ip.2(t, k = 30)

Arguments

`t`	training sample size.
`k`	number of grid points for the alleles frequencies at which the prior is evaluated.

Value

No object is returned by this function.

Note

This function provides the plot of the (unconditional) intrinsic prior only for two alleles data.

Author(s)

Sergio Venturini [email protected]

References

Examples

## Not run: 
# ATTENTION: the following code may take a long time to run! #

### The following code reproduces Figure 3 in Consonni et al. (2011) ###
t <- 30
ip.2(t, 40)

## End(Not run)
## Not run: 
# ATTENTION: the following code may take a long time to run! #

### The following code reproduces Figure 3 in Consonni et al. (2011) ###
t <- 30
ip.2(t, 40)

## End(Not run)

Utility Function

Description

This function provides the calculations needed to plot the (unconditional) intrinsic prior based for a two alleles dataset as described in Consonni et al. (2011).

Usage

  ip.tmp(p11, p21, t)
ip.tmp(p11, p21, t)

Arguments

`p11`	gentotype proportion for the the pair of alleles $A_{1}A_{1}$ .
`p21`	gentotype proportion for the the pair of alleles $A_{2}A_{1}$ .
`t`	training sample size.

Value

ip.tmp returns the value of the (unconditional) intrinsic prior evaluated at the arguments values.

Author(s)

Sergio Venturini [email protected]

References

Utility Function

Description

This function provides the value of the likelihood function of the full (unrestricted) model, as described in Consonni et al. (2011).

Usage

lik.multin(y, p11, p21)
lik.multin(y, p11, p21)

Arguments

`y`	an object of class "HWEdata".
`p11`	gentotype proportion for the alleles pair $A_{1}A_{1}$ .
`p21`	gentotype proportion for the alleles pair $A_{2}A_{1}$ .

Details

This function has been used to generate the likelihood contours that appear in Figure 4 of Consonni et al. (2011) (see the example below).

Value

This function returns the numerical value of the likelihood in correspondence of the argument values.

Note

This function provides the likelihood function value only for the two alleles case.

Author(s)

Sergio Venturini [email protected]

References

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([email protected], 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] <- [email protected]
tables[, 2] <- [email protected]
tables[, 3] <- [email protected]
tables[, 4] <- [email protected]
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)
# 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)

Lindley (1988) Two Alleles Data

Description

Four samples of $n = 100$ genotype counts previously discussed in previously analyzed by Lindley (1988). For the first three sets, the classical "exact" test rejects the null hypothesis of Hardy-Weinberg equilibrium with significance level below 3.4%, whereas for the last data set the Hardy-Weinberg model is not rejected, its p-value being around 20%.

Usage

data(Lindley)data(Lindley)

Format

Four objects of class HWEdata.

Source

Consonni, G., Gutierrez-Pena, E. and Veronese, P. (2008), "Compatible priors for Bayesian model comparison with an application to the Hardy-Weinberg equilibrium model". Test, Vol. 17, No. 3, 585–605.

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 Guo, S.W. and Thompson, E.A. (1992), "Performing the Exact Test of Hardy-Weinberg Proportion for Multiple Alleles". Biometrics, 49, 361–372. Lindley D.V. (1988), "Statistical inference concerning Hardy-Weinberg equilibrium". In: Bernardo, J.M., DeGroot, M.H., Lindley, D.V. and Smith, A.F.M. (eds.), "Bayesian statistics 3". Oxford University Press, 307–326.

Examples

## Not run: 
# ATTENTION: the following code may take a long time to run! #

data(Lindley)
hwe.ibf.exact <- Vectorize(hwe.ibf, "t")
f <- seq(.05, 1, .05)
n <- sum([email protected], na.rm = TRUE)

# Dataset 1 #
plot(dataL1)
npp.exact <- 1/(1 + hwe.ibf.exact(round(f*n), y = dataL1))
npp.std <- 1/(1 + hwe.bf(dataL1))
plot(f, npp.exact, type="l", lwd = 2, xlab = "f = t/n",
	ylab = "Null posterior probability")
abline(h = npp.std, col = gray(.5), lty = "longdash")

# Dataset 2 #
plot(dataL2)
npp.exact <- 1/(1 + hwe.ibf.exact(round(f*n), y = dataL2))
npp.std <- 1/(1 + hwe.bf(dataL2))
plot(f, npp.exact, type="l", lwd = 2, xlab = "f = t/n",
	ylab = "Null posterior probability")
abline(h = npp.std, col = gray(.5), lty = "longdash")

# Dataset 3 #
plot(dataL3)
npp.exact <- 1/(1 + hwe.ibf.exact(round(f*n), y = dataL3))
npp.std <- 1/(1 + hwe.bf(dataL3))
plot(f, npp.exact, type="l", lwd = 2, xlab = "f = t/n",
	ylab = "Null posterior probability")
abline(h = npp.std, col = gray(.5), lty = "longdash")

# Dataset 4 #
plot(dataL4)
npp.exact <- 1/(1 + hwe.ibf.exact(round(f*n), y = dataL4))
npp.std <- 1/(1 + hwe.bf(dataL4))
plot(f, npp.exact, type="l", lwd = 2, xlab = "f = t/n",
	ylab = "Null posterior probability")
abline(h = npp.std, col = gray(.5), lty = "longdash")

## End(Not run)
## Not run: 
# ATTENTION: the following code may take a long time to run! #

data(Lindley)
hwe.ibf.exact <- Vectorize(hwe.ibf, "t")
f <- seq(.05, 1, .05)
n <- sum(dataL1@data.vec, na.rm = TRUE)

# Dataset 1 #
plot(dataL1)
npp.exact <- 1/(1 + hwe.ibf.exact(round(f*n), y = dataL1))
npp.std <- 1/(1 + hwe.bf(dataL1))
plot(f, npp.exact, type="l", lwd = 2, xlab = "f = t/n",
	ylab = "Null posterior probability")
abline(h = npp.std, col = gray(.5), lty = "longdash")

# Dataset 2 #
plot(dataL2)
npp.exact <- 1/(1 + hwe.ibf.exact(round(f*n), y = dataL2))
npp.std <- 1/(1 + hwe.bf(dataL2))
plot(f, npp.exact, type="l", lwd = 2, xlab = "f = t/n",
	ylab = "Null posterior probability")
abline(h = npp.std, col = gray(.5), lty = "longdash")

# Dataset 3 #
plot(dataL3)
npp.exact <- 1/(1 + hwe.ibf.exact(round(f*n), y = dataL3))
npp.std <- 1/(1 + hwe.bf(dataL3))
plot(f, npp.exact, type="l", lwd = 2, xlab = "f = t/n",
	ylab = "Null posterior probability")
abline(h = npp.std, col = gray(.5), lty = "longdash")

# Dataset 4 #
plot(dataL4)
npp.exact <- 1/(1 + hwe.ibf.exact(round(f*n), y = dataL4))
npp.std <- 1/(1 + hwe.bf(dataL4))
plot(f, npp.exact, type="l", lwd = 2, xlab = "f = t/n",
	ylab = "Null posterior probability")
abline(h = npp.std, col = gray(.5), lty = "longdash")

## End(Not run)

Louis and Dempster (1987) Four Alleles Data

Description

Sample of $n = 45$ genotype counts previously discussed in Louis and Dempster (1987) and Guo and Thompson (1992, Example 1). These data are described in Thomson et al. (1986) and concern the antigen class of 45 French type 1 diabetes patients, with the classes being DR1, DR3, DR4, and Y, a fourth class corresponding to all other antigens.

Usage

data(LouisDempster)data(LouisDempster)

Format

An object of class HWEdata.

Source

Louis, E. and Dempster, E. (1987), "An Exact Test for Hardy-Weinberg and Multiple Alleles". Biometrics Vol. 43, No. 4, 805–811.

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 Guo, S.W. and Thompson, E.A. (1992), "Performing the Exact Test of Hardy-Weinberg Proportion for Multiple Alleles". Biometrics, Vol. 49, No. 2, 361–372. Thomson, G., Klitz, W., Louis, E., Lo, S., Bertrams, L., Baur, M., and Neugebauer, M. (1986), "HLA and IDDM predisposition: New aspects". Genetic Epidemiology, Vol. 1, No. 2, 363–368.

Examples

# Example 1 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #


data(LouisDempster)
plot(LouisDempster)
n <- sum([email protected], na.rm = TRUE)
out <- hwe.ibf.mc(LouisDempster, t = n/2, M = 100000, verbose = TRUE)
summary(out, plot = TRUE)

## End(Not run)

# Example 2 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #

M <- 300000
f <- seq(.1, 1, .05)
n <- sum([email protected], na.rm = TRUE)
out <- hwe.ibf.plot(y = LouisDempster, t.vec = round(f*n), M = M)

## End(Not run)
# Example 1 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #


data(LouisDempster)
plot(LouisDempster)
n <- sum(LouisDempster@data.vec, na.rm = TRUE)
out <- hwe.ibf.mc(LouisDempster, t = n/2, M = 100000, verbose = TRUE)
summary(out, plot = TRUE)

## End(Not run)

# Example 2 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #

M <- 300000
f <- seq(.1, 1, .05)
n <- sum(LouisDempster@data.vec, na.rm = TRUE)
out <- hwe.ibf.plot(y = LouisDempster, t.vec = round(f*n), M = M)

## End(Not run)

Plot of Hardy-Weinberg Data and Analysis

Description

Methods for function plot in Package ‘graphics’ to be used with "HWEdata" and "HWEintr" objects.

Methods

signature(x = "HWEdata", y = "missing"): Pictorial representation for a "HWEdata" object.
signature(x = "HWEintr", y = "missing"): Graphical representation of Monte Carlo sums for a "HWEintr" object.

Author(s)

Sergio Venturini [email protected]

References

Two Alleles Simulated Data

Description

Four different samples of $n = 20$ genotype counts simulated under the Hardy-Weinberg equilibrium model.

Usage

data(simdata)data(simdata)

Format

Four objects of class HWEdata.

References

Examples

data(simdata)
summary(dataset1)
plot(dataset1)
summary(dataset2)
plot(dataset2)
summary(dataset3)
plot(dataset3)
summary(dataset4)
plot(dataset4)

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

n <- sum([email protected], 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] <- [email protected]
tables[, 2] <- [email protected]
tables[, 3] <- [email protected]
tables[, 4] <- [email protected]
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)
data(simdata)
summary(dataset1)
plot(dataset1)
summary(dataset2)
plot(dataset2)
summary(dataset3)
plot(dataset3)
summary(dataset4)
plot(dataset4)

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

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)

Summary of Hardy-Weinberg Data and Analysis

Description

Methods for function summary in Package ‘base’ to be used with "HWEdata" and "HWEintr" objects.

Methods

signature(object = "HWEdata"): Extracts the slots of a "HWEdata" object.
signature(object = "HWEintr"): Extracts the slots of a "HWEintr" object.

Author(s)

Sergio Venturini [email protected]

References

Wordsworth et al. (1992) Four Alleles Data

Description

Sample of $n = 230$ genotype counts discussed in Lauretto et al. (2009, Example 3). These data come from a rheumatoid arthritis (RA) study performed by Wordsworth et al. (1992), where two hundred and thirty RA patients were genotyped for the HLA-DR locus. The DR4 allele was subdivided into Dw4, Dw14 and other subtypes. DRX represents all non-DR1, non-Dw4, non-Dw14 alleles.

Usage

data(Wordsworth)data(Wordsworth)

Format

An object of class HWEdata.

Source

Lauretto, M.S., Nakano, F., Faria, S.R., Pereira, C.A.B. and Stern, J.M. (2009), "A straightforward multiallelic significance test for the Hardy-Weinberg equilibrium law". Genetics and Molecular Biology, Vol. 32, No. 3, 619–625.

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 Wordsworth, P., Pile, K.D., Buckley, J.D., Lanchbury, J.S.S., Ollier, B., Lathrop, M. and Bell, J.I. (1992), "HLA heterozygosity contributes to susceptibility to rheumatoid arthritis". American Journal of Human Genetics, 51, 3, 585–591.

Examples

# Example 1 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #


data(Wordsworth)
plot(Wordsworth)
n <- sum([email protected], na.rm = TRUE)
out <- hwe.ibf.mc(Wordsworth, t = n/2, M = 100000, verbose = TRUE)
summary(out, plot = TRUE)

## End(Not run)

# Example 2 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #

data(Wordsworth)
n <- sum([email protected], na.rm = TRUE)
M <- 300000
f <- seq(.1, 1, .05)
out <- hwe.ibf.plot(y = Wordsworth, t.vec = round(f*n), M = M)

## End(Not run)
# Example 1 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #


data(Wordsworth)
plot(Wordsworth)
n <- sum(Wordsworth@data.vec, na.rm = TRUE)
out <- hwe.ibf.mc(Wordsworth, t = n/2, M = 100000, verbose = TRUE)
summary(out, plot = TRUE)

## End(Not run)

# Example 2 #
## Not run: 
# ATTENTION: the following code may take a long time to run! #

data(Wordsworth)
n <- sum(Wordsworth@data.vec, na.rm = TRUE)
M <- 300000
f <- seq(.1, 1, .05)
out <- hwe.ibf.plot(y = Wordsworth, t.vec = round(f*n), M = M)

## End(Not run)

Package 'HWEintrinsic'

Help Index

Objective Bayesian Testing for the Hardy-Weinberg Equilibrium Problem

Description

Details

Author(s)

References

See Also

Graph of the Conditional Intrinsic Prior for a Two Alleles Dataset

Description

Usage

Arguments

Value

Note

Author(s)

References

See Also

Examples

Utility Function

Description

Usage

Arguments

Value

Author(s)

References

See Also

Guo and Thompson (1992) Eight Alleles Simulated Data

Description

Usage

Format

Source

References

Examples

Guo and Thompson (1992) Nine Alleles Data

Description

Usage

Format

Source

References

Examples

Standard Bayes Factor for the Hardy-Weinebrg Testing Problem

Description

Usage

Arguments

Value

Note

Author(s)

References

See Also

Examples

Testing Hardy-Weinberg Equilibrium Using an Intrinsic Prior Approach

Description

Usage

Arguments

Details

Value

Note

Author(s)

References

See Also

Examples

Testing Hardy-Weinberg Equilibrium Using an Intrinsic Prior Approach

Description

Usage

Arguments

Details

Value

Note

Author(s)

References

See Also

Examples

Plot of an Hardy-Weinberg Testing Analysis

Description

Usage

Arguments

Details

Value

Note

Author(s)