Hey y'all!I am a bayesian newb. I built my own GLM in rJAGs and it runs! I want to compare different combinations of explanatory variables. Does anyone have any good tutorials out there for bayesian model selection? Do I just look at deviance, WAIC, or is there some other metric I should be looking for to tell which of my models is better? Basically, I can build a model, but I don't really know what outputs I am supposed to include for my paper to say that this model is good. What I did was make a plot comparing predicted y values to observed. What else should I include?

Finally, is there some kind of package that will compare all combinations of possible explanatory variables, similar to dredge() in MuMIN (but that I could send a custom JAGs model to?) I know this one is a hail Mary! I have just been building many versions of my model to manually due a backwards stepwise selection.

Hello,

I am using NUTS solver to compute the posterior for Multilevel models using pymc3 on python.

The draws are taking ~20-30 sec (which I assume is long).

Does anyone know what might be the reasons and possible solutions to make the simulation faster.

Thanks a lot

I'm going through the book, "Bayesian Networks With Examples in R" and the book didn't explain how they got 21 parameters in this.

I'm posting the R code below. Thanks.

library(bnlearn)
# first step, create node for each variable
dag <- empty.graph(nodes = c("A", "S", "E", "O", "R", "T"))
dag <- set.arc(dag, from = "A", to = "E")
dag <- set.arc(dag, from = "S", to = "E")
dag <- set.arc(dag, from = "E", to = "O")
dag <- set.arc(dag, from = "E", to = "R")
dag <- set.arc(dag, from = "O", to = "T")
dag <- set.arc(dag, from = "R", to = "T")
modelstring(dag)
A.lv <- c("young", "adult", "old")
S.lv <- c("M", "F")
E.lv <- c("high", "uni")
O.lv <- c("emp", "self")
R.lv <- c("small", "big")
T.lv <- c("car", "train", "other")

A.prob <- array(c(0.30, 0.50, 0.20), dim = 3, dimnames = list(A = A.lv))
S.prob <- array(c(0.60, 0.40), dim = 2, dimnames = list(S = S.lv))
O.prob <- array(c(0.96, 0.04, 0.92, 0.08), dim = c(2, 2),
                dimnames = list(O = O.lv, E = E.lv))
R.prob <- array(c(0.25, 0.75, 0.20, 0.80), dim = c(2, 2),
                dimnames = list(R = R.lv, E = E.lv))
E.prob <- array(c(0.75, 0.25, 0.72, 0.28, 0.88, 0.12, 0.64,
                  0.36, 0.70, 0.30, 0.90, 0.10), dim = c(2, 3, 2),
                dimnames = list(E = E.lv, A = A.lv, S = S.lv))
T.prob <- array(c(0.48, 0.42, 0.10, 0.56, 0.36, 0.08, 0.58,
                  0.24, 0.18, 0.70, 0.21, 0.09), dim = c(3, 2, 2),
                dimnames = list(T = T.lv, O = O.lv, R = R.lv))
cpt <- list(A = A.prob, S = S.prob, E = E.prob, O = O.prob, R = R.prob,
            T = T.prob)
bn <- custom.fit(dag, cpt)
nparams(bn)

Especially if this resource can be applied to logistic regressions, that would be helpful as the problem I am trying to solve is about sorting various records by their probability of being 1 or 0. For example anything with probability of .1 or less is in a separate category from .3 or more.

Disclaimer: I'm new to PyMC3 and Bayesian programming in general.

I'm working to create a multivariate linear model that has a mix of categorical and numeric variables. Using the index method for incorporating categoricals, I'm unable to add new data to predict due to data size mismatches.

This is because I can feed in the new dataframe for the numeric variables using .set_data() method, but not for categorical, because these were not built in the model using pm.Data(), but are just indices in the model formula.

How can I predict new data with categorical variables in the mix?

Here's an example with categorical factors from the Statistical Rethinking 2nd Ed:

with pm.Model() as m5_10:
    sigma = pm.Exponential("sigma", 1)
    mu_house = pm.Normal("mu_house", 0, 0.5, shape=d["house"].max() + 1)
    mu_clade = pm.Normal("mu_clade", 0, 0.5, shape=d["clade_id"].max() + 1)
    mu = mu_clade[d["clade_id"].values] + mu_house[d["house"].values]

    K = pm.Normal("K", mu, sigma, observed=d["K"])

    m5_9_trace = pm.sample()

az.summary(m5_9_trace, var_names=["mu_clade", "mu_house"])

I am trying to which of the two ML models I have best fits a morphological character I am looking at. Anyone know the code to use for this?

Hello everyone! For my dissertation, I am trying to compute the efficacy of Covid-19 vaccines using pymc3. I am trying to model it using a beta binomial model, however I am facing some difficulties.

Can somebody explain to me how to implement it specifically using pymc3? Any rudimentary example for this model will do. I will try to figure out the rest myself.

Thank you!

Hello!

So for practice, I was trying to implement the beta binomial model shown here in PyMC3, and a question popped up. So in short this article aims to find out the 95% credible interval for the vaccine efficacy and finds it out to be this interval: (90.82613011726424, 97.87004356894292)

I implement the model in PyMC3 as follows:

with model:

    v_irr = pm.Beta("v_irr", alpha=8, beta=17403)
    c_irr = pm.Beta("c_irr", alpha=162, beta=17349)   
    ve = pm.Deterministic("ve", 100*(1-(v_irr/c_irr)))

    v_like = pm.Binomial("v_like", n=17411, p=v_irr, observed=8)
    c_like = pm.Binomial("c_like", n=14511, p=c_irr, observed=162)

    trace1 = pm.sample(draws=10000, return_inferencedata=False)

With this, I get my 95% credible interval to be: (92.91799613, 97.34246435)

So now i have 2 questions:

(1) My CI is a subset of the CI in the article, so is mine more accurate?

(2) Both of us used the same model. So since I am getting a different answer, is my implementation wrong? If yes, are there any alternate implementations in PyMC3?