Hi all. I would like to start with my appreciation for your work on Bambi and formulae! I think both of these are very valuable additions to the scientific python stack. I saw a thread on twitter that made me think formulae might become the new patsy, so I thought I would describe an affordance in patsy that would be really cool to add to formulae: splines: https://patsy.readthedocs.io/en/latest/spline-regression.html

I would love to be able to include bs(x) in one line, e.g.

model = bmb.Model("test_positive ~ vaccinated + bs(age,df=3) + (vaccinated+vaccinated:day|state)",
                  df, family="bernoulli")

Since this is already implemented in patsy maybe it will not be too hard to add here, as well! Thanks again for all of your work!

I would like to fit a fully specified model. In lme I would write this as:

f = "choice ~ 1 + stimulus*reward + (1+stimulus*reward|subject_id)"`

In formulae this gives the following error:

import formulae as fm
fm.model_description(f)

    891                 return NotImplemented
    892         else:
--> 893             raise ValueError("LHS of group specific term cannot have more than one term.")
    894 
    895     def __repr__(self):

The following does seem to work:

f = "choice ~ 1 + stimulus*reward + (1|subject_id) + (stimulus|subject_id) + (reward|subject_id) + (stimulus:reward|subject_id)"
fm.model_description(f)

Is this the syntax for a fully specified model in formulae? Output:

ModelTerms(
  ResponseTerm(
    Term(
      name= choice
      variable= choice
      kind= None
    )
  ),
  InterceptTerm(),
  Term(
    name= stimulus
    variable= stimulus
    kind= None
  ),
  Term(
    name= reward
    variable= reward
    kind= None
  ),
  InteractionTerm(
    name= stimulus:reward
    variables= {'stimulus', 'reward'}
  ),
  GroupSpecTerm(
    expr= InterceptTerm(),
    factor= Term(
      name= subject_id
      variable= subject_id
      kind= None
    )
  ),
  GroupSpecTerm(
    expr= Term(
      name= stimulus
      variable= stimulus
      kind= None
    ),
    factor= Term(
      name= subject_id
      variable= subject_id
      kind= None
    )
  ),
  GroupSpecTerm(
    expr= Term(
      name= reward
      variable= reward
      kind= None
    ),
    factor= Term(
      name= subject_id
      variable= subject_id
      kind= None
    )
  ),
  GroupSpecTerm(
    expr= InteractionTerm(
      name= stimulus:reward
      variables= {'stimulus', 'reward'}
    ),
    factor= Term(
      name= subject_id
      variable= subject_id
      kind= None
    )
  )
)

Hi.

In Ch3 of the lme4 book manuscript, Bates showed two cases of random effects, one is correlated and the other is uncorrelated (independent):

The correlated one is from formula Reaction ~ 1 + Days + (1 + Days | Subject), the uncorrelated on is from Reaction ~ 1 + Days + (1|Subject) + (0 + Days|Subject).

I tried the same formulae with formulae, the two design matrices are the same:

from formulae import design_matrices
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

data = pd.read_csv('../datasets/sleepstudy.csv')

f1 = 'Reaction ~ 1 + Days + (1 + Days | Subject)'
f2 = 'Reaction ~ 1 + Days + (1|Subject) + (0 + Days|Subject)'

fig, ax = plt.subplots(2, 1, figsize=(12, 6))
for i, formula in enumerate([f1, f2]):
    
    dm = design_matrices(formula, data)
    terms = dm.group.terms_info.keys()

    Zs = [dm.group[key] for key in dm.group.terms_info.keys()]
    Z = np.hstack(Zs)
    Z[Z ==0] = np.nan
    ax[i].imshow(Z.T)
    ax[i].set_title(formula)

Is this an intended behavior for formulae or is it a bug?

The source tarball is broken b/c a required file in setup.py is not shipped with the source. This PR adds a MANIFEST.in file to ensure the necessary files will be packaged.

See https://github.com/conda-forge/staged-recipes/pull/14193 for an example of the problem with the current tarball.

This PR attempts to improve existing online documentation. Many classes, methods, and functions are currently documented, but this documentation is not shown on the webpage. The general aim is to have three broad sections.

• Examples: A collection of examples using formulae to obtain design matrices and fit statistical models.
• API Reference: The documentation for the objects that we expect the users to interact with.
• Internals: A documentation made for people wanting to use formulae for their own libraries.

These tasks may be split into several PRs. The goal of this PR is to improve the Internals reference. Examples will come later.

This PR adds a transformation for survival type responses. It is similar to R's Surv.

Surv(time, status) ~ returns a 2xn matrix with the time and the event status. These can then be used in the model to handle censored and uncensored data differently.

#79

• Add general contrasts so we can use other encodings apart from reference encoding.
• Stop reshaping 1d arrays to column vectors. It's not really necessary.

Adds an equivalent of the bs() function in the splines package in R. This is a stateful transformation by the way (it remember knots, bounds, and whether it has intercept or not).

Using the engine dataset from the R package gamair, we have:

import pandas as pd
from formulae import  design_matrices

data = pd.DataFrame({
    "size": [
        2.78, 2.43, 2.32, 2.43, 2.98, 2.32, 2.32, 2.32, 2.32, 2.13, 
        2.13, 2.13, 2.98, 1.99, 1.99, 1.99, 1.78, 1.58, 1.42
    ],
    "wear": [
        3, 3.1, 4.8, 3.3, 2.8, 2.9, 3.8, 3, 2.7, 3.2, 2.4, 2.6, 1.7,
        2.6, 2.8, 2.4, 2.5, 4.2, 4
    ]
})
design_matrices("wear ~ bs(size, 4)", data).common.design_matrix

array([[1.        , 0.00498077, 0.09481583, 0.56163869, 0.33856471],
       [1.        , 0.10358454, 0.429405  , 0.46238083, 0.00462963],
       [1.        , 0.17899408, 0.48816568, 0.33284024, 0.        ],
       [1.        , 0.10358454, 0.429405  , 0.46238083, 0.00462963],
       [1.        , 0.        , 0.        , 0.        , 1.        ],
       [1.        , 0.17899408, 0.48816568, 0.33284024, 0.        ],
       [1.        , 0.17899408, 0.48816568, 0.33284024, 0.        ],
       [1.        , 0.17899408, 0.48816568, 0.33284024, 0.        ],
       [1.        , 0.17899408, 0.48816568, 0.33284024, 0.        ],
       [1.        , 0.36011331, 0.46706614, 0.16341177, 0.        ],
       [1.        , 0.36011331, 0.46706614, 0.16341177, 0.        ],
       [1.        , 0.36011331, 0.46706614, 0.16341177, 0.        ],
       [1.        , 0.        , 0.        , 0.        , 1.        ],
       [1.        , 0.48758651, 0.37856345, 0.08455375, 0.        ],
       [1.        , 0.48758651, 0.37856345, 0.08455375, 0.        ],
       [1.        , 0.48758651, 0.37856345, 0.08455375, 0.        ],
       [1.        , 0.56530178, 0.19739645, 0.02130178, 0.        ],
       [1.        , 0.39454797, 0.04771909, 0.00187011, 0.        ],
       [1.        , 0.        , 0.        , 0.        , 0.        ]])

EDIT

By mistake, I ended up including here commits related to other topics. Now, it is also possible to pass response levels with index notation without having to use an explicit string. This doesn't work when the level contains spaces (must use a string there) or when the level is a number (use binary() function there). Examples:

This works:

• "y['A'] ~ x1 + x2"
• "y[A] ~ x1 + x2"
• "y['My level contains spaces'] ~ x1 + x2"
• "binary(y, 2)" ~ x1 + x2"

This does not work

• "y[My level contains spaces] ~ x1 + x2"
• "y[1] ~ x1 + x2"

Originally, the group of a group-specific term could only be a single term, not an interaction. This PR aims to enable groups to be made of interactions. It looks like so far it is working, but there are some edge cases with the evaluation with new data that should be addressed.

This feature (problem) was requested (reported) here.

We fail at parsing and resolving terms that are nested function calls. All the following alternatives work without any problem

design_matrices("y ~ f(x1)", data)
design_matrices("y ~ f(arg = x1)", data)
design_matrices("y ~ f(arg = 'blabla')", data)
design_matrices("y ~ f(arg1= a, arg2 = "blabla")", data)

But nested function calls, like the following, don't work

design_matrices("y ~ f(g(x1))", data)
design_matrices("y ~ f(arg = g(x1))", data)
# ...etc

To solve this we need to update the Resolver.visitCallExpr() method as well as how the call term handles Expr.Call() instances, and maybe something else that I'm missing now.

It looks like you're doing a lot of low-level parser generator work by hand. Have you considered using an existing pure Python parser package like rply? We've had a lot of success with that package in hy.

Hello !

I tried to deepcopy a bmb.Model and find out that we face a cannot pickle '_thread.lock' object. After investigation, this is due to saving of thread locked objects within formulae.Environment objects.

Indeed many, locked variables are fetched while calling Environment.capture(). My dummy workaround was to just replace the return cls([frame.f_locals, frame.f_globals])

        if isinstance(env, cls):
            return env
        elif isinstance(env, numbers.Integral):
            depth = env + reference
        else:
            raise TypeError("'env' must be either an integer or an instance of Environment.")
        frame = inspect.currentframe()
        try:
            for _ in range(depth + 1):
                if frame is None:
                    raise ValueError("call-stack is not that deep!")
                frame = frame.f_back
            # return cls([frame.f_locals, frame.f_globals])
            return cls([])
        finally:
            del frame

I'm not experienced enough on Bambi formulae, a good workaround would be to keep only what is necessary for formulae to work in output of the capture() function

I hope it would be helpful

I just tried a model with a different train/test split for the fit() and predict() functionality and hit the following error in the predict() call in bambi and hit the following error:

https://github.com/bambinos/formulae/blob/b926883a4ab744367953095a8540c8bf64ec2449/formulae/terms/variable.py#L223-L226

For a lot of real world use cases it's quite common for unknown categories to show up at inference time (e.g. new users, new customers etc) and it would be great to try and handle this natively in the predict() call such that a 'default' prediction can be made based on the mean across all groups rather than a specific group.

Hi all,

Thank you for the work on this awesome package!

I am using bambi to fit a fairly complex model on a large dataset with a single group effect, which has many individuals. There are around 183,000 rows and 30,542 unique groups. Version 0.3.3 crashes Jupyter reliably when instantiating this design matrix. Interestingly, version 0.2.0 will instantiate it, with the caveat I can't include an interaction term I can in 0.3.3 within bambi (see bambi issue 495).

I've tried a few different approaches, and have found it will instantiate with about 25-50% of the data (using DataFrame.sample(frac=.25)), so it seems more an issue of sheer scale than anything else. I've also tried with Spyder, getting the same issue.

The code below will grab a modified dataset of the same size and structure I am using and set up the model design, which kills my kernel after a minute or two.

import pandas as pd
from formulae import design_matrices

trouble = pd.read_feather('https://osf.io/kw2xh/download')
md = design_matrices('response ~ bs(time, degree=2, df=10) * state + state * level * trait + (1|pid)', data=trouble)

0.2.0 will return this after a little while, however. Any help is greatly appreciated, hopefully this issue isn't localised to my machine!

Tests are quite complete, but they are messy and redundant as well. Also, there are a couple of things we're missing yet (see https://app.codecov.io/gh/bambinos/formulae/blob/master/formulae/terms/terms.py).