diff --git a/Project.toml b/Project.toml index b485404..462f1ef 100644 --- a/Project.toml +++ b/Project.toml @@ -3,7 +3,7 @@ uuid = "a1dec852-9fe5-11e9-361f-8d9fde67cfa2" keywords = ["lenearmodel", "mixedmodel"] desc = "Mixed-effects models with flexible covariance structure." authors = ["Vladimir Arnautov "] -version = "0.16.1" +version = "0.16.2" [deps] DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" diff --git a/docs/src/api.md b/docs/src/api.md index 2976404..d48b79a 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -230,10 +230,28 @@ Metida.thetalength ### Metida.vmatrix +```@docs +Metida.vmatrix +``` + +### Metida.vmatrix! + ```@docs Metida.vmatrix! ``` +### Metida.m2logreml + +```@docs +Metida.m2logreml +``` + +### Metida.logreml + +```@docs +Metida.logreml +``` + ## StatsAPI ### Metida.aic @@ -463,3 +481,14 @@ Metida.raneflenv Metida.edistance ``` +### Metida.m2logml + +```@docs +Metida.m2logml +``` + +### Metida.logml + +```@docs +Metida.logml +``` diff --git a/src/Metida.jl b/src/Metida.jl index e8950eb..c6438b6 100644 --- a/src/Metida.jl +++ b/src/Metida.jl @@ -36,7 +36,7 @@ SPPOW, SpatialPower, SPGAU, SpatialGaussian, RawCoding, UN, Unstructured, CovarianceType, -fit, fit!, LMM, VarEffect, theta, logreml, m2logreml, thetalength, dof_satter, dof_contain, rankx, caic, lcontrast, typeiii, estimate, contrast, +fit, fit!, LMM, VarEffect, theta, logreml, m2logreml, m2logml, logml, thetalength, dof_satter, dof_contain, rankx, caic, lcontrast, typeiii, estimate, contrast, gmatrix, rmatrix, vmatrix!, responsename, nblocks, raneff, AbstractCovarianceType, AbstractCovmatMethod, MetidaModel, getlog, rand, rand!, diff --git a/src/utils.jl b/src/utils.jl index b714afa..b7443e7 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -195,15 +195,50 @@ function varlinkvecapply(v, p; varlinkf = :exp, rholinkf = :sigm) return s end ################################################################################ -function m2logreml(lmm) - return lmm.result.reml + +""" + m2logreml(lmm::LMM, θ = theta(lmm); maxthreads::Int = num_cores()) + +-2 logREML +""" +function m2logreml(lmm::LMM, θ = theta(lmm); maxthreads::Int = num_cores()) + if θ == theta(lmm) + return lmm.result.reml + else + return reml_sweep_β(lmm, LMMDataViews(lmm), θ; maxthreads = maxthreads)[1] + end +end +""" + logreml(lmm::LMM, θ = theta(lmm); maxthreads::Int = num_cores()) + +logREML +""" +function logreml(lmm::LMM, θ = theta(lmm); maxthreads::Int = num_cores()) + return -m2logreml(lmm, θ; maxthreads = maxthreads)/2 end -function logreml(lmm) - return -m2logreml(lmm)/2 + + +""" + m2logml(lmm::LMM, β = coef(lmm), θ = theta(lmm); maxthreads::Int = num_cores()) + +-2 logML +""" +function m2logml(lmm::LMM, β = coef(lmm), θ = theta(lmm); maxthreads::Int = num_cores()) + c = length(lmm.data.yv)*log(2π) + θ₁, θ₂, θ₃, noerror = core_sweep_β(lmm, LMMDataViews(lmm), θ, β, length(lmm.covstr.vcovblock); maxthreads = maxthreads) + if !noerror @warn "There was probably an error during the calculations. The result may be incorrect." end + return θ₁ + θ₃ + c end -function m2logreml(lmm, theta; maxthreads::Int = num_cores()) - return reml_sweep_β(lmm, LMMDataViews(lmm), theta; maxthreads = maxthreads)[1] +""" + logml(lmm::LMM, beta = coef(lmm), θ = theta(lmm); maxthreads::Int = num_cores()) + +logML +""" +function logml(lmm::LMM, beta = coef(lmm), θ = theta(lmm); maxthreads::Int = num_cores()) + return -m2logml(lmm, beta, θ; maxthreads = maxthreads)/2 end + + ################################################################################ function optim_callback(os) diff --git a/test/test.jl b/test/test.jl index 1572367..3ddab10 100644 --- a/test/test.jl +++ b/test/test.jl @@ -51,7 +51,44 @@ include("testdata.jl") random = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG), ) @test Metida.m2logreml(lmm) ≈ 16.241112644506067 atol=1E-6 + @test Metida.logreml(lmm) ≈ -0.5*16.241112644506067 atol=1E-6 + @test Metida.m2logreml(lmm, [0.447322, 0.367367, 0.185068]) ≈ 16.241112644506067 atol=1E-6 + @test Metida.m2logreml(lmm, [0.5, 0.3, 0.2]) ≈ 16.5746217198294 atol=1E-6 + + # core_sweep_β REML, ML estimate test; + reml, θs₂, remlθ₃, noerror = Metida.reml_sweep_β(lmm, Metida.LMMDataViews(lmm), Metida.theta(lmm), Metida.coef(lmm)) + @test reml ≈ 16.241112644506067 atol=1E-6 + @test θs₂ ≈ [55.8946 33.5368 13.4807 14.4666 13.4807 32.8767 + 33.5368 33.5368 6.90537 9.86301 6.90537 19.726 + 13.4807 6.90537 79.7331 0.0 -66.2523 6.57534 + 14.4666 9.86301 0.0 80.226 0.0 9.86301 + 13.4807 6.90537 -66.2523 0.0 79.7331 6.57534 + 32.8767 19.726 6.57534 9.86301 6.57534 32.8767] atol=1E-3 + @test remlθ₃ ≈ 13.999999919958947 atol=1E-6 + + θ₁, θ₂, θ₃, noerror = Metida.core_sweep_β(lmm, Metida.LMMDataViews(lmm), Metida.theta(lmm), Metida.coef(lmm), length(lmm.covstr.vcovblock)) + @test noerror == true + c = (length(lmm.data.yv) - lmm.rankx)*log(2π) + @test θ₁ ≈ -43.860020472151916 atol=1E-6 + @test θ₂ ≈ [55.8946 33.5368 13.4807 14.4666 13.4807 32.8767 + 0.0 33.5368 6.90537 9.86301 6.90537 19.726 + 0.0 0.0 79.7331 0.0 -66.2523 6.57534 + 0.0 0.0 0.0 80.226 0.0 9.86301 + 0.0 0.0 0.0 0.0 79.7331 6.57534 + 0.0 0.0 0.0 0.0 0.0 32.8767] atol=1E-3 + @test θ₃ ≈ 13.999999919958947 atol=1E-6 + logdetθ₂ = logdet(Symmetric(θ₂)) + @test logdetθ₂ ≈ 20.370854266968205 atol=1E-6 + @test θ₁ + logdetθ₂ + θ₃ + c ≈ 16.241112644506067 atol=1E-6 + # ML test + c = length(lmm.data.yv)*log(2π) + @test θ₁ + θ₃ + c ≈ 6.897520775993932 atol=1E-6 + @test Metida.m2logml(lmm, coef(lmm), Metida.theta(lmm); maxthreads = 8) ≈ 6.897520775993932 + @test Metida.m2logml(lmm, coef(lmm)) ≈ 6.897520775993932 + @test Metida.m2logml(lmm) ≈ 6.897520775993932 + @test Metida.logml(lmm) ≈ -0.5*6.897520775993932 + # lmm = Metida.fit(Metida.LMM, Metida.@lmmformula(var~0+sequence+period+formulation, random = formulation|subject:Metida.DIAG), df0) @test Metida.fixedeffn(lmm) == 3