Fix bug in ranef!, change isempty(offset₀) branch in setβθ! (#63)

CI failures are on nightlies only and not related to this change.

src/PIRLS.jl

@@ -211,8 +211,9 @@ function setβθ!{T}(m::GeneralizedLinearMixedModel{T}, v::Vector{T})
     β, lm, offset₀, X, rr = m.β, m.LMM, m.offset₀, m.X, m.resp
     lb = length(β)
     copy!(β, view(v, 1 : lb))
+    isempty(offset₀) ? BLAS.gemv!('N', one(T), X, β, zero(T), rr.offset) :
+        BLAS.gemv!('N', one(T), X, β, one(T), copy!(rr.offset, offset₀))
     setθ!(m.LMM, copy!(m.θ, view(v, (lb + 1) : length(v))))
-    BLAS.gemv!('N', one(T), X, β, one(T), isempty(offset₀) ? fill!(rr.offset, 0) : copy!(rr.offset, offset₀))
     m
 end
 

src/mixedmodel.jl

@@ -105,23 +105,21 @@ function grplevels(m::MixedModel)
 end
 
 """
-    ranef!{T}(v::Vector{Matrix{T}}, m::MixedModel{T}, uscale::Bool = false)
+    ranef!{T}(v::Vector{Matrix{T}}, m::MixedModel{T}, β, uscale::Bool)
 
 Overwrites `v` with the conditional modes of the random effects for `m`.
 
 If `uscale` is `true` the random effects are on the spherical (i.e. `u`) scale, otherwise on the
 original scale
 """
-function ranef!(v::Vector, m::MixedModel, uscale)
+function ranef!{T}(v::Vector, m::LinearMixedModel{T}, β::AbstractArray{T}, uscale::Bool)
     R, Λ = m.R, m.Λ
     k = length(Λ)        # number of random-effects terms
     for j in 1:k
         copy!(v[j], R[j, end])
     end
-    rβ = R[k + 1, end]
-    if !isempty(rβ)      #  in the pirls! function for GLMMs want to skip this
-        β = vec(feR(m) \ rβ)
-        kp1 = k + 1
+    kp1 = k + 1
+    if !isempty(β)       #  in the pirls! function for GLMMs want to skip this
         for j in 1:k     # subtract the fixed-effects contribution
             BLAS.gemv!('N', -1.0, R[j, kp1], β, 1.0, vec(v[j]))
         end
@@ -130,9 +128,13 @@ function ranef!(v::Vector, m::MixedModel, uscale)
         Rjj = R[j, j]
         uj = vec(v[j])
         LinAlg.A_ldiv_B!(isa(Rjj, Diagonal) ? Rjj : UpperTriangular(Rjj), uj)
-        for i in 1:j - 1
-            ui = vec(v[i])
-            ui -= R[i, j] * uj
+        for i in 1:(j - 1)
+            Rij = R[i, j]
+            if isa(Rij, StridedMatrix{T})
+                BLAS.gemv!('N', -one(T), R[i, j], uj, one(T), vec(v[i]))
+            else
+                A_mul_B!(-one(T), R[i, j], uj, one(T), vec(v[i]))
+            end
         end
     end
     if !uscale
@@ -143,6 +145,10 @@ function ranef!(v::Vector, m::MixedModel, uscale)
     v
 end
 
+function ranef!(v::Vector, m::LinearMixedModel, uscale::Bool)
+    ranef!(v, m, feR(m) \ vec(m.R[end-1, end]), uscale)
+end
+
 """
     ranef{T}(m::MixedModel{T}, uscale=false)
 

test/pls.jl

@@ -186,10 +186,10 @@ end
     @test isapprox(objective(fm1), 237585.5534151694, atol = 1e-3)
     ftd1 = fitted(fm1);
     @test size(ftd1) == (73421,)
-    @test isapprox(ftd1[1], 3.1783062136367723, atol = 1e-5)
+    @test isapprox(ftd1[1], 3.19745, atol = 1e-4)
     resid1 = residuals(fm1);
     @test size(resid1) == (73421,)
-    @test isapprox(resid1[1], 1.8216937863632277, atol = 1e-5)
+    @test isapprox(resid1[1], 1.8025537831086234, atol = 1e-5)
 
     fm2 = fit!(lmm(y ~ 1 + service + (1 | s) + (1 | d) + (1 | dept), InstEval));
     @test size(fm2) == (73421,2,4114,3)