pymare.estimators.StanMetaRegression

class StanMetaRegression(tau_prior_scale=None, **sampling_kwargs)[source]

Bases: BaseEstimator

Bayesian meta-regression estimator using Stan.

Parameters:

  • tau_prior_scale (None or float, optional) – Scale of the half-normal prior on tau, the between-group standard deviation. When None (default), it is set to max(std(y), sqrt(mean(v))), the larger of the observed spread of the estimates and the typical sampling standard deviation.

  • **sampling_kwargs – Optional keyword arguments to pass on to CmdStanPy’s sampler (e.g., iter_sampling for the number of post-warmup draws per chain, chains, seed, adapt_delta).

Notes

The model is

y_i &\sim \mathcal{N}(x_i' \beta + \theta_{g(i)}, \sigma_i) \\ \theta_g &\sim \mathcal{N}(0, \tau)

where \sigma_i = \sqrt{v_i} is the known sampling standard deviation of observation i and g(i) is its group. This is the random-effects meta-analysis model of the Stan User’s Guide [1] with that guide’s stated extension to observation-level predictors. The reported tau2 is \tau^2, the between-group variance, matching what every other PyMARE estimator reports under that name.

theta is given a non-centered parameterization (theta = tau * theta_raw with theta_raw standard normal). The centered form produces the funnel geometry that dominates divergences in hierarchical models with few groups, which is this estimator’s principal use case.

\tau gets a half-normal prior, weakly informative per Stan’s recommendations [2] for models with few groups. Its scale is derived from the data rather than fixed, since a fixed scale would be crushingly informative on data measured in thousands and vacuous on data measured in thousandths.

The default scale is max(std(y), sqrt(mean(v))). \tau cannot plausibly exceed the spread of the estimates, and should not be presumed smaller than a typical standard error, so the larger of the two never asserts that \tau is small when either quantity says otherwise. Erring large is deliberate: too small a scale costs coverage, too large costs only precision in \tau^2. validation/stan records the measurements behind that choice. Pass tau_prior_scale to override it.

\beta keeps Stan’s implicit improper uniform prior, so under a diffuse prior on \tau the posterior means agree with VarianceBasedLikelihoodEstimator at method="ML".

The Stan program is compiled on first use and cached beside the installed source, so the cost is paid once per installation rather than per fit.

References

Changed in version 0.0.11:

  • The backend moved from PyStan 3 to CmdStanPy. PyStan’s sampler argument names (num_samples, num_warmup, num_chains, num_thin) are rejected with a message naming their replacements.

  • tau2 is now the between-group variance rather than its square root, and sampling variances are converted to standard deviations before being passed to Stan. Both were wrong before, so posterior estimates change.

  • groups accepts scalar labels of any type, not only integers in 1..k, and fit_dataset() now passes dataset.g rather than dropping it.

  • ci now sets the width of the reported credible interval. It was previously accepted and ignored.

compile(force=False)[source]

Compile the Stan model.

Parameters:

force (bool, optional) – Whether to recompile even when an up-to-date executable already exists. Default = False.

Returns:

The instance, so that compile() can be chained.

Return type:

StanMetaRegression

Notes

Called by fit() when needed, so it never has to be called directly. Calling it in advance is worthwhile when the same estimator will be fitted to several datasets, because the compiled executable does not depend on the data.

The executable is written beside the installed .stan file, where CmdStanPy finds and reuses it on later runs. If that directory is not writable – a read-only site-packages, for instance – it falls back to ~/.pymare/stan and warns once.

fit(y, v, X, groups=None)[source]

Run the Stan sampler and return results.

Parameters:

  • y (numpy.ndarray of shape (K,)) – 1d array of observation-level estimates

  • v (numpy.ndarray of shape (K,)) – 1d array of observation-level variances

  • X (numpy.ndarray of shape (K[, P])) – 1d or 2d array containing observation-level predictors (including intercept); has dimensions K x P, where K is the number of observations and P is the number of predictor variables.

  • groups (None or array-like of shape (K,), optional) – One scalar label per observation, identifying the groups of observations in the y/v/X inputs. Labels may be strings, integers or any other hashable that numpy stores as a single element, and need not be consecutive; they are encoded internally in order of first occurrence by encode_groups(). Composite labels such as tuples are not accepted, because numpy reads a sequence of them as a 2-dimensional array. When None (default), each observation in the inputs is treated as a separate group.

Returns:

The fitted instance.

Return type:

StanMetaRegression

Warns:

UserWarning – If the sampler reported divergent transitions. Divergences mean the sampler could not explore part of the posterior, so the reported means and intervals may be biased; refitting with a larger adapt_delta is the usual remedy.

Notes

This estimator supports (simple) hierarchical models. When multiple observations belong to at least one common sampling unit, the groups argument can specify the nesting structure (i.e., which rows in y, v, and X belong to each group).

The raw CmdStanPy fit is kept on self.result_, so its diagnostics remain reachable – est.result_.diagnose() reports R-hat, effective sample size, E-BFMI and treedepth alongside divergences.

Changed in version 0.0.5: groups accepts arbitrary hashable labels, and passing a numpy array no longer raises.

fit_dataset(dataset, *args, **kwargs)[source]

Apply the current estimator to the passed Dataset container.

A convenience interface that wraps fit() and automatically aligns the variables held in a Dataset with the required arguments.

Parameters:

  • dataset (Dataset) – A PyMARE Dataset instance holding the data.

  • *args – Optional positional arguments to pass onto the fit() method.

  • **kwargs – Optional keyword arguments to pass onto the fit() method.

get_v(dataset)[source]

Get the variances, or an estimate thereof, from the given Dataset.

Parameters:

dataset (Dataset) – The dataset to use to retrieve/estimate v.

Returns:

2-dimensional array of variances/variance estimates.

Return type:

numpy.ndarray

Notes

This is equivalent to directly accessing dataset.v when variances are present, but affords a way of estimating v from sample size (n) for any estimator that implicitly estimates a sigma^2 parameter.

summary(ci=95)[source]

Generate a BayesianMetaRegressionResults object from the fitted estimator.

Parameters:

ci (float, optional) – Desired width of the credible interval, as a percentage. Default = 95.0 (95%).

Return type:

BayesianMetaRegressionResults

Examples using pymare.estimators.StanMetaRegression

Run Estimators on a simulated dataset

Run Estimators on a simulated dataset