I am a beginner and trying to run Heritability analysis. I am reading the tutorial available and hoping that I could follow along with the tutorial. I am wondering where can I get the cg10k data to follow along with the tutorial.

Thanks,

Hi,

Is it possible to fit models where Y is a n x 1 vector and V contains multiple covariance matrices? If not, would it be possible to add such functionality?

My application is that I want to model the data s a function several covariance matrices encoding both gentic and environmental similarity.

All the best,
Espen

chol_gradient: Test Failed at D:\a\VarianceComponentModels.jl\VarianceComponentModels.jl\test\multivariate_calculus_test.jl:87
  Expression: norm(dL1 - dL2) ≈ 0.0
   Evaluated: 5.551115123125783e-17 ≈ 0.0
Stacktrace:
 [1] macro expansion
   @ C:\hostedtoolcache\windows\julia\1.7.0\x64\share\julia\stdlib\v1.7\Test\src\Test.jl:445 [inlined]
 [2] macro expansion
   @ D:\a\VarianceComponentModels.jl\VarianceComponentModels.jl\test\multivariate_calculus_test.jl:87 [inlined]
 [3] macro expansion
   @ C:\hostedtoolcache\windows\julia\1.7.0\x64\share\julia\stdlib\v1.7\Test\src\Test.jl:1283 [inlined]
 [4] top-level scope
   @ D:\a\VarianceComponentModels.jl\VarianceComponentModels.jl\test\multivariate_calculus_test.jl:79

Looks like we would better change a ≈ b to something like isapprox(a, b; atol=eps()).

Hi, Does VarianceComponentModels.jl, as it currently exists, permit the user to specify a model like:

vec(Y) ~ N(XB, Sigma)

where Y is a n by 2 matrix of (bivariate) phenotypes (and vec(Y) is a 2n-length vector), X is a known 2n by 2p matrix, B is a 2p-long vector, and Sigma is a 2n by 2n covariance matrix? I'm assuming that the covariance matrix Sigma = Sigma1 \otimes V1 + Sigma2 \otimes V2, where

Sigma1 & Sigma2 are known n by n matrices and V1 and V2 are unknown 2 by 2 matrices. I'd like to estimate V1, V2, and B.

If this isn't currently implemented, can you provide any suggestions as to how I might build on your existing code to implement it?

Much thanks!

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Hi all,
As far as I know, the current package can't be used in Julia 1.0.1. However, I am wondering whether you might upgrade the package for Julia 1.0.1 soon. Could you please let me know about it? Thanks in advance!

Should be able to deal with this corner case.

I think the format you list in your documentation under installation
Pkg.clone("[email protected]:OpenMendel/VarianceComponentModels.jl.git")
requires a public/private key set-up with github.

I think the alternate format
Pkg.clone("https://github.com/OpenMendel/VarianceComponentModels.jl.git")
would be better. It has the same basic functionality and can be run by anyone. Of course, since it does not require a public/private key, it is not as secure.

Trying to increase the dimension of the problem from 2 to 3 as follows:

generate data from a d-variate response variance component model
using Random, LinearAlgebra
Random.seed!(123)
n = 1000 # no. observations
d = 3 # dimension of responses
m = 3 # no. variance components
p = 3 # no. covariates

X = randn(n, p)

Step 3 model fit:

vcmodel_mle = deepcopy(vcmodel)
@time logl, vcmodel_mle, Σse, Σcov, Bse, Bcov = fit_mle!(vcmodel_mle, vcdata; algo = :MM);

fails:

MethodError: no method matching TwoVarCompVariateRotate(::VarianceComponentVariate{Float64,3,Array{Float64,2},Array{Float64,2},Array{Float64,2}})
Closest candidates are:
TwoVarCompVariateRotate(!Matched::Array{T<:VarianceComponentVariate,N} where N) where T<:VarianceComponentVariate at C:\Users\Denis.julia\packages\VarianceComponentModels\6qDHu\src\VarianceComponentModels.jl:261
TwoVarCompVariateRotate(!Matched::YT<:(Union{AbstractArray{T,2}, AbstractArray{T,1}} where T), !Matched::XT<:(Union{AbstractArray{T,2}, AbstractArray{T,1}} where T), !Matched::Array{T<:AbstractFloat,1}, !Matched::Array{T<:AbstractFloat,2}, !Matched::T<:AbstractFloat) where {T<:AbstractFloat, YT<:(Union{AbstractArray{T,2}, AbstractArray{T,1}} where T), XT<:(Union{AbstractArray{T,2}, AbstractArray{T,1}} where T)} at C:\Users\Denis.julia\packages\VarianceComponentModels\6qDHu\src\VarianceComponentModels.jl:189
TwoVarCompVariateRotate(!Matched::Union{AbstractArray{T,2}, AbstractArray{T,1}} where T, !Matched::Union{AbstractArray{T,2}, AbstractArray{T,1}} where T, !Matched::Array{T,1} where T, !Matched::Array{T,2} where T, !Matched::Real) at C:\Users\Denis.julia\packages\VarianceComponentModels\6qDHu\src\VarianceComponentModels.jl:209
...

Stacktrace:
[1] #fit_mle!#31(::Symbol, ::Function, ::VarianceComponentModel{Float64,3,Array{Float64,2},Array{Float64,2}}, ::VarianceComponentVariate{Float64,3,Array{Float64,2},Array{Float64,2},Array{Float64,2}}) at C:\Users\Denis.julia\packages\VarianceComponentModels\6qDHu\src\two_variance_component.jl:1207
[2] (::getfield(VarianceComponentModels, Symbol("#kw##fit_mle!")))(::NamedTuple{(:algo,),Tuple{Symbol}}, ::typeof(fit_mle!), ::VarianceComponentModel{Float64,3,Array{Float64,2},Array{Float64,2}}, ::VarianceComponentVariate{Float64,3,Array{Float64,2},Array{Float64,2},Array{Float64,2}}) at .\none:0
[3] top-level scope at util.jl:156
[4] top-level scope at In[10]:2

Thank you for the wonderful package. Things seem to run nicely. I have a couple questions.

1. It appears that the package only handles up to two variance components, whether be it univariate or multivariate? Is there a plan to expand to handle more than 3 variance components?

2. What would be the appropriate way of doing inference for these variance components? I imagine for univariate case, you can fit a null model without the variance component of interest, and running some variation of likelihood ratio test as done in GCTA, but would you recommend a different way? Also for multivariate cases, where you have the covariance parameters, is there a way to easily fit the null model with such a term set to zero, and run LRT?

I found that only one random component can be used in this packages, because of the second component is the identity matrix. And I want to use multiple random components in a model. But how can I do ? Thanks.

When using LRT for GWAS, eigen vectors of (co)variance matrices needed to be stored under the null model so that we can rotate covariates including SNPs under all alternative models and to avoid to call 'TwoVarCompModelRotate' many times.