Skip to content

Commit

Permalink
update
Browse files Browse the repository at this point in the history
  • Loading branch information
PharmCat committed Jan 10, 2021
1 parent b979ca6 commit 9d3bcb0
Show file tree
Hide file tree
Showing 13 changed files with 96 additions and 38 deletions.
5 changes: 4 additions & 1 deletion src/gmat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@ function gmat_switch!(G, θ, covstr, i)
G
end
################################################################################
#=
function gmat_base_z!(mx, θ::Vector{T}, covstr) where T
q = sum(length.(covstr.block[1]))
for r = 1:length(covstr.random)
Expand All @@ -44,8 +45,9 @@ function gmat_base_z!(mx, θ::Vector{T}, covstr) where T
end
mx
end
=#
################################################################################

#=
function gmat_base_z2!(mx, θ::Vector{T}, covstr, block) where T
q = sum(length.(covstr.block[1]))
for r = 1:length(covstr.random)
Expand All @@ -62,6 +64,7 @@ function gmat_base_z2!(mx, θ::Vector{T}, covstr, block) where T
end
mx
end
=#
function gmat_base_z2!(mx, θ::Vector{T}, covstr, block, sblock) where T
q = sum(length.(covstr.block[1]))
for r = 1:length(covstr.random)
Expand Down
18 changes: 18 additions & 0 deletions src/linearalgebra.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
"""
A * B * A' + C
"""
#=
function mulαβαtc(A, B, C)
q = size(B, 1)
p = size(A, 1)
Expand All @@ -23,9 +24,11 @@ function mulαβαtc(A, B, C)
end
Symmetric(mx)
end
=#
"""
A * B * A'
"""
#=
function mulαβαt(A, B)
q = size(B, 1)
p = size(A, 1)
Expand Down Expand Up @@ -92,6 +95,7 @@ function mulαβαtc3(A, B, C, X)
mx[p+1:end, 1:p] = X'
mx
end
=#
function mulαβαt3(A, B, X)
q = size(B, 1)
p = size(A, 1)
Expand All @@ -118,6 +122,7 @@ end
A' * B * A -> +θ
A' * B * C -> +β
"""
#=
function mulθβinc!(θ, β, A::AbstractMatrix, B::AbstractMatrix, C::AbstractVector)
q = size(B, 1)
p = size(A, 2)
Expand All @@ -138,10 +143,12 @@ function mulθβinc!(θ, β, A::AbstractMatrix, B::AbstractMatrix, C::AbstractVe
end
θ, β
end
=#
#-------------------------------------------------------------------------------
"""
A' * B * A -> + θ
"""
#=
function mulαtβαinc!(θ, A, B)
q = size(B, 1)
p = size(A, 2)
Expand All @@ -156,6 +163,7 @@ function mulαtβαinc!(θ, A, B)
end
end
end
=#
"""
A * B * A' -> + θ
"""
Expand All @@ -179,6 +187,7 @@ end
"""
A' * B * A -> θ
"""
#=
function mulαtβα!(θ, A, B)
q = size(B, 1)
p = size(A, 2)
Expand All @@ -195,9 +204,11 @@ function mulαtβα!(θ, A, B)
end
θ
end
=#
"""
A * B * A -> θ
"""
#=
function mulαβα!(θ, A, B)
q = size(B, 1)
p = size(A, 2)
Expand All @@ -214,21 +225,26 @@ function mulαβα!(θ, A, B)
end
θ
end
=#
"""
tr(A * B)
"""
#=
function trmulαβ(A, B)
c = 0
@inbounds for n = 1:size(A,1), m = 1:size(B, 1)
c += A[n,m] * B[m, n]
end
c
end
=#
"""
tr(H * A' * B * A)
"""
#=
function trmulhαtβα(H, A, B)
end
=#
"""
(y - X * β)' * V * (y - X * β)
"""
Expand All @@ -252,6 +268,7 @@ end
"""
(y - X * β)
"""
#=
function mulr!(v::AbstractVector, y::AbstractVector, X::AbstractMatrix, β::AbstractVector)
fill!(v, zero(eltype(v)))
q = length(y)
Expand All @@ -276,6 +293,7 @@ function mulr(y::AbstractVector, X::AbstractMatrix, β::AbstractVector)
end
return v
end
=#
"""
A' * B -> + θ
"""
Expand Down
20 changes: 8 additions & 12 deletions src/reml.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
"""
-2 log Restricted Maximum Likelihood;
"""
#=
function reml_sweep(lmm, β, θ::Vector{T})::T where T <: Number
n = length(lmm.data.yv)
N = sum(length.(lmm.data.yv))
Expand All @@ -27,14 +28,16 @@ function reml_sweep(lmm, β, θ::Vector{T})::T where T <: Number
θ₂ = logdet(θ2m)
return -(θ₁ + θ₂ + θ₃ + c)
end

=#
"""
-2 log Restricted Maximum Likelihood; β calculation inside
"""
#=
function reml_sweep_β(lmm, f::Function, θ::Vector{T}) where T <: Number
f(θ)
reml_sweep_β(lmm, θ)
end
=#
"""
-2 log Restricted Maximum Likelihood; β calculation inside
"""
Expand Down Expand Up @@ -69,9 +72,7 @@ function reml_sweep_β(lmm::LMM{T2}, @nospecialize θ::Vector{T}) where T <: Num
catch
θ₁ += Inf
end

#println("logdet(V) $(i) : ", θ₁)

sweep!(Vp, 1:q)
V⁻¹[i] = Symmetric(utriaply!(x -> -x, V))
#-----------------------------------------------------------------------
Expand All @@ -93,7 +94,6 @@ function reml_sweep_β(lmm::LMM{T2}, @nospecialize θ::Vector{T}) where T <: Num
@inbounds θ₃ += mulθ₃(view(lmm.data.yv, lmm.data.block[i]), view(lmm.data.xv, lmm.data.block[i],:), β, V⁻¹[i])
#println("θ₃ : ", θ₃)
end

#logdetθ₂ = logdet(θ₂)
try
logdetθ₂ = logdet(θ₂)
Expand All @@ -102,7 +102,7 @@ function reml_sweep_β(lmm::LMM{T2}, @nospecialize θ::Vector{T}) where T <: Num
end
return θ₁ + logdetθ₂ + θ₃ + c, β, θ₂, θ₃
end

#=
function reml_sweep_β2(lmm::LMM{T2}, @nospecialize θ::Vector{T}) where T <: Number where T2 <: Number
n::Int = length(lmm.data.block)
N::Int = length(lmm.data.yv)
Expand Down Expand Up @@ -157,7 +157,7 @@ function reml_sweep_β2(lmm::LMM{T2}, @nospecialize θ::Vector{T}) where T <: Nu
return θ₁ + logdetθ₂ + θ₃ + c, β, θ₂, θ₃
end

=#
function reml_sweep_β3(lmm::LMM{T2}, @nospecialize θ::Vector{T}) where T <: Number where T2 <: Number
n::Int = length(lmm.data.block)
N::Int = length(lmm.data.yv)
Expand Down Expand Up @@ -187,11 +187,8 @@ function reml_sweep_β3(lmm::LMM{T2}, @nospecialize θ::Vector{T}) where T <: Nu
Vx .= view(lmm.data.xv, lmm.data.block[i],:)
gmat_base_z2!(V, θ, lmm.covstr, lmm.data.block[i], lmm.covstr.sblock[i])
rmat_basep_z2!(V, θ[lmm.covstr.tr[end]], lmm.covstr, lmm.data.block[i], lmm.covstr.sblock[i])

#θ₁ += logdet(V)

θ₁ += logdet(V)

sweep!(Vp, 1:q)
V⁻¹[i] = Symmetric(utriaply!(x -> -x, V))
#-----------------------------------------------------------------------
Expand All @@ -206,9 +203,7 @@ function reml_sweep_β3(lmm::LMM{T2}, @nospecialize θ::Vector{T}) where T <: Nu
#θ₃ += r' * V⁻¹[i] * r
@inbounds θ₃ += mulθ₃(view(lmm.data.yv, lmm.data.block[i]), view(lmm.data.xv, lmm.data.block[i],:), β, V⁻¹[i])
end

logdetθ₂ = logdet(θ₂)

return θ₁ + logdetθ₂ + θ₃ + c, β, θ₂, θ₃
end
#=
Expand Down Expand Up @@ -258,6 +253,7 @@ function reml_sweep_β_ub(lmm::LMM{T2}, θ::Vector{T})::Tuple{T, Vector{T}, Matr
end
=#
################################################################################
#=
function reml_inv_β(lmm::LMM{T2}, θ::Vector{T})::Tuple{T, Vector{T}, Matrix{T}} where T <: Number where T2 <: Number
n::Int = length(lmm.data.block)
N::Int = length(lmm.data.yv)
Expand Down Expand Up @@ -308,5 +304,5 @@ function reml_inv_β(lmm::LMM{T2}, θ::Vector{T})::Tuple{T, Vector{T}, Matrix{T}
end
return θ₁ + logdetθ₂ + θ₃ + c, β, θ₂
end

=#
################################################################################
7 changes: 6 additions & 1 deletion src/reml_gh_deprecated.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
"""
2 log Restricted Maximum Likelihood gradient vector
"""
#=
function reml_grad(yv, Zv, p, Xv, θvec, β)
n = length(yv)
G = gmat(θvec[3:5])
Expand Down Expand Up @@ -40,10 +41,11 @@ function reml_grad(yv, Zv, p, Xv, θvec, β)
end
return - (θ1 .+ θ2 .+ θ3)
end

=#
"""
2 log Restricted Maximum Likelihood hessian matrix
"""
#=
function reml_hessian(yv, Zv, p, Xv, θvec, β)
n = length(yv)
G = gmat(θvec[3:5])
Expand Down Expand Up @@ -91,7 +93,9 @@ function reml_hessian(yv, Zv, p, Xv, θvec, β)
end
return - (θ1 .+ θ2 .+ θ3)
end
=#
################################################################################
#=
function reml_grad2(yv, Zv, p, Xv, θvec, β)
n = length(yv)
for i = 1:n
Expand All @@ -102,3 +106,4 @@ function reml_grad3(yv, Zv, p, Xv, θvec, β)
n = length(yv)
covmat_grad2(vmatvec, Zv, θvec)
end
=#
4 changes: 4 additions & 0 deletions src/rmat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@ function rmat_basep!(mx, θ::AbstractVector{T}, zrv, covstr::CovStructure{T2}) w
end
end
################################################################################
#=
function rmat_basep_z!(mx, θ::AbstractVector{T}, zrv, covstr) where T
for i = 1:length(covstr.block[end])
if covstr.repeated.covtype.s == :SI
Expand All @@ -41,7 +42,9 @@ function rmat_basep_z!(mx, θ::AbstractVector{T}, zrv, covstr) where T
end
mx
end
=#
################################################################################
#=
function rmat_basep_z2!(mx, θ::AbstractVector{T}, covstr, block) where T
subjblock = view(covstr.subjz[end], block, :)
zblock = view(covstr.rz, block, :)
Expand All @@ -53,6 +56,7 @@ function rmat_basep_z2!(mx, θ::AbstractVector{T}, covstr, block) where T
end
mx
end
=#
function rmat_basep_z2!(mx, θ::AbstractVector{T}, covstr, block, sblock) where T
zblock = view(covstr.rz, block, :)
for i = 1:length(sblock[end])
Expand Down
5 changes: 3 additions & 2 deletions src/statsbase.jl
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,5 @@



#=
StatsBase.coef(model::LMM) = error("coef is not defined for $(typeof(model)).")
StatsBase.coefnames(model::LMM) = error("coefnames is not defined for $(typeof(model)).")
Expand Down Expand Up @@ -122,3 +121,5 @@ StatsBase.predict(model::LMM) = error("predict is not defined for $(typeof(model
StatsBase.predict!(model::LMM) = error("predict! is not defined for $(typeof(model)).")
StatsBase.dof_residual(model::LMM) = error("dof_residual is not defined for $(typeof(model)).")
=#
15 changes: 8 additions & 7 deletions src/utils.jl
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
"""
Make X, Z matrices and vector y for each subject;
"""
#=
function subjblocks(df, sbj::Symbol, x::Matrix{T}, z::Matrix{T}, y::Vector{T}, rz::Matrix{T}) where T
u = unique(df[!, sbj])
xa = Vector{Matrix{T}}(undef, length(u))
Expand Down Expand Up @@ -44,6 +45,7 @@ function subjblocks(df, sbj)
r[1] = collect(1:size(df, 1))
r
end
=#
"""
Intersect dataframe.
"""
Expand Down Expand Up @@ -73,7 +75,7 @@ function intersectdf(df, s)::Vector
end
res
end

#=
function intersectsubj(covstr)
a = Vector{Vector{Symbol}}(undef, length(covstr.random)+1)
eq = true
Expand All @@ -89,6 +91,7 @@ function intersectsubj(covstr)
end
intersect(a...), eq
end
=#
function intersectsubj(random, repeated)
a = Vector{Vector{Symbol}}(undef, length(random)+1)
eq = true
Expand All @@ -104,12 +107,12 @@ function intersectsubj(random, repeated)
end
intersect(a...), eq
end

#=
function diffsubj!(a, subj)
push!(a, subj)
symdiff(a...)
end

=#
"""
Variance estimate via OLS and QR decomposition.
"""
Expand Down Expand Up @@ -208,7 +211,7 @@ function varlinkrvecapply2!(v, p; varlinkf = :exp, rholinkf = :sigm)
end

################################################################################

#=
function vmatr(lmm, i)
θ = lmm.result.theta
G = gmat_base(θ, lmm.covstr)
Expand All @@ -225,7 +228,7 @@ function gmatr(lmm, i)
θ = lmm.result.theta
gmat_base(θ, lmm.covstr)
end

=#
################################################################################

function m2logreml(lmm)
Expand All @@ -236,8 +239,6 @@ function logreml(lmm)
end
################################################################################



function optim_callback(os)
false
end
Loading

2 comments on commit 9d3bcb0

@PharmCat
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/27692

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.1.11 -m "<description of version>" 9d3bcb0d617a1b3a1bf01394e2c30ac6369c4b22
git push origin v0.1.11

Please sign in to comment.