[!] logit() and the VGAM package

[jarrod-dalton]

I think there is an error in jagsFunctions

VGAM::logit computes the logit, i.e., log(p/(1-p)).

jagsFunctions has an entry for the JAGS function ilogit() that is mapped to VGAM::logit(), when it should actually be mapped to an inverse logit function.

That function should be base::plogis().

[jarrod-dalton]

This is fixed. I changed the function references for logit and ilogit to base::qlogis() and base::plogis(), respectively.

See the commit

[jarrod-dalton]

Actually, it's not entirely fixed.

It looks like rToJags() will need to be modified, but I'll let you do that so I don't mess it up.

Here's an example where we seem to be getting erroneous results for the PE node:

library(HydeNet)
data(PE)

autoNet <- HydeNetwork(~ wells
                         + pe | wells
                         + d.dimer | pregnant*pe
                         + angio | pe
                         + treat | d.dimer*angio
                         + death | pe*treat,
                       data = PE)

observedNodes <- list(pregnant = 1,
                      wells = 5,
                      death = "Yes")
compiledNet2 <- compileJagsModel(net,
                                data = observedNodes)

post2 <- HydePosterior(compiledNet2,
                       variable.names = monitored,
                       n.iter = 30000)

writeNetworkModel(autoNet, pretty=TRUE)

model{
   wells ~ dnorm(3.7942, 0.6305)
   pe ~ dbern(ilogit(-3.90355 + 0.5757 * wells))
   d.dimer ~ dnorm(210.24251 + 68.37938 * (pe == 2) + 29.29496 * (pregnant == 2), 0.03347)
   pi.pregnant[1] <- 0.9014; pi.pregnant[2] <- 0.0986
   pregnant ~ dcat(pi.pregnant)
   pi.angio <- cpt.angio[(pe+1), ]
   angio ~ dcat(pi.angio)
   treat ~ dbern(ilogit(-5.89316 + 1.73354 * (angio == 2) + 0.01994 * d.dimer))
   pi.death <- cpt.death[(pe+1), (treat+1), ]
   death ~ dcat(pi.death)
}

summary(post2)

    d.dimer      death         pe           pregnant     wells    chain_index   obs_index    
 Min.   :215.5   No :    0   No :30000   Min.   :1   Min.   :5   Min.   :1    Min.   :    1  
 1st Qu.:235.4   Yes:30000   Yes:    0   1st Qu.:1   1st Qu.:5   1st Qu.:1    1st Qu.: 7501  
 Median :239.2                           Median :1   Median :5   Median :1    Median :15000  
 Mean   :239.1                           Mean   :1   Mean   :5   Mean   :1    Mean   :15000  
 3rd Qu.:242.9                           3rd Qu.:1   3rd Qu.:5   3rd Qu.:1    3rd Qu.:22500  
 Max.   :264.3                           Max.   :1   Max.   :5   Max.   :1    Max.   :30000 

With wells == 5, the equation for node pe should evaluate to a probability of 0.2640449. But all our posterior samples have pe = 'No'!

[nutterb]

So....the one line I changed shouldn't affect anything, really (it looks like all of the translations are hard coded, and the jagsFns dataset is more of a reference than anything...maybe I use it somewhere else, but I don't use it in rToJags), but I'm getting reasonable posteriors. So yay? for now?

Here's the code I'm testing on. Is it giving you reasonable posteriors?

library(HydeNet)
data(PE)

autoNet <- HydeNetwork(~ wells
                       + pe | wells
                       + d.dimer | pregnant*pe
                       + angio | pe
                       + treat | d.dimer*angio
                       + death | pe*treat,
                       data = PE)

observedNodes <- list(pregnant = "Yes",
                      wells = 5,
                      death = "Yes")
compiledNet2 <- compileJagsModel(autoNet,
                                 data = observedNodes)

post2 <- HydePosterior(compiledNet2,
                       variable.names = c("wells", "pe", "d.dimer"),
                       n.iter = 30000)
summary(post2)

compiledNet3 <- compileJagsModel(autoNet,
                                 data = list(wells =5))
post3 <- HydePosterior(compiledNet3,
                       data = list(wells = 5), 
                       variable.names = c("wells", "pe", "d.dimer"),
                       n.iter = 30000)
summary(post3)

[jarrod-dalton]

Everything appears to be working here. Maybe I was just hallucinating.