19566: Normalize pca values before glm.

Arvados-DCO-1.1-Signed-off-by: Tom Clegg <tom@curii.com>

diff --git a/glm.go b/glm.go
index 6a829ca547288fd6e82c4ed754520a2acdda389f..6ae98df928ac02a64d5381ff0a4f5bd8bfb0bfe5 100644 (file)
--- a/glm.go
+++ b/glm.go
@@ -12,6 +12,7 @@ import (
 
        "github.com/kshedden/statmodel/glm"
        "github.com/kshedden/statmodel/statmodel"
+       "gonum.org/v1/gonum/stat"
        "gonum.org/v1/gonum/stat/distuv"
 )
 
@@ -22,6 +23,13 @@ var glmConfig = &glm.Config{
        Log:            log.New(io.Discard, "", 0),
 }
 
+func normalize(a []float64) {
+       mean, std := stat.MeanStdDev(a, nil)
+       for i, x := range a {
+               a[i] = (x - mean) / std
+       }
+}
+
 // Logistic regression.
 //
 // onehot is the observed outcome, in same order as sampleInfo, but
@@ -37,6 +45,7 @@ func pvalueGLM(sampleInfo []sampleInfo, onehot []bool, nPCA int) (p float64) {
                                series = append(series, si.pcaComponents[pca])
                        }
                }
+               normalize(series)
                data = append(data, series)
                pcaNames = append(pcaNames, fmt.Sprintf("pca%d", pca))
        }

diff --git a/glm_test.go b/glm_test.go
index 816b74b2c9466d84d0a02dc2341047cace819eaa..b8dbdb489949f10c631775c8b882d8c3a5af7a77 100644 (file)
--- a/glm_test.go
+++ b/glm_test.go
@@ -164,7 +164,7 @@ func (s *glmSuite) TestPvalue(c *check.C) {
 1, 1, 1, 1.23, 2.36
 1, 1, 2, 1.22, 2.32
 1, 1, 3, 1.21, 2.31
-`)), check.Equals, 0.0027896654350661053)
+`)), check.Equals, 0.002789665435066107)
 }
 
 var benchSamples, benchOnehot = func() ([]sampleInfo, []bool) {