Add F test for comparing two models #182

LewisHein · 2017-05-17T17:27:43Z

Here is a rough first draft of what I'm thinking for the F-test for model comparison. If you all think I am on the right track, I'll write up some docs. The addition to the test suite tests a single ANOVA run against the p-value from R's implementation. This is for addressing #181

…redModel...)

nalimilan · 2017-05-17T19:35:34Z

Thanks. I think the function should return a table similar to CoefTable rather than a list of objects when more than two models are passed. That would be more convenient for printing as well as to extract values.

nalimilan · 2017-05-17T19:36:47Z

src/ftest.jl

@@ -0,0 +1,75 @@
+type FTestResult
+    mod1::RegressionModel


Should probably not store the original models, that can take a lot of space and can be tracked separately if needed.

nalimilan · 2017-05-17T19:38:11Z

src/ftest.jl

+    end
+end
+
+export ftest


Exports should go to StatsBase.jl.

Not sure what you mean here -- I can't export ftest until I've defined it.

Sure you can.

Bindings in export statements aren't resolved immediately. As an example,

julia> baremodule Test export A # has no definition end Test julia> using Test # doesn't error even though A is exported and undefined julia> Test.A ERROR: UndefVarError: A not defined

Julia doesn't care until you try to call it.

Yup, you're right. It appears that I made a typo when I tried it the first time

nalimilan · 2017-05-17T19:44:07Z

src/ftest.jl

+end
+
+function ftest(mod1::RegressionModel, mod2::RegressionModel)
+    SSR1 = deviance(mod1.model.rr)


Should call deviance on mod1 to be fully generic.

nalimilan · 2017-05-17T19:44:25Z

src/ftest.jl

+    SSR1 = deviance(mod1.model.rr)
+    SSR2 = deviance(mod2.model.rr)
+
+    nparams1 = length(mod1.model.pp.beta0)


Use dof(mod1).

nalimilan · 2017-05-17T19:45:17Z

src/ftest.jl

+    results::Array{FTestResult, 1}
+end
+
+function ftest(mod1::RegressionModel, mod2::RegressionModel)


The F test is only valid for linear regressions, so the signature should be stricter.

nalimilan · 2017-05-17T19:46:13Z

src/ftest.jl

+
+function ftest(mods::RegressionModel...)
+    nmodels = length(mods)
+    results = Array{FTestResult, 1}((nmodels^2)-nmodels)


What R's anova and Stata's nestreg do is to compare each model to the previous one. Computing all possible comparisons quickly gets very messy.

nalimilan · 2017-05-17T19:49:34Z

It would also be nice to check that the models were fitted on the same data, at least by checking that the response vectors are equal (and maybe even that all common variables are equal).

LewisHein · 2017-05-17T21:19:11Z

It would also be nice to check that the models were fitted on the same data, at least by checking that the response vectors are equal (and maybe even that all common variables are equal).

This would be great; however, I am not sure how to get the common variables of the models. I am pretty new to the internals of GLM.jl, so please forgive my ignorance.

ararslan · 2017-05-17T21:19:20Z

src/GLM.jl

@@ -53,7 +53,8 @@ module GLM
        nobs,           # total number of observations
        predict,        # make predictions
        updateμ!,      # update the response type from the linear predictor
-        wrkresp         # working response
+        wrkresp,        # working response
+	ftest		# compare two models with an F test


Please indent with spaces instead of tabs

ararslan · 2017-05-17T21:19:43Z

src/ftest.jl

+	    MSR1s[i] = result.MSR1
+	    MSR2s[i] = result.MSR2
+	    fstats[i] = result.fstat
+	    pvals[i] = result.pval


Same for these. Indents should be 4 spaces each.

ararslan · 2017-05-17T21:20:04Z

src/ftest.jl

    end

-    return MultiFTestResult(results)
+    return CoefTable([SSR1s, SSR2s, df1s, df2s, MSR1s, MSR2s, fstats, pvals], ["Model 1 SSR", "Model 2 SSR", "Model 1 df", "Model 2 df", "Model 1 MSR", "Model 2 MSR", "F statistic", "p-value"], ["Model $(i-1):$i" for i in 2:nmodels])


It would be good to break this long line into multiple lines for readability.

LewisHein · 2017-05-17T21:29:25Z

Also, switching to LinPredModel has caused my tests to fail, with an error message about having no method matching ftest(::DataFrames.DataFrameRegressionModel{GLM.LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredQR{Float64}},Array{Float64,2}}, ::DataFrames.DataFrameRegressionModel{GLM.LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredQR{Float64}},Array{Float64,2}}).

~~How~~~~ ~~~~do~~~~ ~~~~I~~~~ ~~~~fix~~~~ ~~~~this?~~

EDIT: Got it. and it was simple, too

nalimilan · 2017-05-18T12:21:22Z

src/ftest.jl

+        pvals[i] = result.pval
+    end
+
+    return CoefTable([SSR1s, SSR2s, df1s, df2s, MSR1s, MSR2s, fstats, pvals],


Sorry if I wasn't clear, I didn't mean to (ab)use CoefTable for this, but to take inspiration from it to define the MultiFTestResult type (could be called FTestTable BTW). That type would have the same fields as FTestResult, but those would be vectors rather than scalars (one column for each info).

I'm not even sure we need FTestResult for the special case of a two-model comparison. Opinions?

This isn't an opinion, but rather an explanation. I defined FTestResult purely so that I could define show() for it and have the result nicely printed.

Is there some more standard way to do this?

EDIT: @nalimilan I finally see your point: define just FTestTable, and for comparing two models, return an FTestTable with only one test in it. Correct?

Not totally "standard", but CoefTable (in StatsBase) is a very similar case (lightweight structure to be used where R would use a data frame, which is overkill in Julia).

nalimilan · 2017-05-18T12:28:56Z

This would be great; however, I am not sure how to get the common variables of the models. I am pretty new to the internals of GLM.jl, so please forgive my ignorance.

Unfortunately, you can't get the names of the coefficients at this level, they are only available for DataFrameRegressionModel objects (or other models with a model frame). We could add a specific method which would check this later in DataFrames/StatsModels, but that's beyond the scope of this PR.

What can be done in GLM.jl is to check that the fields of the LmResp object are equal, which should catch most problems.

LewisHein · 2017-05-18T12:32:42Z

Well, this morning with a fresh perspective, I tackled the submodel test. Tell me what you think of my solution.

LewisHein · 2017-05-18T12:45:22Z

But can't you just do return false here?

Where?

nalimilan

Thanks, this looks quite good. Would be nice to check whether any allocation happens when calling ftest (though that can be fixed later).

We also need to solve the naming issue.

nalimilan · 2017-06-23T14:36:32Z

src/ftest.jl

+_diff{N, T}(t::NTuple{N, T})::NTuple{N, T} =  ntuple(i->t[i+1]-t[i], N-1)
+
+import Base: ./
+./{N, T1, T2}(t1::NTuple{N, T1}, t2::NTuple{N, T2}) = ntuple(i->t1[i]/t2[i], N)


Unfortunately you can't do that as that will affect other packages. You can just use ntuple directly, or define a custom function.

nalimilan · 2017-06-23T14:37:30Z

src/ftest.jl

+    pval::Tuple{Vararg{PValue}}
+end
+
+ #function FTestResult{N}(SSR1::Array{Float64, 1}, fstati


Remove this comment (and others), and fix indentation below.

nalimilan · 2017-06-23T14:38:20Z

src/ftest.jl

+
+    nc = 10
+    nr = N
+    outrows = Array{String, 2}(nr+1, nc)


Matrix{String} is slightly nicer to read.

nalimilan · 2017-06-23T14:39:34Z

src/ftest.jl

+            cur_cell_len = length(cur_cell)
+
+            print(io, cur_cell)
+            print(io, RepString(" ", max_colwidths[c]-cur_cell_len+1))


Use " "^(max_colwidths[c]-cur_cell_len+1) as RepString is likely to be deprecated at some point.

nalimilan · 2017-06-23T14:40:15Z

src/ftest.jl

+
+    outrows[1, :] = ["", "Res. DOF",  "DOF",  "ΔDOF",  "SSR",
+                    "ΔSSR",  "R²",  "ΔR²",  "F*",  "p(>F)"]
+    outrows[2, :] = ["Model 1", @sprintf("%.4f", ftr.dof_resid[1]),


Please use a more regular organization, it's very hard to read as-is. Same below.

Also, DOF and number of parameters should be printed as integers, not floats.

nalimilan · 2017-06-23T14:44:49Z

test/runtests.jl

+    @test sprint(show, ftest(mod, nullmod)) == 
+        """
+                Res. DOF DOF    ΔDOF    SSR    ΔSSR    R²      ΔR²    F*       p(>F) 
+        Model 1 10.0000  3.0000         0.1283         0.9603                        


Not a big deal, but I wonder whether it would be possible to force @printf to align number on the decimal separator (see e.g. the R² column).

Maybe, but I'd have to play with it a good bit to figure out how.

Sure, that was just an idea. I know CoefTable has the same issue.

nalimilan · 2017-06-23T14:46:25Z

src/ftest.jl

@@ -0,0 +1,155 @@
+type FTestResult{N}
+    ssr::NTuple{N, Float64}
+    nparams::NTuple{N, Int}


Probably better call this dof, as it's the result of calling dof(model)? Also params is going to have a slightly different meanings from coef (JuliaStats/StatsBase.jl#274).

nalimilan

OK, with the next round we should be ready if @andreasnoack validates and we agree on a name.

nalimilan · 2017-06-23T16:51:42Z

src/ftest.jl


 _diffn{N, T}(t::NTuple{N, T})::NTuple{N, T} =  ntuple(i->t[i]-t[i+1], N-1)

 _diff{N, T}(t::NTuple{N, T})::NTuple{N, T} =  ntuple(i->t[i+1]-t[i], N-1)

-import Base: ./
-./{N, T1, T2}(t1::NTuple{N, T1}, t2::NTuple{N, T2}) = ntuple(i->t1[i]/t2[i], N)
+dividetuple{N, T1, T2}(t1::NTuple{N, T1}, t2::NTuple{N, T2}) = ntuple(i->t1[i]/t2[i], N)


T1 and T2 are not needed. Maybe call this dividetuples?

nalimilan · 2017-06-23T16:52:41Z

src/ftest.jl

-    outrows[2, :] = ["Model 1", @sprintf("%.4f", ftr.dof_resid[1]),
-                     @sprintf("%.4f", ftr.nparams[1]), " ", @sprintf("%.4f", ftr.ssr[1]),
-                     " ", @sprintf("%.4f", ftr.r2[1]), " ", " ", " "]
+    outrows[1, :] = [


I'd group these a bit more, e.g. on two or three lines.

nalimilan · 2017-06-23T16:54:48Z

test/runtests.jl

-        Model 1 10.0000  3.0000         0.1283         0.9603                        
-        Model 2 11.0000  2.0000 -1.0000 3.2292 -3.1008 -0.0000 0.9603 241.6234 <1e-7 
+                Res. DOF DOF ΔDOF SSR    ΔSSR    R²      ΔR²    F*       p(>F) 
+        Model 1 10       3        0.1283         0.9603                        


Can you align integers on the right? That makes it easier to read by columns since digits of the same rank are aligned. (BTW, if you do the same for floats, they will be automatically aligned on the decimal separator since their are fixed to four digits after it.)

nalimilan · 2017-06-23T16:56:38Z

src/ftest.jl

+                     "p(>F)"
+                     ]
+
+    outrows[2, :] = [


Here too you can put an item on the first line, and the closing bracket with the last element. Same below.

nalimilan · 2017-06-25T18:29:15Z

src/ftest.jl

+                     "R²", "ΔR²", "F*", "p(>F)"]
+
+    outrows[2, :] = ["Model 1", @sprintf("%.0d", ftr.dof_resid[1]),
+                    @sprintf("%.0d", ftr.dof[1]), " ",


Should be aligned with " (as it's part of the array, so inside [).

…nto add_anova

ararslan · 2017-06-27T19:52:27Z

src/ftest.jl

+    nmodels = length(mods)
+    for i in 2:nmodels
+        issubmodel(mods[i], mods[i-1]) || 
+        throw(ArgumentError("F test $i is only valid if model $i is nested in model $i-1"))


Should this be $(i-1)?

Yes, good catch

nalimilan · 2017-06-27T19:53:37Z

src/ftest.jl

+        print(io, "\n")
+    end
+end
+


Since it looks like we'll need another round of fixes, it would be nice to remove these two empty lines in your next (and hopefully last) commit.

ararslan · 2017-06-27T19:59:06Z

src/ftest.jl

+Model 1       10   3      0.1283          0.9603                      
+Model 2       11   2   -1 3.2292 -3.1008 -0.0000 0.9603 241.6234 <1e-7
+```
+"""


I'd format it like this:

""" ftest(mod::LinearModel...) For each sequential pair of linear predictors in `mod`, perform an F-test to determine if the first one fits significantly better than the next. A table is returned containing residual degrees of freedom (DOF), degrees of freedom, difference in DOF from the preceding model, sum of squared residuals (SSR), difference in SSR from the preceding model, R², difference in R² from the preceding model, and F-statistic and p-value for the comparison between the two models. !!! note This function can be used to perform an ANOVA by testing the relative fit of two models to the data # Examples Suppose we want to compare the effects of two or more treatments on some result. Because this is an ANOVA, our null hypothesis is that `Result~1` fits the data as well as `Result~Treatment`. ```jldoctest julia> dat = DataFrame(Treatment=[1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2.], Result=[1.1, 1.2, 1, 2.2, 1.9, 2, .9, 1, 1, 2.2, 2, 2]); julia> mod = lm(@formula(Result~Treatment), dat); julia> nullmod = lm(@formula(Result~1), dat); julia> ft = ftest(mod.model, nullmod.model) Res. DOF DOF ΔDOF SSR ΔSSR R² ΔR² F* p(>F) Model 1 10 3 0.1283 0.9603 Model 2 11 2 -1 3.2292 -3.1008 -0.0000 0.9603 241.6234 <1e-7 """

There are a few modifications in here, too lazy to enumerate them rather than just copy-pasting and editing :P

nalimilan · 2017-06-27T20:03:30Z

src/ftest.jl

+    return true
+end
+
+_diffn{N, T}(t::NTuple{N, T})::NTuple{N, T} =  ntuple(i->t[i]-t[i+1], N-1)


::NTuple{N, T} is actually wrong here, since the result has N-1 elements. Just drop the type assertion (same below).

And BTW, that's the cause of the tests failure.

There's also a double space.

nalimilan · 2017-06-27T20:57:42Z

Note that I resolved conflicts with master directly via GitHub, so you'll have either to pull from the branch, or to rebase locally against master and then force-push.

LewisHein · 2017-06-27T21:05:31Z

OK -- fingers crossed...

LewisHein · 2017-06-27T22:27:04Z

I have no idea what happened. I ran the tests before submitting, but now runtests.jl is all messed up.

LewisHein · 2017-06-27T23:28:22Z

Now travis says it's failing on nightly Julia. What now?

nalimilan · 2017-06-28T06:59:39Z

Don't worry about the nightlies, the failure was already there. I'll let @ararslan merge if he's OK with the new version.

nalimilan · 2017-06-29T07:57:41Z

Thanks @LewisHein! Sorry if the process was a bit long. If you want to implement the likelihood ratio test at some point, now that we have settled the API, that would be great!

pbastide · 2017-09-19T17:18:19Z

Hi @LewisHein, @nalimilan and @ararslan,

Thank you for all your work on the GLM package !

We would need the ftest function in the PhyloNetwoks package, that already depends on GLM. Right now we use our own naive implementation of this function, but we would be happy to rely on your release in the future.

However, if I'm correct, right now v0.7.0 does not contain the changes you made here. Are you planing on making this feature available in a release soon ? That would be great for us, so that we could require this new version, and use your function.

Sorry if this is not the right place to ask, or if the answer is already somewhere else. If so, I'd be grateful for you to point me in the right direction.

Thank you !

nalimilan · 2017-09-19T18:17:50Z

Sure. I've just tagged a new release, see JuliaLang/METADATA.jl#11275.

pbastide · 2017-09-20T06:43:34Z

Thank you, that was fast !

Lewis Hein added 4 commits May 17, 2017 11:24

Added F test for comparing two models

1c0bd1e

Fixed most changes from comments by nalimilan

c2e7b7c

Changed spaces to tabs in GLM.jl

cd201c6

Spaces->tabs in ftest.jl and broke long line at the end of ftest(LinP…

f6d7fc4

…redModel...)

nalimilan reviewed May 17, 2017

View reviewed changes

ararslan reviewed May 17, 2017

View reviewed changes

Lewis Hein added 9 commits May 18, 2017 06:32

Changed tests and added check for nesting in ftest

e05b10f

Reformatted comments

885ed1b

Renamed variables in issubmodel

8133560

Performance enhancements and tests for issubmodel

5ff9e05

Remove extra @inbounds in issubmodel

7641032

Switched to always returning FTestTable

6d2c4c1

First draft of documentation

da837d5

Reformat docs and re-order method defs

8970a82

Bugfixes and tests for multi-argument ftest()

5dfd20a

nalimilan reviewed May 18, 2017

View reviewed changes

nalimilan reviewed Jun 23, 2017

View reviewed changes

Lewis Hein added 2 commits June 25, 2017 11:52

Addressed latest comments and fixed doctest

c213f6b

Fixed indentation in ftest.jl

256d848

nalimilan reviewed Jun 25, 2017

View reviewed changes

Lewis Hein and others added 4 commits June 27, 2017 14:59

Docs cleanup and other misc. fixes

81b6baf

Fixed runtests.jl

c94e4ee

Merge branch 'master' into pull-request/1c0bd1e2

78a8a8e

Merge branch 'pull-request/1c0bd1e2' of github.com:LewisHein/GLM.jl i…

dcc5b0b

…nto add_anova

nalimilan approved these changes Jun 27, 2017

View reviewed changes

ararslan reviewed Jun 27, 2017

View reviewed changes

nalimilan reviewed Jun 27, 2017

View reviewed changes

ararslan reviewed Jun 27, 2017

View reviewed changes

nalimilan reviewed Jun 27, 2017

View reviewed changes

Slight docstring adjustments

5e165cd

ararslan merged commit 5392a29 into JuliaStats:master Jun 28, 2017

This was referenced Jul 29, 2017

add effects and ANOVAtest function for lm #70

Closed

Implement ANOVAs for lms #181

Closed

johnmyleswhite mentioned this pull request Jul 29, 2017

Add additional methods #7

Open

35 tasks

nalimilan mentioned this pull request Sep 11, 2017

Additional Methods for RegressionModel JuliaStats/StatsBase.jl#300

Open

nalimilan mentioned this pull request Dec 1, 2017

Tests for RegressionModel JuliaStats/HypothesisTests.jl#121

Open

Add F test for comparing two models #182

Add F test for comparing two models #182

Conversation

LewisHein commented May 17, 2017 • edited Loading

nalimilan commented May 17, 2017

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

nalimilan commented May 17, 2017

LewisHein commented May 17, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

LewisHein commented May 17, 2017 • edited Loading

Choose a reason for hiding this comment

LewisHein May 18, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nalimilan commented May 18, 2017

LewisHein commented May 18, 2017

LewisHein commented May 18, 2017 • edited Loading

nalimilan 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

nalimilan 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

nalimilan commented Jun 27, 2017

LewisHein commented Jun 27, 2017

LewisHein commented Jun 27, 2017

LewisHein commented Jun 27, 2017

nalimilan commented Jun 28, 2017

nalimilan commented Jun 29, 2017

pbastide commented Sep 19, 2017

nalimilan commented Sep 19, 2017 • edited Loading

pbastide commented Sep 20, 2017

LewisHein commented May 17, 2017 •

edited

Loading

LewisHein commented May 17, 2017 •

edited

Loading

LewisHein May 18, 2017 •

edited

Loading

LewisHein commented May 18, 2017 •

edited

Loading

nalimilan commented Sep 19, 2017 •

edited

Loading