Search code examples
rbayesianbayesian-networksrjags

Define conditional linear Gaussian network using rjags

I am struggling to define a conditional linear Gaussian Bayesian network using rjags. (A clg BN is defined by a continuous child node (outcome) having both a continuous normal and a discrete parent (predictors))

For the net below, A is discrete, D and E continuous:

enter image description here

For the rjags model, I supose what I want is for the parameters of node E to be defined on the value node A takes: pseudo code

model {  
  A ~ dcat(c(0.0948, 0.9052 ))
  D ~ dnorm(11.87054, 1/1.503111^2)

  if A==a then E ~ dnorm(6.558366 + 1.180965*D, 1/2.960002^2) 
  if A==b then E ~ dnorm(3.370021 + 1.532289*D, 1/6.554402^2)   
}

I can get something working by using the code below but it would get confusing quickly with more predictor and categorical levels.

library(rjags)

model <- textConnection("model {  
  A ~ dcat(c(0.0948, 0.9052 ))
  D ~ dnorm(11.87054, 1/1.503111^2)

  int = 6.558366 - (A==2)*(6.558366 - 3.370021) 
  slope = 1.180965 - (A==2)*(1.180965 - 1.532289)
  sig = 2.960002 - (A==2)*(2.960002 - 6.554402)

  E ~ dnorm(int + slope*D, 1/sig^2) 
}")

jg <- jags.model(model, n.adapt = 1000

My question: How can i define this model succinctly please?

The data came from

library(bnlearn)
net = model2network("[A][D][E|A:D]")
ft = bn.fit(net, clgaussian.test[c("A", "D", "E")])

coef(ft) 
structure(list(A = structure(c(0.0948, 0.9052), class = "table", .Dim = 2L, .Dimnames = list(
    c("a", "b"))), D = structure(11.8705363469396, .Names = "(Intercept)"), 
    E = structure(c(6.55836552742708, 1.18096500477159, 3.37002124328838, 
    1.53228891423418), .Dim = c(2L, 2L), .Dimnames = list(c("(Intercept)", 
    "D"), c("0", "1")))), .Names = c("A", "D", "E"))

sigma(ft)
structure(list(A = NA, D = 1.50311121682603, E = structure(c(2.96000206596326, 
6.55440224877698), .Names = c("0", "1"))), .Names = c("A", "D", 
"E"))
Solution

  • You just need to use your variable A as an indexing parameter:

    library('rjags')

model <- "
model {  
  A ~ dcat(c(0.0948, 0.9052 ))
  D ~ dnorm(11.87054, 1/1.503111^2)

  ints <- c(6.558366, 3.370021)
  int <- ints[A]
  slopes <- c(1.180965, 1.532289)
  slope <- slopes[A]
  sigs <- c(2.960002, 6.554402)
  sig <- sigs[A]

  E ~ dnorm(int + slope*D, 1/sig^2) 
}
"

jg <- jags.model(textConnection(model), n.adapt = 1000)

    By the way, as you have lots of fixed quantities in the model, it may make more sense to define these in R and then pass them as data to JAGS. This way you can adjust the values and lengths of the vectors (as long as the lengths of catprobs, ints, slopes and sigs match) without having to modify the JAGS code. For example (using runjags for convenience although also possible with jags):

    library("runjags")

model <- "
model {  
  A ~ dcat(catprobs)
  D ~ dnorm(Dmu, Dprec)

  int <- ints[A]
  slope <- slopes[A]
  sig <- sigs[A]

  E ~ dnorm(int + slope*D, 1/sig^2) 

  #data# catprobs, Dmu, Dprec, ints, slopes, sigs
  #monitor# A, D, E
}
"

catprobs <- c(0.0948, 0.9052)
Dmu <- 11.87054
Dprec <- 1/1.503111^2
ints <- c(6.558366, 3.370021)
slopes <- c(1.180965, 1.532289)
sigs <- c(2.960002, 6.554402)

results <- run.jags(model)
results

    Matt