Skip to content

VIPRSBMA

VIPRSBMA

Bases: VIPRSGrid

The VIPRSBMA class is an extension of the VIPRSGrid class that implements Bayesian model averaging for the VIPRS models in the grid. Bayesian model averaging is a technique that allows us to combine the results of multiple models by weighting them according to their evidence. In this context, we weigh the model by their final ELBO values.

For more details on the BMA procedure implemented here, refer to the Supplementary material of:

Zabad S, Gravel S, Li Y. Fast and accurate Bayesian polygenic risk modeling with variational inference. Am J Hum Genet. 2023 May 4;110(5):741-761. doi: 10.1016/j.ajhg.2023.03.009. Epub 2023 Apr 7. PMID: 37030289; PMCID: PMC10183379.

Source code in viprs/model/gridsearch/VIPRSBMA.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
class VIPRSBMA(VIPRSGrid):
    """
    The `VIPRSBMA` class is an extension of the `VIPRSGrid` class that
    implements Bayesian model averaging for the `VIPRS` models in the grid.
    Bayesian model averaging is a technique that allows us to combine the
    results of multiple models by weighting them according to their evidence.
    In this context, we weigh the model by their final ELBO values.

    For more details on the BMA procedure implemented here, refer to the
    Supplementary material of:

    > Zabad S, Gravel S, Li Y. Fast and accurate Bayesian polygenic risk modeling with variational inference.
    Am J Hum Genet. 2023 May 4;110(5):741-761. doi: 10.1016/j.ajhg.2023.03.009.
    Epub 2023 Apr 7. PMID: 37030289; PMCID: PMC10183379.

    """

    def __init__(self,
                 gdl,
                 grid,
                 **kwargs):
        """

        Initialize the `VIPRSBMA` model.

        :param gdl: An instance of `GWADataLoader`
        :param grid: An instance of `HyperparameterGrid`
        :param kwargs: Additional keyword arguments for the VIPRS model
        """

        super().__init__(gdl, grid=grid, **kwargs)

    def average_models(self, normalization='softmax'):
        """
        Use Bayesian model averaging (BMA) to obtain final weights for each parameter.
        We average the weights by using the final ELBO for each model.

        :param normalization: The normalization scheme for the final ELBOs.
        Options are (`softmax`, `sum`).
        :raises KeyError: If the normalization scheme is not recognized.
        """

        if self.n_models < 2:
            return self

        # Extract the models that converged successfully:
        models_to_keep = np.where(self.models_to_keep)[0]

        elbos = self.history['ELBO'][-1][models_to_keep]

        if normalization == 'softmax':
            from scipy.special import softmax
            weights = np.array(softmax(elbos))
        elif normalization == 'sum':
            weights = np.array(elbos)

            # Correction for negative ELBOs:
            weights = weights - weights.min() + 1.
            weights /= weights.sum()
        else:
            raise KeyError("Normalization scheme not recognized. "
                           "Valid options are: `softmax`, `sum`. "
                           "Got: {}".format(normalization))

        if int(self.verbose) > 1:
            logging.info("Averaging PRS models with weights:", weights)

        # Average the model parameters:
        for param in (self.pip, self.post_mean_beta, self.post_var_beta,
                      self.var_gamma, self.var_mu, self.var_tau,
                      self.eta, self.zeta, self.q):

            for c in param:
                param[c] = (param[c][:, models_to_keep]*weights).sum(axis=1)

        # Set the number of models to 1:
        self.n_models = 1

        return self

__init__(gdl, grid, **kwargs)

Initialize the VIPRSBMA model.

Parameters:

Name Type Description Default
gdl

An instance of GWADataLoader

 required
grid

An instance of HyperparameterGrid

 required
kwargs

Additional keyword arguments for the VIPRS model

 {}
Source code in viprs/model/gridsearch/VIPRSBMA.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
def __init__(self,
             gdl,
             grid,
             **kwargs):
    """

    Initialize the `VIPRSBMA` model.

    :param gdl: An instance of `GWADataLoader`
    :param grid: An instance of `HyperparameterGrid`
    :param kwargs: Additional keyword arguments for the VIPRS model
    """

    super().__init__(gdl, grid=grid, **kwargs)

average_models(normalization='softmax')

Use Bayesian model averaging (BMA) to obtain final weights for each parameter. We average the weights by using the final ELBO for each model.

Parameters:

Name Type Description Default
normalization

The normalization scheme for the final ELBOs. Options are (softmax, sum).

 'softmax'

Raises:

Type Description
KeyError

If the normalization scheme is not recognized.
Source code in viprs/model/gridsearch/VIPRSBMA.py 
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
def average_models(self, normalization='softmax'):
    """
    Use Bayesian model averaging (BMA) to obtain final weights for each parameter.
    We average the weights by using the final ELBO for each model.

    :param normalization: The normalization scheme for the final ELBOs.
    Options are (`softmax`, `sum`).
    :raises KeyError: If the normalization scheme is not recognized.
    """

    if self.n_models < 2:
        return self

    # Extract the models that converged successfully:
    models_to_keep = np.where(self.models_to_keep)[0]

    elbos = self.history['ELBO'][-1][models_to_keep]

    if normalization == 'softmax':
        from scipy.special import softmax
        weights = np.array(softmax(elbos))
    elif normalization == 'sum':
        weights = np.array(elbos)

        # Correction for negative ELBOs:
        weights = weights - weights.min() + 1.
        weights /= weights.sum()
    else:
        raise KeyError("Normalization scheme not recognized. "
                       "Valid options are: `softmax`, `sum`. "
                       "Got: {}".format(normalization))

    if int(self.verbose) > 1:
        logging.info("Averaging PRS models with weights:", weights)

    # Average the model parameters:
    for param in (self.pip, self.post_mean_beta, self.post_var_beta,
                  self.var_gamma, self.var_mu, self.var_tau,
                  self.eta, self.zeta, self.q):

        for c in param:
            param[c] = (param[c][:, models_to_keep]*weights).sum(axis=1)

    # Set the number of models to 1:
    self.n_models = 1

    return self