## ----include=FALSE------------------------------------------------------------
##### _Advanced Data Analysis from an Elementary Point of View_ #####
# R code for the chapter "Causal Discovery"
# Please do not re-distribute or use without attribution
# http://www.stat.cmu.edu/~cshalizi/ADAfaEPoV/




## ----results="hide"-----------------------------------------------------------
library(pcalg)
library(SMPracticals)
data(mathmarks)
suffStat <- list(C=cor(mathmarks),n=nrow(mathmarks))
pc.fit <- pc(suffStat, indepTest=gaussCItest, p=ncol(mathmarks),alpha=0.005)


## ----mathmarks-dag,echo=FALSE,results="hide",out.width="0.5\\textwidth"-------
library(Rgraphviz)
plot(pc.fit,labels=colnames(mathmarks),main="Inferred DAG for mathmarks")

## ----eval=FALSE---------------------------------------------------------------
## library(Rgraphviz)
## plot(pc.fit,labels=colnames(mathmarks),main="Inferred DAG for mathmarks")


## -----------------------------------------------------------------------------
summary(pc.fit)


## ----mathmarks-dag2,echo=FALSE,out.width="0.5\\textwidth"---------------------
plot(pc(suffStat, indepTest=gaussCItest, p=ncol(mathmarks),alpha=0.05),
  labels=colnames(mathmarks),main="")

## ----eval=FALSE---------------------------------------------------------------
## plot(pc(suffStat, indepTest=gaussCItest, p=ncol(mathmarks),alpha=0.05),
##   labels=colnames(mathmarks),main="")


## ----include=FALSE------------------------------------------------------------
psychs <- matrix(c(1.0, 0.62, 0.25, 0.16, -0.10, 0.29, 0.18,
               0.62, 1.0, 0.09, 0.28, 0.00, 0.25, 0.15,
               0.25, 0.09, 1.0, 0.07, 0.03, 0.34, 0.19,
               0.16, 0.28, 0.07, 1.0, 0.10, 0.37, 0.41,
               -0.10, 0.00, 0.03, 0.10, 1.0, 0.13, 0.43,
               0.29, 0.25, 0.34, 0.37, 0.13, 1.0, 0.55,
               0.18, 0.15, 0.19, 0.41, 0.43, 0.55, 1.0), ncol=7)
colnames(psychs) <- c("ability","grad","prod","first","sex","cites","pubs")
rownames(psychs) <- colnames(psychs)

## -----------------------------------------------------------------------------
psychs


## ----modeling-psychologists, echo=FALSE, out.width="0.5\\textwidth"-----------
plot(pc(list(C=psychs,n=162),indepTest=gaussCItest,p=7,alpha=0.01),
  labels=colnames(psychs),main="")

## ----eval=FALSE---------------------------------------------------------------
## plot(pc(list(C=psychs,n=162),indepTest=gaussCItest,p=7,alpha=0.01),
##   labels=colnames(psychs),main="")