Two Alleles Simulated Data

Description

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

Usage

data(simdata)

Format

Four objects of class HWEdata.

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

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(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)

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(HWEintrinsic)
Package HWEintrinsic (1.2.2) loaded.
To cite, see citation("HWEintrinsic")

> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/HWEintrinsic/simdata.Rd_%03d_medium.png", width=480, height=480)
> ### Name: simdata
> ### Title: Two Alleles Simulated Data
> ### Aliases: simdata dataset1 dataset2 dataset3 dataset4
> ### Keywords: datasets HWE
> 
> ### ** Examples
> 
> data(simdata)
> summary(dataset1)
$n
[1] 20

$size
[1] 2

$data
     [,1] [,2]
[1,]    3   NA
[2,]    9    8

> plot(dataset1)
> summary(dataset2)
$n
[1] 20

$size
[1] 2

$data
     [,1] [,2]
[1,]    8   NA
[2,]    2   10

> plot(dataset2)
> summary(dataset3)
$n
[1] 20

$size
[1] 2

$data
     [,1] [,2]
[1,]   12   NA
[2,]    5    3

> plot(dataset3)
> summary(dataset4)
$n
[1] 20

$size
[1] 2

$data
     [,1] [,2]
[1,]    2   NA
[2,]   13    5

> plot(dataset4)
> 
> # The following code reproduces Figure 4 in Consonni et al. (2011) #
> ## Not run: 
> ##D # ATTENTION: it may take a long time to run! #
> ##D 
> ##D n <- sum(dataset1@data.vec, na.rm = TRUE)
> ##D f <- c(.1,.5,1)
> ##D t <- round(f*n)
> ##D p11 <- p21 <- seq(0,1,length.out=100)
> ##D ip <- array(NA,c(length(f),length(p11),length(p21)))
> ##D for (k in 1:length(f)) {
> ##D 	ip[k,,] <- outer(X = p11, Y = p21, FUN = Vectorize(ip.tmp), t[k])
> ##D 	print(paste(k," / ",length(f),sep=""), quote = FALSE)
> ##D }
> ##D 
> ##D r <- 2
> ##D R <- r*(r + 1)/2
> ##D l <- 4
> ##D tables <- matrix(NA, nrow = R, ncol = l)
> ##D tables[, 1] <- dataset1@data.vec
> ##D tables[, 2] <- dataset2@data.vec
> ##D tables[, 3] <- dataset3@data.vec
> ##D tables[, 4] <- dataset4@data.vec
> ##D lik <- array(NA, c(l, length(p11), length(p21)))
> ##D M <- 300000
> ##D par(mfrow = c(4, 4))
> ##D for (k in 1:l) {
> ##D 	y <- new("HWEdata", data = tables[, k])
> ##D 	lik[k,,] <- lik.multin(y, p11, p21)
> ##D 	
> ##D 	nlev <- 10
> ##D 	for (q in 1:length(f)) {
> ##D 		contour(p11, p21, ip[q,,], xlab = expression(p[11]),
> ##D 				ylab = expression(p[21]), nlevels = nlev, col = gray(0),
> ##D 				main = "", cex.axis = 1.75, cex.lab = 1.75, labcex = 1.4)
> ##D 		lines(p11^2, 2*p11*(1 - p11), lty = "longdash", col = gray(0), lwd = 2)
> ##D 		contour(p11, p21, lik[k,,], nlevels = nlev, add = TRUE,
> ##D 				col = gray(.7), labcex = 1.2)
> ##D 		abline(a = 1, b = -1, lty = 3, col = gray(.8))
> ##D 	}
> ##D 	hwe.ibf.plot(y = y, t.vec = seq(1,n,1), M = M)
> ##D }
> ## End(Not run)
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>