# Logistic regression burden tests

**URL:** <https://discuss.hail.is/t/logistic-regression-burden-tests/206>\
**Category:** Updates\
**Created:** [May 9, 2017, 1:27pm UTC](https://discuss.hail.is/t/logistic-regression-burden-tests/206 "2017-05-09T13:27:39Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![jbloom](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/jbloom/32/109_2.png) [@jbloom](https://discuss.hail.is/u/jbloom)\
**Post date:** [May 9, 2017, 1:27pm UTC](https://discuss.hail.is/t/logistic-regression-burden-tests/206/1 "2017-05-09T13:27:39Z")

</div>

We’ve added the method [logreg\_burden](https://hail.is/hail/hail.VariantDataset.html#hail.VariantDataset.logreg_burden) to provide a very general framework for aggregating (collapsing) genotypes and then applying logistic regression (Wald, LRT, Firth, or Score). This method returns both a key table of logistic regression results and a key table of collapsed sample scores for those samples with non-missing phenotype and covariates. The latter can then be further analyzed in PySpark or locally in tools like Python and R. Click the method link above for more details.

The performance is very similar to that of [linreg\_burden](https://hail.is/hail/hail.VariantDataset.html#hail.VariantDataset.linreg_burden) since there are typically far fewer keys for regression than variants for aggregation; it takes just a [few minutes](http://discuss.hail.is/t/linear-regression-burden-tests-collapsing-genotypes-by-variant-key/194) to run logistic gene burden starting from the annotated 2353 whole genomes in the 1000 Genomes Project.

---

<div class="post-metadata">

**Author:** ![enriquea](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/enriquea/32/237_2.png) [@enriquea](https://discuss.hail.is/u/enriquea)\
**Post date:** [September 4, 2018, 3:01pm UTC](https://discuss.hail.is/t/logistic-regression-burden-tests/206/2 "2018-09-04T15:01:52Z")

</div>

Hi @jbloom,

Could you please point me out to a similar implementation of logreg\_burden in Hail devel (0.2) version?

Thanks in advance!

---

<div class="post-metadata">

**Author:** ![jbloom](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/jbloom/32/109_2.png) [@jbloom](https://discuss.hail.is/u/jbloom)\
**Post date:** [September 7, 2018, 10:50pm UTC](https://discuss.hail.is/t/logistic-regression-burden-tests/206/3 "2018-09-07T22:50:07Z")

</div>

Hi Enrique,

Regressions in Hail 0.2 are generic, so there are no longer special command for burden tests. Rather, to do logistic gene burden, annotate and group variants by gene to form a gene-by-sample matrix of scores, and then apply logistic regression per gene. Here’s an example that uses the total number minor alleles as a sample’s gene score (for a biallelic dataset with random phenotype and covariates):

```auto
import hail as hl

interval_ht = hl.import_locus_intervals('genes.37.interval_list')

mt = hl.import_vcf('profile.vcf.bgz')
mt = mt.annotate_cols(y = hl.rand_bool(0.5), x1 = hl.rand_norm(), x2 = hl.rand_norm())
mt = hl.variant_qc(mt)
mt = mt.annotate_rows(gene = interval_ht[mt.locus].target)
mt = mt.filter_rows(hl.is_defined(mt.gene))

gene_mt = mt.group_rows_by(mt.gene).aggregate(
    mac = hl.agg.sum(
        hl.cond(mt.variant_qc.AF[1] <= 0.5,
                mt.GT.n_alt_alleles(),
                2 - mt.GT.n_alt_alleles())))

gene_mt = hl.logistic_regression(y=gene_mt.y, x=gene_mt.mac, covariates=[1, gene_mt.x1, gene_mt.x2], test='wald')
gene_mt.logreg.show()

```

Here’s what’s shown:

```auto
+----------+--------------+-----------------------+---------------+----------------+
| gene | logreg.beta | logreg.standard_error | logreg.z_stat | logreg.p_value |
+----------+--------------+-----------------------+---------------+----------------+
| str | float64 | float64 | float64 | float64 |
+----------+--------------+-----------------------+---------------+----------------+
| ATP6V1E1 | 1.69214e-03 | 1.28280e-03 | 1.31910e+00 | 1.87134e-01 |
| CCT8L2 | 4.03789e-02 | 3.46446e-02 | 1.16552e+00 | 2.43808e-01 |
| CECR1 | 1.29431e-03 | 1.82113e-03 | 7.10716e-01 | 4.77260e-01 |
| CECR2 | 4.90592e-04 | 3.88126e-04 | 1.26400e+00 | 2.06230e-01 |
| CECR5 | -2.90466e-04 | 3.41414e-03 | -8.50774e-02 | 9.32200e-01 |
| CECR6 | 7.20520e-03 | 7.78411e-03 | 9.25630e-01 | 3.54638e-01 |
| GAB4 | 3.14315e-03 | 2.69332e-03 | 1.16701e+00 | 2.43205e-01 |
| IL17RA | 3.40054e-03 | 4.98370e-03 | 6.82333e-01 | 4.95028e-01 |
| OR11H1 | -2.85557e-02 | 5.40733e-02 | -5.28093e-01 | 5.97435e-01 |
| POTEH | -6.64950e-03 | 3.76138e-03 | -1.76784e+00 | 7.70884e-02 |
+----------+--------------+-----------------------+---------------+----------------+

+-------------------------+----------------------+---------------------+
| logreg.fit.n_iterations | logreg.fit.converged | logreg.fit.exploded |
+-------------------------+----------------------+---------------------+
| int32 | bool | bool |
+-------------------------+----------------------+---------------------+
| 3 | true | false |
| 3 | true | false |
| 3 | true | false |
| 3 | true | false |
| 2 | true | false |
| 3 | true | false |
| 3 | true | false |
| 3 | true | false |
| 3 | true | false |
| 3 | true | false |
+-------------------------+----------------------+---------------------+
showing top 10 rows

```

---

<div class="post-metadata">

**Author:** ![enriquea](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/enriquea/32/237_2.png) [@enriquea](https://discuss.hail.is/u/enriquea)\
**Post date:** [September 11, 2018, 10:22am UTC](https://discuss.hail.is/t/logistic-regression-burden-tests/206/4 "2018-09-11T10:22:51Z")

</div>

Hi @jbloom

Thank for your quick reply. Could I aggregate by variant type (e.g. silent, damaging, etc…) and perform a logistic regression based on this same example? It looks like is possible to adapt it for performing a more global burden test using log\_regression.

Thanks in advance

---

<div class="post-metadata">

**Author:** ![jbloom](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/jbloom/32/109_2.png) [@jbloom](https://discuss.hail.is/u/jbloom)\
**Post date:** [September 12, 2018, 3:56pm UTC](https://discuss.hail.is/t/logistic-regression-burden-tests/206/5 "2018-09-12T15:56:01Z")

</div>

Absolutely, once you’ve annotated variants with their type, you can use `group_rows_by` and `aggregate` to create a type by sample matrix table (for example, the entries could count the number of variants of that type per sample) and then run logistic regression to test variant type.