Блог

Допустим, у меня есть простая DAG, представляющая смешивающую переменную X = Курение, лечение T и результат Y = Смерть, так что:

T ~ X
Y ~ T X

Возможно ли создать синтетический набор данных, состоящий, скажем, из 1 млн наблюдений, который соответствует некоторым заданным условным вероятностям:

 # Pr(smoking):  
smoking <- data.frame(
  smoking = c(0, 1),
  proba = c(0.7, 0.3)
)

# Pr(treatment | smoking):
treatment <- expand.grid(
  smoking = c(0, 1),
  treatment = c(0, 1)
) %>% arrange(smoking, treatment)
treatment$proba <- c(0.8,  0.2, 0.45, 0.55)

# Pr(death | treatment, smoking):
death <- expand.grid(
  treatment = c(0, 1),
  smoking = c(0,1),
  dead = c(0,1)
) %>% 
  arrange(treatment, smoking, dead)
death$proba <- c(0.9, 0.1, 0.2, 0.8, 0.89, 0.11, 0.5, 0.5)
 

Я могу сделать это вручную здесь, потому что это очень простой DAG, но мне было интересно, можно ли это сделать другим, более масштабируемым способом, используя что-то вроде bnlearn .

Текущее решение:

 db <- data.frame(
  smoking = rbinom(n = 1000000, size = 1, prob = 0.3)
  ) 

db$treatment[db$smoking == 0] <- rbinom(n = sum(db$smoking == 0), size = 1, prob = 0.2)
db$treatment[db$smoking == 1] <- rbinom(n = sum(db$smoking == 1), size = 1, prob = 0.55)

db$dead[db$treatment == 0 amp; db$smoking == 0] <- rbinom(
  n = sum(db$treatment == 0 amp; db$smoking == 0), 
  size = 1, prob = 0.1
  )

db$dead[db$treatment == 0 amp; db$smoking == 1] <- rbinom(
  n = sum(db$treatment == 0 amp; db$smoking == 1), 
  size = 1, prob = 0.8
  )

db$dead[db$treatment == 1 amp; db$smoking == 0] <- rbinom(
  n = sum(db$treatment == 1 amp; db$smoking == 0), 
  size = 1, prob = 0.11
  )

db$dead[db$treatment == 1 amp; db$smoking == 1] <- rbinom(
  n = sum(db$treatment == 1 amp; db$smoking == 1), 
  size = 1, prob = 0.5
  )
 

Будет проще позволить существующим пакетам сделать это за вас; например bnlearn . Вы можете использовать custom.fit для указания DAG и CPTS, а затем использовать rbn для извлечения образцов из него.

Пример

 library(bnlearn)

# Specify DAG
net <- model2network("[treatment|smoking][smoking][death|treatment:smoking]")
graphviz.plot(net)

# Define CPTs
smoking <- matrix(c(0.7, 0.3), ncol = 2, dimnames = list(NULL, c("no", "yes")))
treatment <- matrix(c(0.8,  0.2, 0.45, 0.55), ncol = 2, dimnames = list(c("no", "yes"), c("no", "yes")))
death <- array(c(0.9, 0.1, 0.2, 0.8, 0.89, 0.11, 0.5, 0.5), c(2,2,2), dimnames=list(c("no", "yes"), c("no", "yes"), c("no", "yes")))

# Build BN
fit <- custom.fit(net, dist = list(smoking = smoking, treatment = treatment, death = death))

# Draw samples
set.seed(69395642)
samples <- rbn(fit, n=1e6)