[貝葉斯專著]Bayesian Networks in R: with Applications in Systems Biology

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

Bayesian Networks in R: with Applications in Systems Biology

Radhakrishnan Nagarajan (Author), Marco Scutari (Author), Sophie Lèbre (Author)

Book Description[size=0.86em]Publication Date: May 7 2013 | ISBN-10: 1461464455 | ISBN-13: 978-1461464457 | Edition: 2013
Bayesian Networks in R with Applications in Systems Biology is unique as it introduces the reader to the essential concepts in Bayesian network modeling and inference in conjunction with examples in the open-source statistical environment R. The level of sophistication is also gradually increased across the chapters with exercises and solutions for enhanced understanding for hands-on experimentation of the theory and concepts. The application focuses on systems biology with emphasis on modeling pathways and signaling mechanisms from high-throughput molecular data. Bayesian networks have proven to be especially useful abstractions in this regard. Their usefulness is especially exemplified by their ability to discover new associations in addition to validating known ones across the molecules of interest. It is also expected that the prevalence of publicly available high-throughput biological data sets may encourage the audience to explore investigating novel paradigms using the approaches presented in the book.

Product Details

• Paperback: 157 pages
• Publisher: Springer; 2013 edition (May 7 2013)
• Language: English
• ISBN-10: 1461464455
• ISBN-13: 978-1461464457
• Product Dimensions: 0.9 x 15.4 x 23 cm
• Shipping Weight: 222 g
  

本帖隐藏的内容

Bayesian Networks in R.with Applications in Systems Biology Springer 2013.pdf (3.21 MB, 需要: 100 个论坛币)

2014-7-14 02:29:40 上传

需要: 100 个论坛币  [购买]

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享1 收藏12 回帖

关键词：Applications Application Networks Bayesian Biology essential concepts examples network

1. install.packages("bnlearn")
   
2. library(bnlearn)
   
3. data(lizards)
   
4. lizards = read.table("lizards.txt", header = TRUE)
   
5. 
   
6. #-----------------------------------------------------------------------------#
   
7. 
   
8. str(lizards)
   
9. summary(lizards)
   
10. dim(lizards)
    
11. levels(lizards[, "Species"])
    
12. levels(lizards[, "Height"])
    
13. levels(lizards[, "Diameter"])
    
14. table(lizards[, c(3,2,1)])
    
15. 
    
16. #-----------------------------------------------------------------------------#
    
17. 
    
18. Sagrei.lizards = lizards[lizards$Species == "Sagrei", ]
    
19. Distichus.lizards = lizards[lizards$Species == "Distichus", ]
    
20. par(mfrow = c(2, 2))
    
21. plot(Sagrei.lizards[, "Height"], main = "Perch Height (Sagrei)")
    
22. plot(Distichus.lizards[, "Height"], main = "Perch Height (Distichus)")
    
23. plot(Sagrei.lizards[, "Diameter"], main = "Perch Diameter (Sagrei)")
    
24. plot(Distichus.lizards[, "Diameter"], main = "Perch Diameter (Distichus)")
    
25. 
    
26. #-----------------------------------------------------------------------------#
    
27. 
    
28. set.seed(42)
    
29. diam = numeric(length = nrow(lizards))
    
30. narrow = (lizards$Diameter == "narrow")
    
31. wide = (lizards$Diameter == "wide")
    
32. diam[narrow] = runif(n = 252, min = 2, max = 4)
    
33. diam[wide] = runif(n = 157, min = 4, max = 6)
    
34. new.data = data.frame(
    
35.              Species = lizards[, "Species"],
    
36.              Sim.Diameter = diam)
    
37. 
    
38. 
    
39. #-----------------------------------------------------------------------------#
    
40. 
    
41. summary(new.data[, "Sim.Diameter"])
    
42. is.sagrei = (new.data$Species == "Sagrei")
    
43. summary(new.data[is.sagrei, "Sim.Diameter"])
    
44. summary(new.data[!is.sagrei, "Sim.Diameter"])
    
45. var(new.data[is.sagrei, "Sim.Diameter"])
    
46. var(new.data[!is.sagrei, "Sim.Diameter"])
    
47. 
    
48. #-----------------------------------------------------------------------------#
    
49. 
    
50. boxplot(Sim.Diameter  ìƒ Species, data = new.data, ylab = "Diameter (inches)")
    
51. abline(h = 4, lty = "dashed")

复制代码

Very helpful, I have been looking for similar reference for a while.
THanks.

1. library(vars)
   
2. data(Canada)
   
3. VAR(Canada, p = 2)
   
4. summary(VAR(Canada, p = 2))
   
5. VAR(Canada, p = 2, type = "none")
   
6. VAR(Canada, p = 2, type = "trend")
   
7. VAR(Canada, p = 2, type = "both")
   
8. VAR(Canada, lag.max = 4, ic = "AIC")
   
9. VAR(Canada, lag.max = 4, ic = "SC")
   
10. var.2c = VAR(Canada, p = 2, type = "const")
    
11. stab = stability(var.2c, type = "OLS-CUSUM")
    
12. plot(stab)
    
13. normality.test(var.2c)
    
14. normality.test(var.2c, multivariate.only = FALSE)
    
15. serial.test(var.2c, lags.pt = 16, type = "PT.adjusted")
    
16. arch.test(var.2c)
    
17. arch.test(var.2c, multivariate.only = FALSE)
    
18. 
    
19. #-----------------------------------------------------------------------------#
    
20. 
    
21. library(lars)
    
22. library(GeneNet)
    
23. data(arth800)
    
24. subset = c(60, 141, 260, 333, 365, 424, 441, 512, 521, 578, 789, 799)
    
25. arth12 = arth800.expr[, subset]
    
26. x = arth12[1:(nrow(arth12) - 2), ]     
    
27. y = arth12[-(1:2), "265768_at"]      
    
28. lasso.fit = lars(y = y, x = x, type = "lasso")
    
29. fit.all = lapply(colnames(arth12),
    
30.    function(gene) {
    
31.      y = arth12[-(1:2), gene]
    
32.      lars(y = y, x = x, type = "lasso")
    
33.    })
    
34. plot(lasso.fit)
    
35. coef(lasso.fit)
    
36. lasso.cv = cv.lars(y = y, x = x, mode = "fraction")
    
37. frac = lasso.cv$index[which.min(lasso.cv$cv)]
    
38. predict(lasso.fit, s = frac, type = "coef", mode = "fraction")
    
39. predict(lasso.fit, s = 3, type = "coef", mode = "step")$coefficients
    
40. predict(lasso.fit, s = 0.2, type = "coef", mode = "lambda")$coefficients
    
41. lar.fit = lars(y = y, x = x, type = "lar")
    
42. plot(lar.fit)
    
43. lar.cv = cv.lars(y = y, x = x, type = "lar")
    
44. step.fit = lars(y = y, x = x, type = "stepwise")
    
45. plot(step.fit)
    
46. step.cv = cv.lars(y = y, x = x, type = "stepwise")     
    
47. 
    
48. #-----------------------------------------------------------------------------#
    
49. 
    
50. library(simone)
    
51. simone(arth12, type = "time-course")
    
52. ctrl = setOptions(clusters.crit = "BIC")
    
53. simone(arth12, type = "time-course", clustering = TRUE, control = ctrl)
    
54. plot(simone(arth12, type = "time-course",
    
55.        clustering = TRUE, control = ctrl), output = "BIC")
    
56. plot(simone(arth12, type = "time-course",
    
57.        clustering = TRUE, control = ctrl))
    
58. 
    
59. #-----------------------------------------------------------------------------#
    
60. 
    
61. library(GeneNet)
    
62. dyn = ggm.estimate.pcor(arth800.expr, method = "dynamic")
    
63. arth.arcs = network.test.edges(dyn)
    
64. arth.net = extract.network(arth.arcs, method.ggm = "number", cutoff.ggm = 10)
    
65. arth.net = extract.network(arth.arcs, method.ggm = "prob", cutoff.ggm = 0.05)
    
66. 
    
67. #-----------------------------------------------------------------------------#
    
68. 
    
69. library(G1DBN)
    
70. data(arth800line)
    
71. subset = c(60, 141, 260, 333, 365, 424, 441, 512, 521, 578, 789, 799)
    
72. arth12 = as.matrix(arth800line[, subset])
    
73. step1 = DBNScoreStep1(arth12, method = "ls")
    
74. round(step1$S1ls, 2)[1:6, 1:6]
    
75. edgesG1 = BuildEdges(score = step1$S1ls, threshold = 0.50,  prec = 6)
    
76. nrow(edgesG1)
    
77. step2 = DBNScoreStep2(step1$S1ls, data = arth12,
    
78.        method = "ls", alpha1 = 0.50)
    
79. edgesG = BuildEdges(score = step2, threshold = 0.05, prec = 6)
    
80. 
    
81. #-----------------------------------------------------------------------------#
    
82. 
    
83. library(ARTIVA)
    
84. data(simulatedProfiles)
    
85. targets = c("1", "10", "20", "TF3", "45", "50")
    
86. parents = c("TF1", "TF2", "TF3", "TF4", "TF5")
    
87. DBN = ARTIVAnet(
    
88.         targetData = simulatedProfiles[targets, ],
    
89.         parentData = simulatedProfiles[parents, ],
    
90.         targetNames = targets,
    
91.         parentNames = parents,
    
92.         niter = 50000,
    
93.         savePictures = FALSE)
    
94. head(ARTIVAtest1[, -7])

复制代码

1. isachs = read.table("sachs.interventional.txt",
   
2.            header = TRUE, colClasses = "factor")
   
3. library(gRain)
   
4. library(bnlearn)
   
5. val.str = paste("[PKC][PKA|PKC][praf|PKC:PKA]",
   
6.            "[pmek|PKC:PKA:praf][p44.42|pmek:PKA]",
   
7.            "[pakts473|p44.42:PKA][P38|PKC:PKA]",
   
8.            "[pjnk|PKC:PKA][plcg][PIP3|plcg]",
   
9.            "[PIP2|plcg:PIP3]")
   
10. val = model2network(val.str)
    
11. isachs = isachs[, 1:11]
    
12. for (i in names(isachs))
    
13.   levels(isachs[, i]) = c("LOW", "AVG", "HIGH")
    
14. fitted = bn.fit(val, isachs, method = "bayes")
    
15. jtree = compile(as.grain(fitted))
    
16. jprop = setFinding(jtree, nodes = "p44.42",
    
17.           states = "LOW")
    
18. querygrain(jtree, nodes = "pakts473")$pakts473
    
19. querygrain(jprop, nodes = "pakts473")$pakts473
    
20. querygrain(jtree, nodes = "PKA")$PKA
    
21. querygrain(jprop, nodes = "PKA")$PKA
    
22. names(which.max(querygrain(jprop,
    
23.                     nodes = c("PKA"))$PKA))
    
24. particles = cpdist(fitted, nodes = "pakts473",
    
25.               evidence = (p44.42 == "LOW"))
    
26. prop.table(table(particles))
    
27. particles = cpdist(fitted, nodes = "PKA",
    
28.               evidence = (p44.42 == "LOW"))
    
29. prop.table(table(particles))
    
30. cpquery(fitted,
    
31.   event = (pakts473 == "LOW") & (PKA != "HIGH"),
    
32.   evidence = (p44.42 == "LOW") | (praf == "LOW"))
    
33. 
    
34. #-----------------------------------------------------------------------------#
    
35. 
    
36. x = arth12[1:(nrow(arth12) - 2), ]
    
37. y = arth12[-(1:2), "265768_at"]
    
38. lasso.fit = lars(y = y, x = x, type = "lasso")
    
39. lasso.cv = cv.lars(y = y, x = x, mode = "fraction")
    
40. frac = lasso.cv$index[which.min(lasso.cv$cv)]
    
41. lasso.est = predict(lasso.fit, type = "fit",
    
42.               newx = x, s = frac,
    
43.               mode = "fraction")$fit
    
44. lasso.est
    
45. lasso.pred = predict(lasso.fit, type = "fit",
    
46.                newx = arth12[c("24-1", "24-2"), ],
    
47.                s = frac, mode = "fraction")$fit
    
48. lasso.pred
    
49. library(penalized)
    
50. lambda = optL1(response = y, penalized = x)$lambda
    
51. lasso.t = penalized(response = y, penalized = x,
    
52.             lambda1 = lambda)
    
53. coef(lasso.t)
    
54. dbn1 =
    
55.   model2network("[245094_at][265768_at|245094_at]")
    
56. xp.mean = mean(x[, "245094_at"])
    
57. xp.sd = sd(x[, "245094_at"])
    
58. dbn1.fit =
    
59.   custom.fit(dbn1,
    
60.     dist = list("245094_at" = list(coef = xp.mean,
    
61.              sd = xp.sd), "265768_at" = lasso.t))
    
62. cpquery(dbn1.fit, event = (`265768_at` > 8),
    
63.                   evidence = (`245094_at` > 8))
    
64. cpquery(dbn1.fit, event = (`265768_at` > 8),
    
65.                   evidence = (`245094_at` < 8))
    
66. dist.low = cpdist(dbn1.fit, node = "265768_at",
    
67.              evidence = (`245094_at` < 8))
    
68. dist.high = cpdist(dbn1.fit, node = "265768_at",
    
69.               evidence = (`245094_at` > 8))
    
70. y = arth12[-(1:2), "245094_at"]
    
71. colnames(x)[12] = "245094_at1"
    
72. lambda = optL1(response = y, penalized = x)$lambda
    
73. lasso.s = penalized(response = y, penalized = x,
    
74.              lambda1 = lambda)
    
75. coef(lasso.s)
    
76. dbn2 = empty.graph(c("265768_at", "245094_at",
    
77.         "258736_at", "257710_at", "255070_at",
    
78.         "245319_at", "245094_at1"))
    
79. dbn2 = set.arc(dbn2, "245094_at", "265768_at")
    
80. for (node in names(coef(lasso.s))[-c(1, 6)])
    
81.   dbn2 = set.arc(dbn2, node, "245094_at")
    
82. dbn2 = set.arc(dbn2, "245094_at1", "245094_at")
    
83. dbn2.data = as.data.frame(x[, nodes(dbn2)[1:6]])
    
84. dbn2.data[, "245094_at"] = y
    
85. dbn2.data[, "245094_at1"] = x[, "245094_at"]
    
86. dbn2.fit = bn.fit(dbn2, dbn2.data)
    
87. dbn2.fit[["265768_at"]] = lasso.t
    
88. dbn2.fit[["245094_at"]] = lasso.s
    
89. cpquery(dbn2.fit, event = (`265768_at` > 8),
    
90.   evidence = (`245094_at` > 8) & (`245094_at1` > 8))
    
91. cpquery(dbn2.fit, event = (`265768_at` > 8),
    
92. evidence = (`245094_at1` > 7) & (`245094_at1` < 8))
    
93. cpquery(dbn2.fit, event = (`265768_at` > 8),
    
94.   evidence = TRUE)
    
95. cpd =
    
96.   cpdist(dbn2.fit, node = "245094_at", evidence =
    
97.     (`245094_at1` > 6.5) & (`245094_at1` < 7.5) &
    
98.     (`265768_at` > 7) & (`265768_at` < 8))
    
99. summary(cpd)

复制代码

1. Examples
   
2. data(learning.test)
   
3. res = gs(learning.test)
   
4. cpdag(res)
   
5. vstructs(res)

复制代码

1. Examples
   
2. ## discrete Bayesian network (it is the same with ordinal nodes).
   
3. data(learning.test)
   
4. fitted = bn.fit(hc(learning.test), learning.test)
   
5. # the result should be around 0.025.
   
6. cpquery(fitted, (B == "b"), (A == "a"))
   
7. # for a single observation, predict the value of a single
   
8. # variable conditional on the others.
   
9. var = names(learning.test)
   
10. obs = 2
    
11. str = paste("(", names(learning.test)[-3], "=='",
    
12. sapply(learning.test[obs,-3], as.character), "')",
    
13. sep = "", collapse = " & ")
    
14. str
    
15. str2 = paste("(", names(learning.test)[3], "=='",
    
16. as.character(learning.test[obs, 3]), "')", sep = "")
    
17. str2
    
18. cpquery(fitted, eval(parse(text = str2)), eval(parse(text = str)))
    
19. # do the same with likelihood weighting
    
20. cpquery(fitted, event = eval(parse(text = str2)),
    
21. evidence = as.list(learning.test[2, -3]), method = "lw")
    
22. # conditional distribution of A given C == "c".
    
23. table(cpdist(fitted, "A", (C == "c")))
    
24. ## Gaussian Bayesian network.
    
25. data(gaussian.test)
    
26. fitted = bn.fit(hc(gaussian.test), gaussian.test)
    
27. # the result should be around 0.04.
    
28. cpquery(fitted,
    
29. event = ((A >= 0) & (A <= 1)) & ((B >= 0) & (B <= 3)),
    
30. evidence = (C + D < 10))

复制代码

1. Examples
   
2. bn = model2network("[A][C|A][B|C]")
   
3. dsep(bn, "A", "B", "C")
   
4. bn = model2network("[A][C][B|A:C]")
   
5. dsep(bn, "A", "B", "C")

复制代码