Implementation of CausalCoxMGM algorithm and scripts for analysis of simulated and real-world biomedical datasets.
To install CoxMGM and CausalCoxMGM, run the following command in the terminal:
R CMD INSTALL rCausalMGM
or alternatively:
R CMD INSTALL rCausalMGM/rCausalMGM_1.0.tar.gz
First, we begin by loading the necessray R packages for this analysis.
library(rCausalMGM)
library(survival)
library(dplyr)
library(survminer)
library(ggplot2)Next, we load the WHAS500 dataset and format the data appropriately for analysis by CausalCoxMGM. First, we create our censored variables, Survival and LengthOfStay, using the Surv function from the survival R package. Then, we exclude features that are not clinical measurements or patient outcomes, and convert discrete variables into factors.
data <- read.csv('cvd/whas500.csv', row.names=1)
data$LengthOfStay <- Surv(data$los, 1-data$dstat)
data$Survival <- Surv(data$lenfol, data$fstat)
data <- data %>% select(-c("admitdate", "disdate", "fdate", "dstat", "lenfol", "fstat", "year", "los"))
data <- data %>% mutate_at(c("cvd", "afb", "sho", "chf", "av3", "gender", "miord", "mitype"), factor)
head(data)## age gender hr sysbp diasbp bmi cvd afb sho chf av3 miord mitype
## 1 83 0 89 152 78 25.54051 1 1 0 0 0 1 0
## 2 49 0 84 120 60 24.02398 1 0 0 0 0 0 1
## 3 70 1 83 147 88 22.14290 0 0 0 0 0 0 1
## 4 70 0 65 123 76 26.63187 1 0 0 1 0 0 1
## 5 70 0 63 135 85 24.41255 1 0 0 0 0 0 1
## 6 70 0 76 83 54 23.24236 1 0 0 0 1 0 0
## LengthOfStay Survival
## 1 5 2178+
## 2 5 2172+
## 3 5 2190+
## 4 10 297
## 5 6 2131+
## 6 1+ 1
Next, we perform stratified 5-fold cross-validation to select model hyperparameters. Cross-validation folds are stratified to have approximately the same number of all-cause mortality and hospital discharge events. We measure model performance based on the total deviance of the censored outcomes. Rather than performing a grid search over the hyperparameters, we perform a random search.
set.seed(43)
idx00 <- which(data$Survival[,2]==0 & data$LengthOfStay[,2]==0)
idx10 <- which(data$Survival[,2]==1 & data$LengthOfStay[,2]==0)
idx01 <- which(data$Survival[,2]==0 & data$LengthOfStay[,2]==1)
idx11 <- which(data$Survival[,2]==1 & data$LengthOfStay[,2]==1)
foldid <- rep(0, 500)
foldid[idx00] <- sample(((1:length(idx00))-1) %% 5 + 1)
foldid[idx10] <- sample(((1:length(idx10))-1) %% 5 + 1)
foldid[idx01] <- sample(((1:length(idx01))-1) %% 5 + 1)
foldid[idx11] <- sample(((1:length(idx11))-1) %% 5 + 1)
table(foldid, data$Survival[,2])##
## foldid 0 1
## 1 57 44
## 2 57 43
## 3 57 43
## 4 57 43
## 5 57 42
table(foldid, data$LengthOfStay[,2])##
## foldid 0 1
## 1 8 93
## 2 8 92
## 3 8 92
## 4 8 92
## 5 7 92
lambdas <- runif(100, 0.05, 0.5)
alphas <- runif(100, 0.01, 0.25)
loglik <- matrix(0, 100, 5)
size <- matrix(0, 100, 5)
for (k in 1:5) {
ig.path <- coxmgmPath(data[foldid!=k,], lambdas=lambdas, rank=F)
idx <- 0
for (ig in ig.path$graphs) {
idx <- idx + 1
g <- fciStable(data[foldid!=k,], initialGraph=ig,
alpha=alphas[idx], orientRule="maxp", rank=F)
mb <- g$markov.blankets$Survival
size[idx,k] <- length(mb)
if (length(mb)==1) {
mb <- c(1)
}
f <- as.formula(paste("Survival ~", paste(mb, collapse=" + ")))
res <- coxph(f, data[foldid!=k,])
test.risk <- predict(res, newdata=data[foldid==k,])
res.test <- coxph(Survival ~ offset(test.risk), data[foldid==k,])
loglik[idx,k] <- -as.numeric(logLik(res.test))
mb <- g$markov.blankets$LengthOfStay
size[idx,k] <- size[idx,k] + length(mb)
if (length(mb)==1) {
mb <- c(1)
}
f <- as.formula(paste("LengthOfStay ~", paste(mb, collapse=" + ")))
res <- coxph(f, data[foldid!=k,])
test.risk <- predict(res, newdata=data[foldid==k,])
res.test <- coxph(LengthOfStay ~ offset(test.risk), data[foldid==k,])
loglik[idx,k] <- loglik[idx,k] + -as.numeric(logLik(res.test))
}
}
sizeMean <- rowMeans(size)
loglikMean <- rowMeans(loglik)
loglikSd <- apply(loglik, 1, sd)
plot(sizeMean, loglikMean, pch=19, col='red')
minIdx <- which.min(loglikMean)Now that we have selected the best set of hyperparameters, we learn the final causal graphical model of all-cause mortality and hospital discharge after hospitilazation with acute myocardial infarction. This is done in two stages: first, we learn the undirected CoxMGM that serves as an initial estimate of the adjacencies in the causal graph. Second, we use FCI-Max to prune adjacencies and orient edges.
ig.path <- coxmgmPath(data, lambda=lambdas)
g <- fciStable(data, initialGraph=ig.path$graphs[[minIdx]],
alpha=alphas[minIdx], verbose=T, orientRule="maxp")## Starting FCI-Stable algorithm...
## Starting FAS Stable...
## Searching at depth 0...
## Searching at depth 1...
## Searching at depth 2...
## Searching at depth 3...
## Searching at depth 4...
## FAS Stable Elapsed Time = 0.029 s
## RFCI adjacency pruning...
## Starting Posssible DSep search
## Starting Conservative Orientations...
## Filling Triple Map...
## Orienting colliders...
## Checking Possible-Dsep sets...
## Starting Final Orientations...
## Filling Triple Map...
## Orienting colliders...
## Orienting implied edges...
## FCI-Stable Elapsed Time = 0.25 s
g## Algorithm: CoxMGM-FCI-Max
## Nodes: 15
## Edges: 17
## Unoriented: 4
## Partially Oriented: 6
## Directed: 6
## Bidirected: 1
## lambda = {0.3173977, 0.3173977, 0.3173977, 0.3173977, 0.3173977}
## alpha = 0.02566852
Finally, we can generate a simple plot of the causal graphical model in R.
plot(g, nodeAttr=list(fontsize=36))
The following data files are under the CC0 public domain:
whas500.csv
metabric.rna.full.csv
metabric.rna.erp.full.csv
metabric.rna.ern.full.csv
meta_cohort_common_genes.rds
While all other files composing CausalCoxMGM are under the GPL-3.0 license.