### R code from vignette source 'HardyWeinberg.Rnw'

###################################################
### code chunk number 1: HardyWeinberg.Rnw:604-605
###################################################
options(prompt = "R> ", continue = "+ ", width = 70, useFancyQuotes = FALSE)


###################################################
### code chunk number 2: HardyWeinberg.Rnw:608-610 (eval = FALSE)
###################################################
## install.packages("HardyWeinberg")
## library("HardyWeinberg")


###################################################
### code chunk number 3: HardyWeinberg.Rnw:618-619 (eval = FALSE)
###################################################
## vignette("HardyWeinberg")


###################################################
### code chunk number 4: HardyWeinberg.Rnw:629-631
###################################################
library("HardyWeinberg")
x <- c(MM = 298, MN = 489, NN = 213)


###################################################
### code chunk number 5: HardyWeinberg.Rnw:638-639
###################################################
HW.test <- HWChisq(x, verbose = TRUE)


###################################################
### code chunk number 6: HardyWeinberg.Rnw:645-646
###################################################
HW.test <- HWChisq(x, cc = 0, verbose = TRUE)


###################################################
### code chunk number 7: HardyWeinberg.Rnw:656-657
###################################################
HW.lrtest <- HWLratio(x, verbose = TRUE)


###################################################
### code chunk number 8: HardyWeinberg.Rnw:668-669
###################################################
HW.exacttest <- HWExact(x, verbose = TRUE)


###################################################
### code chunk number 9: HardyWeinberg.Rnw:679-681
###################################################
set.seed(123)
HW.permutationtest <- HWPerm(x, verbose = TRUE)


###################################################
### code chunk number 10: HardyWeinberg.Rnw:688-690
###################################################
x <- c(MN = 489, NN = 213, MM = 298)
HW.test <- HWChisq(x, verbose = TRUE)


###################################################
### code chunk number 11: HardyWeinberg.Rnw:698-699
###################################################
HW.results <- HWAlltests(x, verbose = TRUE, include.permutation.test = TRUE)


###################################################
### code chunk number 12: HardyWeinberg.Rnw:708-710
###################################################
data(Markers)
Markers[1:12,]


###################################################
### code chunk number 13: HardyWeinberg.Rnw:716-720
###################################################
Xt <- table(Markers[,1])
Xv <- as.vector(Xt)
names(Xv) <- names(Xt)
HW.test <- HWChisq(Xv,cc=0,verbose=TRUE)


###################################################
### code chunk number 14: HardyWeinberg.Rnw:727-729
###################################################
set.seed(123)
Results <- HWMissing(Markers[,1], m = 50, method = "sample", verbose=TRUE)


###################################################
### code chunk number 15: HardyWeinberg.Rnw:736-738
###################################################
set.seed(123)
Results <- HWMissing(Markers[, 1:5], m = 50, verbose = TRUE)


###################################################
### code chunk number 16: HardyWeinberg.Rnw:744-746
###################################################
set.seed(123)
Results <- HWMissing(Markers[, 1:5], m = 50, statistic = "exact", verbose = TRUE)


###################################################
### code chunk number 17: HardyWeinberg.Rnw:755-757
###################################################
data(JPTsnps)
Results <- HWPosterior(males=JPTsnps[1,1:3],females=JPTsnps[1,4:6],x.linked=FALSE,precision=0.05)


###################################################
### code chunk number 18: HardyWeinberg.Rnw:764-767
###################################################
data(JPTsnps)
AICs <- HWAIC(JPTsnps[1,1:3],JPTsnps[1,4:6])
AICs


###################################################
### code chunk number 19: HardyWeinberg.Rnw:777-780
###################################################
g1 <- c(0.034, 0.330, 0.636)
g2 <- c(0.349, 0.493, 0.158)
x  <- c(0.270, 0.453,0.277)


###################################################
### code chunk number 20: HardyWeinberg.Rnw:785-788
###################################################
G <- cbind(g1,g2)
contributions <- HWEM(x,G=G)
contributions


###################################################
### code chunk number 21: HardyWeinberg.Rnw:793-796
###################################################
p <- c(af(g1),af(g2))
contributions <- HWEM(x,p=p)
contributions


###################################################
### code chunk number 22: HardyWeinberg.Rnw:805-807
###################################################
SNP1 <- c(A=399,B=205,AA=230,AB=314,BB=107) 
HWChisq(SNP1,cc=0,x.linked=TRUE,verbose=TRUE)


###################################################
### code chunk number 23: HardyWeinberg.Rnw:812-813
###################################################
HWChisq(SNP1[3:5],cc=0)


###################################################
### code chunk number 24: HardyWeinberg.Rnw:821-822
###################################################
HWExact(SNP1,x.linked=TRUE)


###################################################
### code chunk number 25: HardyWeinberg.Rnw:827-828
###################################################
HWExact(SNP1,x.linked=TRUE,pvaluetype="midp")


###################################################
### code chunk number 26: HardyWeinberg.Rnw:834-835
###################################################
HWExact(SNP1[3:5])


###################################################
### code chunk number 27: HardyWeinberg.Rnw:840-841
###################################################
HWPerm(SNP1,x.linked=TRUE)


###################################################
### code chunk number 28: HardyWeinberg.Rnw:846-847
###################################################
HWLratio(SNP1,x.linked=TRUE)


###################################################
### code chunk number 29: HardyWeinberg.Rnw:852-853
###################################################
HWAlltests(SNP1,x.linked=TRUE,include.permutation.test=TRUE)


###################################################
### code chunk number 30: HardyWeinberg.Rnw:858-859
###################################################
AFtest(SNP1)


###################################################
### code chunk number 31: HardyWeinberg.Rnw:869-870
###################################################
HWPosterior(males=SNP1[1:2],females=SNP1[3:5],x.linked=TRUE)


###################################################
### code chunk number 32: HardyWeinberg.Rnw:889-898
###################################################
x <- c(MM = 298, MN = 489, NN = 213)
n <- sum(x)
nM <- mac(x) 
pw4 <- HWPower(n, nM, alpha = 0.05, test = "exact", theta = 4, 
               pvaluetype = "selome")
print(pw4)
pw8 <- HWPower(n, nM, alpha = 0.05, test = "exact", theta = 8, 
               pvaluetype = "selome")
print(pw8)


###################################################
### code chunk number 33: HardyWeinberg.Rnw:917-940
###################################################
set.seed(123)
n <- 100
m <- 100
X1 <- HWData(m, n, p = rep(0.5, m))
X2 <- HWData(m, n)
X3 <- HWData(m, n, p = rep(0.5, m), f = rep(0.5, m))
X4 <- HWData(m, n, f = rep(0.5, m))
X5 <- HWData(m, n, p = rep(c(0.2, 0.4, 0.6, 0.8), 25), conditional = TRUE)
X6 <- HWData(m, n, exactequilibrium = TRUE)
opar <- par(mfrow = c(3, 2),mar = c(1, 0, 3, 0) + 0.1)
par(mfg = c(1, 1))
HWTernaryPlot(X1, main = "(a)")
par(mfg = c(1, 2))
HWTernaryPlot(X2, main = "(b)")
par(mfg = c(2, 1))
HWTernaryPlot(X3, main = "(c)")
par(mfg = c(2, 2))
HWTernaryPlot(X4, main = "(d)")
par(mfg = c(3, 1))
HWTernaryPlot(X5, main = "(e)")
par(mfg = c(3, 2))
HWTernaryPlot(X6, main = "(f)")
par(opar)


###################################################
### code chunk number 34: HardyWeinberg.Rnw:946-947
###################################################
X <- HWData(nm=100,shape1=1,shape2=10)


###################################################
### code chunk number 35: HardyWeinberg.Rnw:955-958 (eval = FALSE)
###################################################
## data("HapMapCHBChr1", package = "HardyWeinberg")
## HWTernaryPlot(HapMapCHBChr1, region = 1)
## HWTernaryPlot(HapMapCHBChr1, region = 7)


###################################################
### code chunk number 36: HardyWeinberg.Rnw:974-981 (eval = FALSE)
###################################################
## set.seed(123)
## data("HapMapCHBChr1", package = "HardyWeinberg")
## HWQqplot(HapMapCHBChr1)
## dev.off()
## set.seed(123)
## SimulatedData <- HWData(nm = 225, n = 84, p = af(HapMapCHBChr1))
## HWQqplot(SimulatedData)


###################################################
### code chunk number 37: HardyWeinberg.Rnw:1004-1006
###################################################
x <- c(fA=182,fB=60,nAB=17,nOO=176)
al.fre <- HWABO(x)


###################################################
### code chunk number 38: HardyWeinberg.Rnw:1015-1017
###################################################
x <- c(AA=20,AB=31,AC=26,BB=15,BC=12,CC=0)
results <- HWTriExact(x)


###################################################
### code chunk number 39: HardyWeinberg.Rnw:1022-1026
###################################################
x <- c(AA=20,AB=31,AC=26,BB=15,BC=12,CC=0)
x <- toTriangular(x)
m <- c(A=0,B=0,C=0)
#results <- HWNetwork(ma=m,fe=x)


###################################################
### code chunk number 40: HardyWeinberg.Rnw:1031-1034
###################################################
males   <- c(A=1,B=21,C=34) 
females <- c(AA=0,AB=1,AC=0,BB=8,BC=24,CC=15)
results <- HWTriExact(females,males)


###################################################
### code chunk number 41: HardyWeinberg.Rnw:1039-1042
###################################################
males   <- c(A=1,B=21,C=34) 
females <- toTriangular(c(AA=0,AB=1,AC=0,BB=8,BC=24,CC=15))
#results <- HWNetwork(ma=males,fe=females)


###################################################
### code chunk number 42: HardyWeinberg.Rnw:1047-1051
###################################################
set.seed(123)
x <- c(AA=20,AB=31,AC=26,BB=15,BC=12,CC=0)
x <- toTriangular(x)
#results <- HWPerm.mult(x)


###################################################
### code chunk number 43: HardyWeinberg.Rnw:1058-1065
###################################################
set.seed(123)
data(NistSTRs)
A1 <- NistSTRs[,1]
A2 <- NistSTRs[,2]
GenotypeCounts <- AllelesToTriangular(A1,A2)
print(GenotypeCounts)
#out <- HWPerm.mult(GenotypeCounts)


###################################################
### code chunk number 44: HardyWeinberg.Rnw:1070-1071
###################################################
#Results <- HWStr(NistSTRs,test="permutation")