diff --git a/.gitignore b/.gitignore index bee34859..b7dfeb7c 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,9 @@ src/temp.jl test/precompile/ test/devtest.jl +test/dstest.jl +test/snoopy.jl examples/ docs/build/ docs/img/ +etc/ diff --git a/Project.toml b/Project.toml index 15b1f2b4..a240aedb 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Metida" uuid = "a1dec852-9fe5-11e9-361f-8d9fde67cfa2" authors = ["Vladimir Arnautov "] -version = "0.1.14" +version = "0.2.0" [deps] CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597" @@ -29,7 +29,8 @@ julia = "1.1, 1.2, 1.3, 1.4, 1.5, 1.6" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "CSV", "DataFrames", "StatsBase"] +test = ["Test", "CSV", "DataFrames", "StatsBase", "StatsModels"] diff --git a/change.log b/change.log index 00f639a6..18f748fe 100644 --- a/change.log +++ b/change.log @@ -1,3 +1,11 @@ +- 0.2.0 + * many changes in VarianceType, covariance structure + * code cleaning + * change methods names + * rmatrix / gmatrix + * StatsBase methods + * test and other + - 0.1.13 * output fix diff --git a/docs/src/api.md b/docs/src/api.md index eaf2a4c6..128f539f 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -49,3 +49,13 @@ Metida.CompoundSymmetry ```@docs Metida.HeterogeneousCompoundSymmetry ``` + +### Metida.gmatrix +```@docs +Metida.gmatrix +``` + +### Metida.rmatrix +```@docs +Metida.rmatrix +``` diff --git a/docs/src/ref.md b/docs/src/ref.md index 1939dcf5..9a38efd5 100644 --- a/docs/src/ref.md +++ b/docs/src/ref.md @@ -3,18 +3,74 @@ ### Sweep algorithm * Sweep based on [SweepOperator.jl](https://github.com/joshday/SweepOperator.jl). Thanks to @joshday and @Hua-Zhou. + * http://hua-zhou.github.io/teaching/biostatm280-2019spring/slides/12-sweep/sweep.html More: -* Section 7.4-7.6 of Numerical Analysis for Statisticians by Kenneth Lange (2010). +* Section 7.4-7.6 of Numerical Analysis for Statisticians by Kenneth Lange (2010). + * The paper A tutorial on the SWEEP operator by James H. Goodnight (1979). -### REML +### REML & Parameter estimation + +* Henderson, C. R., et al. “The Estimation of Environmental and Genetic Trends from Records Subject to Culling.” Biometrics, vol. 15, no. 2, 1959, pp. 192–218. JSTOR, www.jstor.org/stable/2527669. + +* Lindstrom & J.; Bates, M. (1988). Newton—Raphson and EM Algorithms for Linear Mixed-Effects Models for Repeated-Measures Data. Journal of the American Statistical Association. 83. 1014. 10.1080/01621459.1988.10478693. +* Gurka, Matthew. (2006). Selecting the Best Linear Mixed Model under REML. The American Statistician. 60. 19-26. 10.1198/000313006X90396. + +* Sang Hong Lee, Julius H.J. van der Werf. An efficient variance component approach implementing an average information REML suitable for combined LD and linkage mapping with a general complex pedigree. Genetics Selection Evolution, BioMed Central, 2006, 38 (1), pp.25-43. ⟨hal-00894558⟩ + +* F.N. Gumedze, T.T. Dunne, Parameter estimation and inference in the linear mixed model, Linear Algebra and its Applications, Volume 435, Issue 8, 2011, Pages 1920-1944, ISSN 0024-3795, https://doi.org/10.1016/j.laa.2011.04.015. (http://www.sciencedirect.com/science/article/pii/S002437951100320X) ### AI algorithm +* D.L. Johnson, R. Thompson, Restricted Maximum Likelihood Estimation of Variance Components for Univariate Animal Models Using Sparse Matrix Techniques and Average Information, Journal of Dairy Science, Volume 78, Issue 2, 1995, Pages 449-456, ISSN 0022-0302, https://doi.org/10.3168/jds.S0022-0302(95)76654-1. (http://www.sciencedirect.com/science/article/pii/S0022030295766541) + +* Mishchenko, Kateryna & Holmgren, Sverker & Rönnegård, Lars. (2007). Newton-type Methods for REML Estimation in Genetic Analysis of Quantitative Traits. Journal of Computational Methods in Science and Engineering. 8. 10.3233/JCM-2008-81-203. + +* Matilainen K, Mäntysaari EA, Lidauer MH, Strandén I, Thompson R. Employing a Monte Carlo algorithm in Newton-type methods for restricted maximum likelihood estimation of genetic parameters. PLoS One. 2013;8(12):e80821. Published 2013 Dec 10. doi:10.1371/journal.pone.0080821 ### Covariance structures + +* Wolfinger, Russ. (1993). Covariance structure selection in general mixed models. Communications in Statistics-simulation and Computation - COMMUN STATIST-SIMULAT COMPUT. 22. 1079-1106. 10.1080/03610919308813143. + +* Littell, Ramon & Pendergast, Jane & Natarajan, Ranjini. (2000). Modelling covariance structure in the analysis of repeated measures data. Statistics in Medicine. 19. 1793-1819. 10.1002/1097-0258(20000715)19:13%3C1793::AID-SIM482%3E3.0.CO;2-Q. + +* QUINTAL, SILVANA SILVA RED, VIANA, ALEXANDRE PIO, CAMPOS, BIANCA MACHADO, VIVAS, MARCELO, & AMARAL JÚNIOR, ANTONIO TEIXEIRA DO. (2017). ANALYSIS OF STRUCTURES OF COVARIANCE AND REPEATABILITY IN GUAVA SEGREGANTING POPULATION. Revista Caatinga, 30(4), 885-891. https://doi.org/10.1590/1983-21252017v30n408rc + +* McNeish, D., Harring, J. Covariance pattern mixture models: Eliminating random effects to improve convergence and performance. Behav Res 52, 947–979 (2020). https://doi.org/10.3758/s13428-019-01292-4 + +### And more + +* Laird, Nan M., and James H. Ware. “Random-Effects Models for Longitudinal Data.” Biometrics, vol. 38, no. 4, 1982, pp. 963–974. JSTOR, www.jstor.org/stable/2529876. + +* Giesbrecht, F. G., and Burns, J. C. (1985), "Two-Stage Analysis Based on a Mixed Model: Large-sample Asymptotic Theory and Small-Sample Simulation Results," Biometrics, 41, 853-862. + +* Jennrich, R., & Schluchter, M. (1986). Unbalanced Repeated-Measures Models with Structured Covariance Matrices. Biometrics, 42(4), 805-820. doi:10.2307/2530695 + +* Fletcher, Roger (1987), Practical methods of optimization (2nd ed.), New York: John Wiley & Sons, ISBN 978-0-471-91547-8 + +* Wolfinger et al., (1994) Computing gaussian likelihoods and their derivatives for general linear mixed models doi: 10.1137/0915079 + +* Hrong-Tai Fai & Cornelius (1996) Approximate F-tests of multiple degree of freedom hypotheses in generalized least squares analyses of unbalanced split-plot experiments, Journal of Statistical Computation and Simulation, 54:4, 363-378, DOI: 10.1080/00949659608811740 + +* Schaalje GB, McBride JB, Fellingham GW. Adequacy of approximations to distributions of test statistics in complex mixed linear models. J Agric Biol Environ Stat. 2002;7:512–24. + +* Wright, Stephen, and Jorge Nocedal (2006) "Numerical optimization." Springer + +* Van Peer, A. (2010), Variability and Impact on Design of Bioequivalence Studies. Basic & Clinical Pharmacology & Toxicology, 106: 146-153. doi:10.1111/j.1742-7843.2009.00485.x + +### Julia packages + +* Revels, Jarrett & Lubin, Miles & Papamarkou, Theodore. (2016). Forward-Mode Automatic Differentiation in Julia. + +* Mogensen et al., (2018). Optim: A mathematical optimization package for Julia. Journal of Open Source Software, 3(24), 615,doi: 10.21105/joss.00615 + +### CuSolver & CuBLAS + +* https://docs.nvidia.com/cuda/cusolver/index.html + +* https://docs.nvidia.com/cuda/cublas/index.htm diff --git a/docs/src/validation.md b/docs/src/validation.md index a9daa2c7..7051573d 100644 --- a/docs/src/validation.md +++ b/docs/src/validation.md @@ -1,6 +1,81 @@ ## Validation +### sleepstudy.csv +| Model | REML Metida | REML SPSS| +|--------|--------|-------| +| 1 | 1729.4925602367025 | 1729.492560 | +| 2 | 1904.3265170722132 | 1662.172084 | +| 3 | 1772.0953251997046 | 1772.095 | +| 4 | 1730.1895427398322 | 1730.189543 | + + +#### Model 1 + +``` +lmm = Metida.LMM(@formula(Reaction~Days), df; + random = Metida.VarEffect(Metida.SI), + subject = :Subject + ) + Metida.fit!(lmm) +``` + +``` +``` +#### Model 2 + +``` +lmm = Metida.LMM(@formula(Reaction~1), df; + random = Metida.VarEffect(Metida.@covstr(Days), Metida.CS), + subject = :Subject + ) + Metida.fit!(lmm) +``` + +``` +``` + +#### Model 3 + +``` +lmm = Metida.LMM(@formula(Reaction~1), df; + random = Metida.VarEffect(Metida.@covstr(Days), Metida.CSH, subj = :Subject) + ) + Metida.fit!(lmm) +``` + +``` +``` + +#### Model 4 + +``` +lmm = Metida.LMM(@formula(Reaction~1), df; + random = Metida.VarEffect(Metida.@covstr(Days), Metida.ARH, subj = :Subject) + ) + Metida.fit!(lmm) +``` + +``` +``` + +### Pastes.csv + +| Model | REML Metida | REML SPSS| +|--------|--------|-------| +| 5 | 246.99074585348623 | 246.990746 | + +#### Model 5 + +``` +lmm = Metida.LMM(@formula(strength~1), df; +random = [Metida.VarEffect(Metida.SI, subj = :batch), Metida.VarEffect(Metida.SI, subj = [:batch, :cask])] +) +Metida.fit!(lmm) +``` + +``` +``` ### Reference dataset @@ -20,4 +95,4 @@ Anglin, P.M. and R. Gencay (1996) “Semiparametric estimation of a hedonic pric Anglin, P., and Gencay, R. (1996). Semiparametric Estimation of a Hedonic Price Function. Journal of Applied Econometrics, 11, 633–648. Verbeek, M. (2004). A Guide to Modern Econometrics, 2nd ed. Chichester, UK: John Wiley. -O'Brien, R. G., and Kaiser, M. K. (1985) MANOVA method for analyzing repeated measures designs: An extensive primer. Psychological Bulletin 97, 316–333, Table 7. +O'Brien, R. G., and Kaiser, M. K. (1985) MANOVA method for analyzing repeated measures designs: An extensive primer. Psychological Bulletin 97, 316–333, Table 7. diff --git a/src/Metida.jl b/src/Metida.jl index d31aa2fc..f32e4424 100644 --- a/src/Metida.jl +++ b/src/Metida.jl @@ -25,7 +25,6 @@ fit!, LMM, VarEffect, theta, logreml, m2logreml, thetalength include("abstracttype.jl") include("sweep.jl") -include("utils.jl") include("varstruct.jl") include("gmat.jl") include("rmat.jl") @@ -40,6 +39,7 @@ include("fit.jl") include("showutils.jl") include("statsbase.jl") include("initg.jl") +include("utils.jl") function __init__() @@ -48,16 +48,24 @@ end function _precompile_() ccall(:jl_generating_output, Cint, ()) == 1 || return nothing isdefined(Metida, Symbol("#@covstr")) && precompile(Tuple{getfield(Metida, Symbol("#@covstr")), LineNumberNode, Module, Int}) - precompile(Tuple{typeof(Metida.ffx), Int64}) - precompile(Tuple{typeof(Metida.ffxpone), Int64}) precompile(Tuple{typeof(Metida.reml_sweep_β), Metida.LMM{Float64}, Array{Float64, 1}}) precompile(Tuple{typeof(Metida.reml_sweep_β_b), Metida.LMM{Float64}, Array{Float64, 1}}) precompile(Tuple{typeof(Metida.rholinkpsigmoid), Float64}) + precompile(Tuple{typeof(Metida.rholinkpsigmoidr), Float64}) + + precompile(Tuple{typeof(Metida.varlinkvecapply!), Array{Float64, 1}, Array{Symbol, 1}}) + precompile(Tuple{typeof(Metida.varlinkrvecapply!), Array{Float64, 1}, Array{Symbol, 1}}) + precompile(Tuple{typeof(Metida.varlinkvecapply), Array{Float64, 1}, Array{Symbol, 1}}) - precompile(Tuple{typeof(Metida.varlinkvecapply!), Array{Float64, 1}, Array{Function, 1}}) precompile(Tuple{typeof(Metida.vlink), Float64}) + precompile(Tuple{typeof(Metida.vlinkr), Float64}) + + precompile(Tuple{typeof(Metida.initvar), Array{Float64, 1}, Array{Float64, 2}}) + + precompile(Tuple{typeof(Metida.intersectsubj), Array{VarEffect, 1}, VarEffect}) + end _precompile_() #include(".jl") diff --git a/src/deprecated.jl b/src/deprecated.jl index abc0659e..6b8fb760 100644 --- a/src/deprecated.jl +++ b/src/deprecated.jl @@ -179,3 +179,430 @@ function fgh!(F, G, H, x; remlβcalc::Function, remlcalc::Function) nothing end =# + +""" + Make X, Z matrices and vector y for each subject; +""" +#= +function subjblocks(df, sbj::Symbol, x::Matrix{T}, z::Matrix{T}, y::Vector{T}, rz::Matrix{T}) where T + u = unique(df[!, sbj]) + xa = Vector{Matrix{T}}(undef, length(u)) + za = Vector{Matrix{T}}(undef, length(u)) + ya = Vector{Vector{T}}(undef, length(u)) + rza = Vector{Matrix{T}}(undef, length(u)) + @inbounds @simd for i = 1:length(u) + v = findall(x -> x == u[i], df[!, sbj]) + xa[i] = Matrix(view(x, v, :)) + za[i] = Matrix(view(z, v, :)) + ya[i] = Vector(view(y, v)) + rza[i] = Matrix(view(rz, v, :)) + end + return xa, za, rza, ya +end +function subjblocks(df, sbj::Symbol, x::Matrix{T}, z::Matrix{T}, y::Vector{T}, rz::Nothing) where T + u = unique(df[!, sbj]) + xa = Vector{Matrix{T}}(undef, length(u)) + za = Vector{Matrix{T}}(undef, length(u)) + ya = Vector{Vector{T}}(undef, length(u)) + rza = Vector{Matrix{T}}(undef, length(u)) + @inbounds @simd for i = 1:length(u) + v = findall(x -> x == u[i], df[!, sbj]) + xa[i] = Matrix(view(x, v, :)) + za[i] = Matrix(view(z, v, :)) + ya[i] = Vector(view(y, v)) + rza[i] = Matrix{T}(undef, 0, 0) + end + return xa, za, rza, ya +end +function subjblocks(df, sbj) + if isa(sbj, Symbol) + u = unique(df[!, sbj]) + r = Vector{Vector{Int}}(undef, length(u)) + @inbounds @simd for i = 1:length(u) + r[i] = findall(x -> x == u[i], df[!, sbj]) + end + return r + end + r = Vector{Vector{Int}}(undef, 1) + r[1] = collect(1:size(df, 1)) + r +end +=# + + +#= +function intersectsubj(covstr) + a = Vector{Vector{Symbol}}(undef, length(covstr.random)+1) + eq = true + for i = 1:length(covstr.random) + a[i] = covstr.random[i].subj + end + a[end] = covstr.repeated.subj + for i = 2:length(a) + if !(issetequal(a[1], a[i])) + eq = false + break + end + end + intersect(a...), eq +end +=# +#= +function diffsubj!(a, subj) + push!(a, subj) + symdiff(a...) +end +=# +#= +function varlinkvecapply(v, f) + rv = similar(v) + for i = 1:length(v) + rv[i] = f[i](v[i]) + end + rv +end + +function varlinkvecapply!(v, f) + for i = 1:length(v) + v[i] = f[i](v[i]) + end + v +end +=# +#varlinkvecapply(v, f) = map.(f,v) #Test + +#= +function varlinkrvecapply2(v, p; varlinkf = :exp, rholinkf = :sigm) + s = similar(v) + for i = 1:length(v) + if p[i] == :var + s[i] = vlinkr(v[i]) + else + s[i] = rholinksigmoidr(v[i]) + end + end + s +end +=# +################################################################################ +#= +function vmatr(lmm, i) + θ = lmm.result.theta + G = gmat_base(θ, lmm.covstr) + V = mulαβαt(view(lmm.data.zv, lmm.data.block[i],:), G) + if length(lmm.data.zrv) > 0 + rmat_basep!(V, θ[lmm.covstr.tr[end]], view(lmm.data.zrv, lmm.data.block[i],:), lmm.covstr) + else + rmat_basep!(V, θ[lmm.covstr.tr[end]], lmm.data.zrv, lmm.covstr) + end + V +end + +function gmatr(lmm, i) + θ = lmm.result.theta + gmat_base(θ, lmm.covstr) +end +=# + +#= +function rmat_basep_z!(mx, θ::AbstractVector{T}, zrv, covstr) where T + for i = 1:length(covstr.block[end]) + if covstr.repeated.covtype.s == :SI + rmatp_si!(view(mx, covstr.block[end][i], covstr.block[end][i]), θ, view(covstr.rz, covstr.block[end][i], :), covstr.repeated.covtype) + elseif covstr.repeated.covtype.s == :DIAG + rmatp_diag!(view(mx, covstr.block[end][i], covstr.block[end][i]), θ, view(covstr.rz, covstr.block[end][i], :), covstr.repeated.covtype) + elseif covstr.repeated.covtype.s == :AR + rmatp_ar!(view(mx, covstr.block[end][i], covstr.block[end][i]), θ, view(covstr.rz, covstr.block[end][i], :), covstr.repeated.covtype) + elseif covstr.repeated.covtype.s == :ARH + rmatp_arh!(view(mx, covstr.block[end][i], covstr.block[end][i]), θ, view(covstr.rz, covstr.block[end][i], :), covstr.repeated.covtype) + elseif covstr.repeated.covtype.s == :CSH + rmatp_csh!(view(mx, covstr.block[end][i], covstr.block[end][i]), θ, view(covstr.rz, covstr.block[end][i], :), covstr.repeated.covtype) + elseif covstr.repeated.covtype.s == :CS + rmatp_cs!(view(mx, covstr.block[end][i], covstr.block[end][i]), θ, view(covstr.rz, covstr.block[end][i], :), covstr.repeated.covtype) + else + error("Unknown covariance structure!") + end + end + mx +end +=# + +""" + -2 log Restricted Maximum Likelihood; +""" +#= +function reml_sweep(lmm, β, θ::Vector{T})::T where T <: Number + n = length(lmm.data.yv) + N = sum(length.(lmm.data.yv)) + G = gmat_base(θ, lmm.covstr) + c = (N - lmm.rankx)*log(2π) + θ₁ = zero(eltype(θ)) + θ₂ = zero(eltype(θ)) + θ₃ = zero(eltype(θ)) + θ2m = zeros(eltype(θ), lmm.rankx, lmm.rankx) + @inbounds for i = 1:n + q = length(lmm.data.yv[i]) + Vp = mulαβαt3(lmm.data.zv[i], G, lmm.data.xv[i]) + V = view(Vp, 1:q, 1:q) + rmat_basep!(V, θ[lmm.covstr.tr[end]], lmm.data.zrv[i], lmm.covstr) + + θ₁ += logdet(V) + sweep!(Vp, 1:q) + iV = Symmetric(utriaply!(x -> -x, Vp[1:q, 1:q])) + mulαtβαinc!(θ2m, lmm.data.xv[i], iV) + θ₃ += -Vp[end, end] + end + θ₂ = logdet(θ2m) + return -(θ₁ + θ₂ + θ₃ + c) +end +=# + +#= +function gfunc!(g, x, f) + chunk = ForwardDiff.Chunk{1}() + gcfg = ForwardDiff.GradientConfig(f, x, chunk) + ForwardDiff.gradient!(g, f, x, gcfg) +end +function hfunc!(h, x, f) + chunk = ForwardDiff.Chunk{1}() + hcfg = ForwardDiff.HessianConfig(f, x, chunk) + ForwardDiff.hessian!(h, f, x, hcfg) +end +=# + +""" +A * B * A' + C +""" +#= +function mulαβαtc(A, B, C) + q = size(B, 1) + p = size(A, 1) + c = zeros(eltype(B), q) + mx = zeros(eltype(B), p, p) + for i = 1:p + fill!(c, zero(eltype(c))) + @simd for n = 1:q + @simd for m = 1:q + @inbounds c[n] += A[i, m] * B[n, m] + end + end + @simd for n = i:p + @simd for m = 1:q + @inbounds mx[i, n] += A[n, m] * c[m] + end + @inbounds mx[i, n] += C[i, n] + end + end + Symmetric(mx) +end +=# +""" +A * B * A' +""" +#= +function mulαβαt(A, B) + q = size(B, 1) + p = size(A, 1) + c = zeros(eltype(B), q) + mx = zeros(eltype(B), p, p) + for i = 1:p + fill!(c, zero(eltype(c))) + @simd for n = 1:q + @simd for m = 1:q + @inbounds c[n] += A[i, m] * B[n, m] + end + end + @simd for n = i:p + @simd for m = 1:q + @inbounds mx[i, n] += A[n, m] * c[m] + end + end + end + Symmetric(mx) +end +function mulαβαtc2(A, B, C, r::Vector) + q = size(B, 1) + p = size(A, 1) + c = zeros(eltype(B), q) + mx = zeros(eltype(B), p + 1, p + 1) + for i = 1:p + fill!(c, zero(eltype(c))) + @simd for n = 1:q + @simd for m = 1:q + @inbounds c[n] += A[i, m] * B[n, m] + end + end + @simd for n = 1:p + @simd for m = 1:q + @inbounds mx[i, n] += A[n, m] * c[m] + end + @inbounds mx[i, n] += C[i, n] + end + mx[end, i] = r[i] + mx[i, end] = r[i] + end + mx +end +function mulαβαtc3(A, B, C, X) + q = size(B, 1) + p = size(A, 1) + c = zeros(eltype(B), q) + mx = zeros(eltype(B), p + size(X, 2), p + size(X, 2)) + for i = 1:p + fill!(c, zero(eltype(c))) + @simd for n = 1:q + @simd for m = 1:q + @inbounds c[n] += A[i, m] * B[n, m] + end + end + @simd for n = 1:p + @simd for m = 1:q + @inbounds mx[i, n] += A[n, m] * c[m] + end + @inbounds mx[i, n] += C[i, n] + end + end + mx[1:p, p+1:end] = X + mx[p+1:end, 1:p] = X' + mx +end +=# + +""" +A' * B * A -> +θ +A' * B * C -> +β +""" +#= +function mulθβinc!(θ, β, A::AbstractMatrix, B::AbstractMatrix, C::AbstractVector) + q = size(B, 1) + p = size(A, 2) + c = Vector{eltype(θ)}(undef, q) + for i = 1:p + fill!(c, zero(eltype(θ))) + for n = 1:q + for m = 1:q + @inbounds c[n] += B[m, n] * A[m, i] + end + @inbounds β[i] += c[n] * C[n] + end + for n = 1:p + for m = 1:q + @inbounds θ[i, n] += A[m, n] * c[m] + end + end + end + θ, β +end +=# +#------------------------------------------------------------------------------- +""" +A' * B * A -> + θ +""" +#= +function mulαtβαinc!(θ, A, B) + q = size(B, 1) + p = size(A, 2) + c = zeros(eltype(B), q) + for i = 1:p + fill!(c, zero(eltype(c))) + @inbounds for n = 1:q, m = 1:q + c[n] += B[m, n] * A[m, i] + end + @inbounds for n = 1:p, m = 1:q + θ[i, n] += A[m, n] * c[m] + end + end +end +=# + +""" +A' * B * A -> θ +""" +#= +function mulαtβα!(θ, A, B) + q = size(B, 1) + p = size(A, 2) + c = zeros(eltype(B), q) + fill!(θ, zero(eltype(θ))) + for i = 1:p + fill!(c, zero(eltype(c))) + @inbounds for n = 1:q, m = 1:q + c[n] += B[m, n] * A[m, i] + end + @inbounds for n = 1:p, m = 1:q + θ[i, n] += A[m, n] * c[m] + end + end + θ +end +=# +""" +A * B * A -> θ +""" +#= +function mulαβα!(θ, A, B) + q = size(B, 1) + p = size(A, 2) + c = zeros(eltype(B), q) + fill!(θ, zero(eltype(θ))) + for i = 1:p + fill!(c, zero(eltype(c))) + @inbounds for n = 1:q, m = 1:q + c[n] += B[m, n] * A[i, m] + end + @inbounds for n = 1:p, m = 1:q + θ[i, n] += A[m, n] * c[m] + end + end + θ +end +=# +""" +tr(A * B) +""" +#= +function trmulαβ(A, B) + c = 0 + @inbounds for n = 1:size(A,1), m = 1:size(B, 1) + c += A[n,m] * B[m, n] + end + c +end +=# +""" +tr(H * A' * B * A) +""" +#= +function trmulhαtβα(H, A, B) +end +=# + +""" +(y - X * β) +""" +#= +function mulr!(v::AbstractVector, y::AbstractVector, X::AbstractMatrix, β::AbstractVector) + fill!(v, zero(eltype(v))) + q = length(y) + p = size(X, 2) + @simd for n = 1:q + @simd for m = 1:p + @inbounds v[n] += X[n, m] * β[m] + end + v[n] = y[n] - v[n] + end + return v +end +function mulr(y::AbstractVector, X::AbstractMatrix, β::AbstractVector) + v = zeros(eltype(β), length(y)) + q = length(y) + p = size(X, 2) + @simd for n = 1:q + @simd for m = 1:p + @inbounds v[n] += X[n, m] * β[m] + end + v[n] = y[n] - v[n] + end + return v +end +=# diff --git a/src/fit.jl b/src/fit.jl index 3e829a1c..c29f108c 100644 --- a/src/fit.jl +++ b/src/fit.jl @@ -1,16 +1,5 @@ # fit.jl -#= -function gfunc!(g, x, f) - chunk = ForwardDiff.Chunk{1}() - gcfg = ForwardDiff.GradientConfig(f, x, chunk) - ForwardDiff.gradient!(g, f, x, gcfg) -end -function hfunc!(h, x, f) - chunk = ForwardDiff.Chunk{1}() - hcfg = ForwardDiff.HessianConfig(f, x, chunk) - ForwardDiff.hessian!(h, f, x, hcfg) -end -=# + fit_nlopt!(lmm::MetidaModel; kwargs...) = error("MetidaNLopt not found. \n - Run `using MetidaNLopt` before.") """ @@ -53,9 +42,7 @@ function fit!(lmm::LMM{T}; if verbose == :auto verbose = 1 end - #Make varlink function - fv = varlinkvec(lmm.covstr.ct) - fvr = varlinkrvec(lmm.covstr.ct) + # Optimization function if lmm.blocksolve optfunc = reml_sweep_β_b @@ -91,7 +78,7 @@ function fit!(lmm::LMM{T}; copyto!(θ, init) lmmlog!(io, lmm, verbose, LMMLogMsg(:INFO, "Using provided θ: "*string(θ))) else - error("init length $(length(init)) != θ length $(length(init))") + error("init length $(length(init)) != θ length $(length(θ))") end else initθ = sqrt(initvar(lmm.mf.data[lmm.mf.f.lhs.sym], lmm.mm.m)[1]/4) @@ -133,11 +120,11 @@ function fit!(lmm::LMM{T}; end lmmlog!(io, lmm, verbose, LMMLogMsg(:INFO, "First step with AI-like method, θ: "*string(θ))) end - varlinkrvecapply2!(θ, lmm.covstr.ct) + varlinkrvecapply!(θ, lmm.covstr.ct) #Twice differentiable object - vloptf(x) = optfunc(lmm, varlinkvecapply2(x, lmm.covstr.ct))[1] + vloptf(x) = optfunc(lmm, varlinkvecapply(x, lmm.covstr.ct))[1] chunk = ForwardDiff.Chunk{1}() gcfg = ForwardDiff.GradientConfig(vloptf, θ, chunk) hcfg = ForwardDiff.HessianConfig(vloptf, θ, chunk) @@ -156,7 +143,7 @@ function fit!(lmm::LMM{T}; lmm.result.optim = Optim.optimize(td, θ, optmethod, optoptions) end #Theta (θ) vector - lmm.result.theta = varlinkvecapply2!(deepcopy(Optim.minimizer(lmm.result.optim)), lmm.covstr.ct) + lmm.result.theta = varlinkvecapply!(deepcopy(Optim.minimizer(lmm.result.optim)), lmm.covstr.ct) lmmlog!(io, lmm, verbose, LMMLogMsg(:INFO, "Resulting θ: "*string(lmm.result.theta))) try if hcalck diff --git a/src/gmat.jl b/src/gmat.jl index e232fe29..726c0448 100644 --- a/src/gmat.jl +++ b/src/gmat.jl @@ -12,21 +12,21 @@ function gmat_base(θ::Vector{T}, covstr) where T mx end ################################################################################ -function gmat_switch!(G, θ, covstr, i) - if covstr.random[i].covtype.s == :SI - gmat_si!(G, θ[covstr.tr[i]], covstr.q[i], covstr.random[i].covtype) # i > r - elseif covstr.random[i].covtype.s == :DIAG - gmat_diag!(G, θ[covstr.tr[i]], covstr.q[i], covstr.random[i].covtype) - elseif covstr.random[i].covtype.s == :AR - gmat_ar!(G, θ[covstr.tr[i]], covstr.q[i], covstr.random[i].covtype) - elseif covstr.random[i].covtype.s == :ARH - gmat_arh!(G, θ[covstr.tr[i]], covstr.q[i], covstr.random[i].covtype) - elseif covstr.random[i].covtype.s == :CSH - gmat_csh!(G, θ[covstr.tr[i]], covstr.q[i], covstr.random[i].covtype) - elseif covstr.random[i].covtype.s == :CS - gmat_cs!(G, θ[covstr.tr[i]], covstr.q[i], covstr.random[i].covtype) - elseif covstr.random[i].covtype.s == :ZERO - gmat_zero!(G, similar(θ, 0), covstr.q[i], covstr.random[i].covtype) +function gmat_switch!(G, θ, covstr, r) + if covstr.random[r].covtype.s == :SI + gmat_si!(G, θ[covstr.tr[r]], covstr.q[r], covstr.random[r].covtype) # i > r + elseif covstr.random[r].covtype.s == :DIAG + gmat_diag!(G, θ[covstr.tr[r]], covstr.q[r], covstr.random[r].covtype) + elseif covstr.random[r].covtype.s == :AR + gmat_ar!(G, θ[covstr.tr[r]], covstr.q[r], covstr.random[r].covtype) + elseif covstr.random[r].covtype.s == :ARH + gmat_arh!(G, θ[covstr.tr[r]], covstr.q[r], covstr.random[r].covtype) + elseif covstr.random[r].covtype.s == :CSH + gmat_csh!(G, θ[covstr.tr[r]], covstr.q[r], covstr.random[r].covtype) + elseif covstr.random[r].covtype.s == :CS + gmat_cs!(G, θ[covstr.tr[r]], covstr.q[r], covstr.random[r].covtype) + elseif covstr.random[r].covtype.s == :ZERO + gmat_zero!(G, similar(θ, 0), covstr.q[r], covstr.random[r].covtype) else error("Unknown covariance structure!") end @@ -47,7 +47,7 @@ function gmat_base_z!(mx, θ::Vector{T}, covstr) where T end =# ################################################################################ -function gmat_base_z2!(mx, θ::Vector{T}, covstr, block, sblock) where T +function zgz_base_inc!(mx, θ::Vector{T}, covstr, block, sblock) where T q = sum(length.(covstr.block[1])) for r = 1:length(covstr.random) G = zeros(T, covstr.q[r], covstr.q[r]) @@ -66,7 +66,7 @@ function gmat_zero!(mx, θ::Vector{T}, ::Int, ::CovarianceType) where T mx .= zero(T) nothing end -function gmat_si!(mx, θ::Vector{T}, zn::Int, ::CovarianceType) where T +function gmat_si!(mx, θ::Vector{T}, ::Int, ::CovarianceType) where T val = θ[1] ^ 2 for i = 1:size(mx, 1) mx[i, i] = val @@ -79,14 +79,15 @@ function gmat_diag!(mx, θ::Vector{T}, ::Int, ::CovarianceType) where T end nothing end -function gmat_ar!(mx, θ::Vector{T}, zn::Int, ::CovarianceType) where T +function gmat_ar!(mx, θ::Vector{T}, ::Int, ::CovarianceType) where T de = θ[1] ^ 2 - for i = 1:zn + s = size(mx, 1) + for i = 1:s mx[i, i] = de end - if zn > 1 - for m = 1:zn - 1 - for n = m + 1:zn + if s > 1 + for m = 1:s - 1 + for n = m + 1:s ode = de * θ[2] ^ (n - m) @inbounds mx[m, n] = ode @inbounds mx[n, m] = mx[m, n] diff --git a/src/initg.jl b/src/initg.jl index 78835846..108ec746 100644 --- a/src/initg.jl +++ b/src/initg.jl @@ -1,5 +1,4 @@ - function initgstep(f, θ) g = ForwardDiff.gradient(f, θ) he = ForwardDiff.hessian(f, θ) diff --git a/src/linearalgebra.jl b/src/linearalgebra.jl index d9758cef..8e6d87df 100644 --- a/src/linearalgebra.jl +++ b/src/linearalgebra.jl @@ -1,101 +1,7 @@ #linearalgebra.jl """ -A * B * A' + C + """ -#= -function mulαβαtc(A, B, C) - q = size(B, 1) - p = size(A, 1) - c = zeros(eltype(B), q) - mx = zeros(eltype(B), p, p) - for i = 1:p - fill!(c, zero(eltype(c))) - @simd for n = 1:q - @simd for m = 1:q - @inbounds c[n] += A[i, m] * B[n, m] - end - end - @simd for n = i:p - @simd for m = 1:q - @inbounds mx[i, n] += A[n, m] * c[m] - end - @inbounds mx[i, n] += C[i, n] - end - end - Symmetric(mx) -end -=# -""" -A * B * A' -""" -#= -function mulαβαt(A, B) - q = size(B, 1) - p = size(A, 1) - c = zeros(eltype(B), q) - mx = zeros(eltype(B), p, p) - for i = 1:p - fill!(c, zero(eltype(c))) - @simd for n = 1:q - @simd for m = 1:q - @inbounds c[n] += A[i, m] * B[n, m] - end - end - @simd for n = i:p - @simd for m = 1:q - @inbounds mx[i, n] += A[n, m] * c[m] - end - end - end - Symmetric(mx) -end -function mulαβαtc2(A, B, C, r::Vector) - q = size(B, 1) - p = size(A, 1) - c = zeros(eltype(B), q) - mx = zeros(eltype(B), p + 1, p + 1) - for i = 1:p - fill!(c, zero(eltype(c))) - @simd for n = 1:q - @simd for m = 1:q - @inbounds c[n] += A[i, m] * B[n, m] - end - end - @simd for n = 1:p - @simd for m = 1:q - @inbounds mx[i, n] += A[n, m] * c[m] - end - @inbounds mx[i, n] += C[i, n] - end - mx[end, i] = r[i] - mx[i, end] = r[i] - end - mx -end -function mulαβαtc3(A, B, C, X) - q = size(B, 1) - p = size(A, 1) - c = zeros(eltype(B), q) - mx = zeros(eltype(B), p + size(X, 2), p + size(X, 2)) - for i = 1:p - fill!(c, zero(eltype(c))) - @simd for n = 1:q - @simd for m = 1:q - @inbounds c[n] += A[i, m] * B[n, m] - end - end - @simd for n = 1:p - @simd for m = 1:q - @inbounds mx[i, n] += A[n, m] * c[m] - end - @inbounds mx[i, n] += C[i, n] - end - end - mx[1:p, p+1:end] = X - mx[p+1:end, 1:p] = X' - mx -end -=# function mulαβαt3(A, B, X) q = size(B, 1) p = size(A, 1) @@ -119,52 +25,6 @@ function mulαβαt3(A, B, X) mx end """ -A' * B * A -> +θ -A' * B * C -> +β -""" -#= -function mulθβinc!(θ, β, A::AbstractMatrix, B::AbstractMatrix, C::AbstractVector) - q = size(B, 1) - p = size(A, 2) - c = Vector{eltype(θ)}(undef, q) - for i = 1:p - fill!(c, zero(eltype(θ))) - for n = 1:q - for m = 1:q - @inbounds c[n] += B[m, n] * A[m, i] - end - @inbounds β[i] += c[n] * C[n] - end - for n = 1:p - for m = 1:q - @inbounds θ[i, n] += A[m, n] * c[m] - end - end - end - θ, β -end -=# -#------------------------------------------------------------------------------- -""" -A' * B * A -> + θ -""" -#= -function mulαtβαinc!(θ, A, B) - q = size(B, 1) - p = size(A, 2) - c = zeros(eltype(B), q) - for i = 1:p - fill!(c, zero(eltype(c))) - @inbounds for n = 1:q, m = 1:q - c[n] += B[m, n] * A[m, i] - end - @inbounds for n = 1:p, m = 1:q - θ[i, n] += A[m, n] * c[m] - end - end -end -=# -""" A * B * A' -> + θ """ function mulαβαtinc!(θ, A, B) @@ -182,69 +42,6 @@ function mulαβαtinc!(θ, A, B) end θ end -#------------------------------------------------------------------------------- - -""" -A' * B * A -> θ -""" -#= -function mulαtβα!(θ, A, B) - q = size(B, 1) - p = size(A, 2) - c = zeros(eltype(B), q) - fill!(θ, zero(eltype(θ))) - for i = 1:p - fill!(c, zero(eltype(c))) - @inbounds for n = 1:q, m = 1:q - c[n] += B[m, n] * A[m, i] - end - @inbounds for n = 1:p, m = 1:q - θ[i, n] += A[m, n] * c[m] - end - end - θ -end -=# -""" -A * B * A -> θ -""" -#= -function mulαβα!(θ, A, B) - q = size(B, 1) - p = size(A, 2) - c = zeros(eltype(B), q) - fill!(θ, zero(eltype(θ))) - for i = 1:p - fill!(c, zero(eltype(c))) - @inbounds for n = 1:q, m = 1:q - c[n] += B[m, n] * A[i, m] - end - @inbounds for n = 1:p, m = 1:q - θ[i, n] += A[m, n] * c[m] - end - end - θ -end -=# -""" -tr(A * B) -""" -#= -function trmulαβ(A, B) - c = 0 - @inbounds for n = 1:size(A,1), m = 1:size(B, 1) - c += A[n,m] * B[m, n] - end - c -end -=# -""" -tr(H * A' * B * A) -""" -#= -function trmulhαtβα(H, A, B) -end -=# """ (y - X * β)' * V * (y - X * β) """ @@ -266,35 +63,6 @@ function mulθ₃(y::AbstractVector, X::AbstractMatrix, β::AbstractVector, V::A return θ end """ -(y - X * β) -""" -#= -function mulr!(v::AbstractVector, y::AbstractVector, X::AbstractMatrix, β::AbstractVector) - fill!(v, zero(eltype(v))) - q = length(y) - p = size(X, 2) - @simd for n = 1:q - @simd for m = 1:p - @inbounds v[n] += X[n, m] * β[m] - end - v[n] = y[n] - v[n] - end - return v -end -function mulr(y::AbstractVector, X::AbstractMatrix, β::AbstractVector) - v = zeros(eltype(β), length(y)) - q = length(y) - p = size(X, 2) - @simd for n = 1:q - @simd for m = 1:p - @inbounds v[n] += X[n, m] * β[m] - end - v[n] = y[n] - v[n] - end - return v -end -=# -""" A' * B -> + θ """ function mulαtβinc!(θ, A::AbstractMatrix, B::AbstractVector) @@ -310,7 +78,6 @@ function mulαtβinc!(θ, A::AbstractMatrix, B::AbstractVector) end ################################################################################ - function utriaply!(f, m) if size(m, 1) != size(m, 2) error() end @simd for p = 1:size(m, 1) @@ -320,14 +87,3 @@ function utriaply!(f, m) end m end - - -#= -d = 0.4 -akk = [1.0 2.0 4.0 ; 6.0 1.0 3.0; 3 4 5] -b = 0.3 -A = [2.0 5.0 2.3; 6.0 6.0 4.5 ; 1 2 3] -BLAS.syrk!('U', 'N', -d, akk, b, A) --d*akk*akk'+b*A -nsyrk!(-d, akk, b, A) -=# diff --git a/src/lmm.jl b/src/lmm.jl index c8ff744b..540dca41 100644 --- a/src/lmm.jl +++ b/src/lmm.jl @@ -27,6 +27,7 @@ struct LMM{T} <: MetidaModel mm::ModelMatrix covstr::CovStructure{T} data::LMMData{T} + nfixed::Int rankx::UInt32 result::ModelResult blocksolve::Bool @@ -37,20 +38,16 @@ struct LMM{T} <: MetidaModel subject = Vector{Symbol}(undef, 0) elseif isa(subject, Symbol) subject = [subject] - elseif isa(subject, AbstractVector{Symbol}) - # - else - throw(ArgumentError("subject type should be Symbol or Vector{tymbol}")) end - lmmlog = Vector{LMMLogMsg}(undef, 0) - mf = ModelFrame(model, data; contrasts = contrasts) - mm = ModelMatrix(mf) + lmmlog = Vector{LMMLogMsg}(undef, 0) + mf = ModelFrame(model, data; contrasts = contrasts) + mm = ModelMatrix(mf) + nfixed = nterms(mf.f.rhs.terms) if repeated === nothing repeated = VarEffect() end if random === nothing random = VarEffect(Metida.@covstr(0), Metida.RZero(), subj = repeated.subj) - #random = RZero() end if !isa(random, Vector) random = [random] end #blocks @@ -67,7 +64,7 @@ struct LMM{T} <: MetidaModel block = intersectdf(data, subject) lmmdata = LMMData(mm.m, mf.data[mf.f.lhs.sym], block, subject) covstr = CovStructure(random, repeated, data, block) - new{eltype(mm.m)}(model, mf, mm, covstr, lmmdata, rank(mm.m), ModelResult(), blocksolve, lmmlog) + new{eltype(mm.m)}(model, mf, mm, covstr, lmmdata, nfixed, rank(mm.m), ModelResult(), blocksolve, lmmlog) end end ################################################################################ diff --git a/src/reml.jl b/src/reml.jl index 099de21b..7e21b292 100644 --- a/src/reml.jl +++ b/src/reml.jl @@ -24,7 +24,7 @@ function reml_sweep_β_b(lmm, θ::Vector{T}) where T <: Number q = length(lmm.data.block[i]) Vp = mulαβαt3(view(lmm.covstr.z, lmm.data.block[i],:), G, view(lmm.data.xv, lmm.data.block[i],:)) V = view(Vp, 1:q, 1:q) - rmat_basep!(V, θ[lmm.covstr.tr[end]], view(lmm.covstr.rz, lmm.data.block[i],:), lmm.covstr) + rmat_base_inc_b!(V, θ[lmm.covstr.tr[end]], view(lmm.covstr.rz, lmm.data.block[i],:), lmm.covstr) try θ₁ += logdet(cholesky(V)) catch @@ -76,8 +76,8 @@ function reml_sweep_β(lmm, θ::Vector{T}) where T <: Number V = view(Vp, 1:q, 1:q) Vx = view(Vp, 1:q, q+1:q+lmm.rankx) Vx .= view(lmm.data.xv, lmm.data.block[i],:) - gmat_base_z2!(V, θ, lmm.covstr, lmm.data.block[i], lmm.covstr.sblock[i]) - rmat_basep_z2!(V, θ[lmm.covstr.tr[end]], lmm.covstr, lmm.data.block[i], lmm.covstr.sblock[i]) + zgz_base_inc!(V, θ, lmm.covstr, lmm.data.block[i], lmm.covstr.sblock[i]) + rmat_base_inc!(V, θ[lmm.covstr.tr[end]], lmm.covstr, lmm.data.block[i], lmm.covstr.sblock[i]) #----------------------------------------------------------------------- θ₁ += logdet(cholesky(V)) sweep!(Vp, 1:q) @@ -97,32 +97,3 @@ function reml_sweep_β(lmm, θ::Vector{T}) where T <: Number end ################################################################################ ################################################################################ -""" - -2 log Restricted Maximum Likelihood; -""" -#= -function reml_sweep(lmm, β, θ::Vector{T})::T where T <: Number - n = length(lmm.data.yv) - N = sum(length.(lmm.data.yv)) - G = gmat_base(θ, lmm.covstr) - c = (N - lmm.rankx)*log(2π) - θ₁ = zero(eltype(θ)) - θ₂ = zero(eltype(θ)) - θ₃ = zero(eltype(θ)) - θ2m = zeros(eltype(θ), lmm.rankx, lmm.rankx) - @inbounds for i = 1:n - q = length(lmm.data.yv[i]) - Vp = mulαβαt3(lmm.data.zv[i], G, lmm.data.xv[i]) - V = view(Vp, 1:q, 1:q) - rmat_basep!(V, θ[lmm.covstr.tr[end]], lmm.data.zrv[i], lmm.covstr) - - θ₁ += logdet(V) - sweep!(Vp, 1:q) - iV = Symmetric(utriaply!(x -> -x, Vp[1:q, 1:q])) - mulαtβαinc!(θ2m, lmm.data.xv[i], iV) - θ₃ += -Vp[end, end] - end - θ₂ = logdet(θ2m) - return -(θ₁ + θ₂ + θ₃ + c) -end -=# diff --git a/src/rmat.jl b/src/rmat.jl index bd146787..b5183786 100644 --- a/src/rmat.jl +++ b/src/rmat.jl @@ -3,7 +3,7 @@ ################################################################################ ################################################################################ -function rmat_basep!(mx, θ::AbstractVector{T}, zrv, covstr::CovStructure{T2}) where T where T2 +function rmat_base_inc_b!(mx, θ::AbstractVector{T}, zrv, covstr::CovStructure{T2}) where T where T2 if covstr.repeated.covtype.s == :SI rmatp_si!(mx, θ, zrv, covstr.repeated.covtype) elseif covstr.repeated.covtype.s == :DIAG @@ -21,33 +21,11 @@ function rmat_basep!(mx, θ::AbstractVector{T}, zrv, covstr::CovStructure{T2}) w end end ################################################################################ -#= -function rmat_basep_z!(mx, θ::AbstractVector{T}, zrv, covstr) where T - for i = 1:length(covstr.block[end]) - if covstr.repeated.covtype.s == :SI - rmatp_si!(view(mx, covstr.block[end][i], covstr.block[end][i]), θ, view(covstr.rz, covstr.block[end][i], :), covstr.repeated.covtype) - elseif covstr.repeated.covtype.s == :DIAG - rmatp_diag!(view(mx, covstr.block[end][i], covstr.block[end][i]), θ, view(covstr.rz, covstr.block[end][i], :), covstr.repeated.covtype) - elseif covstr.repeated.covtype.s == :AR - rmatp_ar!(view(mx, covstr.block[end][i], covstr.block[end][i]), θ, view(covstr.rz, covstr.block[end][i], :), covstr.repeated.covtype) - elseif covstr.repeated.covtype.s == :ARH - rmatp_arh!(view(mx, covstr.block[end][i], covstr.block[end][i]), θ, view(covstr.rz, covstr.block[end][i], :), covstr.repeated.covtype) - elseif covstr.repeated.covtype.s == :CSH - rmatp_csh!(view(mx, covstr.block[end][i], covstr.block[end][i]), θ, view(covstr.rz, covstr.block[end][i], :), covstr.repeated.covtype) - elseif covstr.repeated.covtype.s == :CS - rmatp_cs!(view(mx, covstr.block[end][i], covstr.block[end][i]), θ, view(covstr.rz, covstr.block[end][i], :), covstr.repeated.covtype) - else - error("Unknown covariance structure!") - end - end - mx -end -=# ################################################################################ -function rmat_basep_z2!(mx, θ::AbstractVector{T}, covstr, block, sblock) where T +function rmat_base_inc!(mx, θ::AbstractVector{T}, covstr, block, sblock) where T zblock = view(covstr.rz, block, :) for i = 1:length(sblock[end]) - rmat_basep!(view(mx, sblock[end][i], sblock[end][i]), θ, view(zblock, sblock[end][i], :), covstr) + rmat_base_inc_b!(view(mx, sblock[end][i], sblock[end][i]), θ, view(zblock, sblock[end][i], :), covstr) end mx end diff --git a/src/showutils.jl b/src/showutils.jl index d1c14158..53a696fc 100644 --- a/src/showutils.jl +++ b/src/showutils.jl @@ -28,25 +28,28 @@ function printmatrix(io::IO, m::Matrix) end end -function rcoefnames(s, t, ::Val) - v = Vector{String}(undef, t) - v .= "─" -end -function rcoefnames(s, t, ::Val{:CSH}) - cn = coefnames(s) - if isa(cn, Vector) - l = length(cn) +function rcoefnames(s, t, ve) + + if ve == :SI + return ["Var"] + elseif ve == :DIAG + return string.(coefnames(s)) + elseif ve == :CS || ve == :AR + return ["Var", "Rho"] + elseif ve == :CSH || ve == :ARH + cn = coefnames(s) + if isa(cn, Vector) + l = length(cn) + else + l = 1 + end + v = Vector{String}(undef, t) + view(v, 1:l) .= string.(cn) + v[end] = "Rho" + return v else - l = 1 + v = Vector{String}(undef, t) + v .= "─" + return v end - v = Vector{String}(undef, t) - view(v, 1:l) .= string.(cn) - v[end] = "Rho" - v -end -function rcoefnames(s, t, ::Val{:SI}) - ["Var"] -end -function rcoefnames(s, t, ::Val{:DIAG}) - string.(coefnames(s)) end diff --git a/src/statsbase.jl b/src/statsbase.jl index f1055abc..68051c24 100644 --- a/src/statsbase.jl +++ b/src/statsbase.jl @@ -15,18 +15,65 @@ StatsBase.islinear(model::LMM) = error("islinear is not defined for $(typeof(mod StatsBase.nulldeviance(model::LMM) = error("nulldeviance is not defined for $(typeof(model)).") -StatsBase.loglikelihood(model::LMM) = - error("loglikelihood is not defined for $(typeof(model)).") - StatsBase.nullloglikelihood(model::LMM) = error("nullloglikelihood is not defined for $(typeof(model)).") StatsBase.score(model::LMM) = error("score is not defined for $(typeof(model)).") -StatsBase.nobs(model::LMM) = error("nobs is not defined for $(typeof(model)).") +=# + +#= +REML: n = total number of observation - number fixed effect parameters; d = number of covariance parameters +ML:, n = total number of observation; d = number of fixed effect parameters + number of covariance parameters. +=# +function StatsBase.nobs(lmm::LMM) + return length(lmm.data.yv) +end + +function StatsBase.dof_residual(lmm::LMM) + nobs(lmm) - lmm.rankx +end + +function StatsBase.dof(lmm::LMM) + lmm.nfixed + lmm.covstr.tl +end -StatsBase.dof(model::LMM) = error("dof is not defined for $(typeof(model)).") +function StatsBase.loglikelihood(lmm::LMM) + -lmm.result.reml/2 +end +function StatsBase.aic(lmm::LMM) + l = loglikelihood(lmm) + d = lmm.covstr.tl + -2l + 2d +end + +function StatsBase.bic(lmm::LMM) + l = loglikelihood(lmm) + d = lmm.covstr.tl + n = nobs(lmm) - lmm.nfixed + -2l + d * log(n) +end + +function StatsBase.aicc(lmm::LMM) + l = loglikelihood(lmm) + d = lmm.covstr.tl + n = nobs(lmm) - lmm.nfixed + -2l + (2d * n) / (n - d - 1.0) +end + +function caic(lmm::LMM) + l = loglikelihood(lmm) + d = lmm.covstr.tl + n = nobs(lmm) - lmm.nfixed + -2l + d * (log(n) + 1.0) +end + +function StatsBase.isfitted(lmm::LMM) + lmm.result.fit +end + +#= StatsBase.mss(model::LMM) = error("mss is not defined for $(typeof(model)).") StatsBase.rss(model::LMM) = error("rss is not defined for $(typeof(model)).") @@ -42,16 +89,6 @@ StatsBase.weights(model::LMM) = error("weights is not defined for $(typeof(model StatsBase.isfitted(model::LMM) = error("isfitted is not defined for $(typeof(model)).") -StatsBase.aic(model::LMM) = -2loglikelihood(model) + 2dof(model) - -function StatsBase.aicc(model::LMM) - k = dof(model) - n = nobs(model) - -2loglikelihood(model) + 2k + 2k*(k+1)/(n-k-1) -end - -StatsBase.bic(model::LMM) = -2loglikelihood(model) + dof(model)*log(nobs(model)) - function StatsBase.r2(model::LMM) Base.depwarn("The default r² method for linear models is deprecated. " * "Packages should define their own methods.", :r2) @@ -120,6 +157,5 @@ StatsBase.predict(model::LMM) = error("predict is not defined for $(typeof(model StatsBase.predict!(model::LMM) = error("predict! is not defined for $(typeof(model)).") -StatsBase.dof_residual(model::LMM) = error("dof_residual is not defined for $(typeof(model)).") =# diff --git a/src/utils.jl b/src/utils.jl index 6679d24d..8d3b0080 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,53 +1,5 @@ """ - Make X, Z matrices and vector y for each subject; -""" -#= -function subjblocks(df, sbj::Symbol, x::Matrix{T}, z::Matrix{T}, y::Vector{T}, rz::Matrix{T}) where T - u = unique(df[!, sbj]) - xa = Vector{Matrix{T}}(undef, length(u)) - za = Vector{Matrix{T}}(undef, length(u)) - ya = Vector{Vector{T}}(undef, length(u)) - rza = Vector{Matrix{T}}(undef, length(u)) - @inbounds @simd for i = 1:length(u) - v = findall(x -> x == u[i], df[!, sbj]) - xa[i] = Matrix(view(x, v, :)) - za[i] = Matrix(view(z, v, :)) - ya[i] = Vector(view(y, v)) - rza[i] = Matrix(view(rz, v, :)) - end - return xa, za, rza, ya -end -function subjblocks(df, sbj::Symbol, x::Matrix{T}, z::Matrix{T}, y::Vector{T}, rz::Nothing) where T - u = unique(df[!, sbj]) - xa = Vector{Matrix{T}}(undef, length(u)) - za = Vector{Matrix{T}}(undef, length(u)) - ya = Vector{Vector{T}}(undef, length(u)) - rza = Vector{Matrix{T}}(undef, length(u)) - @inbounds @simd for i = 1:length(u) - v = findall(x -> x == u[i], df[!, sbj]) - xa[i] = Matrix(view(x, v, :)) - za[i] = Matrix(view(z, v, :)) - ya[i] = Vector(view(y, v)) - rza[i] = Matrix{T}(undef, 0, 0) - end - return xa, za, rza, ya -end -function subjblocks(df, sbj) - if isa(sbj, Symbol) - u = unique(df[!, sbj]) - r = Vector{Vector{Int}}(undef, length(u)) - @inbounds @simd for i = 1:length(u) - r[i] = findall(x -> x == u[i], df[!, sbj]) - end - return r - end - r = Vector{Vector{Int}}(undef, 1) - r[1] = collect(1:size(df, 1)) - r -end -=# -""" - Intersect dataframe. +Intersect dataframe. """ function intersectdf(df, s)::Vector if isa(s, Nothing) return [collect(1:size(df, 1))] end @@ -75,23 +27,9 @@ function intersectdf(df, s)::Vector end res end -#= -function intersectsubj(covstr) - a = Vector{Vector{Symbol}}(undef, length(covstr.random)+1) - eq = true - for i = 1:length(covstr.random) - a[i] = covstr.random[i].subj - end - a[end] = covstr.repeated.subj - for i = 2:length(a) - if !(issetequal(a[1], a[i])) - eq = false - break - end - end - intersect(a...), eq -end -=# +""" +Intersect subject set in effects. +""" function intersectsubj(random, repeated) a = Vector{Vector{Symbol}}(undef, length(random)+1) eq = true @@ -107,12 +45,6 @@ function intersectsubj(random, repeated) end intersect(a...), eq end -#= -function diffsubj!(a, subj) - push!(a, subj) - symdiff(a...) -end -=# """ Variance estimate via OLS and QR decomposition. """ @@ -122,7 +54,17 @@ function initvar(y::Vector, X::Matrix{T}) where T r = y .- X * β sum(x -> x * x, r)/(length(r) - size(X, 2)), β end - +################################################################################ +function nterms(rhs) + if isa(rhs, Term) + p = 1 + elseif isa(rhs, Tuple) + p = length(rhs) + else + p = 0 + end + p +end ################################################################################ # VAR LINK ################################################################################ @@ -157,7 +99,7 @@ function rholinksigmoid2r(ρ::T) where T <: Real end ################################################################################ - +#= function varlinkvec(v) fv = Vector{Function}(undef, length(v)) for i = 1:length(v) @@ -172,24 +114,9 @@ function varlinkrvec(v) end fv end -#= -function varlinkvecapply(v, f) - rv = similar(v) - for i = 1:length(v) - rv[i] = f[i](v[i]) - end - rv -end =# -function varlinkvecapply!(v, f) - for i = 1:length(v) - v[i] = f[i](v[i]) - end - v -end -varlinkvecapply(v, f) = map.(f,v) #Test ################################################################################ -function varlinkvecapply2!(v, p; varlinkf = :exp, rholinkf = :sigm) +function varlinkvecapply!(v, p; varlinkf = :exp, rholinkf = :sigm) for i = 1:length(v) if p[i] == :var v[i] = vlink(v[i]) @@ -199,7 +126,7 @@ function varlinkvecapply2!(v, p; varlinkf = :exp, rholinkf = :sigm) end v end -function varlinkrvecapply2!(v, p; varlinkf = :exp, rholinkf = :sigm) +function varlinkrvecapply!(v, p; varlinkf = :exp, rholinkf = :sigm) for i = 1:length(v) if p[i] == :var v[i] = vlinkr(v[i]) @@ -209,7 +136,7 @@ function varlinkrvecapply2!(v, p; varlinkf = :exp, rholinkf = :sigm) end v end -function varlinkvecapply2(v, p; varlinkf = :exp, rholinkf = :sigm) +function varlinkvecapply(v, p; varlinkf = :exp, rholinkf = :sigm) s = similar(v) for i = 1:length(v) if p[i] == :var @@ -220,38 +147,7 @@ function varlinkvecapply2(v, p; varlinkf = :exp, rholinkf = :sigm) end s end -function varlinkrvecapply2(v, p; varlinkf = :exp, rholinkf = :sigm) - s = similar(v) - for i = 1:length(v) - if p[i] == :var - s[i] = vlinkr(v[i]) - else - s[i] = rholinksigmoidr(v[i]) - end - end - s -end ################################################################################ -#= -function vmatr(lmm, i) - θ = lmm.result.theta - G = gmat_base(θ, lmm.covstr) - V = mulαβαt(view(lmm.data.zv, lmm.data.block[i],:), G) - if length(lmm.data.zrv) > 0 - rmat_basep!(V, θ[lmm.covstr.tr[end]], view(lmm.data.zrv, lmm.data.block[i],:), lmm.covstr) - else - rmat_basep!(V, θ[lmm.covstr.tr[end]], lmm.data.zrv, lmm.covstr) - end - V -end - -function gmatr(lmm, i) - θ = lmm.result.theta - gmat_base(θ, lmm.covstr) -end -=# -################################################################################ - function m2logreml(lmm) lmm.result.reml end @@ -263,3 +159,32 @@ end function optim_callback(os) false end +################################################################################ +""" + gmatrix(lmm::LMM{T}, r::Int) where T +""" +function gmatrix(lmm::LMM{T}, r::Int) where T + if !lmm.result.fit error("Model not fitted!") end + if r > length(lmm.covstr.random) error("Invalid random effect number: $(r)!") end + G = zeros(T, lmm.covstr.q[r], lmm.covstr.q[r]) + gmat_switch!(G, lmm.result.theta, lmm.covstr, r) +end +""" + rmatrix(lmm::LMM{T}, i::Int) where T +""" +function rmatrix(lmm::LMM{T}, i::Int) where T + if !lmm.result.fit error("Model not fitted!") end + if i > length(lmm.data.block) error("Invalid block number: $(i)!") end + q = length(lmm.data.block[i]) + R = zeros(T, q, q) + rmat_base_inc!(R, lmm.result.theta[lmm.covstr.tr[end]], lmm.covstr, lmm.data.block[i], lmm.covstr.sblock[i]) + +end +""" + Return variance-covariance matrix V +""" +#= +function vmatrix() + +end +=# diff --git a/src/varstruct.jl b/src/varstruct.jl index a94a0d74..1e7836b8 100644 --- a/src/varstruct.jl +++ b/src/varstruct.jl @@ -10,36 +10,17 @@ macro covstr(ex) return :(@formula(nothing ~ $ex).rhs) end ################################################################################ -# SIMPLE FUNCTIONS -################################################################################ -function ffx(x::T)::T where T - x -end -function ffxzero(x::T)::T where T - zero(T) -end -function ffxone(x::T)::T where T - one(T) -end -function ffxpone(x::T)::T where T - x + one(T) -end -#= -function ffxmone(x::T)::T where T - x - one(T) -end -function ff2xmone(x::T)::T where T - 2x - one(T) -end -=# -################################################################################ # COVARIANCE TYPE ################################################################################ -struct CovarianceType <: AbstractCovarianceType +struct CovarianceType{T} <: AbstractCovarianceType s::Symbol #Covtype name - f::Function #number of parameters for Z size 2 - v::Function #number of variance parameters for Z size 2 - rho::Function #number of rho parameters for Z size 2 + t::T + function CovarianceType(s, t) + new{typeof(t)}(s, t) + end + function CovarianceType(s) + CovarianceType(s, 0) + end end ################################################################################ """ @@ -50,49 +31,49 @@ Scaled identity covariance type. SI = ScaledIdentity() """ function ScaledIdentity() - CovarianceType(:SI, ffxone, ffxone, ffxzero) + CovarianceType(:SI) end const SI = ScaledIdentity() """ Diag() """ function Diag() - CovarianceType(:DIAG, ffx, ffx, ffxzero) + CovarianceType(:DIAG) end const DIAG = Diag() """ Autoregressive() """ function Autoregressive() - CovarianceType(:AR, x -> 2, ffxone, ffxone) + CovarianceType(:AR) end const AR = Autoregressive() """ HeterogeneousAutoregressive() """ function HeterogeneousAutoregressive() - CovarianceType(:ARH, ffxpone, ffx, ffxone) + CovarianceType(:ARH) end const ARH = HeterogeneousAutoregressive() """ CompoundSymmetry() """ function CompoundSymmetry() - CovarianceType(:CS, x -> 2, ffxone, ffxone) + CovarianceType(:CS) end const CS = CompoundSymmetry() """ HeterogeneousCompoundSymmetry() """ function HeterogeneousCompoundSymmetry() - CovarianceType(:CSH, ffxpone, ffx, ffxone) + CovarianceType(:CSH) end const CSH = HeterogeneousCompoundSymmetry() """ RZero() """ function RZero() - CovarianceType(:ZERO, x -> 0, x -> 0, x -> 0) + CovarianceType(:ZERO) end #TOE @@ -103,6 +84,42 @@ end #ARMA +function covstrparam(ct::CovarianceType, q::Int, p::Int)::Tuple{Int, Int, Int} + if ct.s == :SI + return (1, 0, 1) + elseif ct.s == :DIAG + return (q, 0, q) + elseif ct.s == :VC + return (p, 0, p) + elseif ct.s == :AR + return (1, 1, 2) + elseif ct.s == :ARH + return (q, 1, q + 1) + elseif ct.s == :ARMA + return (1, 2, 3) + elseif ct.s == :CS + return (1, 1, 2) + elseif ct.s == :CSH + return (q, 1, q + 1) + elseif ct.s == :TOEP + return (1, q - 1, q) + elseif ct.s == :TOEPH + return (q, q - 1, 2 * q - 1) + elseif ct.s == :TOEPB + return (1, ct.t - 1, ct.t) + elseif ct.s == :TOEPHB + return (q, ct.t - 1, q + ct.t - 1) + elseif ct.s == :UN + return (q, q * (q + 1) / 2 - q, q * (q + 1) / 2) + elseif ct.s == :ZERO + return (0, 0, 0) + elseif ct.s == :FUNC + error("Not implemented!") + else + error("Unknown covariance type!") + end +end + ################################################################################ # EFFECT ################################################################################ @@ -116,8 +133,8 @@ struct VarEffect covtype::CovarianceType coding::Dict{Symbol, AbstractContrasts} subj::Vector{Symbol} + p::Int function VarEffect(model, covtype::T, coding; subj = nothing) where T <: AbstractCovarianceType - #if isa(model, AbstractTerm) model = tuple(model) end if isa(subj, Nothing) subj = Vector{Symbol}(undef, 0) elseif isa(subj, Symbol) @@ -127,14 +144,22 @@ struct VarEffect else throw(ArgumentError("subj type should be Symbol or Vector{tymbol}")) end - + p = nterms(model) + #= + if isa(model, Term) + p = 1 + elseif isa(model, Tuple) + p = length(model) + else + p = 0 + end + =# if coding === nothing && model !== nothing coding = Dict{Symbol, AbstractContrasts}() elseif coding === nothing && model === nothing coding = Dict{Symbol, AbstractContrasts}() end - #if isa(model, AbstractTerm) model = tuple(model) end - new(model, covtype, coding, subj) + new(model, covtype, coding, subj, p) end function VarEffect(model, covtype::T; coding = nothing, subj = nothing) where T <: AbstractCovarianceType VarEffect(model, covtype, coding; subj = subj) @@ -172,7 +197,7 @@ struct CovStructure{T} <: AbstractCovarianceStructure rz::Matrix{T} # size 2 of z/rz matrix q::Vector{Int} - # number of parametert in each effect + # total number of parameters in each effect t::Vector{Int} # range of each parameters in θ vector tr::Vector{UnitRange{UInt32}} @@ -194,6 +219,10 @@ struct CovStructure{T} <: AbstractCovarianceStructure sblock = Vector{Vector{Vector{Vector{UInt32}}}}(undef, length(blocks)) zrndur = Vector{UnitRange{Int}}(undef, alleffl - 1) rz = Matrix{Float64}(undef, size(data, 1), 0) + #Theta parameter type + ct = Vector{Symbol}(undef, 0) + # Names + rcnames = Vector{String}(undef, 0) # # RANDOM EFFECTS for i = 1:length(random) @@ -208,11 +237,13 @@ struct CovStructure{T} <: AbstractCovarianceStructure schema[i] = apply_schema(random[i].model, StatsModels.schema(data, random[i].coding)) ztemp = modelcols(MatrixTerm(schema[i]), data) q[i] = size(ztemp, 2) - t[i] = random[i].covtype.f(q[i]) + csp = covstrparam(random[i].covtype, q[i], random[i].p) + t[i] = csp[3] z = hcat(z, ztemp) fillur!(zrndur, i, q) fillur!(tr, i, t) subjmatrix!(random[i].subj, data, subjz, i) + updatenametype!(ct, rcnames, csp, schema[i], random[i].covtype.s) end # REPEATED EFFECTS if length(repeated.coding) == 0 @@ -223,51 +254,23 @@ struct CovStructure{T} <: AbstractCovarianceStructure rz = modelcols(MatrixTerm(schema[end]), data) subjmatrix!(repeated.subj, data, subjz, length(subjz)) q[end] = size(rz, 2) - t[end] = repeated.covtype.f(q[end]) + csp = covstrparam(repeated.covtype, q[end], repeated.p) + t[end] = csp[3] tr[end] = UnitRange(sum(t[1:end-1]) + 1, sum(t[1:end-1]) + t[end]) + updatenametype!(ct, rcnames, csp, schema[end], repeated.covtype.s) #Theta length tl = sum(t) - #Theta parameter type - ct = Vector{Symbol}(undef, tl) - # Names - rcnames = Vector{String}(undef, tl) - ctn = 1 - for i = 1:length(random) - if random[i].covtype.s == :ZERO - continue - end - for i2 = 1:random[i].covtype.v(q[i]) - ct[ctn] = :var - ctn +=1 - end - if random[i].covtype.rho(q[i]) > 0 - for i2 = 1:random[i].covtype.rho(q[i]) - ct[ctn] = :rho - ctn +=1 + ######################################################################## + ######################################################################## + for i = 1:length(blocks) + sblock[i] = Vector{Vector{Vector{UInt32}}}(undef, alleffl) + for s = 1:alleffl + sblock[i][s] = Vector{Vector{UInt32}}(undef, 0) + for col in eachcol(view(subjz[s], blocks[i], :)) + if any(col) push!(sblock[i][s], sort!(findall(x->x==true, col))) end end end - view(rcnames, tr[i]) .= rcoefnames(schema[i], t[i], Val{random[i].covtype.s}()) - end - for i2 = 1:repeated.covtype.v(q[end]) - ct[ctn] = :var - ctn +=1 - end - if repeated.covtype.rho(q[end]) > 0 - for i2 = 1:repeated.covtype.rho(q[end]) - ct[ctn] = :rho - ctn +=1 - end end - view(rcnames, tr[end]) .= rcoefnames(schema[end], t[end], Val{repeated.covtype.s}()) - for i = 1:length(blocks) - sblock[i] = Vector{Vector{Vector{UInt32}}}(undef, alleffl) - for s = 1:alleffl - sblock[i][s] = Vector{Vector{UInt32}}(undef, 0) - for col in eachcol(view(subjz[s], blocks[i], :)) - if any(col) push!(sblock[i][s], sort!(findall(x->x==true, col))) end - end - end - end # new{eltype(z)}(random, repeated, schema, rcnames, block, z, sblock, zrndur, rz, q, t, tr, tl, ct) end @@ -285,17 +288,12 @@ function fillur!(ur, i, v) end end ################################################################################ -#= -function schemalength(s) - if isa(s, Tuple) - return length(s) - else - return 1 - end +function updatenametype!(ct, rcnames, csp, schema, s) + append!(ct, fill!(Vector{Symbol}(undef, csp[1]), :var)) + append!(ct, fill!(Vector{Symbol}(undef, csp[2]), :rho)) + append!(rcnames, rcoefnames(schema, csp[3], s)) end -=# -# - +################################################################################ function subjmatrix!(subj, data, subjz, i) if length(subj) > 0 if length(subj) == 1 @@ -347,14 +345,7 @@ function fill_coding_dict!(t::T, d::Dict, data) where T <: Tuple end end ################################################################################ -""" - Return variance-covariance matrix V -""" -#= -function vmat()::AbstractMatrix -end -=# ################################################################################ function Base.show(io::IO, e::VarEffect) println(io, "Effect") diff --git a/test/bench.jl b/test/bench.jl deleted file mode 100644 index c3d2e710..00000000 --- a/test/bench.jl +++ /dev/null @@ -1,124 +0,0 @@ -using DataFrames, Plots, ForwardDiff, StatsBase, StatsModels, TimerOutputs - -function gendata(subjn, sobsn, f, f2) - slist = collect(1:subjn) - sdat = Vector{Int}(undef, 0) - vdat = Vector{eltype(f)}(undef, 0) - for i = 1:length(slist) - v = Vector{Int}(undef, sobsn) - v .= slist[i] - append!(sdat, v) - vf = Vector{eltype(f)}(undef, sobsn) - vf .= rand(f) - append!(vdat, vf) - end - rf = rand(f2, subjn*sobsn) - rv = rand(subjn*sobsn) - DataFrame(subject = sdat, fac = vdat, rfac = rf, var = rv) -end - -function bench1(s, e) - ne = e - ns = s - num = ne - ns + 1 - mem = Vector{Int}(undef, num) - mem2 = Vector{Int}(undef, num) - for i = ns:ne - df = gendata(10, i, ["a", "b"], ["c", "d"]) - categorical!(df, :subject); - categorical!(df, :fac); - categorical!(df, :rfac); - lmm = Metida.LMM(@formula(var~fac), df; - random = Metida.VarEffect(Metida.@covstr(rfac), Metida.VC; subj = :subject), - repeated = Metida.VarEffect(Metida.@covstr(rfac), Metida.VC; subj = :subject), - subject = :subject) - theta = rand(lmm.covstr.tl) - mem[i-ns+1] = @allocated Metida.reml_sweep_β2(lmm, theta) - mem2[i-ns+1] = @allocated ForwardDiff.hessian(x -> Metida.reml_sweep_β2(lmm, x)[1], theta) - end - mem, mem2 -end - -function bench2(s, e) - ne = e - ns = s - num = ne - ns + 1 - mem = Vector{Int}(undef, num) - mem2 = Vector{Int}(undef, num) - for i = ns:ne - df = gendata(10, i, ["a", "b"], ["c", "d"]) - categorical!(df, :subject); - categorical!(df, :fac); - categorical!(df, :rfac); - lmm = Metida.LMM(@formula(var~fac), df; - random = Metida.VarEffect(Metida.@covstr(rfac), Metida.VC; subj = :subject), - repeated = Metida.VarEffect(Metida.@covstr(rfac), Metida.VC; subj = :subject), - subject = :subject) - theta = rand(lmm.covstr.tl) - mem[i-ns+1] = @allocated Metida.reml_sweep_β(lmm, theta) - mem2[i-ns+1] = @allocated ForwardDiff.hessian(x -> Metida.reml_sweep_β(lmm, x)[1], theta) - end - mem, mem2 -end - - -n = 120 -y01, y02 = bench1(3, n) -y11, y12 = bench2(3, n) -x = collect(3:n) - -plot(x, y01) -plot!(x, y11) -plot(x, y02) -p = plot!(x, y12) -png("plot120-VC-VC.png", p) -maxmem = y02[end]/2^20 - -function bench3(s, e) - ne = e - ns = s - num = ne - ns + 1 - mem = Vector{Int}(undef, num) - mem2 = Vector{Int}(undef, num) - for i = ns:ne - df = gendata(10, i, ["a", "b"], ["c", "d"]) - categorical!(df, :subject); - categorical!(df, :fac); - categorical!(df, :rfac); - lmm = Metida.LMM(@formula(var~fac), df; - random = Metida.VarEffect(Metida.@covstr(rfac), Metida.VC; subj = :subject), - repeated = Metida.VarEffect(Metida.@covstr(rfac), Metida.VC; subj = :subject), - subject = :subject) - theta = rand(lmm.covstr.tl) - t1 = @elapsed mem[i-ns+1] = @allocated Metida.fit!(lmm) - t2 = @elapsed mem2[i-ns+1] = @allocated Metida.fit2!(lmm) - println("T1: ", t1, " | T2: ", t2 ) - end - mem, mem2 -end - -n = 100 -y01, y02 = bench3(4, n) -x = collect(4:n) -plot(x, y01) -p = plot!(x, y02) - -png("plot100-VC-VC-NLopt-Optim.png") - -#[1.49136, 2.07002, 0.16827] - - Metida.reml_sweep_β(lmm, lmm.result.theta) - Metida.reml_sweep_β(lmm, [1.49136, 2.07002, 0.16827]) - -v = [0.018569469816145896; 0.015444785130183025; 0.015444785130183028; 0.015444785130182997; 0.015444785130183044; 0.015444785130183032; 0.015444785130183039; 0.015444785130182985; 0.015444785130183034; 0.015444785130183014; 0.01544478513018298; 0.015444785130183032; 0.015444785130183023; 0.01544478513018301; 0.015444785130183039; 0.015444785130183002; 0.015444785130183034; 0.015444785130183032; 0.015444785130183032; 0.01544478513018304; 0.015444785130183039; 0.015444785130183021; 0.015444785130183053; 0.01544478513018301; 0.01544478513018301; 0.018569469816145892] -y = [78.0, 77.0, 78.0, 77.0, 76.0, 74.0, 76.0, 74.0, 76.0, 77.0, 77.0, 77.0, 75.0, 73.0, 72.0, 76.0, 77.0, 76.0, 78.0, 72.0, 77.0, 71.0, 72.0, 77.0, 77.0, 78.0] -m = reshape(v, length(v), 1) -Metida.mulαtβinc!(zeros(1,1), m, y) - - -theta = lmm.result.theta .- 0.1 -ai = ForwardDiff.hessian(x -> Metida.reml_sweep_β2(lmm, x)[4], theta) -g = ForwardDiff.gradient(x -> Metida.reml_sweep_β2(lmm, x)[1], theta) -he = ForwardDiff.hessian(x -> Metida.reml_sweep_β2(lmm, x)[1], theta) - -theta - he*g diff --git a/test/berds.jl b/test/berds.jl index ac4d33dd..299a7c86 100644 --- a/test/berds.jl +++ b/test/berds.jl @@ -25,37 +25,33 @@ REML G - SI, R - SI 30 - 32.418048 =# -df = CSV.File(path*"/csv/berds/rds29.csv", types = Dict(:PK => Float64)) |> DataFrame for i = 1:30 - df = CSV.File(path*"/csv/berds/rds"*string(i)*".csv", types = Dict(:PK => Float64)) |> DataFrame - dropmissing!(df) - transform!(df, :subject => categorical, renamecols=false) - transform!(df, :period => categorical, renamecols=false) - transform!(df, :sequence => categorical, renamecols=false) - transform!(df, :treatment => categorical, renamecols=false) + dfrds = CSV.File(path*"/csv/berds/rds"*string(i)*".csv", types = Dict(:PK => Float64)) |> DataFrame + dropmissing!(dfrds) + transform!(dfrds, :subject => categorical, renamecols=false) + transform!(dfrds, :period => categorical, renamecols=false) + transform!(dfrds, :sequence => categorical, renamecols=false) + transform!(dfrds, :treatment => categorical, renamecols=false) - df.lnpk = log.(df.PK) + dfrds.lnpk = log.(dfrds.PK) @testset " RDS Test $(i) " begin atol=1E-6 - lmm = Metida.LMM(@formula(lnpk~sequence+period+treatment), df; + lmm = Metida.LMM(@formula(lnpk~sequence+period+treatment), dfrds; random = Metida.VarEffect(Metida.@covstr(1), Metida.SI), subject = :subject ) Metida.fit!(lmm) @test lmm.result.reml ≈ remlsb[i] atol=atol - if i == 13 || i == 15 atol = 1E-4 end - lmm = Metida.LMM(@formula(lnpk~sequence+period+treatment), df; + lmm = Metida.LMM(@formula(lnpk~sequence+period+treatment), dfrds; random = Metida.VarEffect(Metida.@covstr(treatment), Metida.CSH), repeated = Metida.VarEffect(Metida.@covstr(treatment), Metida.DIAG), subject = :subject ) Metida.fit!(lmm) @test lmm.result.reml ≈ remls[i] atol=atol - - end end diff --git a/test/lme4.jl b/test/lme4.jl index dc554fbb..3c9c67b7 100644 --- a/test/lme4.jl +++ b/test/lme4.jl @@ -1,14 +1,19 @@ - +################################################################################ +# sleepstudy.csv +################################################################################ df = CSV.File(path*"/csv/lme4/sleepstudy.csv") |> DataFrame transform!(df, :Subject => categorical, renamecols=false) transform!(df, :Days => categorical, renamecols=false) -#= -SPSS -REML 1729.492560 -=# @testset " sleepstudy.csv " begin + #= + SPSS + REML 1729.492560 + R 1375.465318 + Res 987.588488 + =# + #Model 1 lmm = Metida.LMM(@formula(Reaction~Days), df; random = Metida.VarEffect(Metida.SI), subject = :Subject @@ -18,6 +23,10 @@ REML 1729.492560 end @testset " CS sleepstudy.csv " begin + #REML 1903.327 + # 1662.172084 + # 296.693108 + #Model 2 lmm = Metida.LMM(@formula(Reaction~1), df; random = Metida.VarEffect(Metida.@covstr(Days), Metida.CS), subject = :Subject @@ -30,14 +39,41 @@ end ) Metida.fit!(lmm) @test lmm.result.reml ≈ 1729.4925602367027 atol=1E-6 + +end + +@testset " CSH sleepstudy.csv " begin + #REML 1772.095 + #Model 3 + lmm = Metida.LMM(@formula(Reaction~1), df; + random = Metida.VarEffect(Metida.@covstr(Days), Metida.CSH, subj = :Subject) + ) + Metida.fit!(lmm) + @test lmm.result.reml ≈ 1772.0953251997046 atol=1E-6 end +@testset " ARH sleepstudy.csv " begin + #REML 1730.189543 + #Model 4 + lmm = Metida.LMM(@formula(Reaction~1), df; + random = Metida.VarEffect(Metida.@covstr(Days), Metida.ARH, subj = :Subject) + ) + Metida.fit!(lmm) + @test lmm.result.reml ≈ 1730.1895427398322 atol=1E-6 +end + + +################################################################################ +# Penicillin.csv +################################################################################ + df = CSV.File(path*"/csv/lme4/Penicillin.csv") |> DataFrame df.diameter = float.(df.diameter) transform!(df, :plate => categorical, renamecols=false) transform!(df, :sample => categorical, renamecols=false) -@testset " Penicillin.csv " begin +@testset " SI + SI Penicillin.csv " begin + #SPSS 330.860589 lmm = Metida.LMM(@formula(diameter~1), df; random = [Metida.VarEffect(Metida.SI, subj = :plate), Metida.VarEffect(Metida.SI, subj = :sample)] ) @@ -45,13 +81,17 @@ transform!(df, :sample => categorical, renamecols=false) @test lmm.result.reml ≈ 330.86058899109184 atol=1E-6 end +################################################################################ +# Pastes.csv +################################################################################ df = CSV.File(path*"/csv/lme4/Pastes.csv") |> DataFrame transform!(df, :batch => categorical, renamecols=false) transform!(df, :sample => categorical, renamecols=false) transform!(df, :cask=> categorical, renamecols=false) -@testset " Pastes.csv " begin +@testset " SI + SI Pastes.csv " begin + #REML 246.990746 lmm = Metida.LMM(@formula(strength~1), df; random = [Metida.VarEffect(Metida.SI, subj = :batch), Metida.VarEffect(Metida.SI, subj = [:batch, :cask])] ) diff --git a/test/nlopt.jl b/test/nlopt.jl deleted file mode 100644 index c754d92c..00000000 --- a/test/nlopt.jl +++ /dev/null @@ -1,35 +0,0 @@ -using NLopt - -lmm = Metida.LMM(@formula(var~sequence+period+formulation), df; -random = [Metida.VarEffect(Metida.@covstr(formulation), Metida.CSH)], -repeated = Metida.VarEffect(Metida.@covstr(formulation), Metida.VC), -subject = :subject) - -lmmr = Metida.fit!(lmm) - -################################################################################ -opt = NLopt.Opt(:LN_BOBYQA, 5) -NLopt.ftol_rel!(opt, 1.0e-10) -NLopt.ftol_abs!(opt, 1.0e-10) -NLopt.xtol_rel!(opt, 1.0e-10) -NLopt.xtol_abs!(opt, 1.0e-10) -NLopt.initial_step!(opt, [0.005, 0.005, 0.005, 0.005, 0.005]) -fv = Metida.varlinkvec(lmm.covstr.ct) -init = deepcopy( lmm.result.optim.initial_x) -obj = (x,y) -> Metida.reml_sweep_β2(lmm, Metida.varlinkvecapply!(x, fv))[1] -NLopt.min_objective!(opt, obj) -#res = NLopt.optimize!(opt, init) -function fx(opt, i) - init = deepcopy(i) - NLopt.optimize!(opt, init) -end -@benchmark fx($opt, $init) -################################################################################ -#:LN_BOBYQA -#:LN_NEWUOA -#:LN_NEWUOA_BOUND -#:LN_PRAXIS -#:LN_NELDERMEAD -#:LN_SBPLX -#:LN_COBYLA -################################################################################ diff --git a/test/optim.jl b/test/optim.jl deleted file mode 100644 index 43787222..00000000 --- a/test/optim.jl +++ /dev/null @@ -1,10 +0,0 @@ - - -optmethod = Optim.Newton() -optoptions = Optim.Options(g_tol = 1e-12, - iterations = 300, - store_trace = true, - show_trace = true, - allow_f_increases = true) -objopt = x -> Metida.reml_sweep_β(lmm, Metida.varlinkvecapply!(x, fv))[1] -O = Optim.optimize(objopt, deepcopy(lmm.result.optim.initial_x), Optim.Newton(); autodiff = :forward) diff --git a/test/plot100-VC-VC-NLopt-Optim.png b/test/plot100-VC-VC-NLopt-Optim.png deleted file mode 100644 index 6359b47e..00000000 Binary files a/test/plot100-VC-VC-NLopt-Optim.png and /dev/null differ diff --git a/test/plot120-VC-VC.png b/test/plot120-VC-VC.png deleted file mode 100644 index 914cf491..00000000 Binary files a/test/plot120-VC-VC.png and /dev/null differ diff --git a/test/test.jl b/test/test.jl index e1896031..6fe4cc97 100644 --- a/test/test.jl +++ b/test/test.jl @@ -1,6 +1,6 @@ # Metida -using Test, CSV, DataFrames, StatsModels +using Test, CSV, DataFrames, StatsModels, StatsBase path = dirname(@__FILE__) include("testdata.jl") @@ -47,7 +47,31 @@ include("testdata.jl") Metida.fit!(lmm; aifirst = true, init = Metida.theta(lmm)) @test Metida.m2logreml(lmm) ≈ 16.241112644506067 atol=1E-6 + + Metida.gmatrix(lmm, 1) + Metida.rmatrix(lmm, 1) + dof(lmm) + + @test nobs(lmm) == 20 + @test bic(lmm) ≈ 24.558878811225412 atol=1E-6 + @test aic(lmm) ≈ 22.241112644506067 atol=1E-6 + @test aicc(lmm) ≈ 24.241112644506067 atol=1E-6 + @test Metida.caic(lmm) ≈ 27.558878811225412 atol=1E-6 + @test dof_residual(lmm) == 14 + @test isfitted(lmm) == true + end +@testset " Errors " begin + @test_throws ArgumentError Metida.VarEffect(Metida.@covstr(formulation), Metida.DIAG, subj = "subj") + @test_throws ErrorException lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0; + random = Metida.VarEffect(Metida.@covstr(formulation), Metida.CovarianceType(:XX)), + ) + lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0; + random = Metida.VarEffect(Metida.@covstr(formulation), Metida.DIAG), + ) + @test_throws ErrorException Metida.fit!(lmm; init = [0.0, 1.0, 0.0, 0.0, 0.0]) +end + @testset " Model: DIAG/subject + nothing " begin lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0; random = Metida.VarEffect(Metida.@covstr(formulation), Metida.DIAG), diff --git a/test/warn.jl b/test/warn.jl deleted file mode 100644 index 41758cdc..00000000 --- a/test/warn.jl +++ /dev/null @@ -1,54 +0,0 @@ -#Code warntype - -@code_warntype Metida.reml_sweep_β3(lmm, lmm.result.theta) - -@code_typed Metida.reml_sweep_β3(lmm, lmm.result.theta) - -@code_warntype Metida.gmat_base(lmm.result.theta, lmm.covstr) - -@code_typed Metida.gmat_base(lmm.result.theta, lmm.covstr) - -V = [1. 2 3 4; -1 2 3 3; -9 8 2 1; -1 2 1 2] - -@code_warntype Metida.rmat_basep!(V, lmm.result.theta[lmm.covstr.tr[end]], view(lmm.covstr.rz, lmm.data.block[1],:), lmm.covstr) - -@code_warntype Metida.gmat_base_z2!(V, lmm.result.theta, lmm.covstr, lmm.data.block[1], lmm.covstr.sblock[1]) - -@code_warntype Metida.rmat_basep_z2!(V, lmm.result.theta[lmm.covstr.tr[end]], lmm.covstr, lmm.data.block[1], lmm.covstr.sblock[1]) - -@code_warntype Metida.rmatp_si!(V, lmm.result.theta[lmm.covstr.tr[end]], view(lmm.covstr.rz, lmm.covstr.block[end][1], :), lmm.covstr.repeated.covtype) - -@code_typed Metida.rmatp_si!(V, lmm.result.theta[lmm.covstr.tr[end]], view(lmm.covstr.rz, lmm.covstr.block[end][1], :), lmm.covstr.repeated.covtype) - -@code_warntype Metida.rmatp_diag!(V, lmm.result.theta[lmm.covstr.tr[end]], view(lmm.covstr.rz, lmm.covstr.block[end][1], :), lmm.covstr.repeated.covtype) - -@code_warntype Metida.rmatp_csh!(V, lmm.result.theta[lmm.covstr.tr[end]], view(lmm.covstr.rz, lmm.covstr.block[end][1], :), lmm.covstr.repeated.covtype) - -@code_warntype Metida.rmatp_cs!(V, lmm.result.theta[lmm.covstr.tr[end]], view(lmm.covstr.rz, lmm.covstr.block[end][1], :), lmm.covstr.repeated.covtype) - -@code_warntype Metida.rmatp_ar!(V, lmm.result.theta[lmm.covstr.tr[end]], view(lmm.covstr.rz, lmm.covstr.block[end][1], :), lmm.covstr.repeated.covtype) - -@code_warntype Metida.rmatp_arh!(V, lmm.result.theta[lmm.covstr.tr[end]], view(lmm.covstr.rz, lmm.covstr.block[end][1], :), lmm.covstr.repeated.covtype) - - -G = [1. 2.; -2. 3.] -@code_warntype Metida.gmat_si!(G, lmm.result.theta[lmm.covstr.tr[1]], lmm.covstr.q[1], lmm.covstr.random[1].covtype) - -@code_warntype Metida.gmat_diag!(G, lmm.result.theta[lmm.covstr.tr[1]], lmm.covstr.q[1], lmm.covstr.random[1].covtype) - -@code_warntype Metida.gmat_csh!(G, lmm.result.theta[lmm.covstr.tr[1]], lmm.covstr.q[1], lmm.covstr.random[1].covtype) - -@code_warntype Metida.gmat_cs!(G, lmm.result.theta[lmm.covstr.tr[1]], lmm.covstr.q[1], lmm.covstr.random[1].covtype) - -@code_warntype Metida.gmat_ar!(G, lmm.result.theta[lmm.covstr.tr[1]], lmm.covstr.q[1], lmm.covstr.random[1].covtype) - - - -V = zeros(4,4) - -Metida.rmatp_ar!(V, lmm.result.theta[lmm.covstr.tr[end]], view(lmm.covstr.rz, lmm.covstr.block[end][1], :), lmm.covstr.repeated.covtype) -Metida.rmatp_ar!(V, [1.49136, 2.07002, 0.16827][lmm.covstr.tr[end]], view(lmm.covstr.rz, lmm.covstr.block[end][1], :), lmm.covstr.repeated.covtype)