Merge pull request #1 from JuliaStats/master

Tighten code, increase test coverage (#301)

src/MixedModels.jl

@@ -56,7 +56,6 @@ export @formula,
        RandomEffectsTerm,
        ReMat,
        SqrtLink,
-       TestData,
        UniformBlockDiagonal,
        VarCorr,
 
@@ -84,6 +83,7 @@ export @formula,
        GHnorm,
        issingular,
        leverage,
+       logdet,
        loglikelihood,
        lowerbd,
        nobs,
@@ -93,6 +93,7 @@ export @formula,
        predict,
        pwrss,
        ranef,
+       rank,
        refit!,
        replicate,
        residuals,

src/arraytypes.jl

@@ -15,6 +15,11 @@ function UniformBlockDiagonal(dat::Array{T,3}) where {T}
     )
 end
 
+function Base.axes(A::UniformBlockDiagonal)
+    m, n, l = size(A.data)
+    (Base.OneTo(m * l), Base.OneTo(n * l))
+end
+
 function Base.copyto!(dest::UniformBlockDiagonal{T}, src::UniformBlockDiagonal{T}) where {T}
     sdat = src.data
     ddat = dest.data
@@ -42,10 +47,9 @@ function Base.copyto!(dest::Matrix{T}, src::UniformBlockDiagonal{T}) where {T}
 end
 
 function Base.getindex(A::UniformBlockDiagonal{T}, i::Int, j::Int) where {T}
+    @boundscheck checkbounds(A, i, j)
     Ad = A.data
     m, n, l = size(Ad)
-    (0 < i ≤ l * m && 0 < j ≤ l * n) ||
-    throw(IndexError("attempt to access $(l*m) × $(l*n) array at index [$i, $j]"))
     iblk, ioffset = divrem(i - 1, m)
     jblk, joffset = divrem(j - 1, n)
     iblk == jblk ? Ad[ioffset+1, joffset+1, iblk+1] : zero(T)

src/blockdescription.jl

@@ -43,6 +43,8 @@ shorttype(::Diagonal,::Diagonal) = "Diagonal"
 shorttype(::Diagonal,::Matrix) = "Diag/Dense"
 shorttype(::Matrix,::Matrix) = "Dense"
 shorttype(::SparseMatrixCSC,::SparseMatrixCSC) = "Sparse"
+shorttype(::SparseMatrixCSC,::Matrix) = "Sparse/Dense"
+
 
 function Base.show(io::IO, ::MIME"text/plain", b::BlockDescription)
     rowwidth = max(maximum(ndigits, b.blkrows) + 1, 5)

src/bootstrap.jl

@@ -180,8 +180,6 @@ end
 
 issingular(bsamp::MixedModelBootstrap) = map(θ -> any(θ .≈ bsamp.lowerbd), bsamp.θ)
 
-ppoints(n::Integer) = inv(2n):inv(n):1
-
 function Base.propertynames(bsamp::MixedModelBootstrap)
     [:allpars, :objective, :σ, :β, :θ, :σs, :λ, :inds, :lowerbd, :bstr, :fcnames]
 end

src/likelihoodratiotest.jl

@@ -91,14 +91,6 @@ function _likelihoodratiotest(m::Vararg{T}) where T <: MixedModel
     )
 end
 
-function _array_union_nothing(arr::Array{T}) where T
-    Array{Union{T,Nothing}}(arr)
-end
-
-function _prepend_0(arr::Array{T}) where T
-    pushfirst!(copy(arr), -zero(T))
-end
-
 function Base.show(io::IO, lrt::LikelihoodRatioTest; digits=2)
     println(io, "Model Formulae")
 

src/linalg.jl

@@ -25,22 +25,6 @@ function LinearAlgebra.mul!(
     C
 end
 
-LinearAlgebra.mul!(
-    C::Matrix{T},
-    A::BlockedSparse{T},
-    adjB::Adjoint{T,<:Matrix{T}},
-    α::Number,
-    β::Number,
-) where {T} = mul!(C, A.cscmat, adjB, α, β)
-
-LinearAlgebra.mul!(
-    C::BlockedSparse{T},
-    A::BlockedSparse{T},
-    adjB::Adjoint{T,<:BlockedSparse{T}},
-    α::Number,
-    β::Number,
-) where {T} = mul!(C.cscmat, A.cscmat, adjB.parent.cscmat', α, β)
-
 LinearAlgebra.mul!(
     C::StridedVecOrMat{T},
     A::StridedVecOrMat{T},

src/linalg/rankUpdate.jl

@@ -127,8 +127,10 @@ function rankUpdate!(
     C
 end
 
+#=
 rankUpdate!(C::HermOrSym{T,Matrix{T}}, A::BlockedSparse{T}, α = true) where {T} =
     rankUpdate!(C, A.cscmat, α)
+=#
 
 function rankUpdate!(
     C::Diagonal{T,S},

src/linearmixedmodel.jl

@@ -303,17 +303,15 @@ objective and the parameters are printed on stdout at each function evaluation.
 function fit!(m::LinearMixedModel{T}; verbose::Bool = false, REML::Bool = false) where {T}
     optsum = m.optsum
     opt = Opt(optsum)
-    feval = 0
     optsum.REML = REML
     function obj(x, g)
-        isempty(g) || error("gradient not defined")
-        feval += 1
+        isempty(g) || throw(ArgumentError("g should be empty for this objective"))
         val = objective(updateL!(setθ!(m, x)))
-        feval == 1 && (optsum.finitial = val)
-        verbose && println("f_", feval, ": ", round(val, digits = 5), " ", x)
+        verbose && println(round(val, digits = 5), " ", x)
         val
     end
     NLopt.min_objective!(opt, obj)
+    optsum.finitial = obj(optsum.initial, T[])
     fmin, xmin, ret = NLopt.optimize!(opt, copyto!(optsum.final, optsum.initial))
     ## check if small non-negative parameter values can be set to zero
     xmin_ = copy(xmin)
@@ -332,12 +330,12 @@ function fit!(m::LinearMixedModel{T}; verbose::Bool = false, REML::Bool = false)
     ## ensure that the parameter values saved in m are xmin
     updateL!(setθ!(m, xmin))
 
-    optsum.feval = feval
+    optsum.feval = opt.numevals
     optsum.final = xmin
     optsum.fmin = fmin
     optsum.returnvalue = ret
     ret == :ROUNDOFF_LIMITED && @warn("NLopt was roundoff limited")
-    if ret ∈ [:FAILURE, :INVALID_ARGS, :OUT_OF_MEMORY, :FORCED_STOP, :MAXFEVAL_REACHED]
+    if ret ∈ [:FAILURE, :INVALID_ARGS, :OUT_OF_MEMORY, :FORCED_STOP, :MAXEVAL_REACHED]
         @warn("NLopt optimization failure: $ret")
     end
     m
@@ -474,8 +472,10 @@ end
     issingular(m::LinearMixedModel, θ=m.θ)
 
 Test whether the model `m` is singular if the parameter vector is `θ`.
+
+Equality comparisons are used b/c small non-negative θ values are replaced by 0 in `fit!`.
 """
-issingular(m::LinearMixedModel, θ=m.θ) = any(isapprox.(lowerbd(m), θ))
+issingular(m::LinearMixedModel, θ=m.θ) = any(lowerbd(m) .== θ)
 
 function StatsBase.leverage(m::LinearMixedModel{T}) where {T}
     # This can be done more efficiently but reusing existing tools is easier.
@@ -506,11 +506,11 @@ end
 lowerbd(m::LinearMixedModel) = m.optsum.lowerbd
 
 function StatsBase.modelmatrix(m::LinearMixedModel)
-    fetrm = first(m.feterms)
-    if fetrm.rank == size(fetrm, 2)
-        fetrm.x
+    fe = fetrm(m)
+    if fe.rank == size(fe, 2)
+        fe.x
     else
-        fetrm.x[:, invperm(fetrm.piv)]
+        fe.x[:, invperm(fe.piv)]
     end
 end
 

test/Project.toml

@@ -3,7 +3,6 @@ BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 NamedArrays = "86f7a689-2022-50b4-a561-43c23ac3c673"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

test/pirls.jl

@@ -1,25 +1,31 @@
-using DataFrames
-using LinearAlgebra
 using MixedModels
 using Test
 
+using MixedModels: dataset
+
+const fms = Dict(
+    :cbpp => [@formula((incid/hsz) ~ 1 + period + (1|herd))],
+    :contra => [@formula(use ~ 1+age+abs2(age)+urban+livch+(1|urbdist))],
+    :grouseticks => [@formula(ticks ~ 1+year+ch+ (1|index) + (1|brood) + (1|location))],
+    :verbagg => [@formula(r2 ~ 1+anger+gender+btype+situ+(1|subj)+(1|item))],
+)
+
 @testset "contra" begin
-    contra = MixedModels.dataset(:contra)
-    contraform = @formula(use ~ 1+age+abs2(age)+urban+livch+(1|urbdist))
-    gm0 = fit(MixedModel, contraform, contra, Bernoulli(), fast=true);
+    contra = dataset(:contra)
+    gm0 = fit(MixedModel, only(fms[:contra]), contra, Bernoulli(), fast=true);
     @test gm0.lowerbd == zeros(1)
     @test isapprox(gm0.θ, [0.5720734451352923], atol=0.001)
-    @test isapprox(deviance(gm0,true), 2361.657188518064, atol=0.001)
-    gm1 = fit(MixedModel, contraform, contra, Bernoulli());
+    @test isapprox(deviance(gm0), 2361.657188518064, atol=0.001)
+    gm1 = fit(MixedModel, only(fms[:contra]), contra, Bernoulli());
     @test isapprox(gm1.θ, [0.573054], atol=0.005)
     @test lowerbd(gm1) == vcat(fill(-Inf, 7), 0.)
-    @test isapprox(deviance(gm1,true), 2361.54575, rtol=0.00001)
+    @test isapprox(deviance(gm1), 2361.54575, rtol=0.00001)
     @test isapprox(loglikelihood(gm1), -1180.77288, rtol=0.00001)
     @test dof(gm0) == length(gm0.β) + length(gm0.θ)
     @test nobs(gm0) == 1934
     fit!(gm0, fast=true, nAGQ=7)
     @test isapprox(deviance(gm0), 2360.9838, atol=0.001)
-    gm1 = fit(MixedModel, contraform, contra, Bernoulli(), nAGQ=7)
+    gm1 = fit(MixedModel, only(fms[:contra]), contra, Bernoulli(), nAGQ=7)
     @test isapprox(deviance(gm1), 2360.8760, atol=0.001)
     @test gm1.β == gm1.beta
     @test gm1.θ == gm1.theta
@@ -31,7 +37,6 @@ using Test
         @test length(MixedModels.rePCA(gm0)) == 1
         @test length(gm0.rePCA) == 1
     end
-    # gm0.βθ = vcat(gm0.β, gm0.theta)
     # the next three values are not well defined in the optimization
     #@test isapprox(logdet(gm1), 75.7217, atol=0.1)
     #@test isapprox(sum(abs2, gm1.u[1]), 48.4747, atol=0.1)
@@ -41,8 +46,8 @@ using Test
 end
 
 @testset "cbpp" begin
-    cbpp = MixedModels.dataset(:cbpp)
-    gm2 = fit(MixedModel, @formula((incid/hsz) ~ 1 + period + (1|herd)), cbpp, Binomial(), wts=float(cbpp.hsz))
+    cbpp = dataset(:cbpp)
+    gm2 = fit(MixedModel, only(fms[:cbpp]), 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
     @test logdet(gm2) ≈ 16.90105378801136 atol=0.0001
@@ -53,22 +58,20 @@ end
 end
 
 @testset "verbagg" begin
-    gm3 = fit(MixedModel, @formula(r2 ~ 1+anger+gender+btype+situ+(1|subj)+(1|item)),
-        MixedModels.dataset(:verbagg), Bernoulli())
+    gm3 = fit(MixedModel, only(fms[:verbagg]), 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)
     # these two values are not well defined at the optimum
-    @test isapprox(sum(x -> sum(abs2, x), gm3.u), 273.29266717430795, rtol=1e-3)
-    @test sum(gm3.resp.devresid) ≈ 7156.547357801238 rtol=1e-4
+    @test isapprox(sum(x -> sum(abs2, x), gm3.u), 273.29646346940785, rtol=1e-3)
+    @test sum(gm3.resp.devresid) ≈ 7156.550941446312 rtol=1e-4
 end
 
 @testset "grouseticks" begin
     center(v::AbstractVector) = v .- (sum(v) / length(v))
-    grouseticks = MixedModels.dataset(:grouseticks)
+    grouseticks = 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
+    gm4 = fit(MixedModel, only(fms[:grouseticks]), grouseticks, Poisson(), fast=true)  # fails in pirls! with fast=false
     @test isapprox(deviance(gm4), 851.4046, atol=0.001)
     # these two values are not well defined at the optimum
     #@test isapprox(sum(x -> sum(abs2, x), gm4.u), 196.8695297987013, atol=0.1)