Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added opt. fname argument to sparseL and condVar #545

Merged
merged 4 commits into from
Aug 24, 2021
Merged

Added opt. fname argument to sparseL and condVar #545

Changes from 1 commit
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
80beb24
Added opt. fname argument to sparseL and condVar
dmbates Jul 22, 2021
ab113b4
Trim dead code. Add tests.
dmbates Aug 24, 2021
526dd95
Update src/linearmixedmodel.jl
dmbates Aug 24, 2021
05e8311
Blue style
palday Aug 24, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 33 additions & 5 deletions src/linearmixedmodel.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -312,7 +312,7 @@ diagonal blocks from the conditional variance-covariance matrix,
s² Λ(Λ'Z'ZΛ + I)⁻¹Λ'
"""
function condVar(m::LinearMixedModel{T}) where {T}
L = m.L
return [condVar(m, fnm) for fnm in fnames(m)]
s = sdest(m)
@static if VERSION < v"1.6.1"
spL = LowerTriangular(SparseMatrixCSC{T,Int}(sparseL(m)))
Expand Down Expand Up @@ -340,6 +340,23 @@ function condVar(m::LinearMixedModel{T}) where {T}
return val
end

function condVar(m::LinearMixedModel{T}, fname) where {T}
Lblk = LowerTriangular(densify(sparseL(m; fname=fname)))
blk = findfirst(isequal(fname), fnames(m))
λt = Array(m.λ[blk]') .* sdest(m)
vsz = size(λt, 2)
ℓ = length(m.reterms[blk].levels)
val = Array{T}(undef, (vsz, vsz, ℓ))
scratch = Matrix{T}(undef, (size(Lblk, 1), vsz))
for b in 1:ℓ
fill!(scratch, zero(T))
copyto!(view(scratch, (b - 1) * vsz .+ (1:vsz), :), λt)
ldiv!(Lblk, scratch)
mul!(view(val, :, :, b), scratch', scratch)
end
val
end

function _cvtbl(arr::Array{T,3}, trm) where {T}
return merge(
NamedTuple{(fname(trm),)}((trm.levels,)),
Expand Down Expand Up @@ -1002,19 +1019,30 @@ function _coord(A::Matrix)
end

"""
sparseL(m::LinearMixedModel{T}; full::Bool=false) where {T}
sparseL(m::LinearMixedModel; fname::Symbol=first(fnames(m)), full::Bool=false)

Return the lower Cholesky factor `L` as a `SparseMatrix{T,Int32}`.

`full` indicates whether the parts of `L` associated with the fixed-effects and response
are to be included.

When `fname` is the name of a grouping factor the blocks to the left of that block are
dropped.
"""
function sparseL(m::LinearMixedModel{T}; full::Bool=false) where {T}
function sparseL(
m::LinearMixedModel{T};
fname::Symbol=first(fnames(m)),
full::Bool=false,
) where {T}
L, reterms = m.L, m.reterms
nt = length(reterms) + full
sblk = findfirst(isequal(fname), fnames(m))
if isnothing(sblk)
throw(ArgumentError("fname = $fname is not the name of a grouping factor"))
end
blks = sblk:(length(reterms) + full)
rowoffset, coloffset = 0, 0
val = (i=Int32[], j=Int32[], v=T[])
for i in 1:nt, j in 1:i
for i in blks, j in first(blks):i
Lblk = L[block(i, j)]
cblk = _coord(Lblk)
append!(val.i, cblk.i .+ Int32(rowoffset))
Expand Down