Skip to content

Commit 27b9b75

Browse files
mmkim1210mmkim1210
committed
v0.3.4
1 parent 0999550 commit 27b9b75

7 files changed

Lines changed: 59 additions & 40 deletions

File tree

‎Project.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "MultiResponseVarianceComponentModels"
22
uuid = "95bd6349-d788-4b5f-978d-e543c2eb4d7d"
33
authors = ["Hua Zhou <huazhou@ucla.edu> and Minsoo Kim <mmkim1210@gmail.com>"]
4-
version = "0.3.3"
4+
version = "0.3.4"
55

66
[deps]
77
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

‎README.md‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
<p align="center"><img width="100%" style="border-radius: 5px;" src="docs/src/assets/logo.svg"></p>
22

33
# MultiResponseVarianceComponentModels
4-
[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](http://hua-zhou.github.io/MultiResponseVarianceComponentModels.jl/dev)
4+
[![Latest](https://img.shields.io/badge/docs-latest-blue.svg)](http://hua-zhou.github.io/MultiResponseVarianceComponentModels.jl/dev)
55
[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](http://hua-zhou.github.io/MultiResponseVarianceComponentModels.jl/stable)
66
[![CI](https://github.com/Hua-Zhou/MultiResponseVarianceComponentModels.jl/workflows/CI/badge.svg)](https://github.com/Hua-Zhou/MultiResponseVarianceComponentModels.jl/actions)
77
[![Codecov](https://codecov.io/gh/Hua-Zhou/MultiResponseVarianceComponentModels.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/Hua-Zhou/MultiResponseVarianceComponentModels.jl)
@@ -16,7 +16,7 @@ julia> ]
1616
pkg> add https://github.com/Hua-Zhou/MultiResponseVarianceComponentModels.jl.git
1717
```
1818
## Documentation
19-
[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](http://hua-zhou.github.io/MultiResponseVarianceComponentModels.jl/stable)
19+
[![Latest](https://img.shields.io/badge/docs-latest-blue.svg)](http://hua-zhou.github.io/MultiResponseVarianceComponentModels.jl/dev)
2020

2121
## Examples
2222
```julia

‎src/MultiResponseVarianceComponentModels.jl‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
"""
22
__MRVCModels__ stands for __M__ultivariate __R__esponse __V__ariance __C__omponents
3-
linear mixed __Models__. `MRVCModels.jl` permits maximum likelihood (ML) or residual
3+
linear mixed __Models__. MRVCModels.jl permits maximum likelihood (ML) or residual
44
maximum likelihood (REML) estimation and inference.
55
"""
66
module MultiResponseVarianceComponentModels
@@ -139,7 +139,7 @@ reml::Bool pursue REML estimation instead of ML estimation; default false
139139
When there are two variance components, computation can be reduced by avoiding large matrix
140140
inversion per iteration, which is achieved with `MRTVCModel` instance. __MRTVCModels__
141141
stands for __M__ultivariate __R__esponse __T__wo __V__ariance __C__omponents
142-
linear mixed __Models__. `MRVCModel` is more general and is not limited to two variance
142+
linear mixed __Models__. `MRVCModel` is more general, since it is not limited to two variance
143143
components case. For `MRTVCModel`, the number of variance components must be two.
144144
"""
145145
function MRVCModel(

‎src/fit.jl‎

Lines changed: 8 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -17,8 +17,8 @@ log::Bool record iterate history or not; default false
1717
1818
# Extended help
1919
MM algorithm is provably faster than EM algorithm in this setting, so recommend trying
20-
MM algorithm first, which is by default, and switching to EM algorithm if there are
21-
convergence issues.
20+
MM algorithm first, which is the default algorithm, and switching to EM algorithm only if
21+
there are convergence issues.
2222
"""
2323
function fit!(
2424
model :: MRVCModel{T};
@@ -75,31 +75,13 @@ function fit!(
7575
else
7676
throw("Cannot recognize initialization method $init")
7777
end
78-
if model.ymissing
79-
# update conditional variance
80-
C = model.storage_n_miss_n_miss_1
81-
PΩPt = model.storage_nd_nd_miss
82-
PΩPt .= @view model.Ω[model.P, model.P]
83-
nd = size(model.Ω, 1)
84-
n_obs = nd - model.n_miss
85-
sweep!(PΩPt, 1:n_obs)
86-
copytri!(PΩPt, 'U')
87-
copy!(model.storage_n_miss_n_obs_1,
88-
view(PΩPt, (n_obs + 1):nd, 1:n_obs)) # for conditional mean
89-
copy!(model.storage_n_miss_n_miss_1,
90-
view(PΩPt, (n_obs + 1):nd, (n_obs + 1):nd)) # conditional variance
91-
logl = NaN
92-
else
93-
logl = loglikelihood!(model)
94-
end
78+
model.ymissing ? logl = loglikelihood_miss!(model) : logl = loglikelihood!(model)
9579
toc = time()
96-
if !model.ymissing
97-
verbose && println("iter = 0, logl = $logl")
98-
IterativeSolvers.nextiter!(history)
99-
push!(history, :iter , 0)
100-
push!(history, :logl , logl)
101-
push!(history, :itertime, toc - tic)
102-
end
80+
verbose && println("iter = 0, logl = $logl")
81+
IterativeSolvers.nextiter!(history)
82+
push!(history, :iter , 0)
83+
push!(history, :logl , logl)
84+
push!(history, :itertime, toc - tic)
10385
# MM loop
10486
for iter in 1:maxiter
10587
IterativeSolvers.nextiter!(history)

‎test/eigen_test.jl‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -92,8 +92,8 @@ model2 = MRTVCModel(Y, X, V)
9292
model = MRVCModel(Y, X, V)
9393

9494
@testset "fit! two component by MLE with EM" begin
95-
MRVCModels.fit!(model2, algo = :EM, maxiter = 500)
96-
MRVCModels.fit!(model, algo = :EM, maxiter = 500)
95+
MRVCModels.fit!(model2, algo = :EM, verbose = false, maxiter = 500)
96+
MRVCModels.fit!(model, algo = :EM, verbose = false, maxiter = 500)
9797
println("||B̂_MRTVCModel - B̂_MRVCModel|| = $(norm(model2.B - model.B))")
9898
for k in 1:m
9999
println("||Σ̂[$k]_MRTVCModel - Σ̂[$k]_MRVCModel|| = $(norm(model2.Σ[k] - model.Σ[k]))")
@@ -117,8 +117,8 @@ model2 = MRTVCModel(Y, X, V, reml = true)
117117
model = MRVCModel(Y, X, V, reml = true)
118118

119119
@testset "fit! two component by REML with EM" begin
120-
MRVCModels.fit!(model2, algo = :EM, maxiter = 500)
121-
MRVCModels.fit!(model, algo = :EM, maxiter = 500)
120+
MRVCModels.fit!(model2, algo = :EM, verbose = false, maxiter = 500)
121+
MRVCModels.fit!(model, algo = :EM, verbose = false, maxiter = 500)
122122
println("||B̂_MRTVCModel - B̂_MRVCModel|| = $(norm(model2.B_reml - model.B_reml))")
123123
for k in 1:m
124124
println("||Σ̂[$k]_MRTVCModel - Σ̂[$k]_MRVCModel|| = $(norm(model2.Σ[k] - model.Σ[k]))")

‎test/fit_test.jl‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -46,7 +46,7 @@ end
4646
model = MRVCModel(Y, X, V)
4747

4848
@testset "fit! by MLE with MM" begin
49-
@time MRVCModels.fit!(model, algo = :MM, maxiter = 150)
49+
@time MRVCModels.fit!(model, algo = :MM, verbose = false, maxiter = 150)
5050
println("||B_true - B̂|| = $(norm(B_true - model.B))")
5151
for k in 1:m
5252
println("||Σ_true[$k] - Σ̂[$k]|| = $(norm(Σ_true[k] - model.Σ[k]))")
@@ -57,7 +57,7 @@ model = MRVCModel(Y, X, V)
5757
end
5858

5959
@testset "fit! by MLE with EM" begin
60-
@time MRVCModels.fit!(model, algo = :EM, maxiter = 150)
60+
@time MRVCModels.fit!(model, algo = :EM, verbose = false, maxiter = 150)
6161
println("||B_true - B̂|| = $(norm(B_true - model.B))")
6262
for k in 1:m
6363
println("||Σ_true[$k] - Σ̂[$k]|| = $(norm(Σ_true[k] - model.Σ[k]))")
@@ -67,7 +67,7 @@ end
6767
model = MRVCModel(Y, X, V; reml = true)
6868

6969
@testset "fit! by REML with MM" begin
70-
@time MRVCModels.fit!(model, algo = :MM, maxiter = 150)
70+
@time MRVCModels.fit!(model, algo = :MM, verbose = false, maxiter = 150)
7171
println("||B_true - B̂|| = $(norm(B_true - model.B_reml))")
7272
for k in 1:m
7373
println("||Σ_true[$k] - Σ̂[$k]|| = $(norm(Σ_true[k] - model.Σ[k]))")
@@ -78,7 +78,7 @@ model = MRVCModel(Y, X, V; reml = true)
7878
end
7979

8080
@testset "fit! by REML with EM" begin
81-
@time MRVCModels.fit!(model, algo = :EM, maxiter = 150)
81+
@time MRVCModels.fit!(model, algo = :EM, verbose = false, maxiter = 150)
8282
println("||B_true - B̂|| = $(norm(B_true - model.B_reml))")
8383
for k in 1:m
8484
println("||Σ_true[$k] - Σ̂[$k]|| = $(norm(Σ_true[k] - model.Σ[k]))")

‎test/missing_test.jl‎

Lines changed: 38 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,39 @@
11
module MissingTest
22

33
using MultiResponseVarianceComponentModels
4-
using Test
4+
using BenchmarkTools, LinearAlgebra, Profile, Random, StatsBase, Test
5+
6+
const MRVCModels = MultiResponseVarianceComponentModels
7+
Random.seed!(789)
8+
9+
n, p, d, m = 855, 3, 4, 3
10+
X = [ones(n) randn(n, p - 1)] # design matrix including intercept
11+
# V[1] is an AR1(ρ) matrix, with entries ρ^|i-j|
12+
# V[2] has entries i * (n - j + 1) for j ≥ i, then scaled to be a correlation matrix
13+
# V[3] is identity
14+
ρ = 0.5
15+
V = Vector{Matrix{Float64}}(undef, m)
16+
V[1] = [ρ^abs(i - j) for i in 1:n, j in 1:n]
17+
V[2] = [j ≥ i ? i * (n - j + 1) : j * (n - i + 1) for i in 1:n, j in 1:n]
18+
StatsBase.cov2cor!(V[2], [sqrt(V[2][i, i]) for i in 1:n])
19+
V[3] = Matrix(UniformScaling(1.0), n, n)
20+
# true parameter values
21+
B_true = 2 * rand(p, d) # uniform on [0, 2]
22+
Σ_true = [
23+
Matrix(UniformScaling(0.2), d, d),
24+
Matrix(UniformScaling(0.2), d, d),
25+
Matrix(UniformScaling(0.6), d, d)
26+
]
27+
Ω_true = zeros(n * d, n * d)
28+
for k in 1:m
29+
Ω_true .+= kron(Σ_true[k], V[k])
30+
end
31+
y = vec(X * B_true) + cholesky(Symmetric(Ω_true)).L * randn(n * d)
32+
Y = reshape(y, n, d)
33+
34+
Y_miss = Matrix{Union{eltype(Y), Missing}}(missing, size(Y))
35+
copy!(Y_miss, Y)
36+
Y_miss[rand(1:length(Y_miss), n)] .= missing
537

638
@testset "permute" begin
739
Y = reshape(1:16, 4, 4)
@@ -15,4 +47,9 @@ using Test
1547
@test Y_imputed == [3.0 7.5 9.0 13.0; 2.0 7.5 10.0 14.0; 3.0 7.0 11.0 13.5; 4.0 8.0 12.0 13.5]
1648
end
1749

50+
@testset "fit! missing response with MM" begin
51+
model = MRVCModel(Y_miss, X, V; se = false)
52+
# @timev MultiResponseVarianceComponentModels.fit!(model)
53+
end
54+
1855
end

0 commit comments

Comments
 (0)