A multivariate LMM sampler

[urchgene]

Hi,

Thanks for this package. i really like its concept but I'm still not clear on some things.

I am working on the following code for a multi-response LMM:

using Distributions
using DynamicHMC
using ContinuousTransformations
using DiffWrappers
using LaTeXStrings
using MCMCDiagnostics
using StatsBase

##### This struct gave so much problems
#struct Model{T}
#    observations::T
#    cholK::T
#    X::T
#end

struct Model
    observations::Array{Float64,2}
    cholK::Array{Float64,2}
    X::Array{Float64,2}
end

function (m::Model)(θ)
    sigmau, Omegau, sigmae, Omegae, u, beta = θ
    (all(isfinite.(sigmau)) && all(isfinite.(Omegau)) && all(isfinite.(sigmae)) && all(isfinite.(Omegae)))|| return -Inf
    Lsigu = Diagonal(sigmau) * Omegau
    Vu = Lsigu*Lsigu'
    Lsige = Diagonal(sigmae) * Omegae
    Ve = Lsige*Lsige'

    varvecG = kron(m.cholK, Vu)
    a = varvecG * u
    a = reshape(a, size(m.cholK, 1), size(Vu,2))
    mu = m.X*beta' + a

    dx = zeros(size(a,1))

   #### A matrix MvNormal without for-loop will be better but I do not know how this works yet in Julia.
    for i = 1:size(a,1)
        dist = MvNormal(mu[i, :], Ve)
        log_likelihood = sum(logpdf(dist, m.observations[i, :]))
        dx[i] = log_likelihood
    end

    samp_loglik = sum(dx)

    log_prior = sum(logpdf.(Cauchy(2.5), sigmau)) + sum(logpdf.(Cauchy(2.5), sigmae)) + lkj_correlation_cholesky_logpdf(Omegau, 2) + lkj_correlation_cholesky_logpdf(Omegae, 2) + sum(logpdf.(Normal(), u)) + sum(logpdf.(Normal(0, 1), beta))

    log_prior + samp_loglik
end

Y = rand(100,2);
K = rand(100,100); K = K*K'/mean(diag(K*K'));
cholK = chol(K) + 0.001*I;
X = ones(100,1);
ntraits = size(Y, 2);
nfixed = size(X, 2) * ntraits;
nrandom = size(K, 2) * ntraits;

### init params

sigmau = [1, 1];
Omegau = [1 0.1; 0.1 1]; Omegau = chol(Omegau);
sigmae = [1, 1]; 
Omegae = [1 0.1; 0.1 1]; Omegae = chol(Omegae);
u = rand(Normal(), 100*2);
beta = rand(Normal(0, 1000), 2);

#### Run like this ....

model = Model(Y, cholK, X)

model((sigmau, Omegau, sigmae, Omegae, u, beta))

### look at the ContinuousTransformations src to find how to map cov matrices, exp for sigma, real line for beta and u... Must be in same order as parameters...
θ_transformation  = TransformationTuple(
ArrayTransformation(bridge(ℝ, ℝ⁺), ntraits), CorrelationCholeskyFactor(ntraits), ArrayTransformation(bridge(ℝ, ℝ⁺), ntraits), CorrelationCholeskyFactor(ntraits), 
ArrayTransformation(bridge(ℝ, Segment(-100000.0, 100000.0)), nrandom), ArrayTransformation(bridge(ℝ, Segment(-100000.0, 100000.0)), nfixed))

modelT = TransformLogLikelihood(model, θ_transformation)

modelT(zeros(length(modelT)))

modelTG = ForwardGradientWrapper(modelT);

modelTG(zeros(length(modelTG)))

Error snippet was:

julia> modelTG(zeros(length(modelTG)))
ERROR: MethodError: Cannot `convert` an object of type ForwardDiff.Dual{ForwardDiff.Tag{ContinuousTransformations.TransformLogLikelihood{Model,ContinuousTransformations.TransformationTuple{Tuple{ContinuousTransformations.ArrayTransformation{ContinuousTransformations.ComposedTransformation{ContinuousTransformations.Affine{Float64},ContinuousTransformations.Exp},(2,)},ContinuousTransformations.CorrelationCholeskyFactor,ContinuousTransformations.ArrayTransformation{ContinuousTransformations.ComposedTransformation{ContinuousTransformations.Affine{Float64},ContinuousTransformations.Exp},(2,)},ContinuousTransformations.CorrelationCholeskyFactor,ContinuousTransformations.ArrayTransformation{ContinuousTransformations.ComposedTransformation{ContinuousTransformations.Affine{Float64},ContinuousTransformations.Logistic},(200,)},ContinuousTransformations.ArrayTransformation{ContinuousTransformations.ComposedTransformation{ContinuousTransformations.Affine{Float64},ContinuousTransformations.Logistic},(2,)}}}},Float64},Float64,10} to an object of type Float64
This may have arisen from a call to the constructor Float64(...),
since type constructors fall back to convert methods.
Stacktrace:
 [1] setindex!(::Array{Float64,1}, ::ForwardDiff.Dual{ForwardDiff.Tag{ContinuousTransformations.TransformLogLikelihood{Model,ContinuousTransformations.TransformationTuple{Tuple{ContinuousTransformations.ArrayTransformation{ContinuousTransformations.ComposedTransformation{ContinuousTransformations.Affine{Float64},ContinuousTransformations.Exp},(2,)},ContinuousTransformations.CorrelationCholeskyFactor,ContinuousTransformations.ArrayTransformation{ContinuousTransformations.ComposedTransformation{ContinuousTransformations.Affine{Float64},ContinuousTransformations.Exp},(2,)},ContinuousTransformations.CorrelationCholeskyFactor,ContinuousTransformations.ArrayTransformation{ContinuousTransformations.ComposedTransformation{ContinuousTransformations.Affine{Float64},ContinuousTransformations.Logistic},(200,)},ContinuousTransformations.ArrayTransformation{ContinuousTransformations.ComposedTransformation{ContinuousTransformations.Affine{Float64},ContinuousTransformations.Logistic},(2,)}}}},Float64},Float64,10}, ::Int64) at ./array.jl:583
 [2] (::Model)(::Tuple{Array{ForwardDiff.Dual{ForwardDiff.Tag{ContinuousTransformations.TransformLogLikelihood{Model,ContinuousTransformations.TransformationTuple{Tuple{ContinuousTransformations.ArrayTransformation{ContinuousTransformations.ComposedTransformation{ContinuousTransformations.Affine{Float64},ContinuousTransformations.Exp},(2,)},ContinuousTransformations.CorrelationCholeskyFactor,ContinuousTransformations.ArrayTransformation{ContinuousTransformations.ComposedTransformation{ContinuousTransformations.Affine{Float64},ContinuousTransformations.Exp},(2,)},ContinuousTransformations.CorrelationCholeskyFactor,ContinuousTransformations.ArrayTransformation{ContinuousTransformations.ComposedTransformation{ContinuousTransformations.Affine{Float64},ContinuousTransformations.Logistic},(200,)},ContinuousTransformations.ArrayTransformation{ContinuousTransformations.ComposedTransformation{ContinuousTransformations.Affine{Float64},ContinuousTransformations.Logistic},(2,)}}}},Float64},Float64,10},1},LowerTriangular{ForwardDiff.Dual{ForwardDiff.Tag{ContinuousTransformations.TransformLogLikelihood{Model,ContinuousTransformations.TransformationTuple{Tuple{ContinuousTransformations.ArrayTransformation{ContinuousTransformations.ComposedTransformation{ContinuousTransformations.Affine{Float64},ContinuousTransformations.Exp},(2,)},ContinuousTransformations.CorrelationCholeskyFactor,ContinuousTransformations.ArrayTransformation{ContinuousTransformations.ComposedTransformation

Any comments or guidance will be greatly appreciated!

Thanks.

Uche.

[tpapp]

The set of packages is undergoing a reorganization, and changes are in the process of trickling through the registry. I will address this problem shortly.

[tpapp]

I started looking at this, but I need more context. Can you please explain the data generating process (the model) and the parameter structure in detail? In particular, you seem to generate the random observations without using the parameters (which one would use for testing).

[urchgene]

You are right with the data generating process. I just wanted to get it working before making sure it does what it should be doing.

Here is a new data generating process for a correlated response matrix:

l=20
n=15
m=100

M = rand(Binomial(2,.2),m,l)
beta1 = rand(m).*exp(rand(Binomial(5,.2),m))
beta2 = rand(m).*exp(rand(Binomial(5,.1),m))
beta3 = rand(m).*exp(rand(Binomial(5,.1),m)) + beta2

g1 = M'*beta1
g2 = M'*beta2
g3 = M'*beta3
e1 = std(g1)*rand(l)
e2 = (-e1*2*std(g2)/std(g1) + 0.25*std(g2)/std(g1) * rand(l))
e3 = 1*(e1*0.25*std(g2)/std(g1) + 0.25*std(g2)/std(g1) * rand(l))

y1 = 10+g1+e1
y2 = -50+g2+e2
y3 = -5+g3+e3

Y = hcat(y1,y3)
cov(Y)

K = cov(M')
K = K/mean(diag(K))

X = ones(l,1)
Z = eye(l)

The parameters are beta and covariances:

beta is a $1 \times 2$ matrix of grand mean of both responses in Y

The covariance follows the LKJ as in your example but here we have 2 covariances Vu and Ve

This follows a classical linear mixed model; y = Xb + Zu + e with variances for u and e which are random terms in the model except that in this case this is multi-variate.

Lemme know if this explanations is Ok and sorry I replied late.

Uche.

[tpapp]

Sorry, still not getting it. Eg y2 seems not to be used at all, and the es appear to be heavily non-normal.

Is the following a fair description:

1. y, X, and Z are data,
2. b and the covariance matrix of [u,e] are parameters (the latter parametrized as, say, half-Cauchy on variances and LKJ for correlations),
3. [u,e] are conditionally IID for each observations,
4. y = X*b + Z*u .+ e is the data generating process

Then what does the first part have to do with it?

[urchgene]

y2 is not used cos in the code I gave an example of 2 responses instead of the 3 i did with data generation. My bad but the whole idea is to generate a correlated response Y with correlations coming from two random parts [u, e].

So each y has its own [u,e].

y is data, X and Z are design matrices of fixed and random parts respectively. In this example each obs in yi reps a subject in the random part u. This applies to all y in Y.

[b, Vu, Ve, sigmau, sigmae] are parameters.

[u] is a generated quantity from [Vu, CholK] called [a] in the code which form part of the conditional mean part in the model likelihood (given as y ~ MVN(mu=Xb + a, Ve).

[u] is not IID due to [Vu, cholK] but [e] is partially IID cos its distributed as e ~ MVN(0, Ve \otimes I)

if Z was used in the model it will come in here: mu = m.X*beta' + Za The reason is cos in this example it is just identity meaning no repeated records in Y for any subject in a or u as we call it here.

N/B: \otimes means kronecker from Latex synthax!!!

[tpapp]

Still not entirely clear:

1. there is no part "called [a] in the code"
2. cholK (which I guess is the Cholesky factor of the correlation of [u,e]) is also a parameter, isn't it?
3. what is Za?

Instead of mixing LaTeX and Julia code, please kindly provide the following:

1. a function which takes parameters and covariates, and returns data (generating parts randomly as needed),
2. meaningful example values of the parameters (optional, but would be helpful)

Eg

function gendata(covariates, b, Vu, Ve, ...)
    X = covariates.X
    Z = covariates.Z
    # magic happens
    Y # return observables
end

N = 100
b = [1.0, 2.0]
X = randn(N, 2)
# ... other covariates generation
covariates = (X = X, ...)

observables = gendata(data, ...)

Please make sure that it is valid, self-contained Julia code (ie it runs), under Julia 1.0.

[tpapp]

Lack of activity, closing.

[urchgene]

I will follow the new examples package and make a fresh attempt.

Thanks.

[urchgene]

Hi again,

Apologies for not following up with this. I have seen your new examples and they are more compact. Thanks for that. I have also tried to do a univariate one so its easier...

y = Xb + Zu + e ... is the model

################### New DynamicHMC code
############################################

using TransformVariables, LogDensityProblems, DynamicHMC, MCMCDiagnostics,
    Parameters, Statistics, Distributions

# Then define a structure to hold the data: observables, covariates, and the degrees of freedom for the prior.

struct LMM{TY <: AbstractVector, TX <: AbstractMatrix, TZ <: AbstractMatrix, TK <: AbstractMatrix}
    y::TY
    X::TX
    Z::TZ
    K::TK
end

# Then make the type callable with the parameters *as a single argument*.

function (problem::LMM)(θ)
    @unpack y, X, Z, K = problem   # extract the data
    @unpack beta, sigu, sige, a = θ            # works on the named tuple too
    u = sigu*K*a
    loglikelihood(Normal(0, sige), y .- X*β .- Z*u) + logpdf.(Cauchy(2.5), sigu) + logpdf.(Cauchy(2.5), sige) + sum(logpdf.(Normal(0., 1.), a))  
end

N = 100
X = hcat(ones(N), randn(N, 2));
β = [1.0, 2.0, -1.0]
sige = 0.7
sigu = 0.3
Z = zeros(N,N) + 1.0I;
m=100;l=20;
M = rand(Binomial(2,.2),m,l);
K = M*M'/20;
Kchol = cholesky(K + 0.0001I)
K = Array(Kchol.L)
a = rand(Normal(0., 1.), N)
u = sigu*K*a
y = X*β .+ Z*u .+ randn(N) .* sige;
p = LMM(y, X, Z, K);

problem_transformation(p::LMM) = as((beta = as(Array, size(p.X, 2)), sigu = asℝ₊, sige = asℝ₊, a = as(Real, -∞, ∞)))
P = TransformedLogDensity(problem_transformation(p), p)
∇P = FluxGradientLogDensity(P);
chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 1000);

It does not sample cos I must have gotten something wrong. I hope this one is much more easier to understand since I modeled it as per your Linear Regression code.

Thanks below is the error:

ERROR: MethodError: -(::ForwardDiff.Dual{Nothing,Float64,1}, ::Flux.Tracker.TrackedReal{Float64}) is ambiguous. Candidates:
  -(a::Real, b::Flux.Tracker.TrackedReal) in Flux.Tracker at /Users/ugokek/.julia/packages/Flux/rcN9D/src/tracker/scalar.jl:77
  -(x::ForwardDiff.Dual{Tx,V,N} where N where V, y::Real) where Tx in ForwardDiff at /Users/ugokek/.julia/packages/ForwardDiff/KHAmA/src/dual.jl:129
Possible fix, define
  -(::ForwardDiff.Dual{Tx,V,N} where N where V, ::Flux.Tracker.TrackedReal)
Stacktrace:
 [1] zval(::Flux.Tracker.TrackedReal{Float64}, ::Flux.Tracker.TrackedReal{Float64}, ::ForwardDiff.Dual{Nothing,Float64,1}) at /Users/ugokek/.julia/packages/StatsFuns/0W2sM/src/distrs/norm.jl:4
 [2] normlogpdf(::Flux.Tracker.TrackedReal{Float64}, ::Flux.Tracker.TrackedReal{Float64}, ::ForwardDiff.Dual{Nothing,Float64,1}) at /Users/ugokek/.julia/packages/StatsFuns/0W2sM/src/distrs/norm.jl:12
 [3] logpdf(::Normal{Flux.Tracker.TrackedReal{Float64}}, ::ForwardDiff.Dual{Nothing,Float64,1}) at /Users/ugokek/.julia/packages/Distributions/WHjOk/src/univariates.jl:526

Uche

[tpapp]

This works (with LogDensityProblems#master, for th ADgradient):

using TransformVariables, LogDensityProblems, DynamicHMC, MCMCDiagnostics,
    Parameters, Statistics, Distributions, LinearAlgebra

# Then define a structure to hold the data: observables, covariates, and the degrees of freedom for the prior.

struct LMM{TY, TX, TZ, TL}
    "observed outcomes"
    y::TY
    "covariates"
    X::TX
    "covariates for `u`"
    Z::TZ
    "Cholesky factor for the *correlation* matrix `K=L*L'` of `u` (assumed known)."
    L::TL
end

function (problem::LMM)(θ)
    @unpack y, X, Z, L = problem # extract the data
    @unpack β, σᵤ, σₑ, a = θ     # extract the parameters
    u = σᵤ .* (L * a)
    loglik = loglikelihood(Normal(0, σₑ), y .- X*β .- Z*u) + loglikelihood(Normal(0.0, 1.0), a)
    logpri = logpdf(Cauchy(2.5), σᵤ) + logpdf(Cauchy(2.5), σₑ)
    loglik + logpri
end

# generate simulated data
N = 100
X = hcat(ones(N), randn(N, 2));
β = [1.0, 2.0, -1.0]
σₑ = 0.7
σᵤ = 0.3
Z = I

function make_covariance_factor(N; l = 20, d = Binomial(2, 0.2))
    M = rand(d, N, l)
    Σ = M*M' ./ l + 0.0001 * I
    σ = Diagonal(.√diag(Σ))
    cholesky(Symmetric(σ \ Σ / σ)).L
end

L = make_covariance_factor(N)
a = rand(Normal(0.0, 1.0), N)
u = σᵤ*L*a
y = X*β .+ Z*u .+ rand(Normal(0, σₑ), N)
p = LMM(y, X, Z, L);

problem_transformation(p::LMM) =
    as((β = as(Array, size(p.X, 2)), σᵤ = asℝ₊, σₑ = asℝ₊, a = as(Array, size(p.L, 1))))
P = TransformedLogDensity(problem_transformation(p), p)
∇P = ADgradient(Val(:ForwardDiff), P)

# test log density
logdensity(LogDensityProblems.Value, P, randn(dimension(P)))
logdensity(LogDensityProblems.ValueGradient, ∇P, randn(dimension(P)))

chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 1000)

NUTS_statistics(chain)

posterior = P.transformation.(get_position.(chain))

# compare known parameters to posterior
mean(y -> y.β, posterior)
β

mean(y -> y.σₑ, posterior)
σₑ

mean(y -> y.σᵤ, posterior)
σᵤ

# mixing -- rather good
sort(vec(mapslices(effective_sample_size, get_position_matrix(chain); dims = 1)))

Let me know if I can help with anything else.

[tpapp]

Also, would you be willing to contribute this example to the examples repo? I can make the PR, but would appreciate a few words about the context, what kind of model this is, etc. This would help others using this package.

[urchgene]

Hi @tpapp , this is awesome. Thanks. Sure you can use this as an example. This model is a very vital model in genetic analysis covering a wide range in the life sciences from plant and animal breeding to discovery in human health. The next step will be a bivariate example. I will try put something together in the coming days.

Thanks again.

Uche.

[urchgene]

Actually @tpapp, this below worked better:

function make_covariance_factor(N; l = 20, d = Binomial(2, 0.2))
           M = rand(d, N, l)
           Σ = M*M' ./ l + 0.001 * I
           cholesky(Σ).L
       end

Thanks.

Uche.

[tpapp]

But then you are not getting the cholesky factor of a correlation matrix, only of a covariance matrix. Is this intentional? I am fine with it either way.

Does this model have a standard name (in some field/literature)?

[urchgene]

It is intentional to get the cholesky factor for a covariance matrix.

The model is commonly called a Mixed Model in literature. For a single response (univariate mixed model, BLUP model, etc) for a matrix response (multivariate mixed model, multi-trait mixed model etc).

[urchgene]

I'm also wondering if lkj_correlation_cholesky_logpdf is still present in DynamicHMC?

Thanks.

[urchgene]

Hi @tpapp,

So I have come up with a toy example for the bivariate case.

############################################
################### New DynamicHMC code
########### mvLMM model...
############################################

using TransformVariables, LogDensityProblems, DynamicHMC, MCMCDiagnostics,
    Parameters, Statistics, Distributions, LinearAlgebra

# Then define a structure to hold the data: observables, covariates, and the degrees of freedom for the prior.

struct mvLMM{TY, TX, TZ, TL}
    "observed outcomes"
    Y::TY
    "covariates"
    X::TX
    "covariates for `u`"
    Z::TZ
    "Cholesky factor for the *correlation* matrix `K=L*L'` of `u` (assumed known)."
    L::TL
end

function (problem::mvLMM)(θ)
    @unpack Y, X, Z, L = problem # extract the data
    @unpack β, Vu, Ve, a = θ     # extract the parameters
    N = size(L, 1); d = size(Y, 1);
    P = zeros(d,d) + 1.0*I;
    u = kron(Vu, L)* a
    u = reshape(u, (N,d))
    loglik = loglikelihood(MvNormal(fill(0., d), Ve), Y .- (X*β')' .- (Z*u')) + loglikelihood(Normal(0.0, 1.0), a)
    logpri = logpdf(InverseWishart(d, P), Vu) + logpdf(InverseWishart(d, P), Ve)
    loglik + logpri
end

# generate simulated data
N = 100
d = 2
X = hcat(ones(N), randn(N, 2));
β = [1.0 2.0 -1.0; 1.0 2.0 -1.0]
P = zeros(2,2) + 1.0*I;
tt = rand(InverseWishart(2, P),2)
Vu = tt[1]; Ve = tt[2]; 
u = kron(Vu, L)* a; u = reshape(u, (N,d));
Z = I
Y = (X*β')' .+ Z*u' .+ rand(MvNormal([0.,0.], Ve), N)

function make_covariance_factor(N; l = 20, d = Binomial(2, 0.2))
    M = rand(d, N, l)
    Σ = M*M' ./ l + 0.15 * I
    cholesky(Σ).L
end

L = make_covariance_factor(N)
#N = size(L, 1); d = size(Y, 1);
a = rand(Normal(0.0, 1.0), N*d); 
p = mvLMM(Y, X, Z, L);

problem_transformation(p::mvLMM) = as((β = as(Array, size(p.Y, 1) , size(p.X, 2)), Vu = as(Array, size(p.Y, 1), size(p.Y, 1)), Ve = as(Array, size(p.Y, 1), size(p.Y, 1)), a = as(Array, size(p.L, 1)* size(p.Y, 1), 1)))
PL = TransformedLogDensity(problem_transformation(p), p)
∇P = ADgradient(Val(:ForwardDiff), PL)

# test log density
logdensity(LogDensityProblems.Value, PL, randn(dimension(PL)))
logdensity(LogDensityProblems.ValueGradient, ∇P, randn(dimension(PL)))

chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 1000)

NUTS_statistics(chain)

posterior = P.transformation.(get_position.(chain))

# compare known parameters to posterior
mean(y -> y.β, posterior)
β

mean(y -> y.σₑ, posterior)
σₑ

mean(y -> y.σᵤ, posterior)
σᵤ

# mixing -- rather good
sort(vec(mapslices(effective_sample_size, get_position_matrix(chain); dims = 1)))

Problem is

ERROR: PosDefException: matrix is not Hermitian; Cholesky factorization failed.
Stacktrace:
 [1] checkpositivedefinite(::Int64) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/factorization.jl:11
 [2] #cholesky!#91(::Bool, ::Function, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformNamedTuple{(:β, :Vu, :Ve, :a),NTuple{4,TransformVariables.ArrayTransform{TransformVariables.Identity,2}}},mvLMM{Array{Float64,2},Array{Float64,2},UniformScaling{Bool},LowerTriangular{Float64,Array{Float64,2}}}}},Float64},Float64,10},2}, ::Val{false}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/cholesky.jl:182
 [3] #cholesky! at ./none:0 [inlined] (repeats 2 times)
 [4] #cholesky#95(::Bool, ::Function, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformNamedTuple{(:β, :Vu, :Ve, :a),NTuple{4,TransformVariables.ArrayTransform{TransformVariables.Identity,2}}},mvLMM{Array{Float64,2},Array{Float64,2},UniformScaling{Bool},LowerTriangular{Float64,Array{Float64,2}}}}},Float64},Float64,10},2}, ::Val{false}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/cholesky.jl:275

I am thinking this is propagated via MvNormal. What do you suggest? So this problem is mainly because there is no longer an LKJ prior in your new build of the package or am I missing where it is?

Thanks.

Uche.

[tpapp]

You can wrap the variance matrix in a Symmetric which should make the above work.

Regarding LKJ, for now you can copy-paste it from https://github.com/tpapp/ContinuousTransformations.jl/blob/f8eb7fcb6b74c144dc39d36666325d4a34bf470c/src/ContinuousTransformations.jl#L873-L900 but I will find a new home for it soon.

[tpapp]

I ended up packaging them up:

https://github.com/tpapp/AltDistributions.jl

[urchgene]

@tpapp, thanks a lot for the AltDistributions package and your efforts towards answering my questions. I have shopped around quite a bit and DynamicHMC is the most adaptable toolkit for this kind of problems and that is why I am very interested in it. Its flexibility is amazing!

So for the issues:

1. wrapping the variance matrix in a Symmetric did not work. Unless I did something wrong.

2. I used the LKJL but encountered a method error:

using TransformVariables, LogDensityProblems, DynamicHMC, MCMCDiagnostics,
    Parameters, Statistics, Distributions, LinearAlgebra, AltDistributions

# Then define a structure to hold the data: observables, covariates, and the degrees of freedom for the prior.

struct mvLMM{TY, TX, TZ, TL}
    "observed outcomes"
    Y::TY
    "covariates"
    X::TX
    "covariates for `u`"
    Z::TZ
    "Cholesky factor for the *correlation* matrix `K=L*L'` of `u` (assumed known)."
    L::TL
end

function (problem::mvLMM)(θ)
    @unpack Y, X, Z, L = problem # extract the data
    @unpack β, Omegau, Omegae, sigmau, sigmae, a = θ     # extract the parameters
    N = size(Y, 2); d = size(Y, 1);

    Lsigu = Diagonal(sigmau) * Omegau
    Vu = Lsigu*Lsigu'
    Lsige = Diagonal(sigmae) * Omegae
    Ve = Lsige*Lsige'

    u = kron(Vu, L)* a
    u = reshape(u, (N,d))
    Ve = Symmetric(Ve + 0.1*I);
    loglik = loglikelihood(MvNormal(fill(0., d), Ve), Y .- (X*β')' .- (Z*u')) + loglikelihood(Normal(0.0, 1.0), a)
    logpri = sum(logpdf.(Cauchy(2.5), sigmau)) + sum(logpdf.(Cauchy(2.5), sigmae)) + logpdf.(LKJL(d), Omegau) + logpdf.(LKJL(d), Omegae)
    loglik + logpri
end

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

# generate simulated data
N = 100
d = 2
X = hcat(ones(N), randn(N, 2));
β = [1.0 2.0 -1.0; 1.0 2.0 -1.0]
P = zeros(2,2) + 1.0*I;
tt = rand(InverseWishart(2, P),2)
Omegau = Symmetric(tt[1]); Omegae = Symmetric(tt[2]); 
sigmau = rand(Cauchy(2.5), d); sigmae = rand(Cauchy(2.5), d);
Lsigu = Diagonal(sigmau) * Omegau; Lsige = Diagonal(sigmae) * Omegae;

function make_covariance_factor(N; l = 20, d = Binomial(2, 0.2))
    M = rand(d, N, l)
    Σ = M*M' ./ l + 0.15 * I
    cholesky(Σ).L
end

L = make_covariance_factor(N)
a = rand(Normal(0.0, 1.0), N*d);
Vu = Lsigu*Lsigu'; Ve = Lsige*Lsige';
u = kron(Vu, L)* a; u = reshape(u, (N,d));
Z = I
Y = (X*β')' .+ Z*u' .+ rand(MvNormal([0.,0.], Ve), N)

p = mvLMM(Y, X, Z, L);

problem_transformation(p::mvLMM) = as((β = as(Array, size(p.Y, 1) , size(p.X, 2)), Omegau = as(Array, size(p.Y, 1), size(p.Y, 1)), Omegae = as(Array, size(p.Y, 1), size(p.Y, 1)), a = as(Array, size(p.L, 1)* size(p.Y, 1), 1), sigmau = as(Array, size(p.Y, 1)),  sigmae = as(Array, size(p.Y, 1))))
PL = TransformedLogDensity(problem_transformation(p), p)
∇P = ADgradient(Val(:ForwardDiff), PL)

# test log density
logdensity(LogDensityProblems.Value, PL, randn(dimension(PL)))
logdensity(LogDensityProblems.ValueGradient, ∇P, randn(dimension(PL)))

Method error below:

julia> logdensity(LogDensityProblems.Value, PL, randn(dimension(PL)))
ERROR: MethodError: no method matching length(::LKJL{Int64})
Closest candidates are:
  length(::Core.SimpleVector) at essentials.jl:582
  length(::Base.MethodList) at reflection.jl:728
  length(::Core.MethodTable) at reflection.jl:802
  ...
Stacktrace:
 [1] _similar_for(::UnitRange{Int64}, ::Type, ::LKJL{Int64}, ::Base.HasLength) at ./array.jl:532
 [2] _collect(::UnitRange{Int64}, ::LKJL{Int64}, ::Base.HasEltype, ::Base.HasLength) at ./array.jl:563
 [3] collect(::LKJL{Int64}) at ./array.jl:557
 [4] broadcastable(::LKJL{Int64}) at ./broadcast.jl:609
 [5] broadcasted at ./broadcast.jl:1163 [inlined]
 [6] (::mvLMM{Array{Float64,2},Array{Float64,2},UniformScaling{Bool},LowerTriangular{Float64,Array{Float64,2}}})(::NamedTuple{(:β, :Omegau, :Omegae, :a, :sigmau, :sigmae),Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,1},Array{Float64,1}}}) at ./REPL[19]:8
 [7] transform_logdensity(::TransformVariables.TransformNamedTuple{(:β, :Omegau, :Omegae, :a, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},TransformVariables.ArrayTransform{TransformVariables.Identity,2},TransformVariables.ArrayTransform{TransformVariables.Identity,2},TransformVariables.ArrayTransform{TransformVariables.Identity,2},TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}, ::mvLMM{Array{Float64,2},Array{Float64,2},UniformScaling{Bool},LowerTriangular{Float64,Array{Float64,2}}}, ::Array{Float64,1}) at /Users/ugokek/.julia/packages/TransformVariables/G0eaF/src/TransformVariables.jl:144

3. I will like to monitor Vu and Ve in the mvLMM function but they do not have priors. Is this doable?

Thanks again.

Uche.

[tpapp]

Two things:

I don't see why you are broadcasting the logpdf over the correlation Cholesky factors. I would simply use

    logpri = sum(logpdf.(Cauchy(2.5), sigmau)) + sum(logpdf.(Cauchy(2.5), sigmae)) +
        logpdf(LKJL(d), Omegau) + logpdf(LKJL(d), Omegae)

Also, you should ensure that you are transforming into said matrices, eg with

function problem_transformation(p::mvLMM)
    nY = size(p.Y, 1)
    as((β = as(Array, nY, size(p.X, 2)),
        Omegau = CorrCholeskyFactor(nY),
        Omegae = CorrCholeskyFactor(nY),
        a = as(Array, size(p.L, 1)*nY, 1),
        sigmau = as(Array, nY),  sigmae = as(Array, nY)))
end

[tpapp]

Sorry @urchgene, I missed this part:

I will like to monitor Vu and Ve in the mvLMM function but they do not have priors. Is this doable?

What I would do is factor out

    Lsigu = Diagonal(sigmau) * Omegau
    Vu = Lsigu*Lsigu'
    Lsige = Diagonal(sigmae) * Omegae
    Ve = Lsige*Lsige'

into a function, which you can call inside the log pdf and then later for posterior analysis.

[urchgene]

Thanks a lot @tpapp.

Will test how fast this will be for real dataset problems.

[urchgene]

@tpapp, what ways do you suggest may be ways to significantly increase the sampling speed for the models?

[tpapp]

Regarding optimization: I should discuss this in the manual, for which I have no time for at the moment, but I started a discussion at #11 with some suggestions.

[urchgene]

@tpapp Happy New Year and sorry for reopening...

I tried to re-parameterize the model for better optimization and accuracy of results...

############################################
################### New DynamicHMC code
########### mvLMM model...
############################################

using TransformVariables, LogDensityProblems, DynamicHMC, MCMCDiagnostics,
    Parameters, Statistics, Distributions, LinearAlgebra, AltDistributions, PositiveFactorizations, PDMats

# Then define a structure to hold the data: observables, covariates, and the degrees of freedom for the prior.

struct mvLMM{TY, TX, TZ, TL}
    "observed outcomes"
    Y::TY
    "covariates"
    X::TX
    "covariates for `u`"
    Z::TZ
    "Cholesky factor for the *correlation* matrix `K=L*L'` of `u` (assumed known)."
    L::TL
    "something else"
end

##############################################################
##### Function to Build model ############
##############################################################

#function getcov(sigma, Omega)
#    Lsigu = Diagonal(sigma) * Omega
#    Vu = Lsigu*Lsigu'
#    Vu
#end

function getcov(sigma, Omega)
    Vu = Diagonal(sigma) * Omega * Diagonal(sigma)
    Vu
end

function (problem::mvLMM)(θ)
    @unpack Y, X, Z, L = problem # extract the data
    @unpack β, Omegau, Omegae, sigmau, sigmae = θ     # extract the parameters
    N = size(L, 1); d = size(Y, 1);

    Vu = getcov(sigmau, Omegau)
    Ve = getcov(sigmae, Omegae)

    Gcov = kron(Vu, L);
    Rcov = kron(Ve, zeros(N,N) + 1.0*I);

    mu = vec(X*β')    
    Rsig = Gcov + Rcov + 0.001*I;
    Rsig = PDMat(Rsig, cholesky(Positive, Rsig))
    loglik = MvNormal(mu, Rsig)
    vecY = vec(Y); 
    loglik = loglikelihood(loglik, reshape(vecY, (size(vecY, 1),1))) ## loglik fn must be 2-dim
    logpri = sum(logpdf.(Cauchy(2.5), sigmau)) + sum(logpdf.(Cauchy(2.5), sigmae)) + logpdf(LKJL(d), Omegau) + logpdf(LKJL(d), 
Omegae)
    loglik + logpri
end

##############################################################
###### Function to transform parameters for sampling #########
##############################################################

function problem_transformation(p::mvLMM)
    nY = size(p.Y, 1)
    as((β = as(Array, nY, size(p.X, 2)),
        Omegau = CorrCholeskyFactor(nY),
        Omegae = CorrCholeskyFactor(nY),
        sigmau = as(Array, nY),  sigmae = as(Array, nY)))
end

#################################
## Run here
#################################

p = mvLMM(Y, X, Z, L);

PL = TransformedLogDensity(problem_transformation(p), p)
∇P = ADgradient(Val(:ForwardDiff), PL)

# test log density
logdensity(LogDensityProblems.Value, PL, randn(dimension(PL)))
logdensity(LogDensityProblems.ValueGradient, ∇P, randn(dimension(PL)))

chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 1000)

NUTS_statistics(chain)

posterior = PL.transformation.(get_position.(chain))

I will use R to simulate data using RCall because R's modelmatrix function is more accurate...

R code in the file gendata_for_mvLMM.R

l=200;
n<-15;
m<-40;

M<-matrix(rbinom(m*l,2,.25),nrow=l);
rownames(M)<-paste("l",1:nrow(M));
beta1<-rnorm(m)*exp(rbinom(m,5,.25));
beta2<-rnorm(m)*exp(rbinom(m,5,.1));
beta3<- rnorm(m)*exp(rbinom(m,5,.1))+beta2;

g1<-M%*%beta1;
g2<-M%*%beta2;
g3<-M%*%beta3;
e1<-sd(g1)*rnorm(l);
e2<-(-e1*2*sd(g2)/sd(g1)+.25*sd(g2)/sd(g1)*rnorm(l));
e3<-1*(e1*.45*sd(g2)/sd(g1)+.25*sd(g2)/sd(g1)*rnorm(l));

y1<-10+g1+e1;
y2<--50+g2+e2;
y3<--5+g3+e3;

Y<-rbind(t(y1),t(y2), t(y3));

colnames(Y)<-rownames(M);
cov(t(Y));
Y[1:3,1:5];

Kmat<-cov(t(M));
Kmat<-Kmat/mean(diag(Kmat));
Kmat <- 0.95*Kmat + 0.05*diag(nrow(Kmat));
rownames(Kmat)<-colnames(Kmat)<-rownames(M);
X<-matrix(1,nrow=1,ncol=l);
colnames(X)<-rownames(M);
Z<-diag(l);

###########################################################
### Simulate Data using RCall cos modelmatrix in R is better
############################################################

R"source('gendata_for_mvLMM.R')"
@rget Y;
@rget Kmat;
@rget X;
@rget Z;

X = X';

using PositiveFactorizations
using LinearAlgebra

L = cholesky(Positive, Kmat).L;

The code runs and fails at the stage:

logdensity(LogDensityProblems.ValueGradient, ∇P, randn(dimension(PL)))

Here is the error: Method error as it says;

MethodError: no method matching floattype(::Type{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##3#4")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}},Float64},Float64,8}})
Closest candidates are:
  floattype(!Matched::Type{T<:AbstractFloat}) where T<:AbstractFloat at /Users/ugokek/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:264
  floattype(!Matched::Type{T<:Integer}) where T<:Integer at /Users/ugokek/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:265

Stacktrace:
 [1] default_δ(::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##3#4")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}},Float64},Float64,8},2}) at /Users/ugokek/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:269
 [2] default_tol(::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##3#4")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}},Float64},Float64,8},2}) at /Users/ugokek/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:270
 [3] cholesky(::Type{Positive}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##3#4")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}},Float64},Float64,8},2}, ::Type) at /Users/ugokek/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:7 (repeats 2 times)
 [4] (::mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}})(::NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##3#4")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}},Float64},Float64,8},2},UpperTriangular{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##3#4")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}},Float64},Float64,8},Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##3#4")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}},Float64},Float64,8},2}},UpperTriangular{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##3#4")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}},Float64},Float64,8},Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##3#4")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}},Float64},Float64,8},2}},Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##3#4")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}},Float64},Float64,8},1},Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##3#4")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}},Float64},Float64,8},1}}}) at ./In[11]:45
 [5] transform_logdensity(::TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}}, ::mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##3#4")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}},Float64},Float64,8},1}) at /Users/ugokek/.julia/packages/TransformVariables/YDIkx/src/generic.jl:135
 [6] macro expansion at /Users/ugokek/.julia/packages/Parameters/35RKe/src/Parameters.jl:733 [inlined]
 [7] logdensity(::Type{LogDensityProblems.Value}, ::TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##3#4")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}},Float64},Float64,8},1}) at /Users/ugokek/.julia/packages/LogDensityProblems/v6mCF/src/LogDensityProblems.jl:164
 [8] #3 at /Users/ugokek/.julia/packages/LogDensityProblems/v6mCF/src/LogDensityProblems.jl:229 [inlined]
 [9] chunk_mode_gradient!(::DiffResults.MutableDiffResult{1,Float64,Tuple{Array{Float64,1}}}, ::getfield(LogDensityProblems, Symbol("##3#4")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}}, ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##3#4")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}},Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##3#4")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}},Float64},Float64,8},1}}) at /Users/ugokek/.julia/packages/ForwardDiff/KHAmA/src/gradient.jl:139
 [10] gradient! at /Users/ugokek/.julia/packages/ForwardDiff/KHAmA/src/gradient.jl:37 [inlined]
 [11] gradient! at /Users/ugokek/.julia/packages/ForwardDiff/KHAmA/src/gradient.jl:33 [inlined]
 [12] macro expansion at /Users/ugokek/.julia/packages/Parameters/35RKe/src/Parameters.jl:733 [inlined]
 [13] logdensity(::Type{LogDensityProblems.ValueGradient}, ::LogDensityProblems.ForwardDiffLogDensity{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##3#4")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}},Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##3#4")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}},Float64},Float64,8},1}}}, ::Array{Float64,1}) at /Users/ugokek/.julia/packages/LogDensityProblems/v6mCF/src/AD_ForwardDiff.jl:41
 [14] top-level scope at In[14]:1

Thanks again.

Uche.

[tpapp]

I am happy to help. I think your problem comes from

Rsig = PDMat(Rsig, cholesky(Positive, Rsig))

since the method in PositiveFactorizations can't handle generic floats. You may want to open an issue there.

That said, I would always try to use cholesky factors directly (which you seem to be doing, so I don't get why you need to decompose again).

[urchgene]

Thanks for prompt reply @tpapp

I use the PositiveFactorizations because MvNormal fails with Rsig.

PosDefException: matrix is not Hermitian; Cholesky factorization failed.

Stacktrace:
 [1] checkpositivedefinite(::Int64) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/factorization.jl:11
 [2] #cholesky!#91(::Bool, ::Function, ::Array{Float64,2}, ::Val{false}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/cholesky.jl:182
 [3] #cholesky! at ./none:0 [inlined] (repeats 2 times)
 [4] #cholesky#95 at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/cholesky.jl:275 [inlined]
 [5] cholesky at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/cholesky.jl:275 [inlined] (repeats 2 times)
 [6] Type at /Users/ugokek/.julia/packages/PDMats/Kouno/src/pdmat.jl:19 [inlined]
 [7] Type at /Users/ugokek/.julia/packages/Distributions/WHjOk/src/multivariate/mvnormal.jl:208 [inlined]
 [8] (::mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}})(::NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{Array{Float64,2},UpperTriangular{Float64,Array{Float64,2}},UpperTriangular{Float64,Array{Float64,2}},Array{Float64,1},Array{Float64,1}}}) at ./In[4]:45
 [9] transform_logdensity(::TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}}, ::mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}, ::Array{Float64,1}) at /Users/ugokek/.julia/packages/TransformVariables/YDIkx/src/generic.jl:135
 [10] macro expansion at /Users/ugokek/.julia/packages/Parameters/35RKe/src/Parameters.jl:733 [inlined]
 [11] logdensity(::Type{LogDensityProblems.Value}, ::TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :Omegau, :Omegae, :sigmau, :sigmae),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,2},CorrCholeskyFactor,CorrCholeskyFactor,TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}},mvLMM{Array{Float64,2},Adjoint{Float64,Array{Float64,2}},Array{Float64,2},LowerTriangular{Float64,Array{Float64,2}}}}, ::Array{Float64,1}) at /Users/ugokek/.julia/packages/LogDensityProblems/v6mCF/src/LogDensityProblems.jl:164
 [12] top-level scope at In[6]:1

Do you have any ideas cos my matrix is not always PD.

Uche.

[tpapp]

I ran into this a lot, too, so I wrote AltMvNormal. Use Symmetric to pass general matrices.

[urchgene]

Awesome. Thanks @tpapp. I will try it out ASAP.

[urchgene]

Haha... trickle down effects..

AltMvNormal is very cool but loglikelihood from Distributions fail

MethodError: no method matching loglikelihood(::AltMvNormal{Array{Float64,1},LowerTriangular{Float64,Array{Float64,2}}}, ::Array{Float64,2})

[tpapp]

Use logpdf.

[urchgene]

Thanks a lot @tpapp. I really do appreciate your advice cos you don't have to!!!

However, the cholesky problem (non PD) persists:

@timholy suggested this for PositiveFactorizations

see below:

PositiveFactorizations.floattype(::Type{MyWeirdNumberType}) = Float64
That should work if convert(Float64, x::MyWeirdNumberType) works.

I am still very new to Julia types. How can I incorporate this?

Thanks.

I am leaning on this now, so I want to incorporate it into this:

function (problem::mvLMM)(θ)

    @unpack Y, X, Z, L = problem # extract the data
    @unpack β, Omegau, Omegae, sigmau, sigmae = θ     # extract the parameters
    N = size(L, 1); d = size(Y, 1);

    Vu = getcov(sigmau, Omegau)
    Ve = getcov(sigmae, Omegae)

    Gcov = kron(Vu, L);
    Rcov = kron(Ve, zeros(N,N) + 1.0*I);

    mu = vec(X*β')    
    Rsig = Gcov + Rcov + 0.01*I;
    xxx = cholesky(Rsig).L
    loglik = AltMvNormal(Val(:L), mu, xxx)
    vecY = vec(Y); 
    loglik = logpdf(loglik, vecY) ## loglik fn must be 2-dim
    logpri = sum(logpdf.(Cauchy(2.5), sigmau)) + sum(logpdf.(Cauchy(2.5), sigmae)) + logpdf(LKJL(d), Omegau) + logpdf(LKJL(d), 
Omegae)
    loglik + logpri
end

The initial problem with generic float continues with this. So I am guessing @timholy's suggestion will fix it if done well.

Thanks.

[urchgene]

Details below

https://github.com/timholy/PositiveFactorizations.jl/issues/19

[joshualeond]

Hi @urchgene, were you able to resolve these issues? I believe I have a similar issue attempting to use the new LKJ distribution in Distributions.jl with Turing: https://github.com/TuringLang/Bijectors.jl/issues/108

I'm constantly seeing the following error whenever I attempt to sample with NUTS:

ERROR: PosDefException: matrix is not positive definite; Cholesky factorization failed.

@tpapp, I know this is unrelated to your repo but was only curious. Thanks.

[urchgene]

Hi @joshualeond It was not resolved. I could also not figure out how to deal with it.

[tpapp]

As I suggested above, I still think that the best way is to parametrize variance matrices is to use the Cholesky factor, especially using the covariance matrix so one can use different priors for the variances. This is described in detail in the Stan manual; TransformVariables.jl has the relevant transformation implemented TransformVariables.CorrCholeskyFactor.

[joshualeond]

Thanks @urchgene and @tpapp, I'd like to see if the CorrCholeskyFactor could possibly fit into the Turing model spec. I'm not sure it will but I am sure there's useful info from your implementation.

[tpapp]

Closing this for lack of activity (feel free to comment here if you want it reopened). If there is a final example that works, a PR would be welcome.