Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

Powered by ReviewNB

I am facing an error in cell 16, would like some help with resolving it.
Also, the results from manual model (cell 11) and using bambi (cell 13) are quite different, so I want some help with understanding that as well.

View / edit / reply to this conversation on ReviewNB

OriolAbril commented on 2021-07-28T04:19:10Z
----------------------------------------------------------------

@aloctavodia for this to work we need to add the transformed variables to the inferencedata object (or to the dict that is generated from the summary df). Does doing that manually make sense? Is there some bambi magic to help with that? Or maybe this can be removed altogether and keep only waic and loo?

chiral-carbon commented on 2021-08-31T16:33:47Z
----------------------------------------------------------------

the error message has changed here from two commits ago, it is to do with the arguments passed to models_lin[nm].backend.model.logp() or models_quad[nm].backend.model.logp()

Sayam753 commented on 2021-09-08T15:46:06Z
----------------------------------------------------------------

I have two questions and a working solution to fix the error above.

Here comes the questions first -

The description above the code cell says - "Evaluate log likelihoods straight from model.logp"

But the code cell does not calculate log likelihood of the model, instead it computes log probability of the entire model.

For context, PyMC3 model's log prob = sum over all dimensions of log prob calculated for all free RVs + log likelihood of observed RVs

Question 1 - For this part of model selection, are we interested to look into log probability of the entire model or just the model's log likelihood?

Question 2 - If we are interested in model's log likelihood, can we use the log_likelihood present in posterior trace? Just sum it over all dimensions and compare?

Coming to solution to fix the error

1. There are shape mismatches, because bambi creates Intercept prior of shape 1 (reference here) and all variables' mean from arviz summary has shape = ()
2. Then we also need transformed samples to compute model.logp.

I created two helper functions that rectify the shapes and to include transformed samples in the trace.

And things worked. Here is the GitHub gist of my experiment - https://nbviewer.jupyter.org/gist/Sayam753/d605220d570be6697b4755a9f83a5dd4

I think these functions can be used by bambi.

View / edit / reply to this conversation on ReviewNB

OriolAbril commented on 2021-07-28T04:19:10Z
----------------------------------------------------------------

In my opinion, this should create two dictionaries of traces, one for lin and the other for quad, instead of combining them all into the same dict. Something like:

lin_trace_dict = dict()
quad_trace_dict = dict()
for nm in ["k1", "k2", "k3", "k4", "k5"]:
    lin_trace_dict.update({f"{nm}": traces_lin[nm]})
    quad_trace_dict.update({f"{nm}": traces_quad[nm]})

Then we also generate two dfs instead of one.

chiral-carbon commented on 2021-07-30T11:40:53Z
----------------------------------------------------------------

sure yes, it makes more sense

Sayam753 commented on 2021-08-29T20:50:12Z
----------------------------------------------------------------

Why are we creating a new copy of traces_lin and traces_quad ? I think they can be used as it is in downstream tasks.

OriolAbril commented on 2021-08-30T14:28:32Z
----------------------------------------------------------------

I thought we were converting a list into a dict, this is not necessary indeed.

View / edit / reply to this conversation on ReviewNB

OriolAbril commented on 2021-07-28T04:19:11Z
----------------------------------------------------------------

And therefore also call plot_compare twice to get two subplots side by side.

_, axes = plt.subplots(1, 2, sharey=True)
az.plot_compare(dfwaic_lin, ax=axes[0])
az.plot_compare(dfwaic_quad, ax=axes[1])
axes[0].set_title("Linear data")

axes[1].set_title("Quadratic data");

sure yes, it makes more sense

View entire conversation on ReviewNB

View / edit / reply to this conversation on ReviewNB

Sayam753 commented on 2021-08-29T20:38:57Z
----------------------------------------------------------------

Line #133.        Create a polynomial modelspec string for patsy

should be

Line #133.        Create a polynomial modelspec string for bambi

Why are we creating a new copy of traces_lin and traces_quad ? I think they can be used as it is in downstream tasks.

View entire conversation on ReviewNB

I thought we were converting a list into a dict, this is not necessary indeed.

View entire conversation on ReviewNB

View / edit / reply to this conversation on ReviewNB

OriolAbril commented on 2021-08-30T14:39:46Z
----------------------------------------------------------------

This has to be removed

View / edit / reply to this conversation on ReviewNB

OriolAbril commented on 2021-08-30T14:39:47Z
----------------------------------------------------------------

waic stands for "widely applicable information criterion" or "Watanabe–Akaike information criterion" the appliable seems like a momentary confusion.

Also remember to update the descriptions of what the code does. We no longer loop over models to compute waic, giving the dictionary to arviz is enough.

chiral-carbon commented on 2021-08-31T13:13:39Z
----------------------------------------------------------------

will change this.

View / edit / reply to this conversation on ReviewNB

OriolAbril commented on 2021-08-30T14:39:47Z
----------------------------------------------------------------

Now model is converging and results are as expected, k1 for linear data and k2 for quadratic data. Also note the change in scale from deviance to log, so we no longer want lower waic (see https://discourse.pymc.io/t/is-model-with-high-or-low-waic-and-loo-better/7874/2 for more details).

chiral-carbon commented on 2021-08-31T13:15:45Z
----------------------------------------------------------------

right. so should I remove the points mentioning importance/requirement of lower waic everywhere and instead say higher waic is best (due to the change in scale from deviance to log)?

chiral-carbon commented on 2021-08-31T16:39:21Z
----------------------------------------------------------------

should I also mention the scale used for clarification?

OriolAbril commented on 2021-09-08T09:57:45Z
----------------------------------------------------------------

I would add a note about the existence and extended usage of multiple scales and recommend to read the doctring of the function in case of doubt. See also other changes about adding references

View / edit / reply to this conversation on ReviewNB

OriolAbril commented on 2021-08-30T14:39:48Z
----------------------------------------------------------------

same comment here.

View / edit / reply to this conversation on ReviewNB

OriolAbril commented on 2021-08-30T14:40:43Z
----------------------------------------------------------------

This has to be removed

View / edit / reply to this conversation on ReviewNB

OriolAbril commented on 2021-08-30T14:40:43Z
----------------------------------------------------------------

waic stands for "widely applicable information criterion" or "Watanabe–Akaike information criterion" the appliable seems like a momentary confusion.

Also remember to update the descriptions of what the code does. We no longer loop over models to compute waic, giving the dictionary to arviz is enough.

View / edit / reply to this conversation on ReviewNB

OriolAbril commented on 2021-08-30T14:40:44Z
----------------------------------------------------------------

Now model is converging and results are as expected, k1 for linear data and k2 for quadratic data. Also note the change in scale from deviance to log, so we no longer want lower waic (see https://discourse.pymc.io/t/is-model-with-high-or-low-waic-and-loo-better/7874/2 for more details).

View / edit / reply to this conversation on ReviewNB

OriolAbril commented on 2021-08-30T14:40:45Z
----------------------------------------------------------------

same comment here.

will change this.

View entire conversation on ReviewNB

right. so should I remove the points mentioning importance/requirement of lower waic everywhere?

View entire conversation on ReviewNB

View / edit / reply to this conversation on ReviewNB

chiral-carbon commented on 2021-08-31T13:33:41Z
----------------------------------------------------------------

hey @tomicapretto I had run this same code block for GLM using bambi in the previous commit for this notebook but did not get this error. I used bambi==0.5.0 till the previous commit and am using bambi==0.6.0 now, I'm not sure if this error is due to some update or not, could you help?

tomicapretto commented on 2021-08-31T13:44:57Z
----------------------------------------------------------------

If you write "Intercept" instead of "intercept", it should work. But this is still a bug i'd like to fix in the future. Let me know if it fixes the problem!

chiral-carbon commented on 2021-08-31T15:50:32Z
----------------------------------------------------------------

Thank you, it did! very silly of me to not realize this 😅

If you write "Intercept" instead of "intercept", it should work. But this is still a bug i'd like to fix in the future. Let me know if it fixes the problem!

View entire conversation on ReviewNB

Thank you, it did! very silly of me to not realize this 😅

View entire conversation on ReviewNB

the error message has changed here from two commits ago, it is to do with the arguments passed to models_lin[nm].backend.model.logp(). I guess logp is not a valid function name in bambi, but it is in pymc backend?

View entire conversation on ReviewNB

should I also mention the scale used for clarification?

View entire conversation on ReviewNB

I would add a note about the existence and extended usage of multiple scales and recommend to read the doctring of the function in case of doubt. See also other changes about adding references

View entire conversation on ReviewNB

View / edit / reply to this conversation on ReviewNB

OriolAbril commented on 2021-09-08T10:04:17Z
----------------------------------------------------------------

Add an extra sentence to the paragraph "WAIC is implemented in {func}arviz.waic which we will use below"

Note: reviewnb renders everything wrong in the comments, keep in mind the right syntax as copy pasting won't work. ref https://arviz-devs.github.io/arviz/contributing/developer_guide.html#hyperlinks

View / edit / reply to this conversation on ReviewNB

OriolAbril commented on 2021-09-08T10:04:18Z
----------------------------------------------------------------

I would add a link to arviz.loo function. For example, changing the last sentence to "One alternative is applying the numerical approach using the posterier trace as suggested in Vehtari et al 2015, which is implemented in {func}arviz.loo. We will use it via {func}arviz.compare to compare our models.

View / edit / reply to this conversation on ReviewNB

Sayam753 commented on 2021-09-08T14:46:45Z
----------------------------------------------------------------

Line #155.                fml, df, priors={"intercept": bmb.Prior("Normal", mu=0, sigma=100)}, family="gaussian"

There was a mismatch between naming of prior either "intercept" or "Intercept" in comments below. Do we need to correct things here we well?

cc @tomicapretto

tomicapretto commented on 2021-09-08T14:52:41Z
----------------------------------------------------------------

Yes, the intercept is always named "Intercept", so "Intercept" is what you need if you want to put a prior on it.

Yes, the intercept is always named "Intercept", so "Intercept" is what you need if you want to put a prior on it.

View entire conversation on ReviewNB

I have two questions and a working solution to fix the error above.

Here comes the questions first -

The description above the code cell says - "Evaluate log likelihoods straight from model.logp"

But the code cell does not calculate log likelihood of the model, instead it computes log probability of the entire model.

For context, PyMC3 model's log prob = sum over all dimensions of log prob calculated for all free RVs + log likelihood of observed RVs

Question 1 - For this part of model selection, are we interested to look into log probability of the entire model or just the model's log likelihood?

Question 2 - If we are interested in model's log likelihood, can we use the log_likelihood present in posterior trace? Just sum it over all dimensions and compare?

Coming to solution to fix the error

1. There are shape mismatches, because bambi creates Intercept prior of shape 1 (reference here) and all variables' mean from arviz summary has shape = ()
2. Then we also need transformed samples to compute model.logp.

I created two helper functions that rectify the shapes and to include transformed samples in the trace.

And things worked. Here is the GitHub gist of my experiment - https://nbviewer.jupyter.org/gist/Sayam753/d605220d570be6697b4755a9f83a5dd4

I think these functions can be used by bambi.

View entire conversation on ReviewNB