diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 3a176fa06..ecd73dcf3 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -4,8 +4,9 @@ using MixedModels: dataset
 const SUITE = BenchmarkGroup()
 
 const global contrasts = Dict(
-    :mrk17_exp1 => merge(Dict(n => HelmertCoding() for n in (:F, :P, :Q, :lQ, :lT)),
-     Dict(n => Grouping() for n in (:item, :subj))),
+    :mrk17_exp1 => merge(Dict(n => HelmertCoding() for n in (:F, :P, :Q, :lQ, :lT)), Dict(n => Grouping() for n in (:item, :subj))),
+    :kb07 => merge(Dict(n => HelmertCoding() for n in (:spkr, :prec, :load)), Dict(n => Grouping() for n in (:item, :subj))),
+    :insteval => Dict(n => Grouping() for n in (:s, :d)),
 )
 const global fms = Dict(
     :dyestuff => [
diff --git a/docs/src/GaussHermite.md b/docs/src/GaussHermite.md
index 50c867e2d..ee2bc863c 100644
--- a/docs/src/GaussHermite.md
+++ b/docs/src/GaussHermite.md
@@ -33,8 +33,10 @@ For a `k`th order normalized Gauss-Hermite rule the tridiagonal matrix has zeros
 using DataFrames, LinearAlgebra, Gadfly
 sym3 = SymTridiagonal(zeros(3), sqrt.(1:2))
 ev = eigen(sym3);
-show(ev.values)
-show(abs2.(ev.vectors[1,:]))
+ev.values
+```
+```@example Main
+abs2.(ev.vectors[1,:])
 ```
 As a function of `k` this can be written as
 ```@example Main
diff --git a/docs/src/bootstrap.md b/docs/src/bootstrap.md
index 65427de2a..6c0cd3a94 100644
--- a/docs/src/bootstrap.md
+++ b/docs/src/bootstrap.md
@@ -41,7 +41,7 @@ To bootstrap the model parameters, first initialize a random number generator th
 
 ```@example Main
 const rng = MersenneTwister(1234321);
-samp = parametricbootstrap(rng, 10_000, m1, use_threads=true);
+samp = parametricbootstrap(rng, 10_000, m1);
 df = DataFrame(samp.allpars);
 first(df, 10)
 ```
diff --git a/docs/src/constructors.md b/docs/src/constructors.md
index b52403f03..76440afab 100644
--- a/docs/src/constructors.md
+++ b/docs/src/constructors.md
@@ -84,9 +84,9 @@ fm1reml = fit(MixedModel, fm, dyestuff, REML=true)
 ```
 ```@setup Main
 @testset "fm1reml" begin
-    @test deviance(fm1reml) ≈ 319.6542768422538
-    @test varest(fm1reml) ≈ 2451.2499
-    @test VarCorr(fm1reml).σρ.batch.σ[1] ≈ 42.000602
+    @test deviance(fm1reml) ≈ 319.6542768422538 atol=0.001
+    @test varest(fm1reml) ≈ 2451.2499 atol=0.001
+    @test VarCorr(fm1reml).σρ.batch.σ[1] ≈ 42.000602 atol=0.001
 end
 ```
 ### Float-point type in the model
diff --git a/docs/src/optimization.md b/docs/src/optimization.md
index ea58bbf13..6c01e3811 100644
--- a/docs/src/optimization.md
+++ b/docs/src/optimization.md
@@ -120,7 +120,7 @@ identity matrix where the multiple is
 t1.λ
 ```
 
-Because there is only one random-effects term in the model, the matrix $\bf Z$ is the indicators matrix shown as the result of `Matrix(t1)`, but stored in a special sparse format.
+Because there is only one random-effects term in the model, the matrix $\bf Z$ is the indicators matrix shown as the result of `Int.(t1)`, but stored in a special sparse format.
 Furthermore, there is only one block in $\Lambda_\theta$.
 
 
@@ -153,19 +153,23 @@ Int.(t31)
 For this model the matrix $\bf Z$ is the same as that of model `fm2` but the diagonal blocks of $\Lambda_\theta$ are themselves diagonal.
 ```@example Main
 t31.λ
+```
+```@example Main
 MixedModels.getθ(t31)
 ```
-Random-effects terms with distinct grouping factors generate distinct elements of the `allterms` field of the `LinearMixedModel` object.
+Random-effects terms with distinct grouping factors generate distinct elements of the `reterms` field of the `LinearMixedModel` object.
 Multiple `ReMat` objects are sorted by decreasing numbers of random effects.
 ```@example Main
 penicillin = MixedModels.dataset(:penicillin)
 fm4 = fit(MixedModel,
-    @formula(diameter ~ 1 + (1|plate) + (1|sample)),
+    @formula(diameter ~ 1 + (1|sample) + (1|plate)),
     penicillin)
 Int.(first(fm4.reterms))
+```
+```@example Main
 Int.(last(fm4.reterms))
 ```
-Note that the first `ReMat` in `fm4.terms` corresponds to grouping factor `G` even though the term `(1|G)` occurs in the formula after `(1|H)`.
+Note that the first `ReMat` in `fm4.reterms` corresponds to grouping factor `plate` even though the term `(1|plate)` occurs in the formula after `(1|sample)`.
 
 ### Progress of the optimization
 
@@ -183,6 +187,7 @@ fm2.optsum
 ## A blocked Cholesky factor
 
 A `LinearMixedModel` object contains two blocked matrices; a symmetric matrix `A` (only the lower triangle is stored) and a lower-triangular `L` which is the lower Cholesky factor of the updated and inflated `A`.
+In versions 4.0.0 and later of `MixedModels` only the blocks in the lower triangle are stored in `A` and `L`, as a `Vector{AbstractMatrix{T}}`
 ```@docs
 BlockDescription
 ```
@@ -191,6 +196,8 @@ shows the structure of the blocks
 BlockDescription(fm2)
 ```
 
+Another change in v4.0.0 and later is that the last row of blocks is constructed from `m.Xymat` which contains the full-rank model matrix `X` with the response `y` concatenated on the right.
+
 The operation of installing a new value of the variance parameters, `θ`, and updating `L`
 ```@docs
 setθ!
@@ -246,7 +253,7 @@ In a [*generalized linear model*](https://en.wikipedia.org/wiki/Generalized_line
 The scalar distributions of individual responses differ only in their means, which are determined by a *linear predictor* expression $\eta=\bf X\beta$, where, as before, $\bf X$ is a model matrix derived from the values of covariates and $\beta$ is a vector of coefficients.
 
 The unconstrained components of $\eta$ are mapped to the, possiby constrained, components of the mean response, $\mu$, via a scalar function, $g^{-1}$, applied to each component of $\eta$.
-For historical reasons, the inverse of this function, taking components of $\mu$ to the corresponding component of $\eta$ is called the *link function* and more frequently used map from $\eta$ to $\mu$ is the *inverse link*.
+For historical reasons, the inverse of this function, taking components of $\mu$ to the corresponding component of $\eta$ is called the *link function* and the more frequently used map from $\eta$ to $\mu$ is the *inverse link*.
 
 A *generalized linear mixed-effects model* (GLMM) is defined, for the purposes of this package, by
 ```math
@@ -298,7 +305,7 @@ In a call to the `pirls!` function the first argument is a `GeneralizedLinearMix
 The second and third arguments are optional logical values indicating if $\beta$ is to be varied and if verbose output is to be printed.
 
 ```@example Main
-pirls!(mdl, true, true)
+pirls!(mdl, true, false);
 ```
 
 ```@example Main
@@ -344,9 +351,3 @@ mdl1 = @btime fit(MixedModel, vaform, verbagg, Bernoulli())
 
 This fit provided slightly better results (Laplace approximation to the deviance of 8151.400 versus 8151.583) but took 6 times as long.
 That is not terribly important when the times involved are a few seconds but can be important when the fit requires many hours or days of computing time.
-
-The comparison of the slow and fast fit is available in the optimization summary after the slow fit.
-
-```@example Main
-mdl1.LMM.optsum
-```
diff --git a/src/MixedModels.jl b/src/MixedModels.jl
index 2f36a6f38..b44b103e8 100644
--- a/src/MixedModels.jl
+++ b/src/MixedModels.jl
@@ -155,7 +155,7 @@ include("pca.jl")
 include("datasets.jl")
 include("arraytypes.jl")
 include("varcorr.jl")
-include("femat.jl")
+include("Xymat.jl")
 include("remat.jl")
 include("optsummary.jl")
 include("schema.jl")
diff --git a/src/Xymat.jl b/src/Xymat.jl
new file mode 100644
index 000000000..e198ae819
--- /dev/null
+++ b/src/Xymat.jl
@@ -0,0 +1,141 @@
+"""
+    FeTerm{T,S}
+
+Term with an explicit, constant matrix representation
+
+Typically, an `FeTerm` represents the model matrix for the fixed effects. 
+
+!!! note
+    `FeTerm` is not the same as [`FeMat`](@ref)!
+
+# Fields
+* `x`: full model matrix
+* `piv`: pivot `Vector{Int}` for moving linearly dependent columns to the right
+* `rank`: computational rank of `x`
+* `cnames`: vector of column names
+"""
+struct FeTerm{T,S<:AbstractMatrix}
+    x::S
+    piv::Vector{Int}
+    rank::Int
+    cnames::Vector{String}
+end
+
+"""
+    FeTerm(X::AbstractMatrix, cnms)
+    
+Convenience constructor for [`FeTerm`](@ref) that computes the rank and pivot with unit weights.
+
+See the vignette "[Rank deficiency in mixed-effects models](@ref)" for more information on the 
+computation of the rank and pivot.
+"""
+function FeTerm(X::AbstractMatrix{T}, cnms) where {T}
+    if iszero(size(X, 2))
+        return FeTerm{T,typeof(X)}(X, Int[], 0, cnms)
+    end
+    st = statsqr(X)
+    pivot = st.p
+    rank = findfirst(<=(0), diff(st.p))
+    rank = isnothing(rank) ? length(pivot) : rank
+    Xp = pivot == collect(1:size(X, 2)) ? X : X[:, pivot]
+        # single-column rank deficiency is the result of a constant column vector
+        # this generally happens when constructing a dummy response, so we don't
+        # warn.
+    if rank < length(pivot) && size(X,2) > 1
+        @warn "Fixed-effects matrix is rank deficient"
+    end
+    FeTerm{T,typeof(X)}(Xp, pivot, rank, cnms[pivot])
+end
+
+"""
+    FeTerm(X::SparseMatrixCSC, cnms)
+    
+Convenience constructor for a sparse [`FeTerm`](@ref) assuming full rank, identity pivot and unit weights.
+
+Note: automatic rank deficiency handling may be added to this method in the future, as discused in
+the vignette "[Rank deficiency in mixed-effects models](@ref)" for general `FeTerm`.
+"""
+function FeTerm(X::SparseMatrixCSC, cnms::AbstractVector{String})
+    @debug "Full rank is assumed for sparse fixed-effect matrices."
+    rank = size(X,2)
+    FeTerm{eltype(X),typeof(X)}(X, collect(1:rank), rank, collect(cnms))
+end
+
+Base.copyto!(A::FeTerm{T}, src::AbstractVecOrMat{T}) where {T} = copyto!(A.x, src)
+
+Base.eltype(::FeTerm{T}) where {T} = T
+
+function fullrankx(A::FeTerm)
+    x, rnk = A.x, A.rank
+    rnk == size(x, 2) ? x : view(x, :, 1:rnk)  # this handles the zero-columns case
+end
+
+LinearAlgebra.rank(A::FeTerm) = A.rank
+
+"""
+    isfullrank(A::FeTerm)
+
+Does `A` have full column rank?
+"""
+isfullrank(A::FeTerm) = A.rank == length(A.piv)
+
+"""
+    FeMat{T,S}
+
+A matrix and a (possibly) weighted copy of itself.
+
+
+Typically, an `FeMat` represents the fixed-effects model matrix with the response (`y`) concatenated as a final column.
+
+!!! note
+    `FeMat` is not the same as [`FeTerm`](@ref).
+
+# Fields
+- `xy`: original matrix, called `xy` b/c in practice this is `hcat(fullrank(X), y)`
+- `wtxy`: (possibly) weighted copy of `xy` (shares storage with `xy` until weights are applied)
+
+Upon construction the `xy` and `wtxy` fields refer to the same matrix
+"""
+mutable struct FeMat{T,S<:AbstractMatrix} <: AbstractMatrix{T}
+    xy::S
+    wtxy::S
+end
+
+function FeMat(A::FeTerm{T}, y::AbstractVector{T}) where {T}
+    xy = hcat(fullrankx(A), y)
+    FeMat{T,typeof(xy)}(xy, xy)
+end
+
+Base.adjoint(A::FeMat) = Adjoint(A)
+
+Base.eltype(::FeMat{T}) where {T} = T
+
+Base.getindex(A::FeMat, i, j) = getindex(A.xy, i, j)
+
+Base.length(A::FeMat) = length(A.xy)
+
+function *(adjA::Adjoint{T,<:FeMat{T}}, B::FeMat{T}) where {T}
+    adjoint(adjA.parent.wtxy) * B.wtxy
+end
+
+function LinearAlgebra.mul!(R::StridedVecOrMat{T}, A::FeMat{T}, B::StridedVecOrMat{T}) where {T}
+    mul!(R, A.wtxy, B)
+end
+
+function LinearAlgebra.mul!(C, adjA::Adjoint{T,<:FeMat{T}}, B::FeMat{T}) where {T}
+    mul!(C, adjoint(adjA.parent.wtxy), B.wtxy)
+end
+
+function reweight!(A::FeMat{T}, sqrtwts::Vector{T}) where {T}
+    if !isempty(sqrtwts)
+        if A.xy === A.wtxy
+            A.wtxy = similar(A.xy)
+        end
+        mul!(A.wtxy, Diagonal(sqrtwts), A.xy)
+    end
+    A
+end
+
+Base.size(A::FeMat) = size(A.xy)
+
+Base.size(A::FeMat, i::Integer) = size(A.xy, i)
diff --git a/src/arraytypes.jl b/src/arraytypes.jl
index 75d3c02b5..cc03303cc 100644
--- a/src/arraytypes.jl
+++ b/src/arraytypes.jl
@@ -111,5 +111,5 @@ function LinearAlgebra.mul!(
     α,
     β
     ) where {T,P,Ti}
-    mul!(C.cscmat, A, adjB.parent.cscmat', α, β)
+    mul!(C.cscmat, A, adjoint(adjB.parent.cscmat), α, β)
 end
diff --git a/src/blockdescription.jl b/src/blockdescription.jl
index 5b1b50530..e65f15604 100644
--- a/src/blockdescription.jl
+++ b/src/blockdescription.jl
@@ -61,6 +61,3 @@ function Base.show(io::IO, ::MIME"text/plain", b::BlockDescription)
 end
 
 Base.show(io::IO, b::BlockDescription) = show(io, MIME"text/plain"(), b)
-
-@deprecate describeblocks(io, m) show(io, BlockDescription(m))
-@deprecate describeblocks(m) show(stdout, BlockDescription(m))
diff --git a/src/femat.jl b/src/femat.jl
deleted file mode 100644
index d118b857c..000000000
--- a/src/femat.jl
+++ /dev/null
@@ -1,111 +0,0 @@
-"""
-    FeMat{T,S}
-
-Term with an explicit, constant matrix representation
-
-# Fields
-* `x`: matrix
-* `wtx`: weighted matrix
-* `piv`: pivot `Vector{Int}` for moving linearly dependent columns to the right
-* `rank`: computational rank of `x`
-* `cnames`: vector of column names
-
-"""
-mutable struct FeMat{T,S<:AbstractMatrix}
-    x::S
-    wtx::S
-    piv::Vector{Int}
-    rank::Int
-    cnames::Vector{String}
-end
-
-"""
-    FeMat(X::AbstractMatrix, cnms)
-    
-Convenience constructor for [`FeMat`](@ref) that computes the rank and pivot with unit weights.
-
-See the vignette "[Rank deficiency in mixed-effects models](@ref)" for more information on the 
-computation of the rank and pivot.
-"""
-function FeMat(X::AbstractMatrix, cnms)
-    T = eltype(X)
-    if size(X,2) > 0
-        st = statsqr(X)
-        pivot = st.p
-        rank = findfirst(<=(0), diff(st.p))
-        rank = isnothing(rank) ? length(pivot) : rank
-        Xp = pivot == collect(1:size(X, 2)) ? X : X[:, pivot]
-        # single-column rank deficiency is the result of a constant column vector
-        # this generally happens when constructing a dummy response, so we don't
-        # warn.
-        if rank < length(pivot) && size(X,2) > 1
-            @warn "Fixed-effects matrix is rank deficient"
-        end
-    else
-        # although it doesn't take long for an empty matrix,
-        # we can still skip the factorization step, which gets the rank
-        # wrong anyway
-        pivot = Int[]
-        Xp = X
-        rank = 0
-    end
-    FeMat{T,typeof(X)}(Xp, Xp, pivot, rank, cnms[pivot])
-end
-
-"""
-    FeMat(X::SparseMatrixCSC, cnms)
-    
-Convenience constructor for a sparse [`FeMat`](@ref) assuming full rank, identity pivot and unit weights.
-
-Note: automatic rank deficiency handling may be added to this method in the future, as discused in
-the vignette "[Rank deficiency in mixed-effects models](@ref)" for general `FeMat`.
-"""
-function FeMat(X::SparseMatrixCSC, cnms)
-    @debug "Full rank is assumed for sparse fixed-effect matrices."
-    rank = size(X,2)
-    FeMat{eltype(X),typeof(X)}(X, X, 1:rank, rank, cnms)
-end
-
-function reweight!(A::FeMat{T}, sqrtwts::Vector{T}) where {T}
-    if !isempty(sqrtwts)
-        if A.x === A.wtx
-            A.wtx = similar(A.x)
-        end
-        mul!(A.wtx, Diagonal(sqrtwts), A.x)
-    end
-    A
-end
-
-Base.adjoint(A::FeMat) = Adjoint(A)
-
-Base.eltype(A::FeMat{T}) where {T} = T
-
-fullrankwtx(A::FeMat) = rank(A) == size(A, 2) ? A.wtx : A.wtx[:, 1:rank(A)]
-
-Base.length(A::FeMat) = length(A.wtx)
-
-LinearAlgebra.rank(A::FeMat) = A.rank
-
-Base.size(A::FeMat) = size(A.wtx)
-
-Base.size(A::Adjoint{T,<:FeMat{T}}) where {T} = reverse(size(A.parent))
-
-Base.size(A::FeMat, i) = size(A.wtx, i)
-
-Base.copyto!(A::FeMat{T}, src::AbstractVecOrMat{T}) where {T} = copyto!(A.x, src)
-
-*(adjA::Adjoint{T,<:FeMat{T}}, B::FeMat{T}) where {T} =
-    fullrankwtx(adjA.parent)' * fullrankwtx(B)
-
-LinearAlgebra.mul!(R::StridedVecOrMat{T}, A::FeMat{T}, B::StridedVecOrMat{T}) where {T} =
-    mul!(R, A.x, B)
-
-LinearAlgebra.mul!(C, adjA::Adjoint{T,<:FeMat{T}}, B::FeMat{T}) where {T} =
-    mul!(C, fullrankwtx(adjA.parent)', fullrankwtx(B))
-
-"""
-    isfullrank(A::FeMat)
-
-Does `A` have full column rank?
-"""
-isfullrank(A::FeMat) = A.rank == length(A.piv)
diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl
index 0815ad2e5..2540eeba4 100644
--- a/src/generalizedlinearmixedmodel.jl
+++ b/src/generalizedlinearmixedmodel.jl
@@ -54,7 +54,7 @@ struct GeneralizedLinearMixedModel{T<:AbstractFloat} <: MixedModel{T}
 end
 
 function StatsBase.coef(m::GeneralizedLinearMixedModel{T}) where {T}
-    piv = first(m.LMM.feterms).piv
+    piv = m.LMM.feterm.piv
     invpermute!(copyto!(fill(T(-0.0), length(piv)), m.β), piv)
 end
 
@@ -87,7 +87,7 @@ function StatsBase.deviance(m::GeneralizedLinearMixedModel{T}, nAGQ = 1) where {
     u = vec(first(m.u))
     u₀ = vec(first(m.u₀))
     copyto!(u₀, u)
-    ra = RaggedArray(m.resp.devresid, first(m.LMM.allterms).refs)
+    ra = RaggedArray(m.resp.devresid, first(m.LMM.reterms).refs)
     devc0 = sum!(map!(abs2, m.devc0, u), ra)  # the deviance components at z = 0
     sd = map!(inv, m.sd, first(m.LMM.L).diag)
     mult = fill!(m.mult, 0)
@@ -120,8 +120,19 @@ end
 objective(m::GeneralizedLinearMixedModel) = deviance(m)
 
 """
-deviance!(m::GeneralizedLinearMixedModel, nAGQ=1)
+    GLM.wrkresp!(v::AbstractVector{T}, resp::GLM.GlmResp{AbstractVector{T}})
 
+A copy of a method from GLM that generalizes the types in the signature
+"""
+function GLM.wrkresp!(v::AbstractVector{T}, r::GLM.GlmResp{Vector{T}}) where {T<:AbstractFloat}
+    v .= r.eta .+ r.wrkresid
+    isempty(r.offset) && return v
+    v .-= r.offset
+end
+
+"""
+    deviance!(m::GeneralizedLinearMixedModel, nAGQ=1)
+    
 Update `m.η`, `m.μ`, etc., install the working response and working weights in
 `m.LMM`, update `m.LMM.A` and `m.LMM.R`, then evaluate the [`deviance`](@ref).
 """
@@ -359,9 +370,9 @@ function GeneralizedLinearMixedModel(
     if isempty(wts)
         LMM = LinearMixedModel(
             LMM.formula,
-            LMM.allterms,
             LMM.reterms,
-            LMM.feterms,
+            LMM.Xymat,
+            LMM.feterm,
             fill!(similar(y), 1),
             LMM.parmap,
             LMM.dims,
@@ -414,7 +425,7 @@ function Base.getproperty(m::GeneralizedLinearMixedModel, s::Symbol)
         σs(m)
     elseif s == :σρs
         σρs(m)
-    elseif s ∈ (:A, :L, :optsum, :allterms, :reterms, :feterms, :formula)
+    elseif s ∈ (:A, :L, :optsum, :reterms, :Xymat, :feterm, :formula)
         getfield(m.LMM, s)
     elseif s ∈ (:dims, :λ, :lowerbd, :corr, :PCA, :rePCA, :X,)
         getproperty(m.LMM, s)
@@ -508,7 +519,12 @@ function pirls!(
     for j in eachindex(u)         # start from u all zeros
         copyto!(u₀[j], fill!(u[j], 0))
     end
-    varyβ && copyto!(β₀, β)
+    if varyβ
+        copyto!(β₀, β)
+        Llast = last(lm.L)
+        pp1 = size(Llast, 1)
+        Ltru = view(Llast, pp1, 1:(pp1 - 1)) # name read as L'u
+    end
     obj₀ = deviance!(m) * 1.0001
     if verbose
         print("varyβ = ", varyβ, ", obj₀ = ", obj₀)
@@ -519,7 +535,7 @@ function pirls!(
         println()
     end
     for iter = 1:maxiter
-        varyβ && ldiv!(adjoint(feL(m)), copyto!(β, lm.L[end-1]))
+        varyβ && ldiv!(adjoint(feL(m)), copyto!(β, Ltru))
         ranef!(u, m.LMM, β, true) # solve for new values of u
         obj = deviance!(m)        # update GLM vecs and evaluate Laplace approx
         verbose && println(lpad(iter, 4), ": ", obj)
@@ -551,7 +567,7 @@ end
 
 ranef(m::GeneralizedLinearMixedModel; uscale::Bool=false) = ranef(m.LMM, uscale=uscale)
 
-LinearAlgebra.rank(m::GeneralizedLinearMixedModel) = first(m.LMM.feterms).rank
+LinearAlgebra.rank(m::GeneralizedLinearMixedModel) = m.LMM.feterm.rank
 
 """
     refit!(m::GeneralizedLinearMixedModel[, y::Vector];
diff --git a/src/linalg.jl b/src/linalg.jl
index d7783db8d..e63429919 100644
--- a/src/linalg.jl
+++ b/src/linalg.jl
@@ -31,7 +31,7 @@ LinearAlgebra.mul!(
     adjB::Adjoint{T,<:BlockedSparse{T}},
     α::Number,
     β::Number,
-) where {T} = mul!(C, A, adjB.parent.cscmat', α, β)
+) where {T} = mul!(C, A, adjoint(adjB.parent.cscmat), α, β)
 
 LinearAlgebra.mul!(
     C::StridedVector{T},
@@ -39,7 +39,7 @@ LinearAlgebra.mul!(
     B::StridedVector{T},
     α::Number,
     β::Number,
-) where {T} = mul!(C, adjA.parent.cscmat', B, α, β)
+) where {T} = mul!(C, adjoint(adjA.parent.cscmat), B, α, β)
 
 @static if VERSION < v"1.6.0-DEV.1468"
     function LinearAlgebra.ldiv!(
diff --git a/src/linalg/logdet.jl b/src/linalg/logdet.jl
index ec6a2f759..a3aadb063 100644
--- a/src/linalg/logdet.jl
+++ b/src/linalg/logdet.jl
@@ -31,7 +31,9 @@ function LinearAlgebra.logdet(m::LinearMixedModel{T}) where {T}
         s += LD(L[kp1choose2(j)])
     end
     if m.optsum.REML
-        s += LD(L[kp1choose2(nre + 1)])
+        lastL = last(L)
+        s += LD(lastL)  # this adds the log of sqrtpwrss
+        s -= log(last(lastL))
     end
     s + s
 end
diff --git a/src/linearmixedmodel.jl b/src/linearmixedmodel.jl
index 78443b002..56d99e9b7 100644
--- a/src/linearmixedmodel.jl
+++ b/src/linearmixedmodel.jl
@@ -6,9 +6,9 @@ Linear mixed-effects model representation
 ## Fields
 
 * `formula`: the formula for the model
-* `allterms`: a vector of random-effects terms, the fixed-effects terms and the response
 * `reterms`: a `Vector{AbstractReMat{T}}` of random-effects terms.
-* `feterms`: a `Vector{FeMat{T}}` of the fixed-effects model matrix and the response
+* `Xymat`: horizontal concatenation of a full-rank fixed-effects model matrix `X` and response `y` as an `FeMat{T}`
+* `feterm`: the fixed-effects model matrix as an `FeTerm{T}`
 * `sqrtwts`: vector of square roots of the case weights.  Can be empty.
 * `parmap` : Vector{NTuple{3,Int}} of (block, row, column) mapping of θ to λ
 * `dims` : NamedTuple{(:n, :p, :nretrms),NTuple{3,Int}} of dimensions.  `p` is the rank of `X`, which may be smaller than `size(X, 2)`.
@@ -30,9 +30,9 @@ Linear mixed-effects model representation
 """
 struct LinearMixedModel{T<:AbstractFloat} <: MixedModel{T}
     formula::FormulaTerm
-    allterms::Vector{Union{AbstractReMat{T}, FeMat{T}}}
     reterms::Vector{AbstractReMat{T}}
-    feterms::Vector{FeMat{T}}
+    Xymat::FeMat{T}
+    feterm::FeTerm{T}
     sqrtwts::Vector{T}
     parmap::Vector{NTuple{3,Int}}
     dims::NamedTuple{(:n, :p, :nretrms),NTuple{3,Int}}
@@ -86,19 +86,17 @@ Everything else in the structure can be derived from these quantities.
     a future release without being considered a breaking change.
 """
 function LinearMixedModel(y::AbstractArray,
-                           Xs::Tuple, # can't be more specific here without stressing the compiler
-                           form::FormulaTerm, wts = [])
-    y = reshape(float(y), (:, 1)) # y as a floating-point matrix
-    T = promote_type(Float64, eltype(y))  # ensure that eltype of model matrices is at least Float64
-    y = convert(Matrix{T}, y)
+                        Xs::Tuple, # can't be more specific here without stressing the compiler
+                        form::FormulaTerm, wts = [])
+    T = promote_type(Float64, float(eltype(y)))  # ensure eltype of model matrices is at least Float64
 
     reterms = AbstractReMat{T}[]
-    feterms = FeMat{T}[]
+    feterms = FeTerm{T}[]
     for (i, x) in enumerate(Xs)
         if isa(x, AbstractReMat{T})
             push!(reterms, x)
         elseif isa(x, ReMat) # this can occur in weird situation where x is a ReMat{U}
-            # avoid keeping a second copy if unweighted
+                             # avoid keeping a second copy if unweighted
             z = convert(Matrix{T}, x.z)
             wtz = x.z === x.wtz ? z : convert(Matrix{T}, x.wtz)
             S = size(z, 1)
@@ -111,16 +109,14 @@ function LinearMixedModel(y::AbstractArray,
             push!(reterms, x)
         else
             cnames = coefnames(form.rhs[i])
-            push!(feterms, FeMat(x, isa(cnames, String) ? [cnames] : collect(cnames)))
+            push!(feterms, FeTerm(x, isa(cnames, String) ? [cnames] : collect(cnames)))
         end
     end
-    push!(feterms, FeMat(y, [""]))
-
-    return LinearMixedModel(feterms, reterms, form, wts)
+    return LinearMixedModel(convert(Array{T}, y), only(feterms), reterms, form, wts)
 end
 
 """
-    LinearMixedModel(feterms, reterms, form, wts=[])
+    LinearMixedModel(feterm, reterms, form, wts=[])
 
 Private constructor for a `LinearMixedModel` given already assembled fixed and random effects.
 
@@ -133,32 +129,31 @@ can be derived from these quantities.
     This method is internal and experimental and so may change or disappear in
     a future release without being considered a breaking change.
 """
-function LinearMixedModel(feterms::Vector{FeMat{T}}, reterms::Vector{AbstractReMat{T}},
+function LinearMixedModel(y::Array{T}, feterm::FeTerm{T}, reterms::Vector{AbstractReMat{T}},
                           form::FormulaTerm, wts=[]) where T
-
     # detect and combine RE terms with the same grouping var
     if length(reterms) > 1
         reterms = amalgamate(reterms)
     end
 
     sort!(reterms, by = nranef, rev = true)
-    allterms = convert(Vector{Union{AbstractReMat{T},FeMat{T}}}, vcat(reterms, feterms))
+    Xy = FeMat(feterm, vec(y))
     sqrtwts = sqrt.(convert(Vector{T}, wts))
-    reweight!.(allterms, Ref(sqrtwts))
-    A, L = createAL(allterms)
+    reweight!.(reterms, Ref(sqrtwts))
+    reweight!(Xy, sqrtwts)
+    A, L = createAL(reterms, Xy)
     lbd = foldl(vcat, lowerbd(c) for c in reterms)
     θ = foldl(vcat, getθ(c) for c in reterms)
-    X = first(feterms)
     optsum = OptSummary(θ, lbd, :LN_BOBYQA, ftol_rel = T(1.0e-12), ftol_abs = T(1.0e-8))
     fill!(optsum.xtol_abs, 1.0e-10)
     LinearMixedModel(
         form,
-        allterms,
         reterms,
-        feterms,
+        Xy,
+        feterm,
         sqrtwts,
         mkparmap(reterms),
-        (n = size(X, 1), p = X.rank, nretrms = length(reterms)),
+        (n = length(y), p = feterm.rank, nretrms = length(reterms)),
         A,
         L,
         optsum,
@@ -239,14 +234,14 @@ fit(
 )
 
 function StatsBase.coef(m::LinearMixedModel{T}) where {T}
-    piv = first(m.feterms).piv
+    piv = m.feterm.piv
     invpermute!(fixef!(similar(piv, T), m), piv)
 end
 
 βs(m::LinearMixedModel) = NamedTuple{(Symbol.(coefnames(m))...,)}(coef(m))
 
 function StatsBase.coefnames(m::LinearMixedModel)
-    Xtrm = first(m.feterms)
+    Xtrm = m.feterm
     invpermute!(copy(Xtrm.cnames), Xtrm.piv)
 end
 
@@ -296,23 +291,29 @@ function condVar(m::LinearMixedModel{T}) where {T}
     Array{T,3}[reshape(abs2.(ll ./ Ld) .* varest(m), (1, 1, length(Ld)))]
 end
 
-function createAL(allterms::Vector{Union{AbstractReMat{T},FeMat{T}}}) where {T}
-    k = length(allterms)
-    sz = [isa(t, ReMat) ? size(t, 2) : rank(t) for t in allterms]
-    A = AbstractMatrix{T}[]
-    L = AbstractMatrix{T}[]
-    for i = 1:k
-        for j = 1:i
-            Lij = densify(allterms[i]' * allterms[j])
-            push!(L, Lij)
-            push!(A, deepcopy(isa(Lij, BlockedSparse) ? Lij.cscmat : Lij))
+function pushALblock!(A, L, blk)
+    push!(L, blk)
+    push!(A, deepcopy(isa(blk, BlockedSparse) ? blk.cscmat : blk))
+end
+
+function createAL(reterms::Vector{AbstractReMat{T}}, Xy::FeMat{T}) where {T}
+    k = length(reterms)
+    vlen = kchoose2(k+1)
+    A = sizehint!(AbstractMatrix{T}[], vlen)
+    L = sizehint!(AbstractMatrix{T}[], vlen)
+    for i in eachindex(reterms)
+        for j in 1:i
+            pushALblock!(A, L, densify(reterms[i]' * reterms[j]))
         end
     end
-    nretrm = sum(Base.Fix2(isa, ReMat), allterms)
-    for i = 2:nretrm      # check for fill-in due to non-nested grouping factors
-        ci = allterms[i]
-        for j = 1:(i-1)
-            cj = allterms[j]
+    for j in eachindex(reterms)   # can't fold this into the previous loop b/c block order
+        pushALblock!(A, L, densify(Xy' * reterms[j]))
+    end
+    pushALblock!(A, L, densify(Xy'Xy))
+    for i in 2:k      # check for fill-in due to non-nested grouping factors
+        ci = reterms[i]
+        for j in 1:(i-1)
+            cj = reterms[j]
             if !isnested(cj, ci)
                 for l = i:k
                     ind = block(l, i)
@@ -350,7 +351,10 @@ Return the lower Cholesky factor for the fixed-effects parameters, as an `LowerT
 `p × p` matrix.
 """
 function feL(m::LinearMixedModel)
-    LowerTriangular(m.L[kp1choose2(m.dims.nretrms + 1)])
+    XyL = m.L[end]
+    k = size(XyL, 1)
+    inds = Base.OneTo(k - 1)
+    LowerTriangular(view(XyL, inds, inds))
 end
 
 """
@@ -406,8 +410,8 @@ end
 
 function fitted!(v::AbstractArray{T}, m::LinearMixedModel{T}) where {T}
     ## FIXME: Create and use `effects(m) -> β, b` w/o calculating β twice
-    Xtrm = first(m.feterms)
-    vv = mul!(vec(v), Xtrm, fixef!(similar(Xtrm.piv, T), m))
+    Xtrm = m.feterm
+    vv = mul!(vec(v), Xtrm.x, fixef!(similar(Xtrm.piv, T), m))
     for (rt, bb) in zip(m.reterms, ranef(m))
         unscaledre!(vv, rt, bb)
     end
@@ -426,15 +430,16 @@ the length of `v` can be the rank of `X` or the number of columns of `X`.  In th
 case the calculated coefficients are padded with -0.0 out to the number of columns.
 """
 function fixef!(v::AbstractVector{T}, m::LinearMixedModel{T}) where {T}
-    Xtrm = m.allterms[m.dims.nretrms + 1]
-    if isfullrank(Xtrm)
-        ldiv!(feL(m)', copyto!(v, m.L[end-1]))
-    else
-        ldiv!(
-            feL(m)',
-            view(copyto!(fill!(v, -zero(T)), m.L[end-1]), 1:(Xtrm.rank)),
-        )
+    Xtrm = m.feterm
+    fill!(v, -zero(T))
+    XyL = m.L[end]
+    L = feL(m)
+    k = size(XyL, 1)
+    r = size(L, 1)
+    for j in 1:r
+        v[j] = XyL[k, j]
     end
+    ldiv!(feL(m)', length(v) == r ? v : view(v, 1:r))
     v
 end
 
@@ -447,7 +452,7 @@ In the rank-deficient case the truncated parameter vector, of length `rank(m)` i
 This is unlike `coef` which always returns a vector whose length matches the number of
 columns in `X`.
 """
-fixef(m::LinearMixedModel{T}) where {T} = fixef!(Vector{T}(undef, first(m.feterms).rank), m)
+fixef(m::LinearMixedModel{T}) where {T} = fixef!(Vector{T}(undef, m.feterm.rank), m)
 
 """
     fixefnames(m::MixedModel)
@@ -455,7 +460,7 @@ fixef(m::LinearMixedModel{T}) where {T} = fixef!(Vector{T}(undef, first(m.feterm
 Return a (permuted and truncated in the rank-deficient case) vector of coefficient names.
 """
 function fixefnames(m::LinearMixedModel{T}) where {T}
-    Xtrm = first(m.feterms)
+    Xtrm = m.feterm
     Xtrm.cnames[1:Xtrm.rank]
 end
 
@@ -479,12 +484,12 @@ function getθ!(v::AbstractVector{T}, m::LinearMixedModel{T}) where {T}
         throw(DimensionMismatch("length(v) = $(length(v)) ≠ length(m.parmap) = $(length(pmap))"))
     end
     reind = 1
-    λ = first(m.allterms).λ
+    λ = first(m.reterms).λ
     for (k, tp) in enumerate(pmap)
         tp1 = first(tp)
         if reind ≠ tp1
             reind = tp1
-            λ = m.allterms[tp1].λ
+            λ = m.reterms[tp1].λ
         end
         v[k] = λ[tp[2], tp[3]]
     end
@@ -527,7 +532,9 @@ function Base.getproperty(m::LinearMixedModel{T}, s::Symbol) where {T}
     elseif s == :X
         modelmatrix(m)
     elseif s == :y
-        vec(last(m.feterms).x)
+        let xy = m.Xymat.xy
+            view(xy, :, size(xy, 2))
+        end
     elseif s == :rePCA
         rePCA(m)
     else
@@ -540,13 +547,23 @@ function StatsBase.leverage(m::LinearMixedModel{T}) where {T}
     # The i'th leverage value is obtained by replacing the response with the i'th
     # basis vector, updating A and L, then taking the sum of squared values of the
     # last row of L, excluding the last position.
-    yorig = copy(m.y)
-    l = length(m.allterms)
+    yview = m.y
+    yorig = copy(yview)
+    kp1 = length(m.reterms) + 1
+    L = m.L
+    lastL = last(L)
+    pp1 = size(lastL, 1)
+    p = pp1 - 1
     value = map(eachindex(yorig)) do i
-        fill!(m.y, zero(T))
-        m.y[i] = one(T)
+        fill!(yview, zero(T))
+        yview[i] = one(T)
         updateL!(reevaluateAend!(m))
-        sum(j -> sum(abs2, m.L[block(l, j)]), 1:(l-1))
+        s = sum(abs2, view(lastL, pp1, Base.OneTo(p)))
+        for j in eachindex(m.reterms)
+            Lblock = m.L[block(kp1, j)]
+            s += sum(abs2, view(Lblock, pp1, :))
+        end
+        s
     end
     copyto!(m.y, yorig)
     updateL!(reevaluateAend!(m))
@@ -574,14 +591,7 @@ function mkparmap(reterms::Vector{AbstractReMat{T}}) where {T}
     parmap
 end
 
-function StatsBase.modelmatrix(m::LinearMixedModel)
-    fe = first(m.feterms)
-    if fe.rank == size(fe, 2)
-        fe.x
-    else
-        fe.x[:, invperm(fe.piv)]
-    end
-end
+StatsBase.modelmatrix(m::LinearMixedModel) = m.feterm.x
 
 nθ(m::LinearMixedModel) = length(m.parmap)
 
@@ -628,8 +638,8 @@ Base.propertynames(m::LinearMixedModel, private::Bool = false) = (
     :PCA,
     :rePCA,
     :reterms,
-    :feterms,
-    :allterms,
+    :feterm,
+    :Xymat,
     :objective,
     :pvalues,
 )
@@ -649,23 +659,18 @@ Overwrite `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::LinearMixedModel{T},
-    β::AbstractArray{T},
-    uscale::Bool,
-) where {T}
-    (k = length(v)) == m.dims.nretrms || throw(DimensionMismatch(""))
+function ranef!(v::Vector, m::LinearMixedModel{T}, β::AbstractArray{T}, uscale::Bool) where {T}
+    (k = length(v)) == length(m.reterms) || throw(DimensionMismatch(""))
     L = m.L
-    kp2 = length(m.allterms)
-    for j = 1:k
-        mul!(
-            vec(copyto!(v[j], L[block(kp2, j)])),
-            L[block(k + 1, j)]',
-            β,
-            -one(T),
-            one(T),
-        )
+    lind = length(L)
+    for j in k:-1:1
+        lind -= 1
+        Ljkp1 = L[lind]
+        vj = v[j]
+        length(vj) == size(Ljkp1, 2) || throw(DimensionMismatch(""))
+        pp1 = size(Ljkp1, 1)
+        copyto!(vj, view(Ljkp1, pp1, :))
+        mul!(vec(vj), view(Ljkp1, 1:(pp1-1), :)', β, -one(T), one(T))
     end
     for i = k:-1:1
         Lii = L[kp1choose2(i)]
@@ -693,7 +698,7 @@ Return, as a `Vector{Matrix{T}}`, the conditional modes of the random effects in
 If `uscale` is `true` the random effects are on the spherical (i.e. `u`) scale, otherwise on
 the original scale.
 
-For a named variant, see [`@raneftables`](@ref).
+For a named variant, see [`raneftables`](@ref).
 """
 function ranef(m::LinearMixedModel{T}; uscale = false) where {T}
     reterms = m.reterms
@@ -701,7 +706,7 @@ function ranef(m::LinearMixedModel{T}; uscale = false) where {T}
     ranef!(v, m, uscale)
 end
 
-LinearAlgebra.rank(m::LinearMixedModel) = first(m.feterms).rank
+LinearAlgebra.rank(m::LinearMixedModel) = m.feterm.rank
 
 """
     rePCA(m::LinearMixedModel; corr::Bool=true)
@@ -733,16 +738,20 @@ end
 """
     reevaluateAend!(m::LinearMixedModel)
 
-Reevaluate the last column of `m.A` from `m.feterms`.  This function should be called
-after updating the response, `m.feterms[end]`.
+Reevaluate the last column of `m.A` from `m.Xymat`.  This function should be called
+after updating the response.
 """
 function reevaluateAend!(m::LinearMixedModel)
     A = m.A
-    trmn = reweight!(last(m.allterms), m.sqrtwts)
-    nblk = length(m.allterms)
-    for (j, trm) in enumerate(m.allterms)
-        mul!(A[block(nblk, j)], trmn', trm)
+    reterms = m.reterms
+    nre = length(reterms)
+    trmn = reweight!(m.Xymat, m.sqrtwts)
+    ind = kp1choose2(nre)
+    for trm in m.reterms
+        ind += 1
+        mul!(A[ind], trmn', trm)
     end
+    mul!(A[end], trmn', trmn)
     m
 end
 
@@ -758,19 +767,20 @@ function refit!(m::LinearMixedModel; REML=m.optsum.REML)
 end
 
 function refit!(m::LinearMixedModel, y; REML=m.optsum.REML)
-    resp = last(m.feterms)
-    length(y) == size(resp, 1) || throw(DimensionMismatch(""))
+    resp = m.y
+    length(y) == length(resp) || throw(DimensionMismatch(""))
     copyto!(resp, y)
     refit!(m; REML=REML)
 end
 
 StatsBase.residuals(m::LinearMixedModel) = response(m) .- fitted(m)
 
-StatsBase.response(m::LinearMixedModel) = vec(last(m.feterms).x)
+StatsBase.response(m::LinearMixedModel) = m.y
 
 function reweight!(m::LinearMixedModel, weights)
     sqrtwts = map!(sqrt, m.sqrtwts, weights)
-    reweight!.(m.allterms, Ref(sqrtwts))
+    reweight!.(m.reterms, Ref(sqrtwts))
+    reweight!(m.Xymat, sqrtwts)
     updateA!(m)
     updateL!(m)
 end
@@ -842,17 +852,6 @@ end
 Base.show(io::IO, m::LinearMixedModel) = Base.show(io, MIME"text/plain"(), m)
 
 
-#=
-"""
-    sqrtpwrss(m::LinearMixedModel)
-
-Return the square root of the penalized, weighted residual sum-of-squares (pwrss).
-
-This is the element in the lower-right of `m.L`
-"""
-sqrtpwrss(m::LinearMixedModel) = last(m.L)
-=#
-
 """
     ssqdenom(m::LinearMixedModel)
 
@@ -900,30 +899,38 @@ function stderror!(v::AbstractVector{T}, m::LinearMixedModel{T}) where {T}
         scr[i] = true
         v[i] = s * norm(ldiv!(L, scr))
     end
-    invpermute!(v, first(m.feterms).piv)
+    invpermute!(v, m.feterm.piv)
     v
 end
 
 function StatsBase.stderror(m::LinearMixedModel{T}) where {T}
-    stderror!(similar(first(m.feterms).piv, T), m)
+    stderror!(similar(m.feterm.piv, T), m)
 end
 
 """
     updateA!(m::LinearMixedModel)
 
-Update the cross-product array, `m.A`, from `m.allterms`
+Update the cross-product array, `m.A`, from `m.reterms` and `m.Xymat`
 
 This is usually done after a reweight! operation.
 """
 function updateA!(m::LinearMixedModel)
-    allterms = m.allterms
-    k = length(allterms)
+    reterms = m.reterms
+    k = length(reterms)
     A = m.A
-    for j = 1:k
-        for i = j:k
-            mul!(A[block(i,j)], allterms[i]', allterms[j])
+    ind = 1
+    for (i, trmi) in enumerate(reterms)
+        for j = 1:i
+            mul!(A[ind], trmi', reterms[j])
+            ind += 1
         end
     end
+    Xymattr = adjoint(m.Xymat)
+    for trm in reterms
+        mul!(A[ind], Xymattr, trm)
+        ind += 1
+    end
+    mul!(A[end], Xymattr, m.Xymat)
     m
 end
 
@@ -948,29 +955,29 @@ Update the blocked lower Cholesky factor, `m.L`, from `m.A` and `m.reterms` (use
 This is the crucial step in evaluating the objective, given a new parameter value.
 """
 function updateL!(m::LinearMixedModel{T}) where {T}
-    A = m.A
-    L = m.L
-    k = length(m.allterms)
+    A, L, reterms = m.A, m.L, m.reterms
+    k = length(reterms)
     for (l, a) in zip(L, A) # copy each element of A to corresponding element of L
         copyto!(l, a)       # For some reason copyto!.(L, A) allocates a lot of memory
     end
-    for (j, cj) in enumerate(m.reterms)  # pre- and post-multiply by Λ, add I to diagonal
+    for j in eachindex(reterms) # pre- and post-multiply by Λ, add I to diagonal
+        cj = reterms[j]
         scaleinflate!(L[kp1choose2(j)], cj)
-        for i = (j+1):k         # postmultiply column by Λ
+        for i = (j+1):(k+1)     # postmultiply column by Λ
             rmulΛ!(L[block(i,j)], cj)
         end
         for jj = 1:(j-1)        # premultiply row by Λ'
             lmulΛ!(cj', L[block(j, jj)])
         end
     end
-    for j = 1:k                         # blocked Cholesky
+    for j = 1:(k+1)             # blocked Cholesky
         Ljj = L[kp1choose2(j)]
         for jj = 1:(j-1)
             rankUpdate!(Hermitian(Ljj, :L), L[block(j, jj)], -one(T), one(T))
         end
         cholUnblocked!(Ljj, Val{:L})
         LjjT = isa(Ljj, Diagonal) ? Ljj : LowerTriangular(Ljj)
-        for i = (j+1):k
+        for i = (j+1):(k+1)
             Lij = L[block(i, j)]
             for jj = 1:(j-1)
                 mul!(Lij, L[block(i, jj)], L[block(j, jj)]', -one(T), one(T))
diff --git a/src/mixedmodel.jl b/src/mixedmodel.jl
index 8d8b42841..c8ee8e1eb 100644
--- a/src/mixedmodel.jl
+++ b/src/mixedmodel.jl
@@ -57,10 +57,10 @@ end
     vcov(m::MixedModel; corr=false)
 
 Returns the variance-covariance matrix of the fixed effects.
-If `corr=true`, then correlation of fixed effects is returned instead.
+If `corr` is `true`, the correlation of the fixed effects is returned instead.
 """
 function StatsBase.vcov(m::MixedModel; corr=false)
-    Xtrm = first(m isa GeneralizedLinearMixedModel ? m.LMM.feterms : m.feterms)
+    Xtrm = m isa GeneralizedLinearMixedModel ? m.LMM.feterm : m.feterm
     iperm = invperm(Xtrm.piv)
     p = length(iperm)
     r = Xtrm.rank
diff --git a/src/remat.jl b/src/remat.jl
index 75f622140..419964c3a 100644
--- a/src/remat.jl
+++ b/src/remat.jl
@@ -252,21 +252,20 @@ end
 
 *(adjA::Adjoint{T,<:ReMat{T}}, B::ReMat{T}) where {T} = adjA.parent.adjA * sparse(B)
 *(adjA::Adjoint{T,<:FeMat{T}}, B::ReMat{T}) where {T} =
-    mul!(Matrix{T}(undef, rank(adjA.parent), size(B, 2)), adjA, B)
+    mul!(Matrix{T}(undef, size(adjA.parent, 2), size(B, 2)), adjA, B)
 
 function LinearAlgebra.mul!(C::Matrix{T}, adjA::Adjoint{T,<:FeMat{T}}, B::ReMat{T,1},
         α::Number, β::Number) where {T}
     A = adjA.parent
-    Awt = A.wtx
+    Awt = A.wtxy
     n, p = size(Awt)
-    r = A.rank
     m, q = size(B)
-    size(C) == (r, q) && m == n || throw(DimensionMismatch())
+    size(C) == (p, q) && m == n || throw(DimensionMismatch())
     isone(β) || rmul!(C, β)
     zz = B.wtz
     @inbounds for (j, rrj) in enumerate(B.refs)
         αzj = α * zz[j]
-        for i in 1:r
+        for i in 1:p
             C[i, rrj] += αzj * Awt[j, i]
         end
     end
@@ -276,8 +275,8 @@ end
 function LinearAlgebra.mul!(C::Matrix{T}, adjA::Adjoint{T,<:FeMat{T}}, B::ReMat{T,S},
     α::Number, β::Number) where {T,S}
     A = adjA.parent
-    Awt = A.wtx
-    r = rank(A)
+    Awt = A.wtxy
+    r = size(Awt, 2)
     rr = B.refs
     scr = B.scratch
     vscr = vec(scr)
diff --git a/src/simulate.jl b/src/simulate.jl
index 5aa2dfb4b..782b69c5d 100644
--- a/src/simulate.jl
+++ b/src/simulate.jl
@@ -79,7 +79,8 @@ function simulate!(
 
     # the noise term is actually in the GLM and not the LMM part so no noise
     # at the LMM level
-    η = fill!(m.LMM.y, zero(T))
+    η = fill!(copy(m.LMM.y), zero(T))  # ensure that η is a vector - needed for GLM.updateμ! below
+                                       # A better approach is to change the signature for updateμ!
     y = m.resp.y
 
     # assemble the linear predictor
diff --git a/test/UniformBlockDiagonal.jl b/test/UniformBlockDiagonal.jl
index 68e944e1d..4ff35d50d 100644
--- a/test/UniformBlockDiagonal.jl
+++ b/test/UniformBlockDiagonal.jl
@@ -65,6 +65,8 @@ const LMM = LinearMixedModel
         L21 = vf2'vf1
         @test isa(L21, BlockedSparse)
         @test L21[1,1] == 30.0
+        @test size(L21) == (344, 9452)
+        @test size(L21, 1) == 344
         MixedModels.lmulΛ!(vf2', MixedModels.rmulΛ!(L21, vf1))
         @test size(Matrix(L21)) == size(sparse(L21))
 #        L21cb1 = copy(L21.colblocks[1])
diff --git a/test/grouping.jl b/test/grouping.jl
index 930746444..6bc31d5a0 100644
--- a/test/grouping.jl
+++ b/test/grouping.jl
@@ -1,6 +1,11 @@
 using Test
 using StatsModels
 
+@testset "Grouping" begin
+    g = Grouping()
+    @test isnothing(g.levels)
+end
+
 @testset "Grouping pseudo-contrasts" begin
     d = (y = rand(2_000_000), grp=string.([1:1_000_000; 1:1_000_000]))
     ## OOM seems to result in the process being killed on Mac so this messes up CI
@@ -10,6 +15,7 @@ using StatsModels
     @test t isa CategoricalTerm{Grouping}
     @test size(t.contrasts.matrix) == (0,0)
     @test length(t.contrasts.levels) == 1_000_000
+    @test_throws ErrorException StatsModels.modelcols(t, (a = 1.,))
 
     levs = sort(string.(1:1_000_000))
 
diff --git a/test/markdown.jl b/test/markdown.jl
index 82cf3ee7b..6b73713de 100644
--- a/test/markdown.jl
+++ b/test/markdown.jl
@@ -92,7 +92,7 @@ include("modelcache.jl")
 |:----|:------------:|:------------:|:------------:|
 |316  |Diagonal      |              |              |
 |24   |Dense         |Diag/Dense    |              |
-|6    |Dense         |Dense         |Dense         |
+|7    |Dense         |Dense         |Dense         |
 """
 
     end
diff --git a/test/matrixterm.jl b/test/matrixterm.jl
index fe3951480..d3f337de6 100644
--- a/test/matrixterm.jl
+++ b/test/matrixterm.jl
@@ -2,74 +2,70 @@ using LinearAlgebra, MixedModels, StableRNGs, Test, SparseArrays
 
 include("modelcache.jl")
 
-@testset "femat" begin
-    trm = MixedModels.FeMat(hcat(ones(30), repeat(0:9, outer = 3)), ["(Intercept)", "U"])
+@testset "Xymat" begin
+    trm = MixedModels.FeTerm(hcat(ones(30), repeat(0:9, outer = 3)), ["(Intercept)", "U"])
     piv = trm.piv
     ipiv = invperm(piv)
-    @test size(trm) == (30, 2)
-    @test length(trm) == 60
-    @test size(trm') == (2, 30)
-    @test eltype(trm) == Float64
-    @test trm.x === trm.wtx
-    prd = trm'trm
+    mat = MixedModels.FeMat(trm, Float64.(collect(axes(trm.x, 1))))
+    @test size(mat) == (30, 3)
+    @test length(mat) == 90
+    @test size(mat') == (3, 30)
+    @test eltype(mat) == Float64
+    @test mat.xy === mat.wtxy
+    prd = mat'mat
     @test typeof(prd) == Matrix{Float64}
-    @test prd == [30.0 135.0; 135.0 855.0][piv, piv]
+    @test prd[ipiv, ipiv] == [30.0 135.0; 135.0 855.0]
     wts = rand(StableRNG(123454321), 30)
-    MixedModels.reweight!(trm, wts)
-    @test mul!(prd, trm', trm)[ipiv[1], ipiv[1]] ≈ sum(abs2, wts)
+    MixedModels.reweight!(mat, wts)
+    @test mul!(prd, mat', mat)[ipiv[1], ipiv[1]] ≈ sum(abs2, wts)
 
     # empty fixed effects
-    trm = MixedModels.FeMat(ones(10,0), String[])
-    @test size(trm) == (10, 0)
-    @test length(trm) == 0
-    @test size(trm') == (0, 10)
+    trm = MixedModels.FeTerm(ones(10,0), String[])
+    #@test size(trm) == (10, 0)  # There no longer are size and length methods for FeTerm
+    #@test length(trm) == 0
+    #@test size(trm') == (0, 10)
     @test eltype(trm) == Float64
     @test trm.rank == 0
 end
 
-@testset "fematSparse" begin
+@testset "XymatSparse" begin
 
     @testset "sparse and dense yield same fit" begin
         # deepcopy because we're going to modify
-        m = deepcopy(last(models(:insteval)))
+        m = last(models(:insteval))
         # this is kinda sparse: 
-        # julia> mean(first(m.feterms).x)
+        # julia> mean(first(m.feterm).x)
         # 0.10040140325753434
         
-        fe = first(m.feterms)
-        X =  MixedModels.FeMat(SparseMatrixCSC(fe.x), fe.cnames)
+        fe = m.feterm
+        X =  MixedModels.FeTerm(SparseMatrixCSC(fe.x), fe.cnames)
         @test typeof(X.x) <: SparseMatrixCSC
-        @test typeof(X.wtx) <: SparseMatrixCSC
         @test X.rank == 28
         @test X.cnames == fe.cnames
-
-        dense_θ = copy(m.θ)
-        dense_feval = m.optsum.feval
-        m.feterms[1] = X    
-        refit!(m)
-
-        @test dense_θ ≈ m.θ
+        m1 = fit!(LinearMixedModel(collect(m.y), X, deepcopy(m.reterms), m.formula))
+        @test isapprox(m1.θ, m.θ, rtol = 1.0e-5)
     end
 
-    @testset "rank defiency in sparse FeMat" begin
-        trm = MixedModels.FeMat(SparseMatrixCSC(hcat(ones(30), 
+    @testset "rank defiency in sparse FeTerm" begin
+        trm = MixedModels.FeTerm(SparseMatrixCSC(hcat(ones(30), 
                                                      repeat(0:9, outer = 3), 
                                                      2repeat(0:9, outer = 3))), 
                                 ["(Intercept)", "U", "V"])
         piv = trm.piv
         ipiv = invperm(piv)
         @test_broken rank(trm) == 2
-        @test size(trm) == (30, 3)
-        @test length(trm) == 90
-        @test size(trm') == (3, 30)
-        @test eltype(trm) == Float64
-        @test trm.x === trm.wtx
-        prd = trm'trm
+        mat = MixedModels.FeMat(trm, Float64.(collect(axes(trm.x, 1))))
+        @test size(mat) == (30, 4)
+        @test length(mat) == 120
+        @test size(mat') == (4, 30)
+        @test eltype(mat) == Float64
+        @test mat.xy === mat.wtxy
+        prd = MixedModels.densify(mat'mat)
         @test typeof(prd) == typeof(prd)
-        @test prd == [30.0 135.0 270.0; 135.0 855.0 1710.0; 270.0 1710.0 3420][piv, piv]
+        @test prd[ipiv, ipiv] == [30.0 135.0 270.0; 135.0 855.0 1710.0; 270.0 1710.0 3420.0]
         wts = rand(StableRNG(123454321), 30)
-        MixedModels.reweight!(trm, wts)
-        @test mul!(prd, trm', trm)[ipiv[1], ipiv[1]] ≈ sum(abs2, wts)
+        MixedModels.reweight!(mat, wts)
+        @test mul!(prd, mat', mat)[ipiv[1], ipiv[1]] ≈ sum(abs2, wts)
     end
     
-end
\ No newline at end of file
+end
diff --git a/test/pirls.jl b/test/pirls.jl
index c19de046e..69e7274be 100644
--- a/test/pirls.jl
+++ b/test/pirls.jl
@@ -61,7 +61,6 @@ include("modelcache.jl")
     #@test isapprox(sum(gm1.resp.devresid), 2237.349, atol=0.1)
     show(IOBuffer(), gm1)
     show(IOBuffer(), BlockDescription(gm0))
-
 end
 
 @testset "cbpp" begin
diff --git a/test/pls.jl b/test/pls.jl
index 0c5aad948..74db63d06 100644
--- a/test/pls.jl
+++ b/test/pls.jl
@@ -8,7 +8,7 @@ using StatsModels
 using Tables
 using Test
 
-using MixedModels: dataset, likelihoodratiotest
+using MixedModels: likelihoodratiotest
 
 const io = IOBuffer()
 
@@ -17,7 +17,7 @@ include("modelcache.jl")
 @testset "Dyestuff" begin
     fm1 = only(models(:dyestuff))
 
-    @test length(fm1.allterms) == 3
+    @test length(fm1.A) == 3
     @test size(fm1.reterms) == (1, )
     @test lowerbd(fm1) == zeros(1)
     @test fm1.lowerbd == zeros(1)
@@ -31,7 +31,6 @@ include("modelcache.jl")
     @test_logs (:warn, "Model has not been fit") show(fm1)
 
     @test objective(updateL!(setθ!(fm1, [0.713]))) ≈ 327.34216280955366
-    @test_deprecated MixedModels.describeblocks(IOBuffer(), fm1)
 
     show(io, BlockDescription(fm1))
     @test countlines(seekstart(io)) == 3
@@ -92,7 +91,7 @@ include("modelcache.jl")
 
     @test logdet(fm1) ≈ 8.06014522999825 atol=0.001
     @test varest(fm1) ≈ 2451.2501089607676 atol=0.001
-    @test pwrss(fm1) ≈ 73537.49947885796 atol=0.001
+    @test pwrss(fm1) ≈ 73537.50152584909 atol=0.01 # this quantity is not precisely estimated
     @test stderror(fm1) ≈ [17.69455188898009] atol=0.0001
 
     vc = VarCorr(fm1)
@@ -149,7 +148,7 @@ end
     @test objective(fm) ≈ 332.18834867227616 atol=0.001
     @test coef(fm) ≈ [22.97222222222222] atol=0.001
     @test fixef(fm) ≈ [22.97222222222222] atol=0.001
-    @test coef(fm)[1] ≈ mean(dataset(:penicillin).diameter)
+    @test coef(fm)[1] ≈ mean(MixedModels.dataset(:penicillin).diameter)
     @test stderror(fm) ≈ [0.7445960346851368] atol=0.0001
     @test fm.θ ≈ [1.5375772376554968, 3.219751321180035] atol=0.001
     @test first(std(fm)) ≈ [0.8455645948223015] atol=0.0001
@@ -422,7 +421,7 @@ end
     θminqa = [1.6455, -0.2430, 1.0160, 0.8955, 2.7054, 0.0898]
     # very loose tolerance for unstable fit
     # but this is a convenient test of rankUpdate!(::UniformBlockDiagonal)
-    @test all(isapprox.(m.θ, θnlopt; atol=1e-2))
+    @test isapprox(m.θ, θnlopt; atol=5e-2)
 end
 
 @testset "Rank deficient" begin
@@ -437,8 +436,8 @@ end
     @test ct.rownms ==  ["(Intercept)", "x", "x2"]
     @test length(fixefnames(model)) == 2
     @test coefnames(model) == ["(Intercept)", "x", "x2"]
-    piv = first(model.feterms).piv
-    r = first(model.feterms).rank
+    piv = model.feterm.piv
+    r = model.feterm.rank
     @test coefnames(model)[piv][1:r] == fixefnames(model)
 end
 
@@ -471,7 +470,7 @@ end
 end
 
 @testset "unifying ReMat eltypes" begin
-    sleepstudy = dataset(:sleepstudy)
+    sleepstudy = MixedModels.dataset(:sleepstudy)
 
     re = LinearMixedModel(@formula(reaction ~ 1 + days + (1|subj) + (days|subj)), sleepstudy).reterms
     # make sure that the eltypes are still correct
diff --git a/test/runtests.jl b/test/runtests.jl
index ea23b65f8..787c3aa57 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,10 +4,10 @@ import LinearAlgebra: BLAS
 
 # there seem to be processor-specific issues and knowing this is helpful
 println(versioninfo())
-@show BLAS.vendor()
-@static if VERSION < v"1.7.0-DEV.610"
-    # this has gone away but it would be good to still
-    # get a lot of BLAS info for debugging purposes
+@static if VERSION ≥ v"1.7.0-DEV.620"
+    @show getproperty.(BLAS.get_config().loaded_libs, :libname)
+else
+    @show BLAS.vendor()
     if startswith(string(BLAS.vendor()), "openblas")
         println(BLAS.openblas_get_config())
     end