.github/workflows/CI.yml

@@ -18,7 +18,6 @@ jobs:
           - "1.9"
           - "1.10"
           - "1.11"
-          - "nightly"
         os:
           - ubuntu-latest
           - macos-latest

LICENSE

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2021 Hua Zhou <huazhou@ucla.edu> and contributors
+Copyright (c) 2021 Hua Zhou <huazhou@ucla.edu> and Minsoo Kim <mmkim1210@gmail.com> and contributors
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

README.md

@@ -1,10 +1,14 @@
 <p align="center"><img width="100%" style="border-radius: 5px;" src="docs/src/assets/logo.svg"></p>
 
 # MultiResponseVarianceComponentModels
-[![Latest](https://img.shields.io/badge/docs-latest-blue.svg)](http://hua-zhou.github.io/MultiResponseVarianceComponentModels.jl/dev)
-[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](http://hua-zhou.github.io/MultiResponseVarianceComponentModels.jl/stable)
+[![Latest](https://img.shields.io/badge/docs-latest-blueviolet.svg)](http://hua-zhou.github.io/MultiResponseVarianceComponentModels.jl/dev)
+[![Stable](https://img.shields.io/badge/docs-stable-mediumorchid.svg)](http://hua-zhou.github.io/MultiResponseVarianceComponentModels.jl/stable)
+[![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/JuliaDiff/BlueStyle)
+[![License](https://img.shields.io/badge/license-MIT-green.svg)](https://github.com/Hua-Zhou/MultiResponseVarianceComponentModels.jl/blob/master/LICENSE)
 [![CI](https://github.com/Hua-Zhou/MultiResponseVarianceComponentModels.jl/workflows/CI/badge.svg)](https://github.com/Hua-Zhou/MultiResponseVarianceComponentModels.jl/actions)
 [![Codecov](https://codecov.io/gh/Hua-Zhou/MultiResponseVarianceComponentModels.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/Hua-Zhou/MultiResponseVarianceComponentModels.jl)
+![Repo Size](https://img.shields.io/github/repo-size/hua-zhou/MultiResponseVarianceComponentModels.jl)
+[![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
 
 MultiResponseVarianceComponentModels.jl is a <a href="https://julialang.org"><img src="https://julialang.org/assets/infra/julia.ico" width="10em"> Julia </a>package that allows fitting and testing multivariate response variance components linear mixed models of form 
 
@@ -16,7 +20,7 @@ julia> ]
 pkg> add https://github.com/Hua-Zhou/MultiResponseVarianceComponentModels.jl.git
 ```
 ## Documentation
-[![Latest](https://img.shields.io/badge/docs-latest-blue.svg)](http://hua-zhou.github.io/MultiResponseVarianceComponentModels.jl/dev)
+[![Latest](https://img.shields.io/badge/docs-latest-blueviolet.svg)](http://hua-zhou.github.io/MultiResponseVarianceComponentModels.jl/dev)
 
 ## Examples
 ```julia
@@ -25,45 +29,33 @@ using MultiResponseVarianceComponentModels, LinearAlgebra, Random
 # simulation
 begin
     Random.seed!(1234)
-    n = 1_000  # n of observations
-    d = 4      # n of responses
-    p = 10     # n of covariates
-    m = 3      # n of variance components
+    n = 1_000  # number of observations
+    d = 4      # number of responses
+    p = 10     # number of covariates
+    m = 3      # number of variance components
     X = rand(n, p)
-    B = rand(p, d) 
+    B = rand(p, d) # fixed effects
     V = [zeros(n, n) for _ in 1:m] # kernel matrices
     Σ = [zeros(d, d) for _ in 1:m] # variance components
+    Ω = zeros(n * d, n * d) # overall nd-by-nd covariance matrix Ω
     for i in 1:m
         Vi = randn(n, n)
-        copy!(V[i], Vi' * Vi)
+        mul!(V[i], transpose(Vi), Vi)
         Σi = randn(d, d)
-        copy!(Σ[i], Σi' * Σi)
-    end
-    Ω = zeros(n * d, n * d) # overall nd-by-nd covariance matrix Ω
-    for i = 1:m
-        Ω += kron(Σ[i], V[i])
+        mul!(Σ[i], transpose(Σi), Σi)
+        kron_axpy!(Σ[i], V[i], Ω) # Ω = Σ[1]⊗V[1] + ... + Σ[m]⊗V[m]
     end
     Ωchol = cholesky(Ω)
     Y = X * B + reshape(Ωchol.L * randn(n * d), n, d)
 end
 
 # maximum likelihood estimation
 model = MRVCModel(Y, X, V)
-@timev fit!(model) # ~ 30 seconds
+fit!(model)
 
 # residual maximum likelihood estimation
 model = MRVCModel(Y, X, V; reml = true)
-@timev fit!(model) # ~ 30 seconds
-
-# variance components estimates
-model.Σ
-# comparison of true values and estimates
-reduce(hcat, [hcat(vech(Σ[i]), vech(model.Σ[i])) for i in 1:m])
-# sampling variance by inverse of Fisher information matrix
-model.Σcov
-diag(model.Σcov) # m * (binomial(d, 2) + d) parameters
-# log-likelihood
-model.logl
+fit!(model)
 ```
 ## References
 - <u>H. Zhou, L. Hu, J. Zhou, and K. Lange: **MM algorithms for variance components models** (2019) ([link](https://doi.org/10.1080/10618600.2018.1529601))</u>
@@ -72,4 +64,9 @@ model.logl
 ## See also
 - J. Kim, J. Shen, A. Wang, D.V. Mehrotra, S. Ko, J.J. Zhou, and H. Zhou: **VCSEL: Prioritizing SNP-set by penalized variance component selection** (2021) ([link](http://doi.org/10.1214/21-aoas1491))
 - L. Hu, W. Lu, J. Zhou, and H. Zhou: **MM algorithms for variance component estimation and selection in logistic linear mixed models** (2019) ([link](http://doi.org/10.5705/ss.202017.0220))
-- J.J. Zhou, T. Hu, D. Qiao, M.H. Cho, and H. Zhou: **Boosting gene mapping power and efficiency with efficient exact variance component tests of single nucleotide polymorphism sets** (2016) ([link](http://doi.org/10.1534/genetics.116.190454))
+- J.J. Zhou, T. Hu, D. Qiao, M.H. Cho, and H. Zhou: **Boosting gene mapping power and efficiency with efficient exact variance component tests of single nucleotide polymorphism sets** (2016) ([link](http://doi.org/10.1534/genetics.116.190454))
+
+## Contributors
+<a href="https://github.com/Hua-Zhou/MultiResponseVarianceComponentModels.jl/graphs/contributors">
+    <img src="https://contrib.rocks/image?repo=Hua-Zhou/MultiResponseVarianceComponentModels.jl" />
+</a>

docs/src/advanced.md

@@ -35,14 +35,14 @@ Standard errors for our estimates are calculated using the Fisher information ma
 where ``\text{vech}\ \boldsymbol{\Gamma}_i`` creates an ``\frac{d(d+1)}{2} \times 1`` vector from ``\boldsymbol{\Gamma}_i`` by stacking its lower triangular part and ``\boldsymbol{U}_i = (\boldsymbol{I}_d \otimes \boldsymbol{K}_{nd} \otimes \boldsymbol{I}_n) (\boldsymbol{I}_{d^2} \otimes \text{vec}\ \boldsymbol{V}_i) \boldsymbol{D}_{d}``. Here, ``\boldsymbol{K}_{nd}`` is the ``nd \times nd`` [commutation matrix](https://en.wikipedia.org/wiki/Commutation_matrix) and ``\boldsymbol{D}_{d}`` the ``d^2 \times \frac{d(d+1)}{2}`` [duplication matrix](https://en.wikipedia.org/wiki/Duplication_and_elimination_matrices).
 
 # Special case: missing response
-In the setting of missing response, the adjusted MM updates in each interation are
+In the setting of missing response, we let ``\boldsymbol{P}`` be the ``nd \times nd`` permutation matrix such that ``\boldsymbol{P} \cdot \text{vec}\ \boldsymbol{Y} = \begin{bmatrix} \boldsymbol{y}_{\text{obs}} \\ \boldsymbol{y}_{\text{mis}} \end{bmatrix}``, where ``\boldsymbol{y}_{\text{obs}}`` and ``\boldsymbol{y}_{\text{mis}}`` are vectors of observed and missing response values, respectively, in column-major order. If we also let ``\boldsymbol{P} \cdot\text{vec}(\boldsymbol{X}\boldsymbol{B}) = \begin{bmatrix} \boldsymbol{\mu}_{1} \\ \boldsymbol{\mu}_{2} \end{bmatrix}`` and ``\boldsymbol{P} \boldsymbol{\Omega} \boldsymbol{P}^T = \begin{bmatrix} \boldsymbol{\Omega}_{11} & \boldsymbol{\Omega}_{12} \\ \boldsymbol{\Omega}_{21} & \boldsymbol{\Omega}_{22} \end{bmatrix}`` such that conditional mean and variance of ``\boldsymbol{y}_{\text{mis}}`` are ``\boldsymbol{\mu}_2 + \boldsymbol{\Omega}_{21}\boldsymbol{\Omega}_{11}^{-1}(\boldsymbol{y}_{\text{obs}}-\boldsymbol{\mu}_1)`` and ``\boldsymbol{\Omega}_{22} - \boldsymbol{\Omega}_{21}\boldsymbol{\Omega}_{11}^{-1}\boldsymbol{\Omega}_{12}``, respectively, then the adjusted MM updates in each interation become
 ```math
 \begin{aligned}
 \text{vec}\ \boldsymbol{B}^{(t)} &= [(\boldsymbol{I}_d \otimes \boldsymbol{X}^T) \boldsymbol{\Omega}^{-(t)} (\boldsymbol{I}_d \otimes \boldsymbol{X})]^{-1} (\boldsymbol{I}_d \otimes \boldsymbol{X}^T) \boldsymbol{\Omega}^{-(t)} \text{vec}\ \boldsymbol{Z}^{(t)} \\
 \boldsymbol{\Gamma}_i^{(t + 1)} &= \boldsymbol{L}_i^{-(t)T}[\boldsymbol{L}_i^{(t)T}\boldsymbol{\Gamma}_i^{(t)}(\boldsymbol{R}^{*(t)T}\boldsymbol{V}_i\boldsymbol{R}^{*(t)} + \boldsymbol{M}_i^{*(t)})\boldsymbol{\Gamma}_i^{(t)}\boldsymbol{L}_i^{(t)}]^{1/2} \boldsymbol{L}_i^{-(t)},
 \end{aligned}
 ```
-where ``\boldsymbol{Z}^{(t)}`` is the completed response matrix from conditional mean, ``\boldsymbol{M}_i^{*(t)} = (\boldsymbol{I}_d \otimes \boldsymbol{1}_n)^T [(\boldsymbol{1}_d \boldsymbol{1}_d^T \otimes \boldsymbol{V}_i) \odot (\boldsymbol{\Omega}^{-(t)} \boldsymbol{P}^T \boldsymbol{C}^{(t)}\boldsymbol{P}\boldsymbol{\Omega}^{-(t)})] (\boldsymbol{I}_d \otimes \boldsymbol{1}_n)``, and ``\boldsymbol{R}^{*(t)}`` is the ``n \times d`` matrix such that ``\text{vec}\ \boldsymbol{R}^{*(t)} = \boldsymbol{\Omega}^{-(t)} \text{vec}(\boldsymbol{Z}^{(t)} - \boldsymbol{X} \boldsymbol{B}^{(t)})``. Additionally, ``\boldsymbol{P}`` is the ``nd \times nd`` permutation matrix such that ``\boldsymbol{P} \cdot \text{vec}\ \boldsymbol{Y} = \begin{bmatrix} \boldsymbol{y}_{\text{obs}} \\ \boldsymbol{y}_{\text{mis}} \end{bmatrix}``, where ``\boldsymbol{y}_{\text{obs}}`` and ``\boldsymbol{y}_{\text{mis}}`` are vectors of observed and missing response values, respectively, in column-major order, and the block matrix ``\boldsymbol{C}^{(t)}`` is ``\boldsymbol{0}`` except for a lower-right block consisting of conditional variance. As seen, the MM updates are of similar form to the non-missing response case.
+where ``\boldsymbol{Z}^{(t)}`` is the ``n \times d`` completed response matrix from conditional mean such that ``\text{vec}\ \boldsymbol{Z}^{(t)} = \boldsymbol{P}^T \cdot \begin{bmatrix} \boldsymbol{y}_{\text{obs}} \\ \boldsymbol{\mu}_2^{(t)} + \boldsymbol{\Omega}_{21}^{(t)}\boldsymbol{\Omega}_{11}^{-(t)}(\boldsymbol{y}_{\text{obs}}-\boldsymbol{\mu}_1^{(t)}) \end{bmatrix}``, ``\boldsymbol{R}^{*(t)}`` is the ``n \times d`` matrix such that ``\text{vec}\ \boldsymbol{R}^{*(t)} = \boldsymbol{\Omega}^{-(t)} \text{vec}(\boldsymbol{Z}^{(t)} - \boldsymbol{X} \boldsymbol{B}^{(t)})``, and ``\boldsymbol{M}_i^{*(t)} = (\boldsymbol{I}_d \otimes \boldsymbol{1}_n)^T [(\boldsymbol{1}_d \boldsymbol{1}_d^T \otimes \boldsymbol{V}_i) \odot (\boldsymbol{\Omega}^{-(t)} \boldsymbol{P}^T \boldsymbol{C}^{(t)}\boldsymbol{P}\boldsymbol{\Omega}^{-(t)})] (\boldsymbol{I}_d \otimes \boldsymbol{1}_n)``. Here, ``\boldsymbol{C}^{(t)}`` is the block matrix that is ``\boldsymbol{0}`` everywhere except for the lower-right block consisting of conditional variance ``\boldsymbol{\Omega}_{22}^{(t)} - \boldsymbol{\Omega}_{21}^{(t)}\boldsymbol{\Omega}_{11}^{-(t)}\boldsymbol{\Omega}_{12}^{(t)}``. As seen, the MM updates are of similar form to the non-missing response case.
 
 # Special case: ``m = 2``
 When there are ``m = 2`` variance components such that ``\boldsymbol{\Omega} = \boldsymbol{\Gamma}_1 \otimes \boldsymbol{V}_1 + \boldsymbol{\Gamma}_2 \otimes \boldsymbol{V}_2``, repeated inversion of the ``nd \times nd`` matrix ``\boldsymbol{\Omega}`` per iteration can be avoided and reduced to one ``d \times d`` generalized eigen-decomposition per iteration. Without loss of generality, if we assume ``\boldsymbol{V}_2`` to be positive definite, the generalized eigen-decomposition of the matrix pair ``(\boldsymbol{V}_1, \boldsymbol{V}_2)`` yields generalized eigenvalues ``\boldsymbol{d} = (d_1, \dots, d_n)^T`` and generalized eigenvectors ``\boldsymbol{U}`` such that ``\boldsymbol{U}^T \boldsymbol{V}_1 \boldsymbol{U} = \boldsymbol{D} = \text{diag}(\boldsymbol{d})`` and ``\boldsymbol{U}^T \boldsymbol{V}_2 \boldsymbol{U} = \boldsymbol{I}_n``. Similarly, if we let the generalized eigen-decomposition of ``(\boldsymbol{\Gamma}_1^{(t)}, \boldsymbol{\Gamma}_2^{(t)})`` be ``(\boldsymbol{\Lambda}^{(t)}, \boldsymbol{\Phi}^{(t)})`` such that ``\boldsymbol{\Phi}^{(t)T} \boldsymbol{\Gamma}_1^{(t)} \boldsymbol{\Phi}^{(t)} = \boldsymbol{\Lambda}^{(t)} = \text{diag}(\boldsymbol{\lambda^{(t)}})`` and ``\boldsymbol{\Phi}^{(t)T} \boldsymbol{\Gamma}_2^{(t)} \boldsymbol{\Phi}^{(t)} = \boldsymbol{I}_d``, then the MM updates in each iteration become
@@ -78,4 +78,8 @@ where ``\boldsymbol{W}_{ij}`` is the ``d \times d`` matrix that has entries
 (\boldsymbol{W}_{22})_{kl} &= \text{tr}((\lambda_k \boldsymbol{D} + \boldsymbol{I}_n)^{-1}(\lambda_l \boldsymbol{D} + \boldsymbol{I}_n)^{-1})
 \end{aligned}
 ```
-for ``1 \leq k, l \leq d``.
+for ``1 \leq k, l \leq d``.
+
+# References
+- H. Zhou, L. Hu, J. Zhou, and K. Lange: **MM algorithms for variance components models** (2019) ([link](https://doi.org/10.1080/10618600.2018.1529601))
+- M. Kim: **Gene regulation in the human brain and the biological mechanisms underlying psychiatric disorders** (2022) ([link](https://escholarship.org/uc/item/9v08q5f7))

docs/src/examples.md

@@ -4,12 +4,12 @@
 ```@repl 1
 using MultiResponseVarianceComponentModels, LinearAlgebra, Random
 Random.seed!(6789)
-n = 1_000;  # n of observations
-d = 4;      # n of responses
-p = 10;     # n of covariates
-m = 5;      # n of variance components
+n = 1_000;  # number of observations
+d = 4;      # number of responses
+p = 10;     # number of covariates
+m = 5;      # number of variance components
 X = rand(n, p);
-B = rand(p, d)
+B = rand(p, d);
 V = [zeros(n, n) for _ in 1:m]; # kernel matrices
 Σ = [zeros(d, d) for _ in 1:m]; # variance components
 for i in 1:m
@@ -20,7 +20,7 @@ for i in 1:m
 end
 Ω = zeros(n * d, n * d); # overall nd-by-nd covariance matrix Ω
 for i = 1:m
-    Ω += kron(Σ[i], V[i])
+    kron_axpy!(Σ[i], V[i], Ω) # Ω = Σ[1]⊗V[1] + ... + Σ[m]⊗V[m]
 end
 Ωchol = cholesky(Ω);
 Y = X * B + reshape(Ωchol.L * randn(n * d), n, d)
@@ -40,8 +40,8 @@ Then variance components and mean effects estimates can be accessed through
 ```@repl 1
 model.Σ
 model.B
-hcat(vec(B), vec(model.B))
-reduce(hcat, [hcat(vech(Σ[i]), vech(model.Σ[i])) for i in 1:m])
+hcat(vec(B), vec(model.B)) # comparison of true values and estimates
+reduce(hcat, [hcat(vech(Σ[i]), vech(model.Σ[i])) for i in 1:m]) # comparison of true values and estimates
 ```
 
 # Standard errors
@@ -52,7 +52,7 @@ model.Bcov
 ```
 Corresponding standard error of these estimates are
 ```@repl 1
-sqrt.(diag(model.Σcov))
+sqrt.(diag(model.Σcov)) # m * (binomial(d, 2) + d) parameters
 sqrt.(diag(model.Bcov))
 ```
 
@@ -66,11 +66,11 @@ model = MRVCModel(Y, X, V; reml = true)
 Variance components and mean effects estimates and their standard errors can be accessed through:
 ```@repl 1
 model.Σ
-model.B_reml
+model.B_reml # note model.B_reml not model.B
 hcat(vec(B), vec(model.B_reml))
 reduce(hcat, [hcat(vech(Σ[i]), vech(model.Σ[i])) for i in 1:m])
 model.Σcov
-model.Bcov_reml
+model.Bcov_reml # note model.Bcov_reml not model.Bcov
 sqrt.(diag(model.Σcov))
 sqrt.(diag(model.Bcov_reml))
 ```
@@ -101,16 +101,21 @@ function simulate(n, d, p, m)
     B = rand(p, d)
     V = [zeros(n, n) for _ in 1:m]
     Σ = [zeros(d, d) for _ in 1:m]
-    Ω = zeros(n * d, n * d)
     for i in 1:m
         Vi = randn(n, n)
         mul!(V[i], transpose(Vi), Vi)
         Σi = randn(d, d)
         mul!(Σ[i], transpose(Σi), Σi)
-        kron_axpy!(Σ[i], V[i], Ω) # Ω = Σ[1]⊗V[1] + ... + Σ[m]⊗V[m]
     end
-    Ωchol = cholesky(Ω)
-    Y = X * B + reshape(Ωchol.L * randn(n * d), n, d)
+    Y = X * B
+    for i in 1:m
+        Z = randn(n, d)
+        Σchol = cholesky(Σ[i])
+        Vchol = cholesky(V[i])
+        rmul!(Z, transpose(Σchol.L))
+        lmul!(Vchol.L, Z)
+        Y .= Y .+ Z
+    end
     Y, X, V, B, Σ
 end
 Y, X, V, B, Σ = simulate(5_000, 4, 10, 2)