- model, err = glm.NewGLM(dataset, "outcome", append([]string{"constants", "variant"}, pcaNames...), glmConfig)
- if err != nil {
- return math.NaN()
+
+ return func(onehot []bool) (p float64) {
+ defer func() {
+ if recover() != nil {
+ // typically "matrix singular or near-singular with condition number +Inf"
+ p = math.NaN()
+ }
+ }()
+
+ variant := make([]statmodel.Dtype, 0, len(sampleInfo))
+ row := 0
+ for _, si := range sampleInfo {
+ if si.isTraining {
+ if onehot[row] {
+ variant = append(variant, 1)
+ } else {
+ variant = append(variant, 0)
+ }
+ row++
+ }
+ }
+
+ data := append([][]statmodel.Dtype{data[0], variant}, data[1:]...)
+ names := append([]string{"outcome", "variant"}, names[1:]...)
+ dataset := statmodel.NewDataset(data, names)
+
+ model, err := glm.NewGLM(dataset, "outcome", names[1:], glmConfig)
+ if err != nil {
+ return math.NaN()
+ }
+ resultComp := model.Fit()
+ logComp := resultComp.LogLike()
+ dist := distuv.ChiSquared{K: 1}
+ return dist.Survival(-2 * (logCov - logComp))