Leave one out cross validation #101

mortonjt · 2017-01-26T05:38:05Z

There are a few things that are added to this PR namely

percent_explained : calculates the percentage of the total variance explained for each balance axis (think about the analogue for principal components).
leave one out cross validation : calculates model error on training dataset and prediction error on testing dataset
leave one variable out cross validation : calculates R^2 and model error when each variable is left out of the model. This gives us an idea about the variable importance for each variable. But more work concerning effect size of each variable will need to be done in the future

In addition, the OLS procedure was slightly optimized to reduce the amount of formula -> matrix conversions that needed to be computed. Will need to propagate this to the remaining methods. Thanks @Josh-Lefler for discovering this!

Remember, this is the last major PR that needs to be merged in. Reviews will be greatly appreciated. Thanks!

coveralls · 2017-01-26T05:50:42Z

Coverage increased (+0.05%) to 98.142% when pulling 129f81c on mortonjt:loo into abee6d8 on biocore:master.

josenavas

Thanks @mortonjt ! This looks really good. Note however that my knowledge on the math behind is rather limited, so I just checked style and small changes that can improve the code. Nothing blocking though.

josenavas · 2017-01-26T21:08:28Z

gneiss/regression/_ols.py

+def _fit_ols(y, x, **kwargs):
+    """ Perform the basic ols regression."""
+    exog_data = x
+    submodels = []


Minor, and it depends on how big len(y.columns) is but this can be written as:

submodels = [smf.OLS(endog=y[b], exog=x, **kwargs) for b in y.columns]

Appending to a list may take more than expected if the list is really long.

I didn't know that. Thanks!

josenavas · 2017-01-26T21:10:07Z

gneiss/regression/_ols.py

@@ -298,3 +318,144 @@ def r2(self):

        sst = sse + ssr
        return 1 - sse / sst
+
+    @property
+    def mse(self):


Is this property described somewhere (i.e. doc)?

I'll add it in. It is a measure of Mean of Sum of Squares Error, which is used for calculating R^2 and the cross validation statistics.

josenavas · 2017-01-26T21:14:11Z

gneiss/regression/_ols.py

+        fit
+        """
+        params = {'fit_regularized': False}
+        for key in params:


Are you expecting to grow params in the future? As it is, this for loop only executes one operation so I don't see the point of having it.

The main reason for having the kwargs argument was to pass these into the fit function, which have a ton of parameters in them. One advantage of doing this is if the fit function in statsmodels changes, we aren't forced to change our APIs to fit their changes.

josenavas · 2017-01-26T21:15:52Z

gneiss/regression/_ols.py

+                            columns=exog_names)
+        results = pd.DataFrame(index=self.balances.index,
+                               columns=['mse', 'pred_err'])
+        i = 0


I think there are 2 more pythonic way of doing this is, I'm unsure which one is the best for your case. Since i is just an index, this should work:

for sample_id, (inidx, outidx) in zip(self.balances.index, cv_iter): etc...

However, If you still need i for something else that I may be missing, then this is the pythonic way of doing the index counting:

for i, (inidx, outidx) in enumerate(cv_iter): ....

And remove the i+=1 from the end of the for loop.

Will use enumerate instead. Thanks!

josenavas · 2017-01-26T21:21:20Z

gneiss/regression/_ols.py

+
+        Parameters
+        ----------
+        regularized : bool


regularized is not named parameter of this function

I'll move this up to the params section. We do actually need this

Scratch that. I'll delete it.

josenavas · 2017-01-26T21:21:45Z

gneiss/regression/_ols.py

+        """
+
+        params = {'regularized': False}
+        for key in params:


Same comment as above regarding the for loop.

josenavas · 2017-01-26T21:22:10Z

gneiss/regression/_ols.py

+        results = pd.DataFrame(index=exog_names,
+                               columns=['mse', 'Rsquared'])
+        i = 0
+        for inidx, outidx in cv_iter:


Same comment as above regarding the usage of the i variable.

josenavas · 2017-01-26T21:26:00Z

gneiss/regression/tests/test_ols.py

+        # The L1 regularization either has numerical issues
+        # or is random ...
+        res_coefs = res.coefficients()
+        pdt.assert_index_equal(exp_coefs.index, res_coefs.index)


Minor and I'm unsure if it applies in the pandas testing functions: in general we use assert*(observed, expected), but I'm unsure if the pandas error messages are expecting this convention.

Going to remove this test since we won't be needing it.

josenavas · 2017-01-26T21:27:07Z

gneiss/regression/tests/test_ols.py

+            'y1': [-2.074816e+10, -2.039271e+09, -1.843143e+09,
+                   2.261170e+10, -7.925699e+07]
+            }, index=['Intercept', 'x1', 'x2', 'x3', 'x4']).T
+        # The L1 regularization either has numerical issues


From my ignorant view, if it is random, can you set a seed? Otherwise, can you test some of the values up to a certain decimal position?

So this is actually depending on statsmodels/statsmodels#3393

Not sure what is going on there - but maybe we shouldn't include any regularization fitting here until that issue gets resolved. I'm going to remove those fits in the next commit.

mortonjt · 2017-01-27T22:13:40Z

gneiss/regression/_ols.py

        for s in self.submodels:
            # assumes that the underlying submodels have implemented `fit`.
-            m = s.fit(**kwargs)
+            if regularized:
+                m = s.fit(**kwargs)


This is backwards. This should be fit_regularized, not fit. I'm removing the fit_regularized functionality for the time being to resolve this. But this actually messed up the test_mse test. I have fixed that in the upcoming commit

@josenavas

Addressed @josenavas's comments

josenavas · 2017-01-27T22:18:32Z

Once tests pass

coveralls · 2017-01-27T22:53:16Z

Coverage increased (+0.1%) to 98.229% when pulling ae53632 on mortonjt:loo into abee6d8 on biocore:master.

mortonjt · 2017-02-08T00:05:21Z

@sjanssen2 @qiyunzhu could you review / merge? Thanks!

qiyunzhu

Good job! Minor comments.

qiyunzhu · 2017-02-08T00:36:01Z

gneiss/regression/_ols.py

+    """ Perform the basic ols regression."""
+    # mixed effects code is obtained here:
+    # http://stackoverflow.com/a/22439820/1167475
+    submodels = [smf.OLS(endog=y[b], exog=x, **kwargs) for b in y.columns]


Maybe these two lines can be merged into one? (Not sure if what I proposed is Pythonic, if not sorry.)

qiyunzhu · 2017-02-08T00:37:10Z

gneiss/regression/_ols.py

-            # assumes that the underlying submodels have implemented `fit`.
-            m = s.fit(**kwargs)
-            self.results.append(m)
+    def fit(self, regularized=False, **kwargs):


So this rewritting is basically an optimization of the old function. Is that right?

qiyunzhu · 2017-02-08T00:38:49Z

gneiss/regression/_ols.py

+                               columns=['mse', 'Rsquared'])
+        for i, (inidx, outidx) in enumerate(cv_iter):
+            feature_id = exog_names[i]
+            res_i = _fit_ols(endog, exog.loc[:, inidx], **kwargs)


This part is bit awkward to me. res_i is defined twice with different data types. Is this Pythonic?

ok, I'll rename it

Basically, all this is doing is

Retrieving the model specifications for the linear regression (model_i)

Actually fitting the model with a numerical optimization (res_i)

Sounds good!

qiyunzhu · 2017-02-08T00:39:34Z

gneiss/regression/_ols.py

+            results.loc[feature_id, 'Rsquared'] = 1 - sse / sst
+            # degrees of freedom for residuals
+            dfe = res_i[0].df_resid
+            results.loc[feature_id, 'mse'] = sse / dfe


Can this dfe be zero?

If that ever happens, would indicate overfitting to the extreme (i.e. more parameters than data points). I think it is safe to say that this will never happen.

qiyunzhu · 2017-02-08T00:40:14Z

gneiss/regression/_ols.py

+        # Using sum of squares error calculation (df=1)
+        # instead of population variance (df=0).
+        axis_vars = np.var(self.balances, ddof=1, axis=0)
+        return axis_vars / axis_vars.sum()


Can this axis_vars.sum() be zero?

That would mean that the variance is zero. For practical reasons, that will not happen.

coveralls · 2017-02-08T07:13:24Z

Coverage increased (+0.1%) to 98.22% when pulling 89d7810 on mortonjt:loo into abee6d8 on biocore:master.

coveralls · 2017-02-08T07:14:34Z

Coverage increased (+0.1%) to 98.22% when pulling 89d7810 on mortonjt:loo into abee6d8 on biocore:master.

coveralls · 2017-02-08T07:14:55Z

Coverage increased (+0.1%) to 98.22% when pulling 89d7810 on mortonjt:loo into abee6d8 on biocore:master.

mortonjt added 5 commits January 25, 2017 21:15

Adding test cases for loo

ddbcf1c

Adding lovo, percent_explained and mse tests

bbb6a02

flake8

c4f47fd

flake8

c440608

fixing dumb numerical issues

129f81c

josenavas reviewed Jan 26, 2017

View reviewed changes

mortonjt commented Jan 27, 2017

View reviewed changes

mortonjt added 2 commits January 27, 2017 14:13

Removing regularization

be4633e

Addressed @josenavas's comments

flake9

823eccd

josenavas approved these changes Jan 27, 2017

View reviewed changes

mortonjt added 2 commits January 27, 2017 14:39

fixing precision in test

cc227a8

Bumping down precision

ae53632

mortonjt mentioned this pull request Jan 28, 2017

Lme tweaks #104

Closed

qiyunzhu reviewed Feb 8, 2017

View reviewed changes

mortonjt added 3 commits February 7, 2017 23:10

Update _ols.py

353f1b7

renaming first res_i to model_i

f2bd9fb

renaming first res_i to model_i

89d7810

qiyunzhu merged commit 173bf10 into biocore:master Feb 8, 2017

Leave one out cross validation #101

Leave one out cross validation #101

Conversation

mortonjt commented Jan 26, 2017 • edited Loading

coveralls commented Jan 26, 2017

josenavas left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

josenavas commented Jan 27, 2017

coveralls commented Jan 27, 2017

mortonjt commented Feb 8, 2017

qiyunzhu left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

coveralls commented Feb 8, 2017

coveralls commented Feb 8, 2017

coveralls commented Feb 8, 2017

mortonjt commented Jan 26, 2017 •

edited

Loading