bss.models package¶
Submodules¶
bss.models.probit module¶
-
class
bss.models.probit.Probit(X, Y, R, target_sparsity=0.01, gamma0_v=1.0, lambda_params=(1e-06, 1e-06), nu_params=(1e-06, 1e-06), xi=0.999999, xi_prior_shape=(1, 1), check_finite=True, min_eigenval=0, jitter=1e-06)¶ The Probit model used for modeling Sparse Regression using a Gaussian field. [EA14].
\[y|X,\beta,\beta_0, \nu \propto \mathcal{N}(\beta_0 1_n + X \beta, \nu^{-1} I_n)\]Parameters: X : ndarray
The predictor matrix of real numbers, n x p in size, where n is the no. of samples (genotypes) and p is the no. of features (SNPs).
Y : ndarray
The response vector of real numbers, n x 1 in size, with each value representing the phenotype value for the sample.
R : ndarray
The covariance matrix for the SNPs, p x p in size. The matrix may not be positive-definite, but is converted to one internally.
target_sparsity : float
The proportion of included predictors. For example, a value of 0.01 indicates that around 1% of total SNPs are expected be included in our model. This value affects the probit threshold gamma_0 of the model.
gamma0_v : float
Variance of the probit threshold gamma_0
lambda_params : tuple
Shape parameter and Inverse-scale parameter of the gamma prior placed on the model parameter lambda, where lambda is the inverse squared global scale parameter for the regression weights.
nu_params : tuple
Shape parameter and Inverse-scale parameter of the gamma prior placed on the model parameter nu, where nu is the residual precision.
xi : float
The shrinkage constant in the interval [0,1] to regularize the covariance matrix towards the identity matrix. This ensures that the covariance matrix is positive definite. A larger xi value biases our estimate towards the supplied R matrix, a lower value biases it towards the identity matrix. If None, then xi is sampled from a beta distribution with shape parameters specified by the tuple xi_prior_shape.
xi_prior_shape : tuple
Shape parameters of the beta prior placed on the model parameter xi, specified as a 2-tuple of real values. This argument is ignored and xi is not sampled, if it is specified explicitly using the xi parameter.
check_finite : bool
Whether to check that the input matrices contain only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs. This parameter is passed on to several linear algebra functions in scipy internally.
min_eigenval : float
Minimum Eigenvalue we can accept in the covariance matrix. Any eigenvalues encountered below this threshold are set to zero, and the resulting covariance matrix normalized to give ones on the diagonal.
jitter : float
A small value to add to the diagonals of the covariance matrix to avoid conditioning issues.
Methods
-
log_joint()¶ The joint log-likelihood value of given model parameters.
Returns: float
A scalar joint log-likelihood value of the entire model, given the model parameters and the true Y values stored in this model object.
-
log_marg_like(gamma, gamma0, lamb, nu)¶ The marginal likelihood log-value of given model parameters.
Parameters: gamma : ndarray
A dx1 vector of probit values, one per SNP
gamma0 : float
The probit threshold, All SNPs with gamma > gamma0 are ‘activated’.
lamb : float
Inverse squared global scale parameter for the regression weights.
nu : float
The residual precision of the model.
Returns: float
A scalar log-likelihood value of the model, given the model parameters (and the true Y values stored in this model object).
-
ppi_distribution(*args, **kwargs)¶ Normal distribution for the posterior probability of inclusion.
Parameters: gamma : ndarray
A dx1 vector of probit values, one per SNP
gamma0 : float
The probit threshold, All SNPs with gamma > gamma0 are ‘activated’.
lamb : float
Inverse squared global scale parameter for the regression weights.
Returns: Mvn
A multivariate normal distribution object representing the posterior probability of inclusion
Notes
- We’re interested in the posterior probability of inclusion:
- \[p(y|X,\gamma,\gamma_0,\nu,\lambda)\]
- marginalizing out the effect size captured by \(\beta\):
- \[\beta | \nu,\lambda,\Gamma \propto \mathcal{N}(0, (\nu\lambda)^{-1}\Gamma)\]
The degenerate Gaussian form of the \(\beta\) prior above (note that the covariance matrix \(\Gamma\) is a diagonal matrix of indicator values) allows us to perform this marginalization in closed form:
\[ \begin{align}\begin{aligned}p(y|X,\gamma,\gamma_0,\nu,\lambda)\\= \int \int \mathcal{N} (y|\beta_01_n + X\beta,\nu^{-1}I_n \; \mathcal{N} (\beta|0, (\nu\lambda)^{-1}\Gamma)) \; \mathcal{N}(\beta_0|0,(\nu\lambda)^{-1}) \; d\beta d\beta_0\\= \int \mathcal{N} (y|\beta_0 1_n, \nu^{-1}(\lambda^{-1}X\Gamma X^T + I_n)) \; \mathcal{N}(\beta_0 | 0, (\nu\lambda)^{-1}) d\beta_0\\= \mathcal{N} (y|0, \nu^{-1} \lambda^{-1}(1_n1_n^T + X\Gamma X^T) + I_n))\end{aligned}\end{align} \]
-
probit_distribution(*args, **kwargs)¶ The probit distribution of the model, driven by a latent Gaussian field.
Parameters: xi : float
The shrinkage constant in the interval [0,1] to regularize the covariance matrix towards the identity matrix. This ensures that the covariance matrix is positive definite. A larger xi value biases our estimate towards the supplied R matrix, a lower value biases it towards the identity matrix.
Returns: The normal distribution of the probits (gamma) of the model.
-
run_mcmc(iters=1000, burn_in=100, detailed=False)¶ Execute a run of MCMC on this model given the current state of model parameters.
Parameters: iters : int
The total no. of iterations (excluding the burn-in) to run the MCMC simulations.
burn_in : int
The total no. of burn-in iterations of the MCMC simulation.
detailed : bool, optional
Whether to return a detailed trace of individual model parameters, or a trace of the joint log-likelihood value.
Returns: list, dictionary
A list of scalars (joint log-likelihood value of the entire model), or a dictionary of lists (with keys inclusion, gamma0, lambda, nu, xi, joint, likelihood), and traces of individual model parameters as values in the dictionary).
-
update_gamma()¶ Apply MCMC transition operator to model parameter \(\gamma\). For updating \(\gamma\), we use elliptical slice sampling (ESS) described in [MAM10]
Returns: On return, the \(\gamma\) parameter of the model has been updated using a new sample. Notes
ESS samples efficiently and robustly from latent Gaussian models when significant covariance structure is imposed by the prior, as in the Gaussian processes and the present structured sparsity model.
ESS generates random elliptical loci using hte Gaussian prior and then searches along these loci to find acceptable points for slice sampling. When the data ar weakly informative and the prior is strong, as is the case here, the elliptical loci effectively captures the dependence between the variables and enable faster mixing. Here, using ESS for \(\gamma\) enables us to avoid directly sampling over the large discrete space of sparsity patterns that makes unstructured spike-and-slab computationally challenging.
-
update_gamma0()¶ Apply MCMC transition operator to the model parameter \(\gamma_0\), the sparsity threshold.
Returns: On return, the \(\gamma_0\) parameter of the model has been updated using a new sample. Notes
The parameter \(\gamma_0\) specifies the probit threshold and, conditioned on \(\gamma\), it determines which entries on the diagonal of \(\Gamma\) are zero and which are one. The conditional density of \(\gamma_0\) does not have a simple closed form, but can be efficiently sampled using the exponential-expansion slice sampling algorithm described in [Nea03]
-
update_lambda()¶ Apply MCMC transition operator to model parameter \(\lambda\), the global weight inverse-scale parameter.
Returns: On return, the \(\lambda\) parameter of the model has been updated using a new sample. Notes
The parameter \(\lambda\) determines the scale of the “slab” portion of the weight prior. The conditional density of \(\lambda\) does not have a simple closed form, but can be efficiently sampled using the exponential-expansion slice sampling algorithm described in [Nea03]
-
update_nu()¶ Apply MCMC transition operator to model parameter \(\nu\), the residual precision.
Returns: On return, the \(\nu\) parameter of the model has been updated using this analytical approach. Notes
The scalar nu determines the precision of the residual Gaussian noise of the response variables. With the choice of a conjugate gamma prior distribution, the conditional posterior is also gamma:
\[ \begin{align}\begin{aligned}p(\nu | y, X, \Gamma, \lambda) \propto \mathcal{N}(y|0,\nu^{-1}(\lambda^{-1}(1_n1_n^T + X \Gamma X^T) + I_n)) \, Gam(\nu|a_\nu, b_\nu)\\= Gam(\nu | a_\nu^{(n)}, b_\nu^{(n)})\\a_\nu^{(n)} = a_\nu + \frac{N}{2}\\b_\nu^{(n)} = b_\nu + \frac{1}{2} y^T (\lambda^{-1} (1_n 1_n^T + X \Gamma X^T) + I_n)^{-1} y\end{aligned}\end{align} \]
-
update_parameters()¶ Update all model parameters in a single MCMC step.
Returns: On return, the gamma, gamma0, lambda, nu and optionally the xi model parameter have been updated.
-
update_xi()¶ Apply MCMC transition operator to the shrinkage factor \(\xi\), used to regularize the covariance matrix towards the identity matrix.
Returns: On return, the \(\xi\) parameter of the model has been updated using a new sample.
-