From b3dcd2b0c8517469d0a0b6d956db1bae6ab69c5a Mon Sep 17 00:00:00 2001 From: Quarto GHA Workflow Runner Date: Sun, 4 Feb 2024 15:22:59 +0000 Subject: [PATCH] Built site for gh-pages --- .nojekyll | 2 +- notebooks/getting_started.html | 33 ++++++++++++++++++++++++++++----- search.json | 2 +- 3 files changed, 30 insertions(+), 7 deletions(-) diff --git a/.nojekyll b/.nojekyll index fc9c6ab79..3b6d16051 100644 --- a/.nojekyll +++ b/.nojekyll @@ -1 +1 @@ -e1947167 \ No newline at end of file +7956d809 \ No newline at end of file diff --git a/notebooks/getting_started.html b/notebooks/getting_started.html index 097d3d1ff..56d25f15a 100644 --- a/notebooks/getting_started.html +++ b/notebooks/getting_started.html @@ -819,16 +819,12 @@

Families

The following table summarizes the currently available families and their associated links:

----- + @@ -836,136 +832,163 @@

Families

+ + + + + + + + + + + + + + + + + + + + + + + + + + +
Family name Response distribution Default linkExample notebook
asymmetriclaplace AsymmetricLaplace identityQuantile Regression
bernoulli Bernoulli logitLogistic Regression
beta Beta logitBeta Regression
beta_binomial BetaBinomial logitTo be added
binomial Binomial logitHierarchical Logistic Regression
categorical Categorical softmaxCategorical Regression
cumulative Cumulative logitOrdinal Models
dirichlet_multinomial DirichletMultinomial logitTo be added
exponential Exponential logSurvival Models
gamma Gamma inverseGamma Regression
gaussian Normal identityMultiple Linear Regression
hurdle_gamma HurdleGamma logTo be added
hurdle_lognormal HurdleLogNormal identityTo be added
hurdle_negativebinomial HurdleNegativeBinomial logTo be added
hurdle_poisson HurdlePoisson logHurdle Poisson Regression
multinomial Multinomial softmaxTo be added
negativebinomial NegativeBinomial logNegative Binomial Regression
laplace Laplace identityTo be added
poisson Poisson logGaussian Processes with a Poisson likelihood
sratio StoppingRatio logitOrdinal Models
t StudentT identityRobust Linear Regression
vonmises VonMises tan(x / 2)Circular Regression
wald InverseGaussian inverse squaredWald Regression
weibull Weibull logTo be added
zero_inflated_binomial ZeroInflatedBinomial logitTo be added
zero_inflated_negativebinomial ZeroInflatedNegativeBinomial logTo be added
zero_inflated_poisson ZeroInflatedPoisson logZero Inflated Poisson Regression
diff --git a/search.json b/search.json index 82de7ee51..18f72fca0 100644 --- a/search.json +++ b/search.json @@ -235,7 +235,7 @@ "href": "notebooks/getting_started.html", "title": "Bambi", "section": "", - "text": "Bambi requires a working Python interpreter (3.8+). We recommend installing Python and key numerical libraries using the Anaconda Distribution, which has one-click installers available on all major platforms.\nAssuming a standard Python environment is installed on your machine (including pip), Bambi itself can be installed in one line using pip:\npip install bambi\nAlternatively, if you want the bleeding edge version of the package, you can install from GitHub:\npip install git+https://github.com/bambinos/bambi.git\n\n\nSuppose we have data for a typical within-subjects psychology experiment with 2 experimental conditions. Stimuli are nested within condition, and subjects are crossed with condition. We want to fit a model predicting reaction time (RT) from the common effect of condition, group specific intercepts for subjects, group specific condition slopes for students, and group specific intercepts for stimuli. Using Bambi we can fit this model and summarize its results as follows:\nimport bambi as bmb\n\n# Assume we already have our data loaded as a pandas DataFrame\nmodel = bmb.Model(\"rt ~ condition + (condition|subject) + (1|stimulus)\", data)\nresults = model.fit(draws=5000, chains=2)\naz.plot_trace(results)\naz.summary(results)\n\n\n\n\n\n\nimport warnings\n\nwarnings.filterwarnings(\"ignore\", category=FutureWarning)\n\n\nimport arviz as az\nimport bambi as bmb\nimport numpy as np\nimport pandas as pd\n\n\naz.style.use(\"arviz-darkgrid\")\n\n\n\n\nCreating a new model in Bambi is simple:\n\n# Read in a tab-delimited file containing our data\ndata = pd.read_table(\"data/my_data.txt\", sep=\"\\t\")\n\n# Initialize the model\nmodel = bmb.Model(\"y ~ x + z\", data)\n\n# Inspect model object\nmodel\n\n Formula: y ~ x + z\n Family: gaussian\n Link: mu = identity\n Observations: 50\n Priors: \n target = mu\n Common-level effects\n Intercept ~ Normal(mu: 0.1852, sigma: 2.5649)\n x ~ Normal(mu: 0.0, sigma: 2.231)\n z ~ Normal(mu: 0.0, sigma: 2.4374)\n \n Auxiliary parameters\n sigma ~ HalfStudentT(nu: 4.0, sigma: 1.013)\n\n\nTypically, we will initialize a Bambi Model by passing it a model formula and a pandas DataFrame. Other arguments such as family, priors, and link are available. By default, it uses family=\"gaussian\" which implies a linear regression with normal error. We get back a model that we can immediately fit by calling model.fit().\n\n\n\nAs with most mixed effect modeling packages, Bambi expects data in “long” format–meaning that each row should reflects a single observation at the most fine-grained level of analysis. For example, given a model where students are nested into classrooms and classrooms are nested into schools, we would want data with the following kind of structure:\n\n\n\n\nstudent\ngender\ngpa\nclass\nschool\n\n\n\n\n1\nF\n3.4\n1\n1\n\n\n2\nF\n3.7\n1\n1\n\n\n3\nM\n2.2\n1\n1\n\n\n4\nF\n3.9\n2\n1\n\n\n5\nM\n3.6\n2\n1\n\n\n6\nM\n3.5\n2\n1\n\n\n7\nF\n2.8\n3\n2\n\n\n8\nM\n3.9\n3\n2\n\n\n9\nF\n4.0\n3\n2\n\n\n\n\n\n\n\n\nModels are specified in Bambi using a formula-based syntax similar to what one might find in R packages like lme4 or brms using the Python formulae library. A couple of examples illustrate the breadth of models that can be easily specified in Bambi:\n\ndata = pd.read_csv(\"data/rrr_long.csv\")\ndata.head(10)\n\n\n\n\n\n \n \n \n uid\n condition\n gender\n age\n study\n self_perf\n stimulus\n value\n \n \n \n \n 0\n 1.0\n 0.0\n 1.0\n 24.0\n 0.0\n 8.0\n rating_c1\n 3.0\n \n \n 1\n 2.0\n 1.0\n 0.0\n 27.0\n 0.0\n 9.0\n rating_c1\n 7.0\n \n \n 2\n 3.0\n 0.0\n 1.0\n 25.0\n 0.0\n 3.0\n rating_c1\n 5.0\n \n \n 3\n 5.0\n 0.0\n 1.0\n 20.0\n 0.0\n 3.0\n rating_c1\n 7.0\n \n \n 4\n 8.0\n 1.0\n 1.0\n 19.0\n 0.0\n 6.0\n rating_c1\n 6.0\n \n \n 5\n 9.0\n 0.0\n 1.0\n 22.0\n 0.0\n 3.0\n rating_c1\n 6.0\n \n \n 6\n 10.0\n 1.0\n 1.0\n 49.0\n 0.0\n 4.0\n rating_c1\n 6.0\n \n \n 7\n 11.0\n 0.0\n 0.0\n 24.0\n 0.0\n 5.0\n rating_c1\n 7.0\n \n \n 8\n 12.0\n 1.0\n 0.0\n 26.0\n 0.0\n 6.0\n rating_c1\n 2.0\n \n \n 9\n 13.0\n 0.0\n 1.0\n 23.0\n 0.0\n 7.0\n rating_c1\n 1.0\n \n \n\n\n\n\n\n# Number of rows with missing values\ndata.isna().any(axis=1).sum()\n\n401\n\n\nWe pass dropna=True to tell Bambi to drop rows containing missing values. The number of rows dropped is different from the number of rows with missing values because Bambi only considers columns involved in the model.\n\n# Common (or fixed) effects only\nbmb.Model(\"value ~ condition + age + gender\", data, dropna=True)\n\nAutomatically removing 33/6940 rows from the dataset.\n\n\n Formula: value ~ condition + age + gender\n Family: gaussian\n Link: mu = identity\n Observations: 6907\n Priors: \n target = mu\n Common-level effects\n Intercept ~ Normal(mu: 4.5457, sigma: 28.4114)\n condition ~ Normal(mu: 0.0, sigma: 12.0966)\n age ~ Normal(mu: 0.0, sigma: 1.3011)\n gender ~ Normal(mu: 0.0, sigma: 13.1286)\n \n Auxiliary parameters\n sigma ~ HalfStudentT(nu: 4.0, sigma: 2.4186)\n\n\n\n# Common effects and group specific (or random) intercepts for subject\nbmb.Model(\"value ~ condition + age + gender + (1|uid)\", data, dropna=True)\n\nAutomatically removing 33/6940 rows from the dataset.\n\n\n Formula: value ~ condition + age + gender + (1|uid)\n Family: gaussian\n Link: mu = identity\n Observations: 6907\n Priors: \n target = mu\n Common-level effects\n Intercept ~ Normal(mu: 4.5457, sigma: 28.4114)\n condition ~ Normal(mu: 0.0, sigma: 12.0966)\n age ~ Normal(mu: 0.0, sigma: 1.3011)\n gender ~ Normal(mu: 0.0, sigma: 13.1286)\n \n Group-level effects\n 1|uid ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 28.4114))\n \n Auxiliary parameters\n sigma ~ HalfStudentT(nu: 4.0, sigma: 2.4186)\n\n\n\n# Multiple, complex group specific effects with both\n# group specific slopes and group specific intercepts\nbmb.Model(\n \"value ~ condition + age + gender + (1|uid) + (condition|study) + (condition|stimulus)\", \n data, \n dropna=True\n)\n\nAutomatically removing 33/6940 rows from the dataset.\n\n\n Formula: value ~ condition + age + gender + (1|uid) + (condition|study) + (condition|stimulus)\n Family: gaussian\n Link: mu = identity\n Observations: 6907\n Priors: \n target = mu\n Common-level effects\n Intercept ~ Normal(mu: 4.5457, sigma: 28.4114)\n condition ~ Normal(mu: 0.0, sigma: 12.0966)\n age ~ Normal(mu: 0.0, sigma: 1.3011)\n gender ~ Normal(mu: 0.0, sigma: 13.1286)\n \n Group-level effects\n 1|uid ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 28.4114))\n 1|study ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 28.4114))\n condition|study ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 12.0966))\n 1|stimulus ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 28.4114))\n condition|stimulus ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 12.0966))\n \n Auxiliary parameters\n sigma ~ HalfStudentT(nu: 4.0, sigma: 2.4186)\n\n\nEach of the above examples specifies a full model that can be fitted using PyMC by doing\nresults = model.fit()\n\n\nWhen a categorical common effect with N levels is added to a model, by default, it is coded by N-1 dummy variables (i.e., reduced-rank coding). For example, suppose we write \"y ~ condition + age + gender\", where condition is a categorical variable with 4 levels, and age and gender are continuous variables. Then our model would contain an intercept term (added to the model by default, as in R), three dummy-coded variables (each contrasting the first level of condition with one of the subsequent levels), and continuous predictors for age and gender. Suppose, however, that we would rather use full-rank coding of conditions. If we explicitly remove the intercept –as in \"y ~ 0 + condition + age + gender\"– then we get the desired effect. Now, the intercept is no longer included, and condition will be coded using 4 dummy indicators, each one coding for the presence or absence of the respective condition without reference to the other conditions.\nGroup specific effects are handled in a comparable way. When adding group specific intercepts, coding is always full-rank (e.g., when adding group specific intercepts for 100 schools, one gets 100 dummy-coded indicators coding each school separately, and not 99 indicators contrasting each school with the very first one). For group specific slopes, coding proceeds the same way as for common effects. The group specific effects specification \"(condition|subject)\" would add an intercept for each subject, plus N-1 condition slopes (each coded with respect to the first, omitted, level as the referent). If we instead specify \"(0+condition|subject)\", we get N condition slopes and no intercepts.\n\n\n\nOnce a model is fully specified, we need to run the PyMC sampler to generate parameter estimates. If we’re using the one-line fit() interface, sampling will begin right away:\n\nmodel = bmb.Model(\"value ~ condition + age + gender + (1|uid)\", data, dropna=True)\nresults = model.fit()\n\nAutomatically removing 33/6940 rows from the dataset.\nAuto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\n\n\nMultiprocess sampling (2 chains in 2 jobs)\nNUTS: [value_sigma, Intercept, condition, age, gender, 1|uid_sigma, 1|uid_offset]\n\n\n\n\n\n\n\n \n \n 100.00% [4000/4000 00:34<00:00 Sampling 2 chains, 0 divergences]\n \n \n\n\nSampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 35 seconds.\nWe recommend running at least 4 chains for robust computation of convergence diagnostics\n\n\nThe above code obtains 1,000 draws (the default value) and return them as an InferenceData instance.\n\n\n\n\n\n\nTip\n\n\n\nInferenceData is a rich data structure to store and manipulate data such as posterior samples, prior/posterior predictive samples, observations, etc. It is based on xarray, a library offering N-dimensional labeled arrays (you can think of it as a generalization of both Numpy arrays and Pandas dataframes). To learn how to perform common operations with InferenceData, like indexing, selection etc please check this and for details of the InferenceData Schema see this specification.\n\n\nIn this case, the fit() method accepts optional keyword arguments to pass onto PyMC’s sample() method, so any methods accepted by sample() can be specified here. We can also explicitly set the number of draws via the draws argument. For example, if we call fit(draws=2000, chains=2), the PyMC sampler will sample two chains in parallel, drawing 2,000 draws for each one. We could also specify starting parameter values, the step function to use, and so on (for full details, see the PyMC documentation).\nAlternatively, we can build a model, but not fit it.\n\nmodel = bmb.Model(\"value ~ condition + age + gender + (1|uid)\", data, dropna=True)\nmodel.build()\n\nAutomatically removing 33/6940 rows from the dataset.\n\n\nBuilding without sampling can be useful if we want to inspect the internal PyMC model before we start the (potentially long) sampling process. Once we’re satisfied, and wish to run the sampler, we can then simply call model.fit(), and the sampler will start running. Another good reason to build a model is to generate plot of the marginal priors using model.plot_priors().\n\nmodel.plot_priors();\n\nSampling: [1|uid_sigma, Intercept, age, condition, gender, value_sigma]\n\n\n\n\n\n\n\n\n\nBayesian inference requires one to specify prior probability distributions that represent the analyst’s belief (in advance of seeing the data) about the likely values of the model parameters. In practice, analysts often lack sufficient information to formulate well-defined priors, and instead opt to use “weakly informative” priors that mainly serve to keep the model from exploring completely pathological parts of the parameter space (e.g., when defining a prior on the distribution of human heights, a value of 3,000 cms should be assigned a probability of exactly 0).\nBy default, Bambi will intelligently generate weakly informative priors for all model terms, by loosely scaling them to the observed data. Currently, Bambi uses a methodology very similar to the one described in the documentation of the R package rstanarm. While the default priors will behave well in most typical settings, there are many cases where an analyst will want to specify their own priors–and in general, when informative priors are available, it’s a good idea to use them.\nFortunately, Bambi is built on top of PyMC, which means that we can seamlessly use any of the over 40 Distribution classes defined in PyMC. We can specify such priors in Bambi using the Prior class, which initializes with a name argument (which must map on exactly to the name of a valid PyMC Distribution) followed by any of the parameters accepted by the corresponding distribution. For example:\n\n# A Laplace prior with mean of 0 and scale given by a HalfNormal with a scale of 1\nmy_favorite_prior = bmb.Prior(\"Laplace\", mu=0, b=bmb.Prior(\"HalfNormal\", sigma=1))\n\n# Set the prior when adding a term to the model; more details on this below.\npriors = {\"1|uid\": my_favorite_prior}\nbmb.Model(\"value ~ condition + (1|uid)\", data, priors=priors, dropna=True)\n\nAutomatically removing 9/6940 rows from the dataset.\n\n\n Formula: value ~ condition + (1|uid)\n Family: gaussian\n Link: mu = identity\n Observations: 6931\n Priors: \n target = mu\n Common-level effects\n Intercept ~ Normal(mu: 4.5516, sigma: 8.4548)\n condition ~ Normal(mu: 0.0, sigma: 12.1019)\n \n Group-level effects\n 1|uid ~ Laplace(mu: 0.0, b: HalfNormal(sigma: 1.0))\n \n Auxiliary parameters\n sigma ~ HalfStudentT(nu: 4.0, sigma: 2.4197)\n\n\nPriors specified using the Prior class can be nested to arbitrary depths–meaning, we can set any of a given prior’s argument to point to another Prior instance. This is particularly useful when specifying hierarchical priors on group specific effects, where the individual group specific slopes or intercepts are constrained to share a common source distribution:\n\nsubject_sd = bmb.Prior(\"HalfCauchy\", beta=5)\nsubject_prior = bmb.Prior(\"Normal\", mu=0, sd=subject_sd)\npriors = {\"1|uid\": subject_prior}\nbmb.Model(\"value ~ condition + (1|uid)\", data, priors=priors, dropna=True)\n\nAutomatically removing 9/6940 rows from the dataset.\n\n\n Formula: value ~ condition + (1|uid)\n Family: gaussian\n Link: mu = identity\n Observations: 6931\n Priors: \n target = mu\n Common-level effects\n Intercept ~ Normal(mu: 4.5516, sigma: 8.4548)\n condition ~ Normal(mu: 0.0, sigma: 12.1019)\n \n Group-level effects\n 1|uid ~ Normal(mu: 0.0, sd: HalfCauchy(beta: 5.0))\n \n Auxiliary parameters\n sigma ~ HalfStudentT(nu: 4.0, sigma: 2.4197)\n\n\nThe above prior specification indicates that the individual subject intercepts are to be treated as if they are randomly sampled from the same underlying normal distribution, where the variance of that normal distribution is parameterized by a separate hyperprior (a half-cauchy with beta = 5).\nIt’s important to note that explicitly setting priors by passing in Prior objects will disable Bambi’s default behavior of scaling priors to the data in order to ensure that they remain weakly informative. This means that if you specify your own prior, you have to be sure not only to specify the distribution you want, but also any relevant scale parameters. For example, the 0.5 in Prior(\"Normal\", mu=0, sd=0.5) will be specified on the scale of the data, not the bounded partial correlation scale that Bambi uses for default priors. This means that if your outcome variable has a mean value of 10,000 and a standard deviation of, say, 1,000, you could potentially have some problems getting the model to produce reasonable estimates, since from the perspective of the data, you’re specifying an extremely strong prior.\n\n\nBambi’s priors are a thin layer on top of PyMC distributions. If you want to ask for a prior distribution by name, it must be the name of a PyMC distribution. But sometimes we want to use more complex distributions as priors. For all those cases, Bambi’s Prior class allow users to pass a function that returns a distribution that will be used as the prior. See the following example:\ndef CustomPrior(name, *args, dims=None, **kwargs):\n return pm.Normal(name, *args, dims=dims, **kwargs)\n\npriors = {\"x\": Prior(\"CustomPrior\", mu=0, sigma=5, dist=CustomPrior)}\nmodel = Model(\"y ~ x\", data, priors=priors)\nThe example above is trival because it’s just a wrapper of the pm.Normal distribution. But we can use this pattern to construct more complex distributions, such as a Truncated Laplace distribution shown below.\ndef TruncatedLaplace(name, mu,b,lower,upper,*args, dims=None, **kwargs):\n lap_dist = pm.Laplace.dist(mu=mu, b=b)\n return pm.Truncated(name, lap_dist, lower=lower, upper=upper, *args, dims=dims, **kwargs)\nIn summary, custom priors allow for greater flexibility by combining existing PyMC distributions in different ways. If you need to use a distribution that’s not implemented in PyMC, please check the link for further details.\n\n\n\n\nBambi supports the construction of mixed models with non-normal response distributions (i.e., generalized linear mixed models, or GLMMs). GLMMs are specified in the same way as LMMs, except that the user must specify the distribution to use for the response, and (optionally) the link function with which to transform the linear model prediction into the desired non-normal response. The easiest way to construct a GLMM is to simple set the family when creating the model:\n\ndata = bmb.load_data(\"admissions\")\nmodel = bmb.Model(\"admit ~ gre + gpa + rank\", data, family=\"bernoulli\")\nresults = model.fit()\n\nModeling the probability that admit==1\nAuto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\nMultiprocess sampling (2 chains in 2 jobs)\nNUTS: [Intercept, gre, gpa, rank]\n\n\n\n\n\n\n\n \n \n 100.00% [4000/4000 00:07<00:00 Sampling 2 chains, 0 divergences]\n \n \n\n\nSampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 7 seconds.\nWe recommend running at least 4 chains for robust computation of convergence diagnostics\n\n\nIf no link argument is explicitly set (see below), the canonical link function (or an otherwise sensible default) will be used.\n\n\n\nFollowing the convention used in many R packages, the response distribution to use for a GLMM is specified in a Family class that indicates how the response variable is distributed, as well as the link function transforming the linear response to a non-linear one.\nThe following table summarizes the currently available families and their associated links:\n\n\n\n\n\n\n\n\n\nFamily name\nResponse distribution\nDefault link\n\n\n\n\nasymmetriclaplace\nAsymmetricLaplace\nidentity\n\n\nbernoulli\nBernoulli\nlogit\n\n\nbeta\nBeta\nlogit\n\n\nbeta_binomial\nBetaBinomial\nlogit\n\n\nbinomial\nBinomial\nlogit\n\n\ncategorical\nCategorical\nsoftmax\n\n\ncumulative\nCumulative\nlogit\n\n\ndirichlet_multinomial\nDirichletMultinomial\nlogit\n\n\nexponential\nExponential\nlog\n\n\ngamma\nGamma\ninverse\n\n\ngaussian\nNormal\nidentity\n\n\nhurdle_gamma\nHurdleGamma\nlog\n\n\nhurdle_lognormal\nHurdleLogNormal\nidentity\n\n\nhurdle_negativebinomial\nHurdleNegativeBinomial\nlog\n\n\nhurdle_poisson\nHurdlePoisson\nlog\n\n\nmultinomial\nMultinomial\nsoftmax\n\n\nnegativebinomial\nNegativeBinomial\nlog\n\n\nlaplace\nLaplace\nidentity\n\n\npoisson\nPoisson\nlog\n\n\nsratio\nStoppingRatio\nlogit\n\n\nt\nStudentT\nidentity\n\n\nvonmises\nVonMises\ntan(x / 2)\n\n\nwald\nInverseGaussian\ninverse squared\n\n\nweibull\nWeibull\nlog\n\n\nzero_inflated_binomial\nZeroInflatedBinomial\nlogit\n\n\nzero_inflated_negativebinomial\nZeroInflatedNegativeBinomial\nlog\n\n\nzero_inflated_poisson\nZeroInflatedPoisson\nlog\n\n\n\n\nAlthough the easiest way to specify a family is by name, using one of the options listed in the table above, users can also create and use their own family, providing enormous flexibility. In the following example, we show how the built-in Bernoulli family could be constructed on-the-fly:\n\nfrom scipy import special\n\n# Construct likelihood distribution ------------------------------\n# This must use a valid PyMC distribution name.\n# 'parent' is the name of the variable that represents the mean of the distribution. \n# The mean of the Bernoulli family is given by 'p'.\nlikelihood = bmb.Likelihood(\"Bernoulli\", parent=\"p\")\n\n# Set link function ----------------------------------------------\n# There are two alternative approaches.\n# 1. Pass a name that is known by Bambi\nlink = bmb.Link(\"logit\")\n\n# 2. Build everything from scratch\n# link: A function that maps the response to the linear predictor\n# linkinv: A function that maps the linear predictor to the response\n# linkinv_backend: A function that maps the linear predictor to the response\n# that works with PyTensor tensors.\n# bmb.math.sigmoid is a PyTensor tensor function wrapped by PyMC and Bambi \nlink = bmb.Link(\n \"my_logit\", \n link=special.logit,\n linkinv=special.expit,\n linkinv_backend=bmb.math.sigmoid\n)\n\n# Construct the family -------------------------------------------\n# Families are defined by a name, a Likelihood and a Link.\nfamily = bmb.Family(\"bernoulli\", likelihood, link)\n\n# Now it's business as usual\n# Note: 'gre' is standardized to mean 0 and variance 1\nmodel = bmb.Model(\"admit ~ scale(gre) + gpa + rank\", data, family=family)\nmodel.build()\nmodel.graph()\n\n\n\n\n\nresults = model.fit()\n\nAuto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\nMultiprocess sampling (2 chains in 2 jobs)\nNUTS: [Intercept, scale(gre), gpa, rank]\n\n\n\n\n\n\n\n \n \n 100.00% [4000/4000 00:03<00:00 Sampling 2 chains, 0 divergences]\n \n \n\n\nSampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 4 seconds.\nWe recommend running at least 4 chains for robust computation of convergence diagnostics\n\n\nThe above example produces results identical to simply setting family='bernoulli'.\nOne complication in specifying a custom Family is that one must pass both a link function and an inverse link function which must be able to operate over PyTensor tensors rather than numpy arrays, so you’ll probably need to rely on tensor operations provided in pytensor.tensor (many of which are also wrapped by PyMC) when defining a new link.\n\n\n\nWhen a model is fitted, it returns a InferenceData object containing data related to the model and the posterior. This object can be passed to many functions in ArviZ to obtain numerical and visuals diagnostics and plot in general.\n\n\n\nTo visualize a plot of the posterior estimates and sample traces for all parameters, simply pass the InferenceData object to the arviz function az._plot_trace:\n\naz.plot_trace(results, compact=False);\n\n\n\n\nMore details on this plot are available in the ArviZ documentation.\n\n\n\nIf you prefer numerical summaries of the posterior estimates, you can use the az.summary() function from ArviZ which provides a pandas DataFrame with some key summary and diagnostics info on the model parameters, such as the 94% highest posterior density intervals\n\naz.summary(results)\n\n\n\n\n\n \n \n \n mean\n sd\n hdi_3%\n hdi_97%\n mcse_mean\n mcse_sd\n ess_bulk\n ess_tail\n r_hat\n \n \n \n \n Intercept\n -2.129\n 1.139\n -4.216\n 0.117\n 0.024\n 0.017\n 2245.0\n 1512.0\n 1.0\n \n \n scale(gre)\n 0.268\n 0.126\n 0.024\n 0.496\n 0.003\n 0.002\n 2321.0\n 1658.0\n 1.0\n \n \n gpa\n 0.787\n 0.323\n 0.131\n 1.335\n 0.007\n 0.005\n 2119.0\n 1439.0\n 1.0\n \n \n rank\n -0.566\n 0.129\n -0.813\n -0.332\n 0.003\n 0.002\n 2584.0\n 1778.0\n 1.0\n \n \n\n\n\n\nIf you want to view summaries or plots for specific parameters, you can pass a list of its names:\n\n# show the names of all variables stored in the InferenceData object\nlist(results.posterior.data_vars)\n\n['Intercept', 'scale(gre)', 'gpa', 'rank']\n\n\nYou can find detailed, worked examples of fitting Bambi models and working with the results in the example notebooks here.\n\n\n\nBambi is just a high-level interface to PyMC. As such, Bambi internally stores virtually all objects generated by PyMC, making it easy for users to retrieve, inspect, and modify those objects. For example, the Model class created by PyMC (as opposed to the Bambi class of the same name) is accessible from model.backend.model.\n\ntype(model.backend.model)\n\npymc.model.core.Model\n\n\n\nmodel.backend.model\n\n\\[\n \\begin{array}{rcl}\n \\text{Intercept} &\\sim & \\operatorname{Normal}(0,~23.4)\\\\\\text{scale(gre)} &\\sim & \\operatorname{Normal}(0,~2.5)\\\\\\text{gpa} &\\sim & \\operatorname{Normal}(0,~6.58)\\\\\\text{rank} &\\sim & \\operatorname{Normal}(0,~2.65)\\\\\\text{admit} &\\sim & \\operatorname{Bernoulli}(f(\\text{Intercept},~\\text{rank},~\\text{gpa},~\\text{scale(gre)}))\n \\end{array}\n \\]\n\n\n\nmodel.backend.model.observed_RVs\n\n[admit ~ Bernoulli(f(Intercept, rank, gpa, scale(gre)))]\n\n\n\nmodel.backend.model.unobserved_RVs\n\n[Intercept ~ Normal(0, 23.4),\n scale(gre) ~ Normal(0, 2.5),\n gpa ~ Normal(0, 6.58),\n rank ~ Normal(0, 2.65)]\n\n\n\n%load_ext watermark\n%watermark -n -u -v -iv -w\n\nLast updated: Sat Nov 11 2023\n\nPython implementation: CPython\nPython version : 3.10.6\nIPython version : 8.5.0\n\narviz : 0.16.1\nscipy : 1.11.2\nbambi : 0.13.1.dev7+gc6e5dbb5.d20231111\nnumpy : 1.23.4\npandas: 2.1.0\n\nWatermark: 2.4.3" + "text": "Bambi requires a working Python interpreter (3.8+). We recommend installing Python and key numerical libraries using the Anaconda Distribution, which has one-click installers available on all major platforms.\nAssuming a standard Python environment is installed on your machine (including pip), Bambi itself can be installed in one line using pip:\npip install bambi\nAlternatively, if you want the bleeding edge version of the package, you can install from GitHub:\npip install git+https://github.com/bambinos/bambi.git\n\n\nSuppose we have data for a typical within-subjects psychology experiment with 2 experimental conditions. Stimuli are nested within condition, and subjects are crossed with condition. We want to fit a model predicting reaction time (RT) from the common effect of condition, group specific intercepts for subjects, group specific condition slopes for students, and group specific intercepts for stimuli. Using Bambi we can fit this model and summarize its results as follows:\nimport bambi as bmb\n\n# Assume we already have our data loaded as a pandas DataFrame\nmodel = bmb.Model(\"rt ~ condition + (condition|subject) + (1|stimulus)\", data)\nresults = model.fit(draws=5000, chains=2)\naz.plot_trace(results)\naz.summary(results)\n\n\n\n\n\n\nimport warnings\n\nwarnings.filterwarnings(\"ignore\", category=FutureWarning)\n\n\nimport arviz as az\nimport bambi as bmb\nimport numpy as np\nimport pandas as pd\n\n\naz.style.use(\"arviz-darkgrid\")\n\n\n\n\nCreating a new model in Bambi is simple:\n\n# Read in a tab-delimited file containing our data\ndata = pd.read_table(\"data/my_data.txt\", sep=\"\\t\")\n\n# Initialize the model\nmodel = bmb.Model(\"y ~ x + z\", data)\n\n# Inspect model object\nmodel\n\n Formula: y ~ x + z\n Family: gaussian\n Link: mu = identity\n Observations: 50\n Priors: \n target = mu\n Common-level effects\n Intercept ~ Normal(mu: 0.1852, sigma: 2.5649)\n x ~ Normal(mu: 0.0, sigma: 2.231)\n z ~ Normal(mu: 0.0, sigma: 2.4374)\n \n Auxiliary parameters\n sigma ~ HalfStudentT(nu: 4.0, sigma: 1.013)\n\n\nTypically, we will initialize a Bambi Model by passing it a model formula and a pandas DataFrame. Other arguments such as family, priors, and link are available. By default, it uses family=\"gaussian\" which implies a linear regression with normal error. We get back a model that we can immediately fit by calling model.fit().\n\n\n\nAs with most mixed effect modeling packages, Bambi expects data in “long” format–meaning that each row should reflects a single observation at the most fine-grained level of analysis. For example, given a model where students are nested into classrooms and classrooms are nested into schools, we would want data with the following kind of structure:\n\n\n\n\nstudent\ngender\ngpa\nclass\nschool\n\n\n\n\n1\nF\n3.4\n1\n1\n\n\n2\nF\n3.7\n1\n1\n\n\n3\nM\n2.2\n1\n1\n\n\n4\nF\n3.9\n2\n1\n\n\n5\nM\n3.6\n2\n1\n\n\n6\nM\n3.5\n2\n1\n\n\n7\nF\n2.8\n3\n2\n\n\n8\nM\n3.9\n3\n2\n\n\n9\nF\n4.0\n3\n2\n\n\n\n\n\n\n\n\nModels are specified in Bambi using a formula-based syntax similar to what one might find in R packages like lme4 or brms using the Python formulae library. A couple of examples illustrate the breadth of models that can be easily specified in Bambi:\n\ndata = pd.read_csv(\"data/rrr_long.csv\")\ndata.head(10)\n\n\n\n\n\n \n \n \n uid\n condition\n gender\n age\n study\n self_perf\n stimulus\n value\n \n \n \n \n 0\n 1.0\n 0.0\n 1.0\n 24.0\n 0.0\n 8.0\n rating_c1\n 3.0\n \n \n 1\n 2.0\n 1.0\n 0.0\n 27.0\n 0.0\n 9.0\n rating_c1\n 7.0\n \n \n 2\n 3.0\n 0.0\n 1.0\n 25.0\n 0.0\n 3.0\n rating_c1\n 5.0\n \n \n 3\n 5.0\n 0.0\n 1.0\n 20.0\n 0.0\n 3.0\n rating_c1\n 7.0\n \n \n 4\n 8.0\n 1.0\n 1.0\n 19.0\n 0.0\n 6.0\n rating_c1\n 6.0\n \n \n 5\n 9.0\n 0.0\n 1.0\n 22.0\n 0.0\n 3.0\n rating_c1\n 6.0\n \n \n 6\n 10.0\n 1.0\n 1.0\n 49.0\n 0.0\n 4.0\n rating_c1\n 6.0\n \n \n 7\n 11.0\n 0.0\n 0.0\n 24.0\n 0.0\n 5.0\n rating_c1\n 7.0\n \n \n 8\n 12.0\n 1.0\n 0.0\n 26.0\n 0.0\n 6.0\n rating_c1\n 2.0\n \n \n 9\n 13.0\n 0.0\n 1.0\n 23.0\n 0.0\n 7.0\n rating_c1\n 1.0\n \n \n\n\n\n\n\n# Number of rows with missing values\ndata.isna().any(axis=1).sum()\n\n401\n\n\nWe pass dropna=True to tell Bambi to drop rows containing missing values. The number of rows dropped is different from the number of rows with missing values because Bambi only considers columns involved in the model.\n\n# Common (or fixed) effects only\nbmb.Model(\"value ~ condition + age + gender\", data, dropna=True)\n\nAutomatically removing 33/6940 rows from the dataset.\n\n\n Formula: value ~ condition + age + gender\n Family: gaussian\n Link: mu = identity\n Observations: 6907\n Priors: \n target = mu\n Common-level effects\n Intercept ~ Normal(mu: 4.5457, sigma: 28.4114)\n condition ~ Normal(mu: 0.0, sigma: 12.0966)\n age ~ Normal(mu: 0.0, sigma: 1.3011)\n gender ~ Normal(mu: 0.0, sigma: 13.1286)\n \n Auxiliary parameters\n sigma ~ HalfStudentT(nu: 4.0, sigma: 2.4186)\n\n\n\n# Common effects and group specific (or random) intercepts for subject\nbmb.Model(\"value ~ condition + age + gender + (1|uid)\", data, dropna=True)\n\nAutomatically removing 33/6940 rows from the dataset.\n\n\n Formula: value ~ condition + age + gender + (1|uid)\n Family: gaussian\n Link: mu = identity\n Observations: 6907\n Priors: \n target = mu\n Common-level effects\n Intercept ~ Normal(mu: 4.5457, sigma: 28.4114)\n condition ~ Normal(mu: 0.0, sigma: 12.0966)\n age ~ Normal(mu: 0.0, sigma: 1.3011)\n gender ~ Normal(mu: 0.0, sigma: 13.1286)\n \n Group-level effects\n 1|uid ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 28.4114))\n \n Auxiliary parameters\n sigma ~ HalfStudentT(nu: 4.0, sigma: 2.4186)\n\n\n\n# Multiple, complex group specific effects with both\n# group specific slopes and group specific intercepts\nbmb.Model(\n \"value ~ condition + age + gender + (1|uid) + (condition|study) + (condition|stimulus)\", \n data, \n dropna=True\n)\n\nAutomatically removing 33/6940 rows from the dataset.\n\n\n Formula: value ~ condition + age + gender + (1|uid) + (condition|study) + (condition|stimulus)\n Family: gaussian\n Link: mu = identity\n Observations: 6907\n Priors: \n target = mu\n Common-level effects\n Intercept ~ Normal(mu: 4.5457, sigma: 28.4114)\n condition ~ Normal(mu: 0.0, sigma: 12.0966)\n age ~ Normal(mu: 0.0, sigma: 1.3011)\n gender ~ Normal(mu: 0.0, sigma: 13.1286)\n \n Group-level effects\n 1|uid ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 28.4114))\n 1|study ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 28.4114))\n condition|study ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 12.0966))\n 1|stimulus ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 28.4114))\n condition|stimulus ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 12.0966))\n \n Auxiliary parameters\n sigma ~ HalfStudentT(nu: 4.0, sigma: 2.4186)\n\n\nEach of the above examples specifies a full model that can be fitted using PyMC by doing\nresults = model.fit()\n\n\nWhen a categorical common effect with N levels is added to a model, by default, it is coded by N-1 dummy variables (i.e., reduced-rank coding). For example, suppose we write \"y ~ condition + age + gender\", where condition is a categorical variable with 4 levels, and age and gender are continuous variables. Then our model would contain an intercept term (added to the model by default, as in R), three dummy-coded variables (each contrasting the first level of condition with one of the subsequent levels), and continuous predictors for age and gender. Suppose, however, that we would rather use full-rank coding of conditions. If we explicitly remove the intercept –as in \"y ~ 0 + condition + age + gender\"– then we get the desired effect. Now, the intercept is no longer included, and condition will be coded using 4 dummy indicators, each one coding for the presence or absence of the respective condition without reference to the other conditions.\nGroup specific effects are handled in a comparable way. When adding group specific intercepts, coding is always full-rank (e.g., when adding group specific intercepts for 100 schools, one gets 100 dummy-coded indicators coding each school separately, and not 99 indicators contrasting each school with the very first one). For group specific slopes, coding proceeds the same way as for common effects. The group specific effects specification \"(condition|subject)\" would add an intercept for each subject, plus N-1 condition slopes (each coded with respect to the first, omitted, level as the referent). If we instead specify \"(0+condition|subject)\", we get N condition slopes and no intercepts.\n\n\n\nOnce a model is fully specified, we need to run the PyMC sampler to generate parameter estimates. If we’re using the one-line fit() interface, sampling will begin right away:\n\nmodel = bmb.Model(\"value ~ condition + age + gender + (1|uid)\", data, dropna=True)\nresults = model.fit()\n\nAutomatically removing 33/6940 rows from the dataset.\nAuto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\n\n\nMultiprocess sampling (2 chains in 2 jobs)\nNUTS: [value_sigma, Intercept, condition, age, gender, 1|uid_sigma, 1|uid_offset]\n\n\n\n\n\n\n\n \n \n 100.00% [4000/4000 00:34<00:00 Sampling 2 chains, 0 divergences]\n \n \n\n\nSampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 35 seconds.\nWe recommend running at least 4 chains for robust computation of convergence diagnostics\n\n\nThe above code obtains 1,000 draws (the default value) and return them as an InferenceData instance.\n\n\n\n\n\n\nTip\n\n\n\nInferenceData is a rich data structure to store and manipulate data such as posterior samples, prior/posterior predictive samples, observations, etc. It is based on xarray, a library offering N-dimensional labeled arrays (you can think of it as a generalization of both Numpy arrays and Pandas dataframes). To learn how to perform common operations with InferenceData, like indexing, selection etc please check this and for details of the InferenceData Schema see this specification.\n\n\nIn this case, the fit() method accepts optional keyword arguments to pass onto PyMC’s sample() method, so any methods accepted by sample() can be specified here. We can also explicitly set the number of draws via the draws argument. For example, if we call fit(draws=2000, chains=2), the PyMC sampler will sample two chains in parallel, drawing 2,000 draws for each one. We could also specify starting parameter values, the step function to use, and so on (for full details, see the PyMC documentation).\nAlternatively, we can build a model, but not fit it.\n\nmodel = bmb.Model(\"value ~ condition + age + gender + (1|uid)\", data, dropna=True)\nmodel.build()\n\nAutomatically removing 33/6940 rows from the dataset.\n\n\nBuilding without sampling can be useful if we want to inspect the internal PyMC model before we start the (potentially long) sampling process. Once we’re satisfied, and wish to run the sampler, we can then simply call model.fit(), and the sampler will start running. Another good reason to build a model is to generate plot of the marginal priors using model.plot_priors().\n\nmodel.plot_priors();\n\nSampling: [1|uid_sigma, Intercept, age, condition, gender, value_sigma]\n\n\n\n\n\n\n\n\n\nBayesian inference requires one to specify prior probability distributions that represent the analyst’s belief (in advance of seeing the data) about the likely values of the model parameters. In practice, analysts often lack sufficient information to formulate well-defined priors, and instead opt to use “weakly informative” priors that mainly serve to keep the model from exploring completely pathological parts of the parameter space (e.g., when defining a prior on the distribution of human heights, a value of 3,000 cms should be assigned a probability of exactly 0).\nBy default, Bambi will intelligently generate weakly informative priors for all model terms, by loosely scaling them to the observed data. Currently, Bambi uses a methodology very similar to the one described in the documentation of the R package rstanarm. While the default priors will behave well in most typical settings, there are many cases where an analyst will want to specify their own priors–and in general, when informative priors are available, it’s a good idea to use them.\nFortunately, Bambi is built on top of PyMC, which means that we can seamlessly use any of the over 40 Distribution classes defined in PyMC. We can specify such priors in Bambi using the Prior class, which initializes with a name argument (which must map on exactly to the name of a valid PyMC Distribution) followed by any of the parameters accepted by the corresponding distribution. For example:\n\n# A Laplace prior with mean of 0 and scale given by a HalfNormal with a scale of 1\nmy_favorite_prior = bmb.Prior(\"Laplace\", mu=0, b=bmb.Prior(\"HalfNormal\", sigma=1))\n\n# Set the prior when adding a term to the model; more details on this below.\npriors = {\"1|uid\": my_favorite_prior}\nbmb.Model(\"value ~ condition + (1|uid)\", data, priors=priors, dropna=True)\n\nAutomatically removing 9/6940 rows from the dataset.\n\n\n Formula: value ~ condition + (1|uid)\n Family: gaussian\n Link: mu = identity\n Observations: 6931\n Priors: \n target = mu\n Common-level effects\n Intercept ~ Normal(mu: 4.5516, sigma: 8.4548)\n condition ~ Normal(mu: 0.0, sigma: 12.1019)\n \n Group-level effects\n 1|uid ~ Laplace(mu: 0.0, b: HalfNormal(sigma: 1.0))\n \n Auxiliary parameters\n sigma ~ HalfStudentT(nu: 4.0, sigma: 2.4197)\n\n\nPriors specified using the Prior class can be nested to arbitrary depths–meaning, we can set any of a given prior’s argument to point to another Prior instance. This is particularly useful when specifying hierarchical priors on group specific effects, where the individual group specific slopes or intercepts are constrained to share a common source distribution:\n\nsubject_sd = bmb.Prior(\"HalfCauchy\", beta=5)\nsubject_prior = bmb.Prior(\"Normal\", mu=0, sd=subject_sd)\npriors = {\"1|uid\": subject_prior}\nbmb.Model(\"value ~ condition + (1|uid)\", data, priors=priors, dropna=True)\n\nAutomatically removing 9/6940 rows from the dataset.\n\n\n Formula: value ~ condition + (1|uid)\n Family: gaussian\n Link: mu = identity\n Observations: 6931\n Priors: \n target = mu\n Common-level effects\n Intercept ~ Normal(mu: 4.5516, sigma: 8.4548)\n condition ~ Normal(mu: 0.0, sigma: 12.1019)\n \n Group-level effects\n 1|uid ~ Normal(mu: 0.0, sd: HalfCauchy(beta: 5.0))\n \n Auxiliary parameters\n sigma ~ HalfStudentT(nu: 4.0, sigma: 2.4197)\n\n\nThe above prior specification indicates that the individual subject intercepts are to be treated as if they are randomly sampled from the same underlying normal distribution, where the variance of that normal distribution is parameterized by a separate hyperprior (a half-cauchy with beta = 5).\nIt’s important to note that explicitly setting priors by passing in Prior objects will disable Bambi’s default behavior of scaling priors to the data in order to ensure that they remain weakly informative. This means that if you specify your own prior, you have to be sure not only to specify the distribution you want, but also any relevant scale parameters. For example, the 0.5 in Prior(\"Normal\", mu=0, sd=0.5) will be specified on the scale of the data, not the bounded partial correlation scale that Bambi uses for default priors. This means that if your outcome variable has a mean value of 10,000 and a standard deviation of, say, 1,000, you could potentially have some problems getting the model to produce reasonable estimates, since from the perspective of the data, you’re specifying an extremely strong prior.\n\n\nBambi’s priors are a thin layer on top of PyMC distributions. If you want to ask for a prior distribution by name, it must be the name of a PyMC distribution. But sometimes we want to use more complex distributions as priors. For all those cases, Bambi’s Prior class allow users to pass a function that returns a distribution that will be used as the prior. See the following example:\ndef CustomPrior(name, *args, dims=None, **kwargs):\n return pm.Normal(name, *args, dims=dims, **kwargs)\n\npriors = {\"x\": Prior(\"CustomPrior\", mu=0, sigma=5, dist=CustomPrior)}\nmodel = Model(\"y ~ x\", data, priors=priors)\nThe example above is trival because it’s just a wrapper of the pm.Normal distribution. But we can use this pattern to construct more complex distributions, such as a Truncated Laplace distribution shown below.\ndef TruncatedLaplace(name, mu,b,lower,upper,*args, dims=None, **kwargs):\n lap_dist = pm.Laplace.dist(mu=mu, b=b)\n return pm.Truncated(name, lap_dist, lower=lower, upper=upper, *args, dims=dims, **kwargs)\nIn summary, custom priors allow for greater flexibility by combining existing PyMC distributions in different ways. If you need to use a distribution that’s not implemented in PyMC, please check the link for further details.\n\n\n\n\nBambi supports the construction of mixed models with non-normal response distributions (i.e., generalized linear mixed models, or GLMMs). GLMMs are specified in the same way as LMMs, except that the user must specify the distribution to use for the response, and (optionally) the link function with which to transform the linear model prediction into the desired non-normal response. The easiest way to construct a GLMM is to simple set the family when creating the model:\n\ndata = bmb.load_data(\"admissions\")\nmodel = bmb.Model(\"admit ~ gre + gpa + rank\", data, family=\"bernoulli\")\nresults = model.fit()\n\nModeling the probability that admit==1\nAuto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\nMultiprocess sampling (2 chains in 2 jobs)\nNUTS: [Intercept, gre, gpa, rank]\n\n\n\n\n\n\n\n \n \n 100.00% [4000/4000 00:07<00:00 Sampling 2 chains, 0 divergences]\n \n \n\n\nSampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 7 seconds.\nWe recommend running at least 4 chains for robust computation of convergence diagnostics\n\n\nIf no link argument is explicitly set (see below), the canonical link function (or an otherwise sensible default) will be used.\n\n\n\nFollowing the convention used in many R packages, the response distribution to use for a GLMM is specified in a Family class that indicates how the response variable is distributed, as well as the link function transforming the linear response to a non-linear one.\nThe following table summarizes the currently available families and their associated links:\n\n\n\n\nFamily name\nResponse distribution\nDefault link\nExample notebook\n\n\n\n\nasymmetriclaplace\nAsymmetricLaplace\nidentity\nQuantile Regression\n\n\nbernoulli\nBernoulli\nlogit\nLogistic Regression\n\n\nbeta\nBeta\nlogit\nBeta Regression\n\n\nbeta_binomial\nBetaBinomial\nlogit\nTo be added\n\n\nbinomial\nBinomial\nlogit\nHierarchical Logistic Regression\n\n\ncategorical\nCategorical\nsoftmax\nCategorical Regression\n\n\ncumulative\nCumulative\nlogit\nOrdinal Models\n\n\ndirichlet_multinomial\nDirichletMultinomial\nlogit\nTo be added\n\n\nexponential\nExponential\nlog\nSurvival Models\n\n\ngamma\nGamma\ninverse\nGamma Regression\n\n\ngaussian\nNormal\nidentity\nMultiple Linear Regression\n\n\nhurdle_gamma\nHurdleGamma\nlog\nTo be added\n\n\nhurdle_lognormal\nHurdleLogNormal\nidentity\nTo be added\n\n\nhurdle_negativebinomial\nHurdleNegativeBinomial\nlog\nTo be added\n\n\nhurdle_poisson\nHurdlePoisson\nlog\nHurdle Poisson Regression\n\n\nmultinomial\nMultinomial\nsoftmax\nTo be added\n\n\nnegativebinomial\nNegativeBinomial\nlog\nNegative Binomial Regression\n\n\nlaplace\nLaplace\nidentity\nTo be added\n\n\npoisson\nPoisson\nlog\nGaussian Processes with a Poisson likelihood\n\n\nsratio\nStoppingRatio\nlogit\nOrdinal Models\n\n\nt\nStudentT\nidentity\nRobust Linear Regression\n\n\nvonmises\nVonMises\ntan(x / 2)\nCircular Regression\n\n\nwald\nInverseGaussian\ninverse squared\nWald Regression\n\n\nweibull\nWeibull\nlog\nTo be added\n\n\nzero_inflated_binomial\nZeroInflatedBinomial\nlogit\nTo be added\n\n\nzero_inflated_negativebinomial\nZeroInflatedNegativeBinomial\nlog\nTo be added\n\n\nzero_inflated_poisson\nZeroInflatedPoisson\nlog\nZero Inflated Poisson Regression\n\n\n\n\nAlthough the easiest way to specify a family is by name, using one of the options listed in the table above, users can also create and use their own family, providing enormous flexibility. In the following example, we show how the built-in Bernoulli family could be constructed on-the-fly:\n\nfrom scipy import special\n\n# Construct likelihood distribution ------------------------------\n# This must use a valid PyMC distribution name.\n# 'parent' is the name of the variable that represents the mean of the distribution. \n# The mean of the Bernoulli family is given by 'p'.\nlikelihood = bmb.Likelihood(\"Bernoulli\", parent=\"p\")\n\n# Set link function ----------------------------------------------\n# There are two alternative approaches.\n# 1. Pass a name that is known by Bambi\nlink = bmb.Link(\"logit\")\n\n# 2. Build everything from scratch\n# link: A function that maps the response to the linear predictor\n# linkinv: A function that maps the linear predictor to the response\n# linkinv_backend: A function that maps the linear predictor to the response\n# that works with PyTensor tensors.\n# bmb.math.sigmoid is a PyTensor tensor function wrapped by PyMC and Bambi \nlink = bmb.Link(\n \"my_logit\", \n link=special.logit,\n linkinv=special.expit,\n linkinv_backend=bmb.math.sigmoid\n)\n\n# Construct the family -------------------------------------------\n# Families are defined by a name, a Likelihood and a Link.\nfamily = bmb.Family(\"bernoulli\", likelihood, link)\n\n# Now it's business as usual\n# Note: 'gre' is standardized to mean 0 and variance 1\nmodel = bmb.Model(\"admit ~ scale(gre) + gpa + rank\", data, family=family)\nmodel.build()\nmodel.graph()\n\n\n\n\n\nresults = model.fit()\n\nAuto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\nMultiprocess sampling (2 chains in 2 jobs)\nNUTS: [Intercept, scale(gre), gpa, rank]\n\n\n\n\n\n\n\n \n \n 100.00% [4000/4000 00:03<00:00 Sampling 2 chains, 0 divergences]\n \n \n\n\nSampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 4 seconds.\nWe recommend running at least 4 chains for robust computation of convergence diagnostics\n\n\nThe above example produces results identical to simply setting family='bernoulli'.\nOne complication in specifying a custom Family is that one must pass both a link function and an inverse link function which must be able to operate over PyTensor tensors rather than numpy arrays, so you’ll probably need to rely on tensor operations provided in pytensor.tensor (many of which are also wrapped by PyMC) when defining a new link.\n\n\n\nWhen a model is fitted, it returns a InferenceData object containing data related to the model and the posterior. This object can be passed to many functions in ArviZ to obtain numerical and visuals diagnostics and plot in general.\n\n\n\nTo visualize a plot of the posterior estimates and sample traces for all parameters, simply pass the InferenceData object to the arviz function az._plot_trace:\n\naz.plot_trace(results, compact=False);\n\n\n\n\nMore details on this plot are available in the ArviZ documentation.\n\n\n\nIf you prefer numerical summaries of the posterior estimates, you can use the az.summary() function from ArviZ which provides a pandas DataFrame with some key summary and diagnostics info on the model parameters, such as the 94% highest posterior density intervals\n\naz.summary(results)\n\n\n\n\n\n \n \n \n mean\n sd\n hdi_3%\n hdi_97%\n mcse_mean\n mcse_sd\n ess_bulk\n ess_tail\n r_hat\n \n \n \n \n Intercept\n -2.129\n 1.139\n -4.216\n 0.117\n 0.024\n 0.017\n 2245.0\n 1512.0\n 1.0\n \n \n scale(gre)\n 0.268\n 0.126\n 0.024\n 0.496\n 0.003\n 0.002\n 2321.0\n 1658.0\n 1.0\n \n \n gpa\n 0.787\n 0.323\n 0.131\n 1.335\n 0.007\n 0.005\n 2119.0\n 1439.0\n 1.0\n \n \n rank\n -0.566\n 0.129\n -0.813\n -0.332\n 0.003\n 0.002\n 2584.0\n 1778.0\n 1.0\n \n \n\n\n\n\nIf you want to view summaries or plots for specific parameters, you can pass a list of its names:\n\n# show the names of all variables stored in the InferenceData object\nlist(results.posterior.data_vars)\n\n['Intercept', 'scale(gre)', 'gpa', 'rank']\n\n\nYou can find detailed, worked examples of fitting Bambi models and working with the results in the example notebooks here.\n\n\n\nBambi is just a high-level interface to PyMC. As such, Bambi internally stores virtually all objects generated by PyMC, making it easy for users to retrieve, inspect, and modify those objects. For example, the Model class created by PyMC (as opposed to the Bambi class of the same name) is accessible from model.backend.model.\n\ntype(model.backend.model)\n\npymc.model.core.Model\n\n\n\nmodel.backend.model\n\n\\[\n \\begin{array}{rcl}\n \\text{Intercept} &\\sim & \\operatorname{Normal}(0,~23.4)\\\\\\text{scale(gre)} &\\sim & \\operatorname{Normal}(0,~2.5)\\\\\\text{gpa} &\\sim & \\operatorname{Normal}(0,~6.58)\\\\\\text{rank} &\\sim & \\operatorname{Normal}(0,~2.65)\\\\\\text{admit} &\\sim & \\operatorname{Bernoulli}(f(\\text{Intercept},~\\text{rank},~\\text{gpa},~\\text{scale(gre)}))\n \\end{array}\n \\]\n\n\n\nmodel.backend.model.observed_RVs\n\n[admit ~ Bernoulli(f(Intercept, rank, gpa, scale(gre)))]\n\n\n\nmodel.backend.model.unobserved_RVs\n\n[Intercept ~ Normal(0, 23.4),\n scale(gre) ~ Normal(0, 2.5),\n gpa ~ Normal(0, 6.58),\n rank ~ Normal(0, 2.65)]\n\n\n\n%load_ext watermark\n%watermark -n -u -v -iv -w\n\nLast updated: Sat Nov 11 2023\n\nPython implementation: CPython\nPython version : 3.10.6\nIPython version : 8.5.0\n\narviz : 0.16.1\nscipy : 1.11.2\nbambi : 0.13.1.dev7+gc6e5dbb5.d20231111\nnumpy : 1.23.4\npandas: 2.1.0\n\nWatermark: 2.4.3" }, { "objectID": "notebooks/hierarchical_binomial_bambi.html",