Skip to content

graphld.likelihood

Gaussian likelihood functions for GWAS summary statistics under an infinitesimal model.

Model

The likelihood of GWAS summary statistics under an infinitesimal model is:

$$\beta \sim N(0, D)$$

$$z|\beta \sim N(n^{1/2}R\beta, R)$$

where:

  • $\beta$ is the effect-size vector in s.d-per-s.d. units
  • $D$ is a diagonal matrix of per-variant heritabilities
  • $z$ is the GWAS summary statistic vector
  • $R$ is the LD correlation matrix
  • $n$ is the sample size

Precision-Premultiplied Statistics

The likelihood functions operate on precision-premultiplied GWAS summary statistics:

$$pz = n^{-1/2} R^{-1}z \sim N(0, M), \quad M = D + n^{-1}R^{-1}$$

For model context and practical interpretation, see the Likelihood Functions guide.

likelihood

Functions for computing likelihoods in the LDGM model.

gaussian_likelihood 

gaussian_likelihood(pz: ndarray, M: PrecisionOperator) -> float

Compute log-likelihood of GWAS summary statistics under a Gaussian model.

The model is

beta ~ MVN(0, D) z|beta ~ MVN(sqrt(n)Rbeta, R) where R is the LD matrix, n the sample size pz = inv(R) * z / sqrt(n) M = cov(pz) = D + inv(R)/n

Parameters:

Name Type Description Default
pz ndarray

Array of precision-premultiplied GWAS effect size estimates

 required
M PrecisionOperator

PrecisionOperator. This should be the covariance of pz.

 required

Returns:

Type Description
float

Log-likelihood value
Source code in src/graphld/likelihood.py 
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
def gaussian_likelihood(
    pz: np.ndarray,
    M: PrecisionOperator,
) -> float:
    """Compute log-likelihood of GWAS summary statistics under a Gaussian model.

    The model is:
        beta ~ MVN(0, D)
        z|beta ~ MVN(sqrt(n)*R*beta, R) where R is the LD matrix, n the sample size
        pz = inv(R) * z / sqrt(n)
        M = cov(pz) = D + inv(R)/n

    Args:
        pz: Array of precision-premultiplied GWAS effect size estimates
        M: PrecisionOperator. This should be the covariance of pz.

    Returns:
        Log-likelihood value

    """
    # Following scipy's convention:
    # log_pdf = -0.5 * (n * log(2π) + log|Σ| + x^T Σ^{-1} x)
    n = len(pz)
    logdet = M.logdet()

    # Compute quadratic form
    b = M.solve(pz)
    quad = np.sum(pz * b)

    # Compute log likelihood
    ll = -0.5 * (n * np.log(2 * np.pi) + logdet + quad)

    return ll

gaussian_likelihood_gradient 

gaussian_likelihood_gradient(pz: ndarray, M: PrecisionOperator, del_M_del_a: Optional[ndarray] = None, n_samples: int = 10, seed: Optional[int] = None, trace_estimator: Optional[str] = 'xdiag') -> np.ndarray

Computes the score under a Gaussian model.

The model is:
beta ~ MVN(0, D)
z|beta ~ MVN(sqrt(n)*R*beta, R) where R is the LD matrix, n the sample size
pz = inv(R) * z / sqrt(n)
M = cov(pz) = D + inv(R)/n

Parameters:

Name Type Description Default
pz ndarray

Array of precision-premultiplied GWAS effect size estimates

 required
M PrecisionOperator

PrecisionOperator. This should be the covariance of pz.

 required
del_M_del_a Optional[ndarray]

Matrix of derivatives of M's diagonal elements wrt parameters a

 None
n_samples int

Number of probe vectors for Hutchinson's method or xdiag

 10
seed Optional[int]

Random seed for generating probe vectors

 None
trace_estimator Optional[str]

Method for computing the trace estimator. Options: "exact", "hutchinson", "xdiag"

 'xdiag'

Returns:

Type Description
ndarray

Array of diagonal elements of the gradient wrt M's diagonal elements,
ndarray

or with respect to parameters a if del_M_del_a is provided
Source code in src/graphld/likelihood.py 
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
86
87
88
89
90
91
92
93
94
95
96
97
def gaussian_likelihood_gradient(
    pz: np.ndarray,
    M: PrecisionOperator,
    del_M_del_a: Optional[np.ndarray] = None,
    n_samples: int = 10,
    seed: Optional[int] = None,
    trace_estimator: Optional[str] = "xdiag",
) -> np.ndarray:
    """Computes the score under a Gaussian model.

        The model is:
        beta ~ MVN(0, D)
        z|beta ~ MVN(sqrt(n)*R*beta, R) where R is the LD matrix, n the sample size
        pz = inv(R) * z / sqrt(n)
        M = cov(pz) = D + inv(R)/n

    Args:
        pz: Array of precision-premultiplied GWAS effect size estimates
        M: PrecisionOperator. This should be the covariance of pz.
        del_M_del_a: Matrix of derivatives of M's diagonal elements wrt parameters a
        n_samples: Number of probe vectors for Hutchinson's method or xdiag
        seed: Random seed for generating probe vectors
        trace_estimator: Method for computing the trace estimator.
            Options: "exact", "hutchinson", "xdiag"

    Returns:
        Array of diagonal elements of the gradient wrt M's diagonal elements,
        or with respect to parameters a if del_M_del_a is provided
    """
    # Compute b = M^(-1) * pz
    b = M.solve(pz)

    # Compute diagonal elements of M^(-1)
    minv_diag = M.inverse_diagonal(method=trace_estimator,
                                    n_samples=n_samples,
                                    seed=seed)

    # Compute gradient diagonal elements
    node_grad = -0.5 * (minv_diag.flatten() - b.flatten()**2)

    if del_M_del_a is None:
        return node_grad
    with np.errstate(divide="ignore", invalid="ignore", over="ignore"):
        return node_grad @ del_M_del_a

gaussian_likelihood_hessian 

gaussian_likelihood_hessian(pz: ndarray, M: PrecisionOperator, del_M_del_a: Optional[ndarray] = None, trace_estimator: str = _TRACE_ESTIMATOR_DEFAULT, n_samples: int = 100, seed: Optional[int] = None, *, diagonal_method: Optional[str] = None) -> np.ndarray

Computes the average information matrix of the Gaussian log-likelihood.

The model is

beta ~ MVN(0, D) z|beta ~ MVN(sqrt(n)Rbeta, R) where R is the LD matrix, n the sample size pz = inv(R) * z / sqrt(n) M = cov(pz) = D + inv(R)/n

Parameters:

Name Type Description Default
pz ndarray

Array of precision-premultiplied GWAS effect size estimates

 required
M PrecisionOperator

PrecisionOperator. This should be the covariance of pz.

 required
del_M_del_a Optional[ndarray]

Matrix of derivatives of M's diagonal elements wrt parameters a. If None, only the diagonal elements are computed.

 None
trace_estimator str

Method for computing the diagonal trace estimator when del_M_del_a is None. Options accepted by PrecisionOperator include "exact", "hutchinson", and "xdiag". Defaults to "xdiag".

 _TRACE_ESTIMATOR_DEFAULT
n_samples int

Number of probe vectors for Hutchinson's method or xdiag

 100
seed Optional[int]

Random seed for generating probe vectors

 None
diagonal_method Optional[str]

Deprecated alias for trace_estimator.

 None

Returns:

Type Description
ndarray

Matrix of second derivatives wrt parameters a, or array of diagonal elements
ndarray

if del_M_del_a is None
Source code in src/graphld/likelihood.py 
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
def gaussian_likelihood_hessian(
    pz: np.ndarray,
    M: PrecisionOperator,
    del_M_del_a: Optional[np.ndarray] = None,
    trace_estimator: str = _TRACE_ESTIMATOR_DEFAULT,
    n_samples: int = 100,
    seed: Optional[int] = None,
    *,
    diagonal_method: Optional[str] = None,
) -> np.ndarray:
    """Computes the average information matrix of the Gaussian log-likelihood.

    The model is:
        beta ~ MVN(0, D)
        z|beta ~ MVN(sqrt(n)*R*beta, R) where R is the LD matrix, n the sample size
        pz = inv(R) * z / sqrt(n)
        M = cov(pz) = D + inv(R)/n

    Args:
        pz: Array of precision-premultiplied GWAS effect size estimates
        M: PrecisionOperator. This should be the covariance of pz.
        del_M_del_a: Matrix of derivatives of M's diagonal elements wrt parameters a.
            If None, only the diagonal elements are computed.
        trace_estimator: Method for computing the diagonal trace estimator when
            del_M_del_a is None. Options accepted by PrecisionOperator include
            "exact", "hutchinson", and "xdiag". Defaults to "xdiag".
        n_samples: Number of probe vectors for Hutchinson's method or xdiag
        seed: Random seed for generating probe vectors
        diagonal_method: Deprecated alias for trace_estimator.

    Returns:
        Matrix of second derivatives wrt parameters a, or array of diagonal elements
        if del_M_del_a is None
    """
    trace_estimator_was_supplied = trace_estimator is not _TRACE_ESTIMATOR_DEFAULT
    if diagonal_method is not None:
        if trace_estimator_was_supplied and trace_estimator != diagonal_method:
            raise ValueError(
                "trace_estimator and diagonal_method specify different Hessian "
                "diagonal estimators; use trace_estimator only."
            )
        warnings.warn(
            "diagonal_method is deprecated; use trace_estimator instead.",
            DeprecationWarning,
            stacklevel=2,
        )
        trace_estimator = diagonal_method
    elif not trace_estimator_was_supplied:
        trace_estimator = "xdiag"

    # Compute b = M^(-1) * pz
    b = M.solve(pz)
    if b.ndim == 1:
        b = b.reshape(-1, 1)

    # If del_M_del_a is None, compute only the diagonal of the Hessian
    if del_M_del_a is None:
        minv_diag = M.inverse_diagonal(
            method=trace_estimator, n_samples=n_samples, seed=seed
        )
        hess_diag = -0.5 * minv_diag.flatten() * b.flatten()**2
        return hess_diag

    # Compute b_scaled = b .* del_sigma_del_a
    b_scaled = b * del_M_del_a

    # Compute M^(-1) * b_scaled
    minv_b_scaled = M.solve(b_scaled)

    # Compute Hessian: -1/2 * b_scaled^T * M^(-1) * b_scaled
    with np.errstate(divide="ignore", invalid="ignore", over="ignore"):
        hess = -0.5 * (b_scaled.T @ minv_b_scaled)

    return hess