2890 字

# 因果结构提取

``````E <- 6*rnorm(1000,10,1)+rnorm(1000)
G <- rnorm(1000,50,2)+rnorm(1000)
M <- E*0.7+G*0.3+rnorm(1000)
Epi <- E*0.1+G*0.9+rnorm(1000)
Asthma <- 0.6*M+0.4*Epi+rnorm(1000)
data <- cbind.data.frame(E,G,M,Epi,Asthma)
cor(data)``````
``````##               E        G      M    Epi Asthma
## E       1.00000 -0.02042 0.9647 0.2565 0.8551
## G      -0.02042  1.00000 0.1314 0.8550 0.3549
## M       0.96474  0.13144 1.0000 0.3733 0.9159
## Epi     0.25647  0.85500 0.3733 1.0000 0.5888
## Asthma  0.85509  0.35488 0.9159 0.5888 1.0000``````
``pairs(data)``

``cor.test(Asthma,E)``
``````##
##  Pearson's product-moment correlation
##
## data:  Asthma and E
## t = 52, df = 998, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8375 0.8709
## sample estimates:
##    cor
## 0.8551``````
``cor.test(Asthma,G)``
``````##
##  Pearson's product-moment correlation
##
## data:  Asthma and G
## t = 12, df = 998, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2995 0.4079
## sample estimates:
##    cor
## 0.3549``````

``````library(pcalg)
suffStat <- list(C=cor(data), n = nrow(data))
pc.fit <- pc(suffStat, gaussCItest, p = ncol(data), alpha = 0.01, verbose = TRUE)``````
``````## Order=0; remaining edges:20
## x= 1  y= 2  S=  : pval = 0.519
## x= 1  y= 3  S=  : pval = 0
## x= 1  y= 4  S=  : pval = 1.2e-16
## x= 1  y= 5  S=  : pval = 0
## x= 2  y= 3  S=  : pval = 2.987e-05
## x= 2  y= 4  S=  : pval = 0
## x= 2  y= 5  S=  : pval = 1.067e-31
## x= 3  y= 1  S=  : pval = 0
## x= 3  y= 2  S=  : pval = 2.987e-05
## x= 3  y= 4  S=  : pval = 3.173e-35
## x= 3  y= 5  S=  : pval = 0
## x= 4  y= 1  S=  : pval = 1.2e-16
## x= 4  y= 2  S=  : pval = 0
## x= 4  y= 3  S=  : pval = 3.173e-35
## x= 4  y= 5  S=  : pval = 4.855e-101
## x= 5  y= 1  S=  : pval = 0
## x= 5  y= 2  S=  : pval = 1.067e-31
## x= 5  y= 3  S=  : pval = 0
## x= 5  y= 4  S=  : pval = 4.855e-101
## Order=1; remaining edges:18
## x= 1  y= 3  S= 4 : pval = 0
## x= 1  y= 3  S= 5 : pval = 0
## x= 1  y= 4  S= 3 : pval = 2.175e-46
## x= 1  y= 4  S= 5 : pval = 3.191e-101
## x= 1  y= 5  S= 3 : pval = 2.84e-18
## x= 1  y= 5  S= 4 : pval = 0
## x= 2  y= 3  S= 4 : pval = 1.207e-38
## x= 2  y= 3  S= 5 : pval = 1.84e-72
## x= 2  y= 4  S= 3 : pval = 0
## x= 2  y= 4  S= 5 : pval = 0
## x= 2  y= 5  S= 3 : pval = 4.245e-101
## x= 2  y= 5  S= 4 : pval = 1.425e-31
## x= 3  y= 1  S= 2 : pval = 0
## x= 3  y= 1  S= 4 : pval = 0
## x= 3  y= 1  S= 5 : pval = 0
## x= 3  y= 2  S= 1 : pval = 1.174e-94
## x= 3  y= 2  S= 4 : pval = 1.207e-38
## x= 3  y= 2  S= 5 : pval = 1.84e-72
## x= 3  y= 4  S= 1 : pval = 1.184e-65
## x= 3  y= 4  S= 2 : pval = 1.006e-69
## x= 3  y= 4  S= 5 : pval = 4.762e-71
## x= 3  y= 5  S= 1 : pval = 5.469e-142
## x= 3  y= 5  S= 2 : pval = 0
## x= 3  y= 5  S= 4 : pval = 0
## x= 4  y= 1  S= 2 : pval = 8.164e-77
## x= 4  y= 1  S= 3 : pval = 2.175e-46
## x= 4  y= 1  S= 5 : pval = 3.191e-101
## x= 4  y= 2  S= 1 : pval = 0
## x= 4  y= 2  S= 3 : pval = 0
## x= 4  y= 2  S= 5 : pval = 0
## x= 4  y= 3  S= 1 : pval = 1.184e-65
## x= 4  y= 3  S= 2 : pval = 1.006e-69
## x= 4  y= 3  S= 5 : pval = 4.762e-71
## x= 4  y= 5  S= 1 : pval = 2.8e-195
## x= 4  y= 5  S= 2 : pval = 7.757e-101
## x= 4  y= 5  S= 3 : pval = 5.666e-140
## x= 5  y= 1  S= 2 : pval = 0
## x= 5  y= 1  S= 3 : pval = 2.84e-18
## x= 5  y= 1  S= 4 : pval = 0
## x= 5  y= 2  S= 1 : pval = 4.455e-179
## x= 5  y= 2  S= 3 : pval = 4.245e-101
## x= 5  y= 2  S= 4 : pval = 1.425e-31
## x= 5  y= 3  S= 1 : pval = 5.469e-142
## x= 5  y= 3  S= 2 : pval = 0
## x= 5  y= 3  S= 4 : pval = 0
## x= 5  y= 4  S= 1 : pval = 2.8e-195
## x= 5  y= 4  S= 2 : pval = 7.757e-101
## x= 5  y= 4  S= 3 : pval = 5.666e-140
## Order=2; remaining edges:18
## x= 1  y= 3  S= 4 5 : pval = 2.984e-295
## x= 1  y= 4  S= 3 5 : pval = 3.844e-29
## x= 1  y= 5  S= 3 4 : pval = 0.577
## x= 2  y= 3  S= 4 5 : pval = 2.05e-08
## x= 2  y= 4  S= 3 5 : pval = 4.093e-267
## x= 2  y= 5  S= 3 4 : pval = 0.4735
## x= 3  y= 1  S= 2 4 : pval = 0
## x= 3  y= 1  S= 2 5 : pval = 4.25e-304
## x= 3  y= 1  S= 4 5 : pval = 2.984e-295
## x= 3  y= 2  S= 1 4 : pval = 1.062e-28
## x= 3  y= 2  S= 1 5 : pval = 3.746e-09
## x= 3  y= 2  S= 4 5 : pval = 2.05e-08
## x= 3  y= 4  S= 1 2 : pval = 0.1595
## x= 3  y= 5  S= 1 2 : pval = 1.519e-51
## x= 3  y= 5  S= 1 4 : pval = 1.218e-71
## x= 3  y= 5  S= 2 4 : pval = 0
## x= 4  y= 1  S= 2 3 : pval = 1.981e-08
## x= 4  y= 1  S= 2 5 : pval = 0.1367
## x= 4  y= 2  S= 1 3 : pval = 0
## x= 4  y= 2  S= 1 5 : pval = 3.064e-224
## x= 4  y= 2  S= 3 5 : pval = 4.093e-267
## x= 4  y= 5  S= 1 2 : pval = 7.358e-24
## x= 4  y= 5  S= 1 3 : pval = 1.651e-120
## x= 4  y= 5  S= 2 3 : pval = 7.409e-36
## x= 5  y= 3  S= 1 2 : pval = 1.519e-51
## x= 5  y= 3  S= 1 4 : pval = 1.218e-71
## x= 5  y= 3  S= 2 4 : pval = 0
## x= 5  y= 4  S= 1 2 : pval = 7.358e-24
## x= 5  y= 4  S= 1 3 : pval = 1.651e-120
## x= 5  y= 4  S= 2 3 : pval = 7.409e-36
## Order=3; remaining edges:10``````
``plot(pc.fit,main='Demo',labels=c("E","G","M","Epi","Asthma"))``
``## Loading required namespace: Rgraphviz``