### R code from vignette source 'vignettes/qpgraph/inst/doc/qpgraphSimulate.Rnw'

###################################################
### code chunk number 1: setup
###################################################
pdf.options(useDingbats=FALSE)
options(width=80)
rm(list=ls())
try(detach("package:mvtnorm"), silent=TRUE)
try(detach("package:qtl"), silent=TRUE)


###################################################
### code chunk number 2: qpgraphSimulate.Rnw:114-120
###################################################
library(qpgraph)

args(erGraphParam)
args(erMarkedGraphParam)
args(dRegularGraphParam)
args(dRegularMarkedGraphParam)


###################################################
### code chunk number 3: qpgraphSimulate.Rnw:125-129
###################################################
erGraphParam()
erMarkedGraphParam()
dRegularGraphParam()
dRegularMarkedGraphParam()


###################################################
### code chunk number 4: qpgraphSimulate.Rnw:134-136
###################################################
showClass("graphParam")
showClass("markedGraphParam")


###################################################
### code chunk number 5: qpgraphSimulate.Rnw:147-148
###################################################
args(rgraphBAM)


###################################################
### code chunk number 6: qpgraphSimulate.Rnw:154-156
###################################################
rgraphBAM(erGraphParam())
rgraphBAM(n=2, dRegularGraphParam())


###################################################
### code chunk number 7: 3regulargraph
###################################################
set.seed(1234)
g <- rgraphBAM(dRegularMarkedGraphParam(pI=2, pY=10, d=3))
plot(g)


###################################################
### code chunk number 8: qpgraphSimulate.Rnw:194-195
###################################################
args(rUGgmm)


###################################################
### code chunk number 9: qpgraphSimulate.Rnw:206-215
###################################################
rUGgmm(dRegularGraphParam(p=4, d=2))
rUGgmm(matrix(c(0, 1, 0, 1,
                1, 0, 1, 0,
                0, 1, 0, 1,
                1, 0, 1, 0), nrow=4, byrow=TRUE))
rUGgmm(matrix(c(1, 2,
                2, 3,
                3, 4,
                4, 1), ncol=2, byrow=TRUE))


###################################################
### code chunk number 10: qpgraphSimulate.Rnw:222-225
###################################################
set.seed(12345)
gmm <- rUGgmm(dRegularGraphParam(p=4, d=2))
summary(gmm)


###################################################
### code chunk number 11: qpgraphSimulate.Rnw:229-236
###################################################
class(gmm)
names(gmm)
gmm$X
gmm$p
gmm$g
gmm$mean
gmm$sigma


###################################################
### code chunk number 12: 4cyclegraph
###################################################
plot(gmm)
round(solve(gmm$sigma), digits=1)


###################################################
### code chunk number 13: qpgraphSimulate.Rnw:271-272
###################################################
rmvnorm(10, gmm)


###################################################
### code chunk number 14: qpgraphSimulate.Rnw:277-281
###################################################
set.seed(123)
X <- rmvnorm(10000, gmm)
round(solve(cov(X)), digits=1)
round(solve(gmm$sigma), digits=1)


###################################################
### code chunk number 15: qpgraphSimulate.Rnw:345-346
###################################################
args(rHMgmm)


###################################################
### code chunk number 16: qpgraphSimulate.Rnw:355-364
###################################################
rHMgmm(dRegularMarkedGraphParam(pI=1, pY=3, d=2))
rHMgmm(matrix(c(0, 1, 0, 1,
                1, 0, 1, 0,
                0, 1, 0, 1,
                1, 0, 1, 0), nrow=4, byrow=TRUE), I=1)
rHMgmm(matrix(c(1, 2,
                2, 3,
                3, 4,
                4, 1), ncol=2, byrow=TRUE), I=1)


###################################################
### code chunk number 17: qpgraphSimulate.Rnw:374-377
###################################################
set.seed(12345)
gmm <- rHMgmm(dRegularMarkedGraphParam(pI=1, pY=3, d=2))
summary(gmm)


###################################################
### code chunk number 18: qpgraphSimulate.Rnw:382-395
###################################################
class(gmm)
names(gmm)
gmm$X
gmm$I
gmm$Y
gmm$p
gmm$pI
gmm$pY
gmm$g
gmm$mean()
gmm$sigma
gmm$a
gmm$eta2


###################################################
### code chunk number 19: hmgmm
###################################################
plot(gmm)


###################################################
### code chunk number 20: qpgraphSimulate.Rnw:421-422
###################################################
rcmvnorm(10, gmm)


###################################################
### code chunk number 21: qpgraphSimulate.Rnw:427-433
###################################################
set.seed(123)
X <- rcmvnorm(10000, gmm)
csigma <- (1/10000)*sum(X[, gmm$I] == 1)*cov(X[X[, gmm$I]==1, gmm$Y]) +
          (1/10000)*sum(X[, gmm$I] == 2)*cov(X[X[, gmm$I]==2, gmm$Y])
round(solve(csigma), digits=1)
round(solve(gmm$sigma), digits=1)


###################################################
### code chunk number 22: qpgraphSimulate.Rnw:438-443
###################################################
smean <- apply(X[, gmm$Y], 2, function(x, i) tapply(x, i, mean), X[, gmm$I])
smean
gmm$mean()
abs(smean[1, ] - smean[2, ])
gmm$a


###################################################
### code chunk number 23: qpgraphSimulate.Rnw:467-469
###################################################
eQTLcrossParam()
args(eQTLcrossParam)


###################################################
### code chunk number 24: qpgraphSimulate.Rnw:475-476
###################################################
reQTLcross(n=2, eQTLcrossParam())


###################################################
### code chunk number 25: qpgraphSimulate.Rnw:489-492
###################################################
data(map10, package="qtl")
class(map10)
head(lapply(map10, head))


###################################################
### code chunk number 26: qpgraphSimulate.Rnw:497-500
###################################################
eqtl <- reQTLcross(eQTLcrossParam(map=map10, genes=100))
class(eqtl)
eqtl


###################################################
### code chunk number 27: qpgraphSimulate.Rnw:507-508
###################################################
eqtl$model


###################################################
### code chunk number 28: geneticmap
###################################################
par(mfrow=c(1,2))
plot(map10)
plot(eqtl, main="eQTL model with cis-associations only")


###################################################
### code chunk number 29: qpgraphSimulate.Rnw:532-535
###################################################
set.seed(123)
eqtl <- reQTLcross(eQTLcrossParam(map=map10, genes=100,
                                  cis=0.5, trans=c(5, 5)), a=5)


###################################################
### code chunk number 30: qpgraphSimulate.Rnw:545-547
###################################################
head(ciseQTL(eqtl))
transeQTL(eqtl)


###################################################
### code chunk number 31: transeqtl
###################################################
plot(eqtl, main="eQTL model with trans-eQTL")


###################################################
### code chunk number 32: qpgraphSimulate.Rnw:567-570
###################################################
set.seed(123)
cross <- sim.cross(map10, eqtl)
cross


###################################################
### code chunk number 33: qpgraphSimulate.Rnw:584-589
###################################################
allcis <- ciseQTL(eqtl)
allcis[allcis$chrom==1, ]
gene <- allcis[2, "gene"]
chrom <- allcis[2, "chrom"]
location <- allcis[2, "location"]


###################################################
### code chunk number 34: qpgraphSimulate.Rnw:592-595
###################################################
connectedGenes <- graph::inEdges(gene, eqtl$model$g)[[1]]
connectedGenes <- connectedGenes[connectedGenes %in% eqtl$model$Y]
connectedGenes


###################################################
### code chunk number 35: qpgraphSimulate.Rnw:606-609
###################################################
library(qtl)

out.mr <- scanone(cross, method="mr", pheno.col=c(gene, connectedGenes))


###################################################
### code chunk number 36: lodprofiles
###################################################
plot(out.mr, chr=chrom, ylab="LOD score", lodcolumn=1:3)
abline(v=allcis[allcis$gene == gene, "location"])


###################################################
### code chunk number 37: qpgraphSimulate.Rnw:630-633
###################################################
out.perm <- scanone(cross, method="mr", pheno.col=c(gene, connectedGenes),
                    n.perm=100, verbose=FALSE)
summary(out.perm, alpha=c(0.05, 0.10))


###################################################
### code chunk number 38: qpgraphSimulate.Rnw:638-639
###################################################
summary(out.mr, perms=out.perm, alpha=0.05)


###################################################
### code chunk number 39: qpgraphSimulate.Rnw:649-652
###################################################
allcis[allcis$gene %in% connectedGenes, ]
alltrans <- transeQTL(eqtl)
alltrans[alltrans$gene %in% connectedGenes, ]


###################################################
### code chunk number 40: info
###################################################
toLatex(sessionInfo())