[query] Lower linear SKAT by danking · Pull Request #12637 · hail-is/hail

danking · 2023-01-30T16:24:25Z

CHANGELOG: Query-on-Batch now supports hl.skat(..., logistic=False).

I also added actual tests for hl.skat, which were lost at some point.

I am somewhat not confident in my documentation and comments, because the SKAT paper is terse and unclear. I would really apprecaiate strong criticism of the documentation and the code comments.

CHANGELOG: Query-on-Batch now supports `hl.skat(..., logistic=False)`. I also added actual tests for `hl.skat`, which were lost at some point. I am somewhat not confident in my documentation and comments, because the SKAT paper is terse and unclear. I would really apprecaiate strong criticism of the documentation and the code comments.

danking · 2023-01-30T16:32:59Z

These are the rendered docs. They are not included in the actual docs because I do not want people using this function yet.

patrick-schultz

Some notes after a first pass

patrick-schultz · 2023-01-30T17:22:26Z

hail/python/hail/methods/statgen.py

+    4. Multiplying an orthogonal matrix by a vector of independent normal variables produces a new
+       vector of independent normal variables.


This is only true if you replace "independent" with "i.i.d."

patrick-schultz · 2023-01-30T17:24:01Z

hail/python/hail/methods/statgen.py

+    .. math::
+
+        \begin{align*}
+        U \Lambda U^T &= G W G^T \quad\quad U \textrm{ orthonormal, } \Lambda \textrm{ diagonal} \\


U orthogonal

patrick-schultz · 2023-01-30T18:11:24Z

hail/python/hail/methods/statgen.py

+    ht = ht.annotate(
+        Q = (((ht.y_residual @ ht.G) * ht.weight) @ ht.G.T) @ ht.y_residual.T
+    )


More efficient:

((ht.y_residual @ ht.G).map(lambda x: x**2) * ht.weight).sum(0)

patrick-schultz · 2023-01-30T18:13:53Z

hail/python/hail/methods/statgen.py

+    #
+    # We avoid the square root in order to avoid complex numbers.
+
+    Q, _ = hl.nd.qr(ht.covmat)


Already computed this

patrick-schultz · 2023-01-30T18:20:35Z

hail/python/hail/methods/statgen.py

+
+    Q, _ = hl.nd.qr(ht.covmat)
+    C0 = Q.T @ ht.G
+    singular_values = hl.nd.svd(ht.G.T @ (ht.G * ht.weight) - C0.T @ (C0 * ht.weight), compute_uv=False)


As discussed in zulip, I think a better way to compute this is:

R = nd.qr(ht.G - Q @ (Q.T @ G), mode="r") singular_values = hl.nd.svd(R, compute_uv=False) eigenvalues = singular_values.map(lambda x: x**2)

(replace singular_values with eigenvalues below)

danking · 2023-01-31T23:57:16Z

Comments addressed. Still gotta work on the verbiage to clearly explain the present of P_0.

x

danking · 2023-02-01T22:56:04Z

@patrick-schultz Alright. This passes tests locally. I'm happy with the docs verbiage. I am eager for your review!

patrick-schultz

A few docs comments. Everything else looks great!

patrick-schultz · 2023-02-02T15:29:06Z

hail/python/hail/methods/statgen.py

+    .. math::
+
+        \begin{align*}
+        X &: R^{N \times K} \\


Might mention that X is covariates. I was momentarily confused since in the args x is the genotypes.

patrick-schultz · 2023-02-02T15:30:36Z

hail/python/hail/methods/statgen.py

+        h &\sim N(0, 1) \\
+        h &= \frac{1}{\widehat{\sigma}} r \\


I don't know what the standard stats idiom is, but these feel backwards to me. I want to read this as "if we define h = ..., then h is distributed as N(0, 1)". The other way sounds like "draw a random normal variate h, then h = this other thing we already computed".

patrick-schultz · 2023-02-02T15:34:53Z

hail/python/hail/methods/statgen.py

+    .. math::
+
+        \begin{align*}
+        U \Lambda U &= B \quad\quad \Lambda \textrm{ orthogonal } U \textrm{ diagonal} \\


Should be U \Lambda U^T, Lambda diagonal, U orthogonal

patrick-schultz · 2023-02-02T15:56:45Z

hail/python/hail/methods/statgen.py

+    #     B = A A.T
+    #     Q = h.T B h
+    #
+    # This is called a "quadratic form". It is a weighted sum of the squares of the entries of h,


Not quite. It's a weighted sum of products of pairs of entries of h (w_{00} h_0^2 + w_{01} h_0 h_1 + ....). The eigendecomposition converts it to a sum of squares of normals.

patrick-schultz · 2023-02-02T16:01:30Z

hail/python/hail/methods/statgen.py

+    # Since B is a real symmetric matrix, U is orthogonal. U and W are not necessarily the same
+    # matrix but their determinants are +-1 so the squared singular values and eigenvalues differ by
+    # at most a sign.


U and W are the same, in the sense that eigenvectors / singular vectors are only defined up to sign anyways. But we don't care about the eigenvectors, do we?

The determinant bit is confusing (their determinants are +-1 because they're orthogonal, but that doesn't tell us anything), and the squared singular values are equal to the eigenvalues; no sign ambiguity there. B is pos. def., it has positive eigenvalues.

I think you can just drop this paragraph.

patrick-schultz · 2023-02-02T16:09:57Z

hail/python/hail/methods/statgen.py

+    ).or_error(hl.format('hl._linear_skat: every weight must be positive, in group %s, the weights were: %s',
+                         ht.group, weights_arr))
+    singular_values = hl.nd.svd(A, compute_uv=False)
+    # SVD(M) = U S V. U and V are unitary, therefore SVD(k M) = U (k S) V.


Not sure what this is about

In the next line, we scale the singular values by \sigma^2 instead of multiplying A by sqrt(\sigma^2).

I added a blank line to clarify to which expression the comment applies.

Oh I see. That doesn't have anything to do with unitarity. Scalars commute with all matrices.

patrick-schultz · 2023-02-02T16:12:20Z

hail/python/hail/methods/statgen.py

+        # I *think* the reasoning for taking the complement of the CDF value is:
+        #
+        # 1. Q is a measure of variance and thus positive.
+        #
+        # 2. We want to know the probability of obtaining a variance even larger ("more extreme")
+        #
+        # Ergo, we want to check the right-tail of the distribution.


I addressed every comment except for the one about the SVD(k M) comment; let me know if that makes sense now.

patrick-schultz · 2023-02-02T21:58:17Z

hail/python/hail/methods/statgen.py

+    .. math::
+
+        \begin{align*}
+        U \Lambda U &= B \quad\quad \Lambda \textrm{ diagonal } U \textrm{ orthogonal} \\


Still missing a transpose: U \Lambda U.T

fixed that one and a couple other ones where the transpose was missing

danking assigned patrick-schultz Jan 30, 2023

danking force-pushed the lowered-linear-skat branch from 2892148 to 6f5e795 Compare January 30, 2023 16:25

patrick-schultz previously requested changes Jan 30, 2023

View reviewed changes

fix

3b26208

fix test

25d1c51

danking force-pushed the lowered-linear-skat branch from 03a8191 to a45a1b7 Compare February 1, 2023 19:40

fixes to tests

8139db1

danking force-pushed the lowered-linear-skat branch from a45a1b7 to ce65d52 Compare February 1, 2023 19:50

Dan King added 6 commits February 1, 2023 14:52

fix import

2fb05d2

address comments

ee5759b

fix

08344ce

fix scoping of covmat_Q

362dbb2

REVERT THIS

170a2c6

fix

5bd1f52

danking force-pushed the lowered-linear-skat branch from ce65d52 to 5bd1f52 Compare February 1, 2023 19:52

Dan King added 8 commits February 1, 2023 15:28

fix calculation of singular values to account for s2 and weight

d83c3bc

passes

d960600

fixes

f7a7e8b

final fixes to tests

c9ad2ba

new test, rework docs, check for negative weights

d402000

reverting docs changes

b9dc6a6

add missing test file

0a24277

lots of docs fixes, notational changes etc.

843a9c1

simplify Q statistic calculation

d4f53c4

danking mentioned this pull request Feb 2, 2023

[query] lower logistic SKAT #12643

Merged

maybe fix docs

3316640

patrick-schultz previously requested changes Feb 2, 2023

View reviewed changes

Dan King added 2 commits February 2, 2023 14:31

fix p-value in docs test

918e1ac

address comments

e733239

patrick-schultz previously requested changes Feb 2, 2023

View reviewed changes

fix U transposes

02d0ef5

patrick-schultz approved these changes Feb 3, 2023

View reviewed changes

danking merged commit 743db89 into hail-is:main Feb 3, 2023

		4. Multiplying an orthogonal matrix by a vector of independent normal variables produces a new
		vector of independent normal variables.

Conversation

danking commented Jan 30, 2023

Uh oh!

danking commented Jan 30, 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!

Choose a reason for hiding this comment

Uh oh!

danking commented Jan 31, 2023

Uh oh!

danking commented Feb 1, 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!

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!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants