[query] lower logistic SKAT by danking · Pull Request #12643 · hail-is/hail

danking · 2023-02-02T02:10:45Z

CHANGELOG: hl.skat(..., logistic=True) now supported in the Batch backend.

CHANGELOG: `hl.skat(..., logistic=True)` now supported in the Batch backend.

danking · 2023-02-06T15:23:46Z

    """
-    if hl.current_backend().requires_lowering and not logistic:
-        ht = hl._linear_skat(key_expr, weight_expr, y, x, covariates, max_size, accuracy, iterations)
+    if hl.current_backend().requires_lowering:


In the second phase, we'll benchmark and address speed problems. Until we're certain it is fast enough, I don't want to enable it for all backends.

danking · 2023-02-06T15:25:31Z

+        [2, 1, 1, 1, 0, 1, 1, 2, 1, 1, 2, 1, 0, 0, 1],
+        [1, 0, 1, 1, 1, 2, 0, 2, 1, 1, 0, 1, 1, 0, 0],
+        [0, 2, 0, 0, 2, 1, 1, 2, 2, 1, 1, 1, 0, 1, 1],
+        [1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 0]]


I needed to create a matrix on which the logistic model converged.

danking · 2023-02-06T15:27:39Z

    assert result.fault == 0


-def test_linear_skat_R_truth():


This test is inadvertently duplicated in main. I just removed it entirely, but the diff is a bit awkward.a

danking · 2023-02-06T15:44:20Z

Hmm. I trust the code now. I test against several R SKAT runs. I'm not sure I understand how we derive that Q is generalized chi-squared distributed. We use the residual phenotypes in the calculation of Q, but those are inverse-logit transformed normal variables. The derivation for the linear case doesn't apply, as far as I can tell. I assume the residuals are Bernoulli distributed? Maybe not. I guess the phenotypes are Bernoulli but the errors aren't? I'm not sure.

danking · 2023-02-13T21:49:26Z

bump @patrick-schultz ideally we would merge this before Wednesday's meeting!

patrick-schultz · 2023-02-06T17:42:10Z

+    # See linear SKAT code comment for an extensive description of the mathematics here.
+
+    sqrtv = hl.sqrt(ht.s2)
+    Q, _ = hl.nd.qr((ht.covmat.T * sqrtv).T)


Should use ht.covmat.T * sqrtv.reshape(-1, 1)

covmat is K by N_samples. sqrtv is N_samples. The code I have above scales each sample by its s2, e.g.:

In [7]: import hail as hl ...: import numpy as np ...: s2 = np.array([1,2,3,4]) ...: covmat = np.array([[1,1],[1,2],[1,3],[1,1]]) ...: ...: hl.eval((hl.literal(covmat).T * hl.literal(s2)).T) Out[7]: array([[1, 1], [2, 4], [3, 9], [4, 4]])

If I use your expression,

In [11]: import hail as hl ...: import numpy as np ...: s2 = np.array([1,2,3,4]) ...: covmat = np.array([[1,1],[1,2],[1,3],[1,1]]) ...: ...: hl.eval((hl.literal(covmat).T * hl.literal(s2).reshape(-1, 1))) ... Error summary: HailException: Incompatible NDArray shapes: [ 2 4 ] vs [ 4 1 ]

But maybe you meant to say this?

In [11]: import hail as hl ...: import numpy as np ...: s2 = np.array([1,2,3,4]) ...: covmat = np.array([[1,1],[1,2],[1,3],[1,1]]) ...: ...: hl.eval((hl.literal(covmat) * hl.literal(s2).reshape(-1, 1))) Out[11]: array([[1, 1], [2, 4], [3, 9], [4, 4]])

I've made the last change since it seems better than all the transposition, but I must admit to not understanding how -1 functions in reshape.

Oops, yes, that's what I meant.

The -1 thing is a numpy behavior we mimic. From their docs:

One shape dimension can be -1. In this case, the value is inferred from the length of the array and remaining dimensions.

So reshape(-1, 1) is a universal "make this a column vector" command.

patrick-schultz · 2023-02-06T17:44:02Z

+    sqrtv = hl.sqrt(ht.s2)
+    Q, _ = hl.nd.qr((ht.covmat.T * sqrtv).T)
+    weights_arr = ht.weight._data_array()
+    G_scaled = (ht.G.T * sqrtv).T


patrick-schultz · 2023-02-06T17:49:58Z

+    1. The residual phenotypes are Bernoulli distributed with mean :math:`p` and variance
+       :math:`\sigma^2 = p(1 - p)` where :math:`p` is the best-fit probability.


Should be "the phenotypes are Bernoulli ..."

I changed this to have a sample index since the sigmas are sample indexed in the logistic case.

patrick-schultz · 2023-02-13T22:06:08Z

+    We can transform the residuals into standard normal variables by normalizing by their
+    variance.


I think everything below this doesn't make much sense. The residuals are certainly not normal -- they can only take two possible values. I tried to look into how to justify the logistic case, and made a little progress ("Pearson residuals" are relevant), but I think it's best to not try to write a complete derivation here. It should be enough to summarize the math from the paper that we implement.

patrick-schultz · 2023-02-13T22:07:31Z

bump @patrick-schultz ideally we would merge this before Wednesday's meeting!

My bad, didn't realize I hadn't submitted my review.

danking · 2023-02-15T23:33:28Z

@patrick-schultz , I've reworked the docs in a less wrong way. Let me know what you think!

x

patrick-schultz

Looks great, just a couple small fixes

patrick-schultz · 2023-02-16T17:10:00Z

+        X &: R^{N \times K} \\
+        G &: \{0, 1, 2\}^{N \times M} \\
+        \\
+        \textrm{logit}(P(y=1)) &= \beta_0 X + \beta_1 G + \varepsilon \quad\quad \varepsilon \sim N(0, \sigma^2)


This is still saying the residuals are normally distributed. I think the right way to express the logistic regression model is

y ~ Bernoulli(logit^{-1}(beta_0 X + beta_1 G))

patrick-schultz · 2023-02-16T17:11:33Z

+
+    .. math::
+
+        \textrm{logit}(P(y=1)) = \beta_\textrm{null} X + \varepsilon \quad\quad \varepsilon \sim N(0, \sigma^2)


patrick-schultz · 2023-02-16T17:18:27Z

+    Recall that the eigenvalues of a symmetric matrix and its transpose are the same, so we can
+    instead consider :math:`Z^T Z` (note that we elide transpositions of symmetric matrices):


This isn't saying what you mean (Z Z^T is symmetric, so its transpose is Z Z^T). The only way I know to say this is that the eigenvalues of both Z Z^T and Z^T Z are the squared singular values of Z.

patrick-schultz · 2023-02-16T17:36:48Z

+    # See linear SKAT code comment for an extensive description of the mathematics here.
+
+    sqrtv = hl.sqrt(ht.s2)
+    Q, _ = hl.nd.qr((ht.covmat.T * sqrtv).T)


Oops, yes, that's what I meant.

The -1 thing is a numpy behavior we mimic. From their docs:

One shape dimension can be -1. In this case, the value is inferred from the length of the array and remaining dimensions.

So reshape(-1, 1) is a universal "make this a column vector" command.

danking · 2023-02-16T23:31:44Z

Thanks for clarifying these points and pointing out the errors!

x

patrick-schultz

🎉

danking added the stacked PR label Feb 2, 2023

danking force-pushed the lowered-logistic-skat branch 4 times, most recently from 121d147 to bce7b7d Compare February 6, 2023 14:13

[query] lower logistic SKAT

6165ffb

CHANGELOG: `hl.skat(..., logistic=True)` now supported in the Batch backend.

danking removed the stacked PR label Feb 6, 2023

danking commented Feb 6, 2023

View reviewed changes

danking force-pushed the lowered-logistic-skat branch from bce7b7d to 6165ffb Compare February 6, 2023 15:26

danking commented Feb 6, 2023

View reviewed changes

danking assigned patrick-schultz Feb 6, 2023

Dan King added 2 commits February 6, 2023 10:33

minor doc fix

bc6b233

remove wrong sentence from docs

c187a78

Dan King added 3 commits February 9, 2023 15:35

covariates are required now

72ad2f0

support max_iter = 0 in logreg_fit

3b9364b

fix syntax

732de03

patrick-schultz previously requested changes Feb 13, 2023

View reviewed changes

Dan King added 3 commits February 15, 2023 14:38

wip docs

735ed6a

docs iwp

1aa181d

final take on the math

042c5e9

patrick-schultz previously requested changes Feb 16, 2023

View reviewed changes

docs fixes

f3bff04

patrick-schultz approved these changes Feb 17, 2023

View reviewed changes

danking merged commit 9c5851c into hail-is:main Feb 17, 2023

		1. The residual phenotypes are Bernoulli distributed with mean :math:`p` and variance
		:math:`\sigma^2 = p(1 - p)` where :math:`p` is the best-fit probability.

		We can transform the residuals into standard normal variables by normalizing by their
		variance.


		.. math::

		\textrm{logit}(P(y=1)) = \beta_\textrm{null} X + \varepsilon \quad\quad \varepsilon \sim N(0, \sigma^2)

		Recall that the eigenvalues of a symmetric matrix and its transpose are the same, so we can
		instead consider :math:`Z^T Z` (note that we elide transpositions of symmetric matrices):

Uh oh!

Conversation

danking commented Feb 2, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

danking commented Feb 6, 2023

Uh oh!

danking commented Feb 13, 2023

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

patrick-schultz commented Feb 13, 2023

Uh oh!

danking commented Feb 15, 2023

Uh oh!

patrick-schultz left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

danking commented Feb 16, 2023

Uh oh!

patrick-schultz left a comment

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

danking commented Feb 2, 2023 •

edited

Loading