Add functions dataset and datasets (#254)

* Add non-exported dataset and datasets functions, depend on Feather

* Change test scripts to use dataset

Author: dmbates

Project.toml

@@ -7,6 +7,7 @@ version = "2.2.0"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+Feather = "becb17da-46f6-5d3c-ad1b-1c5fe96bc73c"
 GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd"
@@ -29,6 +30,7 @@ TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"
 BlockArrays = "0.5, 0.6, 0.7, 0.8, 0.9, 0.10, 0.11"
 BlockDiagonals = "0.1"
 Distributions = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22"
+Feather = "0.5"
 GLM = "1"
 NLopt = "0.5"
 NamedArrays = "0.9"

src/MixedModels.jl

@@ -3,6 +3,7 @@ module MixedModels
 using BlockArrays
 using BlockDiagonals
 using Distributions
+using Feather
 using GLM
 using LinearAlgebra
 using NamedArrays

src/utilities.jl

@@ -115,3 +115,25 @@ function replicate(f::Function, n::Integer; use_threads=false)
     end
     results
 end
+
+"""
+    dataset(nm)
+
+Return the data frame of test data set named `nm`, which can be a `String` or `Symbol`
+"""
+function dataset(nm::AbstractString)
+    path = joinpath(TestData, nm * ".feather")
+    if !isfile(path)
+        throw(ArgumentError(
+            "Dataset \"$nm\" is not available.\nUse MixedModels.datasets() for available names."))
+    end
+    Feather.read(path)
+end
+dataset(nm::Symbol) = dataset(string(nm))
+
+"""
+    datasets()
+
+Return a vector of names of the available test data sets
+"""
+datasets() = first.(Base.Filesystem.splitext.(filter(Base.Fix2(endswith, ".feather"), readdir(TestData))))

test/FactorReTerm.jl

@@ -9,15 +9,12 @@ using Test
 
 const LMM = LinearMixedModel
 
-data(nm::AbstractString) = Feather.read(joinpath(MixedModels.TestData, string(nm, ".feather")))
-data(nm::Symbol) = data(string(nm))
-
 @testset "scalarReMat" begin
-    ds = data("dyestuff")
+    ds = MixedModels.dataset("dyestuff")
     f1 = @formula(yield ~ 1 + (1|batch))
     y1, Xs1 = modelcols(apply_schema(f1, schema(ds), LMM), ds)
     sf = Xs1[2]
-    psts = data("pastes")
+    psts = MixedModels.dataset("pastes")
     f2 = @formula(strength ~ 1 + (1|sample) + (1|batch))
     y2, Xs2 = modelcols(apply_schema(f2, schema(psts), LMM), psts)
     sf1 = Xs2[2]
@@ -83,7 +80,7 @@ data(nm::Symbol) = data(string(nm))
 end
 
 @testset "RandomEffectsTerm" begin
-    slp = data("sleepstudy")
+    slp = MixedModels.dataset("sleepstudy")
     contrasts =  Dict{Symbol,Any}()
 
     @testset "Detect same variable as blocking and experimental" begin
@@ -102,7 +99,7 @@ end
 
 @testset "Categorical Blocking Variable" begin
     # deepcopy because we're going to modify it
-    slp = deepcopy(data("sleepstudy"))
+    slp = deepcopy(MixedModels.dataset("sleepstudy"))
     contrasts =  Dict{Symbol,Any}()
     f = @formula(reaction ~ 1 + (1|subj))
 

test/UniformBlockDiagonal.jl

@@ -7,9 +7,6 @@ using SparseArrays
 using StatsModels
 using Test
 
-dat(nm::AbstractString) = Feather.read(joinpath(MixedModels.TestData, nm * ".feather"))
-dat(nm::Symbol) = dat(string(nm))
-
 const LMM = LinearMixedModel
 
 @testset "UBlk" begin
@@ -59,7 +56,7 @@ const LMM = LinearMixedModel
 
     @testset "updateL" begin
         @test ones(2, 2) == MixedModels.rankUpdate!(Hermitian(zeros(2, 2)), ones(2))
-        d3 = dat(:d3)
+        d3 = MixedModels.dataset(:d3)
         sch = schema(d3)
         vf1 = modelcols(apply_schema(@formula(y ~ 1 + u + (1+u|g)), sch, LMM), d3)[2][2]
         vf2 = modelcols(apply_schema(@formula(y ~ 1 + u + (1+u|h)), sch, LMM), d3)[2][2]

test/fit.jl

@@ -1,31 +1,27 @@
 using MixedModels, Feather, Test
 
-data(nm::AbstractString) = Feather.read(joinpath(MixedModels.TestData, string(nm, ".feather")))
-
-data(nm::Symbol) = data(string(nm))
-
 @testset "linear" begin
-    m1 = fit(MixedModel, @formula(yield ~ 1 + (1|batch)), data(:dyestuff))
+    m1 = fit(MixedModel, @formula(yield ~ 1 + (1|batch)), MixedModels.dataset(:dyestuff))
     @test first(m1.θ) ≈ 0.7525806757718846 rtol=1.0e-5
 end
 
 @testset "generalized" begin
     gm1 = fit(MixedModel, @formula(use ~ 1 + urban + livch + age + abs2(age) + (1|dist)),
-              data(:contra), Bernoulli())
+              MixedModels.dataset(:contra), Bernoulli())
     @test deviance(gm1) ≈ 2372.7286 atol=1.0e-3
 end
 
 @testset "Normal-IdentityLink" begin
-    @test isa(fit(MixedModel, @formula(yield ~ 1 + (1|batch)), data(:dyestuff), Normal()),
+    @test isa(fit(MixedModel, @formula(yield ~ 1 + (1|batch)), MixedModels.dataset(:dyestuff), Normal()),
               LinearMixedModel)
     @test_throws(ArgumentError("use LinearMixedModel for Normal distribution with IdentityLink"),
                  fit(GeneralizedLinearMixedModel,
                      @formula(yield ~ 1 + (1|batch)),
-                     data(:dyestuff)))
+                     MixedModels.dataset(:dyestuff)))
 end
 
 @testset "Normal Distribution GLMM" begin
-    @test_broken(isa(fit(MixedModel, @formula(yield ~ 1 + (1|batch)), data(:dyestuff),
+    @test_broken(isa(fit(MixedModel, @formula(yield ~ 1 + (1|batch)), MixedModels.dataset(:dyestuff),
                          Normal(), LogLink),
                      GeneralizedLinearMixedModel))
 end

test/missing.jl

@@ -1,9 +1,7 @@
 using Feather, MixedModels, Test
 
-dat(nm::AbstractString) = Feather.read(joinpath(MixedModels.TestData, string(nm, ".feather")))
-dat(nm::Symbol) = dat(string(nm))
 # deepcopy because we're going to modify it
-slp = deepcopy(dat(:sleepstudy))
+slp = deepcopy(MixedModels.dataset(:sleepstudy))
 slp[!,:days] = Array{Union{Missing, Float64},1}(slp[!,:days])
 slp[1,:days] = missing
 
@@ -21,13 +19,13 @@ slp[1,:days] = missing
 @testset "Missing Omit" begin
     @testset "Missing from unused variables" begin
         # missing from unused variables should have no impact
-        m1 = fit(MixedModel, @formula(reaction ~ 1 + (1|subj)), dat(:sleepstudy))
+        m1 = fit(MixedModel, @formula(reaction ~ 1 + (1|subj)), MixedModels.dataset(:sleepstudy))
         m1_missing = fit(MixedModel, @formula(reaction ~ 1 + (1|subj)), slp)
         @test isapprox(m1.θ, m1_missing.θ, rtol=1.0e-12)
     end
 
     @testset "Missing from used variables" begin
-        m1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1|subj)), dat(:sleepstudy))
+        m1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1|subj)), MixedModels.dataset(:sleepstudy))
         m1_missing = fit(MixedModel, @formula(reaction ~ 1 + days + (1|subj)), slp)
         @test nobs(m1) - nobs(m1_missing) == 1
     end

test/pirls.jl

@@ -1,10 +1,7 @@
 using DataFrames, Feather, LinearAlgebra, MixedModels, Test
 
-dat(nm::AbstractString) = Feather.read(joinpath(MixedModels.TestData, string(nm, ".feather")))
-dat(nm::Symbol) = dat(string(nm))
-
 @testset "contra" begin
-    contra = dat(:contra)
+    contra = MixedModels.dataset(:contra)
     contraform = @formula(use ~ 1+age+abs2(age)+urban+livch+(1|urbdist))
     gm0 = fit(MixedModel, contraform, contra, Bernoulli(), fast=true);
     @test gm0.lowerbd == zeros(1)
@@ -35,7 +32,7 @@ dat(nm::Symbol) = dat(string(nm))
 end
 
 @testset "cbpp" begin
-    cbpp = dat(:cbpp)
+    cbpp = MixedModels.dataset(:cbpp)
     gm2 = fit(MixedModel, @formula((incid/hsz) ~ 1 + period + (1|herd)), cbpp, Binomial(), wts=float(cbpp.hsz))
     @test deviance(gm2,true) ≈ 100.09585619892968 atol=0.0001
     @test sum(abs2, gm2.u[1]) ≈ 9.723054788538546 atol=0.0001
@@ -48,7 +45,7 @@ end
 
 @testset "verbagg" begin
     gm3 = fit(MixedModel, @formula(r2 ~ 1+anger+gender+btype+situ+(1|subj)+(1|item)),
-        dat(:verbagg), Bernoulli())
+        MixedModels.dataset(:verbagg), Bernoulli())
     @test deviance(gm3) ≈ 8151.40 rtol=1e-5
     @test lowerbd(gm3) == vcat(fill(-Inf, 6), zeros(2))
     @test fitted(gm3) == predict(gm3)
@@ -59,7 +56,7 @@ end
 
 @testset "grouseticks" begin
     center(v::AbstractVector) = v .- (sum(v) / length(v))
-    grouseticks = dat(:grouseticks)
+    grouseticks = MixedModels.dataset(:grouseticks)
     grouseticks.ch = center(grouseticks.height)
     gm4 = fit(MixedModel, @formula(ticks ~ 1+year+ch+ (1|index) + (1|brood) + (1|location)),
         grouseticks, Poisson(), fast=true)  # fails in pirls! with fast=false

test/pls.jl

@@ -1,13 +1,10 @@
 using BlockArrays, DataFrames, Feather, LinearAlgebra, MixedModels, NamedArrays
 using Random, SparseArrays, Statistics, Tables, Test
 
-dat(nm::AbstractString) = Feather.read(joinpath(MixedModels.TestData, nm * ".feather"))
-dat(nm::Symbol) = dat(string(nm))
-
 const LMM = LinearMixedModel
 
 @testset "Dyestuff" begin
-    ds = dat(:dyestuff)
+    ds = MixedModels.dataset(:dyestuff)
     fm1 = LMM(@formula(yield ~ 1 + (1|batch)), ds)
 
     @test length(fm1.allterms) == 3
@@ -86,7 +83,7 @@ const LMM = LinearMixedModel
 end
 
 @testset "Dyestuff2" begin
-    ds2 = dat(:dyestuff2)
+    ds2 = MixedModels.dataset(:dyestuff2)
     fm = fit(MixedModel, @formula(yield ~ 1 + (1|batch)), ds2)
     @test lowerbd(fm) == zeros(1)
     show(IOBuffer(), fm)
@@ -97,12 +94,12 @@ end
     @test stderror(fm) ≈ [0.6669857396443261]
     @test coef(fm) ≈ [5.6656]
     @test logdet(fm) ≈ 0.0
-    refit!(fm, float(dat(:dyestuff)[!, :yield]))
+    refit!(fm, float(MixedModels.dataset(:dyestuff)[!, :yield]))
     @test objective(fm) ≈ 327.3270598811428 atol=0.001
 end
 
 @testset "penicillin" begin
-    pen = dat(:penicillin)
+    pen = MixedModels.dataset(:penicillin)
     fm = LMM(@formula(diameter ~ 1 + (1 | plate) + (1 | sample)), pen);
     @test size(fm) == (144, 1, 30, 2)
     @test fm.θ == ones(2)
@@ -126,7 +123,7 @@ end
 end
 
 @testset "pastes" begin
-    fm = LMM(@formula(strength ~ (1|sample) + (1|batch)), dat(:pastes))
+    fm = LMM(@formula(strength ~ (1|sample) + (1|batch)), MixedModels.dataset(:pastes))
     @test size(fm) == (60, 1, 40, 2)
     @test fm.θ == ones(2)
     @test lowerbd(fm) == zeros(2)
@@ -146,7 +143,7 @@ end
 end
 
 @testset "InstEval" begin
-    insteval = dat(:insteval)
+    insteval = MixedModels.dataset(:insteval)
     fm1 = LMM(@formula(y ~ 1 + service + (1|s) + (1|d) + (1|dept)), insteval)
     @test size(fm1) == (73421, 2, 4114, 3)
     @test fm1.θ == ones(3)
@@ -170,7 +167,7 @@ end
 end
 
 @testset "sleep" begin
-    slp = dat(:sleepstudy)
+    slp = MixedModels.dataset(:sleepstudy)
     fm = LinearMixedModel(@formula(reaction ~ 1 + days + (1+days|subj)), slp);
     @test lowerbd(fm) == [0.0, -Inf, 0.0]
     A11 = fm.A[Block(1,1)]
@@ -294,7 +291,7 @@ end
 end
 
 @testset "d3" begin
-    fm = updateL!(LMM(@formula(y ~ 1 + u + (1+u|g) + (1+u|h) + (1+u|i)), dat(:d3)));
+    fm = updateL!(LMM(@formula(y ~ 1 + u + (1+u|g) + (1+u|h) + (1+u|i)), MixedModels.dataset(:d3)));
     @test pwrss(fm) ≈ 5.1261847180180885e6 rtol = 1e-6
     @test objective(fm) ≈ 901641.2930413672 rtol = 1e-6
     fit!(fm)
@@ -305,7 +302,7 @@ end
 end
 
 @testset "simulate!" begin
-    ds = dat(:dyestuff)
+    ds = MixedModels.dataset(:dyestuff)
     fm = fit(MixedModel, @formula(yield ~ 1 + (1|batch)), ds)
     refit!(simulate!(Random.MersenneTwister(1234321), fm))
     @test deviance(fm) ≈ 339.0218639362958 atol=0.001

test/utilities.jl

@@ -21,3 +21,10 @@ end
 
 	@test all(sort!(single_thread) .== sort!(multi_thread))
 end
+
+@testset "datasets" begin
+	@test isa(MixedModels.datasets(), Vector{String})
+	@test size(MixedModels.dataset(:dyestuff)) == (30, 2)
+	@test size(MixedModels.dataset("dyestuff")) == (30, 2)
+	@test_throws ArgumentError MixedModels.dataset(:foo)
+end