glm.go

   1 // Copyright (C) The Lightning Authors. All rights reserved.
   2 //
   3 // SPDX-License-Identifier: AGPL-3.0
   4
   5 package lightning
   6
   7 import (
   8         "fmt"
   9         "math"
  10
  11         "github.com/kshedden/statmodel/glm"
  12         "github.com/kshedden/statmodel/statmodel"
  13 )
  14
  15 var glmConfig = &glm.Config{
  16         Family:         glm.NewFamily(glm.BinomialFamily),
  17         FitMethod:      "IRLS",
  18         ConcurrentIRLS: 1000,
  19 }
  20
  21 func pvalueGLM(sampleInfo []sampleInfo, onehot []bool) float64 {
  22         nPCA := len(sampleInfo[0].pcaComponents)
  23         pcaNames := make([]string, 0, nPCA)
  24         data := make([][]statmodel.Dtype, 0, nPCA)
  25         for pca := 0; pca < nPCA; pca++ {
  26                 series := make([]statmodel.Dtype, 0, len(sampleInfo))
  27                 for _, si := range sampleInfo {
  28                         if si.isTraining {
  29                                 series = append(series, si.pcaComponents[pca])
  30                         }
  31                 }
  32                 data = append(data, series)
  33                 pcaNames = append(pcaNames, fmt.Sprintf("pca%d", pca))
  34         }
  35
  36         variant := make([]statmodel.Dtype, 0, len(sampleInfo))
  37         outcome := make([]statmodel.Dtype, 0, len(sampleInfo))
  38         for row, si := range sampleInfo {
  39                 if si.isTraining {
  40                         if onehot[row] {
  41                                 variant = append(variant, 1)
  42                         } else {
  43                                 variant = append(variant, 0)
  44                         }
  45                         if si.isCase {
  46                                 outcome = append(outcome, 1)
  47                         } else {
  48                                 outcome = append(outcome, 0)
  49                         }
  50                 }
  51         }
  52         data = append(data, variant, outcome)
  53
  54         dataset := statmodel.NewDataset(data, append(pcaNames, "variant", "outcome"))
  55         model, err := glm.NewGLM(dataset, "outcome", pcaNames, glmConfig)
  56         if err != nil {
  57                 return math.NaN()
  58         }
  59         resultCov := model.Fit()
  60         logCov := resultCov.LogLike()
  61         model, err = glm.NewGLM(dataset, "outcome", append([]string{"variant"}, pcaNames...), glmConfig)
  62         if err != nil {
  63                 return math.NaN()
  64         }
  65         resultComp := model.Fit()
  66         logComp := resultComp.LogLike()
  67         return chisquared.Survival(-2 * (logCov - logComp))
  68 }