Add a formula(model) accessor, loglikelihood(model, :), and update tests (#226)

* add formula accessor method

* bump patch

* bump dataframes compat to 1

* fix tests

* trying to fix version problems on 1.0

* how about this one

* add tests for loglikelihood method

Author: kleinschmidt

Project.toml

@@ -1,6 +1,6 @@
 name = "StatsModels"
 uuid = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
-version = "0.6.22"
+version = "0.6.23"
 
 [deps]
 DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
@@ -14,11 +14,11 @@ StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [compat]
-CategoricalArrays = "0.8"
+CategoricalArrays = "0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.10"
 DataAPI = "1.1"
-DataFrames = "0.21"
-DataStructures = "0.17.0, 0.18"
-ShiftedArrays = "1.0.0"
+DataFrames = "1"
+DataStructures = "0.17, 0.18"
+ShiftedArrays = "1"
 StatsBase = "0.33.5"
 StatsFuns = "0.9"
 Tables = "0.2, 1"

docs/src/api.md

@@ -64,6 +64,7 @@ fit
 response
 modelmatrix
 lrtest
+formula
 ```
 
 ### Traits

src/StatsModels.jl

@@ -35,6 +35,7 @@ export
     coefnames,
     dropterm,
     setcontrasts!,
+    formula,
 
     AbstractTerm,
     ConstantTerm,

src/statsmodel.jl

@@ -96,6 +96,16 @@ for (modeltype, dfmodeltype) in ((:StatisticalModel, TableStatisticalModel),
     end
 end
 
+"""
+    formula(model)
+
+Retrieve formula from a fitted or specified model
+"""
+function formula end
+
+formula(m::TableStatisticalModel) = m.mf.f
+formula(m::TableRegressionModel) = m.mf.f
+
 @doc """
     fit(Mod::Type{<:StatisticalModel}, f::FormulaTerm, data, args...;
         contrasts::Dict{Symbol}, kwargs...)
@@ -132,6 +142,7 @@ StatsBase.r2(mm::TableRegressionModel) = r2(mm.model)
 StatsBase.adjr2(mm::TableRegressionModel) = adjr2(mm.model)
 StatsBase.r2(mm::TableRegressionModel, variant::Symbol) = r2(mm.model, variant)
 StatsBase.adjr2(mm::TableRegressionModel, variant::Symbol) = adjr2(mm.model, variant)
+StatsBase.loglikelihood(mm::TableModels, c::Colon) = loglikelihood(mm.model, c)
 
 function _return_predictions(T, yp::AbstractVector, nonmissings, len)
     out = Vector{Union{eltype(yp),Missing}}(missing, len)

test/contrasts.jl

@@ -227,7 +227,7 @@
         f_sdiff = apply_schema(f, schema(d2, Dict(:x => sdiff_hyp)))
         f_effects = apply_schema(f, schema(d2, Dict(:x => effects_hyp)))
 
-        y_means = by(d2, :x, :y => mean).y_mean
+        y_means = combine(groupby(d2, :x), :y => mean).y_mean
 
         y, X_sdiff = modelcols(f_sdiff, d2)
         @test X_sdiff \ y ≈ [mean(y_means); diff(y_means)]

test/modelmatrix.jl

@@ -159,7 +159,6 @@
         d = DataFrame(x = repeat([:a, :b], outer = 4),
                       y = repeat([:c, :d], inner = 2, outer = 2),
                       z = repeat([:e, :f], inner = 4))
-        categorical!(d)
         cs = Dict([Symbol(name) => EffectsCoding() for name in names(d)])
         d.n = 1.:8
 

test/runtests.jl

@@ -4,6 +4,7 @@ using SparseArrays
 
 using StatsModels
 using DataFrames
+using CategoricalArrays
 using StatsBase
 
 using StatsModels: ContrastsMatrix

test/statsmodel.jl

@@ -47,6 +47,8 @@ StatsBase.deviance(mod::DummyMod) = sum((response(mod) .- predict(mod)).^2)
 # Incorrect but simple definition
 StatsModels.isnested(mod1::DummyMod, mod2::DummyMod; atol::Real=0.0) =
     dof(mod1) <= dof(mod2)
+StatsBase.loglikelihood(mod::DummyMod) = -sum((response(mod) .- predict(mod)).^2)
+StatsBase.loglikelihood(mod::DummyMod, ::Colon) = -(response(mod) .- predict(mod)).^2
 
 # A dummy RegressionModel type that does not support intercept
 struct DummyModNoIntercept <: RegressionModel
@@ -122,6 +124,12 @@ Base.show(io::IO, m::DummyModTwo) = println(io, m.msg)
     ## coefnames delegated to model frame by default
     @test coefnames(m) == coefnames(ModelFrame(f, d)) == ["(Intercept)", "x1", "x2", "x1 & x2"]
 
+    @test formula(m) == m.mf.f
+
+    ## loglikelihood methods from StatsBase
+    @test length(loglikelihood(m, :)) == nrow(d)
+    @test sum(loglikelihood(m, :)) == loglikelihood(m) == -deviance(m)
+
     ## test prediction method
     ## vanilla
     @test predict(m) == [ ones(size(d,1)) Array(d.x1) Array(d.x2) Array(d.x1).*Array(d.x2) ] * collect(1:4)