-
Notifications
You must be signed in to change notification settings - Fork 12
Description
I suspect this is not a bug but rather either something that just can not be done, or I am doing something wrong.
This first code block works fine. Basic problem is to have linear predictors for coefficients of a beta distribution for the probability parameter of some binomial trials under two conditions.
`library(tidyverse)
library(causact)
set.seed(2023)
n_subjects <- 200
df <- tibble(id=1:(n_subjects)) %>%
mutate(
cond = sample(c('ctrl','trtmt'),nrow(.), replace=TRUE) ,
disc_cov = sample(c('a','b','c'),nrow(.), replace=TRUE),
cont_cov = rnorm(nrow(.), 0,1),
trials = sample(5:25, nrow(.), replace=TRUE),
lin_pred = 0.2cont_cov+0.3case_when(disc_cov=='a'~1, disc_cov=='b'0, disc_cov=='c' -1)+if_else(cond=='trtmt',1,0),
prob = 1 / (1+exp(-lin_pred)),
successes = rbinom(nrow(.),trials, prob)
) %>%
select(successes, cond, trials)
graph <- dag_create() %>%
dag_node("Successes", "s", rhs = binomial(trials, prob), data = df$successes) %>%
dag_node("Probability", "prob", child = "s", rhs = beta(lin_pred_a, lin_pred_b)) %>%
dag_node("Linear Predictor A", "lin_pred_a", rhs = c_cond_a, child = "prob") %>%
dag_node("Linear Predictor B", "lin_pred_b", rhs = c_cond_b, child = "prob") %>%
dag_node("Cond Coefficient A", "c_cond_a", rhs = exponential(1), child = c("lin_pred_a")) %>%
dag_node("Cond Coefficient B", "c_cond_b", rhs = exponential(1), child = c("lin_pred_b")) %>%
dag_node("Trials", "trials", child = "s", data = df$trials) %>%
dag_plate("Cond A", "cond_a", nodeLabels = c("c_cond_a"), data = df$cond, addDataNode = TRUE) %>%
dag_plate("Cond B", "cond_b", nodeLabels = c("c_cond_b"), data = df$cond, addDataNode = TRUE) %>%
dag_plate("Observation", "i", nodeLabels = c("s", "trials", "prob", "lin_pred_a", "lin_pred_b","cond_a","cond_b"))
graph %>% dag_render()
drawsDF <- graph %>% dag_numpyro()
drawsDF %>% dagp_plot()`
Plots, etc of the above look ok.
But, this code
graph <- dag_create() %>% dag_node("Successes", "s", rhs = binomial(trials, prob), data = df$successes) %>% dag_node("Probability", "prob", child = "s", rhs = beta(lin_pred_a, lin_pred_b)) %>% dag_node("Linear Predictor A", "lin_pred_a", rhs = c_cond_a, child = "prob") %>% dag_node("Linear Predictor B", "lin_pred_b", rhs = c_cond_b, child = "prob") %>% dag_node("Cond Coefficient A", "c_cond_a", rhs = exponential(1), child = c("lin_pred_a")) %>% dag_node("Cond Coefficient B", "c_cond_b", rhs = exponential(1), child = c("lin_pred_b")) %>% dag_node("Trials", "trials", child = "s", data = df$trials) %>% dag_plate("Cond", "cond", nodeLabels = c("c_cond_a", "c_cond_b"), data = df$cond, addDataNode = TRUE) %>% dag_plate("Observation", "i", nodeLabels = c("s", "trials", "prob", "lin_pred_a", "lin_pred_b","cond"))
gives an error : Error in dag_edge(., from = childrenNewEdgeDF$parentNodeID, to = childrenNewEdgeDF$childID) :
You have attempted to connect 8 to itself, but this would create a cycle and is not a Directed Acyclic Graph. You cannot connect a node to itself in a DAG.
What I am trying to do is merge the two nodes "cond_a" and "cond_b" from the first example into one node "cond" in the second, as both refer to the same data field (df$cond)
I suspect this is not an actual bug - but maybe there is a way to do it?