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