From 9c516fdccbd7e36055612c178f1c0d31d9e27d9e Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 21 Mar 2026 11:03:51 -0700 Subject: [PATCH 01/18] Project.toml: support Symbolics v7 & Utils v4 --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index b445f3db..eabc5b36 100644 --- a/Project.toml +++ b/Project.toml @@ -40,8 +40,8 @@ Optim = "1" PrettyTables = "3" ProximalAlgorithms = "0.7" StatsBase = "0.33, 0.34" -Symbolics = "4, 5, 6" -SymbolicUtils = "1.4 - 1.5, 1.7, 2, 3" +Symbolics = "4, 5, 6, 7" +SymbolicUtils = "1.4 - 1.5, 1.7, 2, 3, 4" StatsAPI = "1" [extras] From 6e1ffaaa76f21bfd44fe8ee506fa370f9d4e22f4 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 21 Mar 2026 11:03:51 -0700 Subject: [PATCH 02/18] prepare_start_params(): tighten type check --- src/optimizer/abstract.jl | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/optimizer/abstract.jl b/src/optimizer/abstract.jl index e9a8c47b..0c7913c4 100644 --- a/src/optimizer/abstract.jl +++ b/src/optimizer/abstract.jl @@ -154,11 +154,24 @@ function prepare_start_params(start_val::AbstractVector, model::AbstractSem; kwa "The length of `start_val` vector ($(length(start_val))) does not match the number of model parameters ($(nparams(model))).", ), ) + (eltype(start_val) <: Number) || throw( + TypeError( + :prepare_start_params, + "start_val elements must be numeric", + Number, + eltype(start_val), + ), + ) return start_val end function prepare_start_params(start_val::AbstractDict, model::AbstractSem; kwargs...) - return [start_val[param] for param in params(model)] # convert to a vector + # convert to a vector + return prepare_start_params( + [start_val[param] for param in params(model)], + model; + kwargs..., + ) end # get from the ParameterTable (potentially from a different model with match param names) From 32068dee990fd5474387aeaef008186f08e2e627 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 21 Mar 2026 11:03:51 -0700 Subject: [PATCH 03/18] SemImplied/SemLossFun: drop meanstructure kwarg - for SemImplied require spec::SemSpec as positional - for SemLossFunction require implied argument --- src/frontend/specification/Sem.jl | 10 +++++- src/implied/RAM/generic.jl | 24 +++++-------- src/implied/RAM/symbolic.jl | 32 ++++++----------- src/implied/abstract.jl | 14 -------- src/loss/ML/FIML.jl | 11 +++--- src/loss/ML/ML.jl | 3 +- src/loss/WLS/WLS.jl | 27 ++++++++------ .../recover_parameters_twofact.jl | 2 +- test/unit_tests/model.jl | 35 ++++--------------- 9 files changed, 59 insertions(+), 99 deletions(-) diff --git a/src/frontend/specification/Sem.jl b/src/frontend/specification/Sem.jl index 53858abd..a47bad4b 100644 --- a/src/frontend/specification/Sem.jl +++ b/src/frontend/specification/Sem.jl @@ -109,7 +109,15 @@ function get_fields!(kwargs, specification, observed, implied, loss) # implied if !isa(implied, SemImplied) - implied = implied(; specification, kwargs...) + # FIXME remove this implicit logic + # SemWLS only accepts vech-ed implied covariance + if isa(loss, Type) && (loss <: SemWLS) && !haskey(kwargs, :vech) + implied_kwargs = copy(kwargs) + implied_kwargs[:vech] = true + else + implied_kwargs = kwargs + end + implied = implied(specification; implied_kwargs...) end kwargs[:implied] = implied diff --git a/src/implied/RAM/generic.jl b/src/implied/RAM/generic.jl index 4c1fa323..d57500a3 100644 --- a/src/implied/RAM/generic.jl +++ b/src/implied/RAM/generic.jl @@ -6,14 +6,10 @@ Model implied covariance and means via RAM notation. # Constructor - RAM(;specification, - meanstructure = false, - gradient = true, - kwargs...) + RAM(; specification, gradient = true, kwargs...) # Arguments - `specification`: either a `RAMMatrices` or `ParameterTable` object -- `meanstructure::Bool`: does the model have a meanstructure? - `gradient::Bool`: is gradient-based optimization used # Extended help @@ -53,9 +49,9 @@ Vector of indices of each parameter in the respective RAM matrix: - `ram.M_indices` Additional interfaces -- `ram.F⨉I_A⁻¹` -> ``F(I-A)^{-1}`` -- `ram.F⨉I_A⁻¹S` -> ``F(I-A)^{-1}S`` -- `ram.I_A` -> ``I-A`` +- `F⨉I_A⁻¹(::RAM)` -> ``F(I-A)^{-1}`` +- `F⨉I_A⁻¹S(::RAM)` -> ``F(I-A)^{-1}S`` +- `I_A(::RAM)` -> ``I-A`` Only available in gradient! calls: - `ram.I_A⁻¹` -> ``(I-A)^{-1}`` @@ -90,15 +86,13 @@ end ### Constructors ############################################################################################ -function RAM(; - specification::SemSpecification, +function RAM( + spec::SemSpecification; + #vech = false, gradient_required = true, - meanstructure = false, kwargs..., ) - ram_matrices = convert(RAMMatrices, specification) - - check_meanstructure_specification(meanstructure, ram_matrices) + ram_matrices = convert(RAMMatrices, spec) # get dimensions of the model n_par = nparams(ram_matrices) @@ -126,7 +120,7 @@ function RAM(; end # μ - if meanstructure + if !isnothing(ram_matrices.M) MS = HasMeanStruct M_pre = materialize(ram_matrices.M, rand_params) ∇M = gradient_required ? sparse_gradient(ram_matrices.M) : nothing diff --git a/src/implied/RAM/symbolic.jl b/src/implied/RAM/symbolic.jl index df7c497a..0f7868da 100644 --- a/src/implied/RAM/symbolic.jl +++ b/src/implied/RAM/symbolic.jl @@ -12,12 +12,10 @@ Subtype of `SemImplied` that implements the RAM notation with symbolic precomput gradient = true, hessian = false, approximate_hessian = false, - meanstructure = false, kwargs...) # Arguments - `specification`: either a `RAMMatrices` or `ParameterTable` object -- `meanstructure::Bool`: does the model have a meanstructure? - `gradient::Bool`: is gradient-based optimization used - `hessian::Bool`: is hessian-based optimization used - `approximate_hessian::Bool`: for hessian based optimization: should the hessian be approximated @@ -79,20 +77,16 @@ end ### Constructors ############################################################################################ -function RAMSymbolic(; - specification::SemSpecification, - loss_types = nothing, - vech = false, - simplify_symbolics = false, - gradient = true, - hessian = false, - meanstructure = false, - approximate_hessian = false, +function RAMSymbolic( + spec::SemSpecification; + vech::Bool = false, + simplify_symbolics::Bool = false, + gradient::Bool = true, + hessian::Bool = false, + approximate_hessian::Bool = false, kwargs..., ) - ram_matrices = convert(RAMMatrices, specification) - - check_meanstructure_specification(meanstructure, ram_matrices) + ram_matrices = convert(RAMMatrices, spec) n_par = nparams(ram_matrices) par = (Symbolics.@variables θ[1:n_par])[1] @@ -102,10 +96,6 @@ function RAMSymbolic(; M = !isnothing(ram_matrices.M) ? materialize(Num, ram_matrices.M, par) : nothing F = ram_matrices.F - if !isnothing(loss_types) && any(T -> T <: SemWLS, loss_types) - vech = true - end - I_A⁻¹ = neumann_series(A) # Σ @@ -146,7 +136,7 @@ function RAMSymbolic(; end # μ - if meanstructure + if !isnothing(ram_matrices.M) MS = HasMeanStruct μ_sym = eval_μ_symbolic(M, I_A⁻¹, F; simplify = simplify_symbolics) μ_eval! = Symbolics.build_function(μ_sym, par, expression = Val{false})[2] @@ -230,10 +220,10 @@ end ############################################################################################ # expected covariations of observed vars -function eval_Σ_symbolic(S, I_A⁻¹, F; vech = false, simplify = false) +function eval_Σ_symbolic(S, I_A⁻¹, F; vech::Bool = false, simplify::Bool = false) Σ = F * I_A⁻¹ * S * permutedims(I_A⁻¹) * permutedims(F) Σ = Array(Σ) - vech && (Σ = Σ[tril(trues(size(F, 1), size(F, 1)))]) + vech && (Σ = SEM.vech(Σ)) if simplify Threads.@threads for i in eachindex(Σ) Σ[i] = Symbolics.simplify(Σ[i]) diff --git a/src/implied/abstract.jl b/src/implied/abstract.jl index 6d298f65..d4868d74 100644 --- a/src/implied/abstract.jl +++ b/src/implied/abstract.jl @@ -31,17 +31,3 @@ function check_acyclic(A::AbstractMatrix; verbose::Bool = false) return A end end - -# Verify that the `meanstructure` argument aligns with the model specification. -function check_meanstructure_specification(meanstructure, ram_matrices) - if meanstructure & isnothing(ram_matrices.M) - throw(ArgumentError( - "You set `meanstructure = true`, but your model specification contains no mean parameters." - )) - end - if !meanstructure & !isnothing(ram_matrices.M) - throw(ArgumentError( - "If your model specification contains mean parameters, you have to set `Sem(..., meanstructure = true)`." - )) - end -end \ No newline at end of file diff --git a/src/loss/ML/FIML.jl b/src/loss/ML/FIML.jl index da5ccb7c..8572b15a 100644 --- a/src/loss/ML/FIML.jl +++ b/src/loss/ML/FIML.jl @@ -75,15 +75,16 @@ Can handle observed data with missing values. # Constructor - SemFIML(; observed::SemObservedMissing, specification, kwargs...) + SemFIML(; observed::SemObservedMissing, implied::SemImplied, kwargs...) # Arguments -- `observed`: the observed data with missing values (see [`SemObservedMissing`](@ref)) -- `specification`: [`SemSpecification`](@ref) object +- `observed::SemObservedMissing`: the observed part of the model + (see [`SemObservedMissing`](@ref)) +- `implied::SemImplied`: the implied part of the model # Examples ```julia -my_fiml = SemFIML(observed = my_observed, specification = my_parameter_table) +my_fiml = SemFIML(observed = my_observed, implied = my_implied) ``` # Interfaces @@ -118,7 +119,7 @@ function SemFIML(; observed::SemObservedMissing, implied, specification, kwargs. ExactHessian(), [SemFIMLPattern(pat) for pat in observed.patterns], zeros(nobserved_vars(observed), nobserved_vars(observed)), - CommutationMatrix(nvars(specification)), + CommutationMatrix(nvars(implied)), nothing, ) end diff --git a/src/loss/ML/ML.jl b/src/loss/ML/ML.jl index ce77ea9c..aae1dada 100644 --- a/src/loss/ML/ML.jl +++ b/src/loss/ML/ML.jl @@ -8,11 +8,10 @@ Maximum likelihood estimation. # Constructor - SemML(;observed, meanstructure = false, approximate_hessian = false, kwargs...) + SemML(; observed, approximate_hessian = false, kwargs...) # Arguments - `observed::SemObserved`: the observed part of the model -- `meanstructure::Bool`: does the model have a meanstructure? - `approximate_hessian::Bool`: if hessian-based optimization is used, should the hessian be swapped for an approximation # Examples diff --git a/src/loss/WLS/WLS.jl b/src/loss/WLS/WLS.jl index b7f66d55..9de011f6 100644 --- a/src/loss/WLS/WLS.jl +++ b/src/loss/WLS/WLS.jl @@ -10,8 +10,7 @@ At the moment only available with the `RAMSymbolic` implied type. # Constructor SemWLS(; - observed, - meanstructure = false, + observed, implied, wls_weight_matrix = nothing, wls_weight_matrix_mean = nothing, approximate_hessian = false, @@ -19,7 +18,7 @@ At the moment only available with the `RAMSymbolic` implied type. # Arguments - `observed`: the `SemObserved` part of the model -- `meanstructure::Bool`: does the model have a meanstructure? +- `implied::SemImplied`: the implied part of the model - `approximate_hessian::Bool`: should the hessian be swapped for an approximation - `wls_weight_matrix`: the weight matrix for weighted least squares. Defaults to GLS estimation (``0.5*(D^T*kron(S,S)*D)`` where D is the duplication matrix @@ -29,7 +28,7 @@ At the moment only available with the `RAMSymbolic` implied type. # Examples ```julia -my_wls = SemWLS(observed = my_observed) +my_wls = SemWLS(observed = my_observed, implied = my_implied) ``` # Interfaces @@ -50,12 +49,11 @@ SemWLS{HE}(args...) where {HE <: HessianEval} = SemWLS{HE, map(typeof, args)...}(HE(), args...) function SemWLS(; - observed, - implied, - wls_weight_matrix = nothing, - wls_weight_matrix_mean = nothing, - approximate_hessian = false, - meanstructure = false, + observed::SemObserved, + implied::SemImplied, + wls_weight_matrix::Union{AbstractMatrix, Nothing} = nothing, + wls_weight_matrix_mean::Union{AbstractMatrix, Nothing} = nothing, + approximate_hessian::Bool = false, kwargs..., ) if observed isa SemObservedMissing @@ -81,6 +79,10 @@ function SemWLS(; nobs_vars = nobserved_vars(observed) tril_ind = filter(x -> (x[1] >= x[2]), CartesianIndices(obs_cov(observed))) s = obs_cov(observed)[tril_ind] + size(s) == size(implied.Σ) || + throw(DimensionMismatch("SemWLS requires implied covariance to be in vech-ed form " * + "(vectorized lower triangular part of Σ matrix): $(size(s)) expected, $(size(implied.Σ)) found.\n" * + "$(nameof(typeof(implied))) must be constructed with vech=true.")) # compute V here if isnothing(wls_weight_matrix) @@ -94,9 +96,12 @@ function SemWLS(; "wls_weight_matrix has to be of size $(length(tril_ind))×$(length(tril_ind))", ) end + size(wls_weight_matrix) == (length(s), length(s)) || + DimensionMismatch("wls_weight_matrix has to be of size $(length(s))×$(length(s))") - if meanstructure + if MeanStruct(implied) == HasMeanStruct if isnothing(wls_weight_matrix_mean) + @warn "Computing WLS weight matrix for the meanstructure using obs_cov()" wls_weight_matrix_mean = inv(obs_cov(observed)) else size(wls_weight_matrix_mean) == (nobs_vars, nobs_vars) || DimensionMismatch( diff --git a/test/examples/recover_parameters/recover_parameters_twofact.jl b/test/examples/recover_parameters/recover_parameters_twofact.jl index 9f9503af..a4bd7d5f 100644 --- a/test/examples/recover_parameters/recover_parameters_twofact.jl +++ b/test/examples/recover_parameters/recover_parameters_twofact.jl @@ -53,7 +53,7 @@ start = [ repeat([0.5], 4) ] -implied_ml = RAMSymbolic(; specification = ram_matrices, start_val = start) +implied_ml = RAMSymbolic(ram_matrices; start_val = start) implied_ml.Σ_eval!(implied_ml.Σ, true_val) diff --git a/test/unit_tests/model.jl b/test/unit_tests/model.jl index d9f9254b..fbe2a937 100644 --- a/test/unit_tests/model.jl +++ b/test/unit_tests/model.jl @@ -46,13 +46,16 @@ function test_params_api(semobj, spec::SemSpecification) @test @inferred(param_labels(semobj)) == param_labels(spec) end -@testset "Sem(implied=$impliedtype, loss=SemML)" for impliedtype in (RAM, RAMSymbolic) - +@testset "Sem(implied=$impliedtype, loss=$losstype)" for (impliedtype, losstype) in [ + (RAM, SemML), + (RAMSymbolic, SemML), + (RAMSymbolic, SemWLS), +] model = Sem( specification = ram_matrices, observed = obs, implied = impliedtype, - loss = SemML, + loss = losstype, ) @test model isa Sem @@ -71,29 +74,3 @@ end @test @inferred(nsamples(model)) == nsamples(obs) end - -@testset "Sem(implied=RAMSymbolic, loss=SemWLS)" begin - - model = Sem( - specification = ram_matrices, - observed = obs, - implied = RAMSymbolic, - loss = SemWLS, - ) - - @test model isa Sem - @test @inferred(implied(model)) isa RAMSymbolic - @test @inferred(observed(model)) isa SemObserved - - test_vars_api(model, ram_matrices) - test_params_api(model, ram_matrices) - - test_vars_api(implied(model), ram_matrices) - test_params_api(implied(model), ram_matrices) - - @test @inferred(loss(model)) isa SemLoss - semloss = loss(model).functions[1] - @test semloss isa SemWLS - - @test @inferred(nsamples(model)) == nsamples(obs) -end \ No newline at end of file From e81cec008cb99ef492841bd1c4e0503ac45cca38 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 21 Mar 2026 11:17:30 -0700 Subject: [PATCH 04/18] refactor Sem, SemEnsemble, SemLoss --- src/StructuralEquationModels.jl | 17 +- src/additional_functions/simulation.jl | 87 ++-- src/additional_functions/start_val/common.jl | 17 + .../start_val/start_fabin3.jl | 16 +- .../start_val/start_simple.jl | 34 +- src/frontend/finite_diff.jl | 35 ++ src/frontend/fit/fitmeasures/RMSEA.jl | 32 +- src/frontend/fit/fitmeasures/chi2.jl | 73 ++- src/frontend/fit/fitmeasures/dof.jl | 11 +- src/frontend/fit/fitmeasures/minus2ll.jl | 34 +- src/frontend/fit/standard_errors/hessian.jl | 25 +- src/frontend/pretty_printing.jl | 8 +- src/frontend/specification/Sem.jl | 439 +++++++++++++----- src/implied/RAM/generic.jl | 13 +- src/implied/RAM/symbolic.jl | 12 +- src/implied/empty.jl | 6 +- src/loss/ML/FIML.jl | 86 ++-- src/loss/ML/ML.jl | 90 ++-- src/loss/WLS/WLS.jl | 94 ++-- src/loss/abstract.jl | 42 ++ src/loss/constant/constant.jl | 28 +- src/loss/regularization/ridge.jl | 9 +- src/objective_gradient_hessian.jl | 251 +++++----- src/optimizer/abstract.jl | 14 +- src/types.jl | 243 +++------- test/examples/multigroup/build_models.jl | 222 +++------ test/examples/political_democracy/by_parts.jl | 284 +++++------ .../political_democracy/constraints.jl | 4 +- .../political_democracy/constructor.jl | 169 +++---- .../recover_parameters_twofact.jl | 24 +- test/unit_tests/model.jl | 5 +- 31 files changed, 1203 insertions(+), 1221 deletions(-) create mode 100644 src/additional_functions/start_val/common.jl create mode 100644 src/frontend/finite_diff.jl create mode 100644 src/loss/abstract.jl diff --git a/src/StructuralEquationModels.jl b/src/StructuralEquationModels.jl index 19dd6f43..d98e7925 100644 --- a/src/StructuralEquationModels.jl +++ b/src/StructuralEquationModels.jl @@ -44,6 +44,7 @@ include("frontend/specification/EnsembleParameterTable.jl") include("frontend/specification/StenoGraphs.jl") include("frontend/fit/summary.jl") include("frontend/StatsAPI.jl") +include("frontend/finite_diff.jl") # pretty printing include("frontend/pretty_printing.jl") # observed @@ -53,26 +54,28 @@ include("observed/covariance.jl") include("observed/missing_pattern.jl") include("observed/missing.jl") include("observed/EM.jl") -# constructor -include("frontend/specification/Sem.jl") -include("frontend/specification/documentation.jl") # implied include("implied/abstract.jl") include("implied/RAM/symbolic.jl") include("implied/RAM/generic.jl") include("implied/empty.jl") # loss +include("loss/abstract.jl") include("loss/ML/ML.jl") include("loss/ML/FIML.jl") include("loss/regularization/ridge.jl") include("loss/WLS/WLS.jl") include("loss/constant/constant.jl") +# constructor +include("frontend/specification/Sem.jl") +include("frontend/specification/documentation.jl") # optimizer include("optimizer/abstract.jl") include("optimizer/Empty.jl") include("optimizer/optim.jl") # helper functions include("additional_functions/helper.jl") +include("additional_functions/start_val/common.jl") include("additional_functions/start_val/start_fabin3.jl") include("additional_functions/start_val/start_simple.jl") include("additional_functions/artifacts.jl") @@ -94,14 +97,11 @@ include("frontend/fit/standard_errors/z_test.jl") include("frontend/fit/standard_errors/confidence_intervals.jl") export AbstractSem, - AbstractSemSingle, - AbstractSemCollection, coef, coefnames, coeftable, Sem, SemFiniteDiff, - SemEnsemble, MeanStruct, NoMeanStruct, HasMeanStruct, @@ -116,8 +116,8 @@ export AbstractSem, start_val, start_fabin3, start_simple, + AbstractLoss, SemLoss, - SemLossFunction, SemML, SemFIML, em_mvn, @@ -125,6 +125,9 @@ export AbstractSem, SemConstant, SemWLS, loss, + nsem_terms, + sem_terms, + sem_term, SemOptimizer, optimizer, optimizer_engine, diff --git a/src/additional_functions/simulation.jl b/src/additional_functions/simulation.jl index 4839bc27..6d694c97 100644 --- a/src/additional_functions/simulation.jl +++ b/src/additional_functions/simulation.jl @@ -43,36 +43,42 @@ function update_observed end # change observed (data) without reconstructing the whole model ############################################################################################ +# don't change non-SEM terms +replace_observed(loss::AbstractLoss; kwargs...) = loss + # use the same observed type as before -replace_observed(model::AbstractSemSingle; kwargs...) = - replace_observed(model, typeof(observed(model)).name.wrapper; kwargs...) +replace_observed(loss::SemLoss; kwargs...) = + replace_observed(loss, typeof(SEM.observed(loss)).name.wrapper; kwargs...) + +# construct a new observed type +replace_observed(loss::SemLoss, observed_type; kwargs...) = + replace_observed(loss, observed_type(; kwargs...); kwargs...) -function replace_observed(model::AbstractSemSingle, observed_type; kwargs...) - new_observed = observed_type(; kwargs...) +function replace_observed(loss::SemLoss, new_observed::SemObserved; kwargs...) kwargs = Dict{Symbol, Any}(kwargs...) + old_observed = SEM.observed(loss) + implied = SEM.implied(loss) # get field types kwargs[:observed_type] = typeof(new_observed) - kwargs[:old_observed_type] = typeof(model.observed) - kwargs[:implied_type] = typeof(model.implied) - kwargs[:loss_types] = [typeof(lossfun) for lossfun in model.loss.functions] + kwargs[:old_observed_type] = typeof(old_observed) # update implied - new_implied = update_observed(model.implied, new_observed; kwargs...) + new_implied = update_observed(implied, new_observed; kwargs...) kwargs[:implied] = new_implied + kwargs[:implied_type] = typeof(new_implied) kwargs[:nparams] = nparams(new_implied) # update loss - new_loss = update_observed(model.loss, new_observed; kwargs...) - - return Sem(new_observed, new_implied, new_loss) + return update_observed(loss, new_observed; kwargs...) end -function update_observed(loss::SemLoss, new_observed; kwargs...) - new_functions = Tuple( - update_observed(lossfun, new_observed; kwargs...) for lossfun in loss.functions - ) - return SemLoss(new_functions, loss.weights) +replace_observed(loss::LossTerm; kwargs...) = + LossTerm(replace_observed(loss.loss; kwargs...), loss.id, loss.weight) + +function replace_observed(sem::Sem; kwargs...) + updated_terms = Tuple(replace_observed(term; kwargs...) for term in loss_terms(sem)) + return Sem(updated_terms...) end function replace_observed( @@ -111,39 +117,38 @@ end # simulate data ############################################################################################ """ - (1) rand(model::AbstractSemSingle, params, n) - - (2) rand(model::AbstractSemSingle, n) + rand(sem::Union{Sem, SemLoss, SemImplied}, [params], n) -Sample normally distributed data from the model-implied covariance matrix and mean vector. +Sample from the multivariate normal distribution implied by the SEM model. # Arguments -- `model::AbstractSemSingle`: model to simulate from. -- `params`: parameter values to simulate from. -- `n::Integer`: Number of samples. +- `sem`: SEM model to use. Ensemble models with multiple SEM terms are not supported. +- `params`: optional SEM model parameters to simulate from, otherwise uses the + current state of implied covariances and means. +- `n::Integer`: Number of samples to draw. # Examples ```julia rand(model, start_simple(model), 100) ``` """ -function Distributions.rand( - model::AbstractSemSingle{O, I, L}, - params, - n::Integer, -) where {O, I <: Union{RAM, RAMSymbolic}, L} - update!(EvaluationTargets{true, false, false}(), model.implied, model, params) - return rand(model, n) -end - -function Distributions.rand( - model::AbstractSemSingle{O, I, L}, - n::Integer, -) where {O, I <: Union{RAM, RAMSymbolic}, L} - if MeanStruct(model.implied) === NoMeanStruct - data = permutedims(rand(MvNormal(Symmetric(model.implied.Σ)), n)) - elseif MeanStruct(model.implied) === HasMeanStruct - data = permutedims(rand(MvNormal(model.implied.μ, Symmetric(model.implied.Σ)), n)) +function Distributions.rand(implied::SemImplied, params, n::Integer) + if !isnothing(params) + # update the implied covariances with the new model params + update!(EvaluationTargets{true, false, false}(), implied, params) + end + Σ = Symmetric(implied.Σ) + if MeanStruct(implied) === NoMeanStruct + return permutedims(rand(MvNormal(Σ), n)) + elseif MeanStruct(implied) === HasMeanStruct + return permutedims(rand(MvNormal(implied.μ, Σ), n)) end - return data end + +Distributions.rand(loss::SemLoss, params, n::Integer) = rand(SEM.implied(loss), params, n) + +Distributions.rand(model::Sem, params, n::Integer) = rand(sem_term(model), params, n) + +# rand() overloads without SEM params +Distributions.rand(implied::Union{SemImplied, SemLoss, Sem}, n::Integer) = + Distributions.rand(implied, nothing, n) diff --git a/src/additional_functions/start_val/common.jl b/src/additional_functions/start_val/common.jl new file mode 100644 index 00000000..92c85d6f --- /dev/null +++ b/src/additional_functions/start_val/common.jl @@ -0,0 +1,17 @@ + +# start values for SEM Models (including ensembles) +function start_values(f, model::AbstractSem; kwargs...) + start_vals = fill(0.0, nparams(model)) + + # initialize parameters using the SEM loss terms + # (first SEM loss term that sets given parameter to nonzero value) + for term in loss_terms(model) + issemloss(term) || continue + term_start_vals = f(loss(term); kwargs...) + for (i, val) in enumerate(term_start_vals) + iszero(val) || (start_vals[i] = val) + end + end + + return start_vals +end diff --git a/src/additional_functions/start_val/start_fabin3.jl b/src/additional_functions/start_val/start_fabin3.jl index 53d3442a..54337028 100644 --- a/src/additional_functions/start_val/start_fabin3.jl +++ b/src/additional_functions/start_val/start_fabin3.jl @@ -7,12 +7,17 @@ Not available for ensemble models. function start_fabin3 end # splice model and loss functions -function start_fabin3(model::AbstractSemSingle; kwargs...) - return start_fabin3(model.observed, model.implied, model.loss.functions..., kwargs...) +function start_fabin3(model::SemLoss; kwargs...) + return start_fabin3(model.observed, model.implied; kwargs...) end -function start_fabin3(observed::SemObserved, implied::SemImplied, args...; kwargs...) - return start_fabin3(implied.ram_matrices, obs_cov(observed), obs_mean(observed)) +function start_fabin3(observed::SemObserved, implied::SemImplied; kwargs...) + return start_fabin3( + implied.ram_matrices, + obs_cov(observed), + # ignore observed means if no meansturcture + !isnothing(implied.ram_matrices.M) ? obs_mean(observed) : nothing, + ) end function start_fabin3( @@ -161,3 +166,6 @@ end function is_in_Λ(ind_vec, F_ind) return any(ind -> !(ind[2] ∈ F_ind) && (ind[1] ∈ F_ind), ind_vec) end + +# ensembles +start_fabin3(model::AbstractSem; kwargs...) = start_values(start_fabin3, model; kwargs...) diff --git a/src/additional_functions/start_val/start_simple.jl b/src/additional_functions/start_val/start_simple.jl index 4fbc8719..afdbf92e 100644 --- a/src/additional_functions/start_val/start_simple.jl +++ b/src/additional_functions/start_val/start_simple.jl @@ -15,34 +15,11 @@ Return a vector of simple starting values. """ function start_simple end -# Single Models ---------------------------------------------------------------------------- -function start_simple(model::AbstractSemSingle; kwargs...) - return start_simple(model.observed, model.implied, model.loss.functions...; kwargs...) -end - -function start_simple(observed, implied, args...; kwargs...) - return start_simple(implied.ram_matrices; kwargs...) -end - -# Ensemble Models -------------------------------------------------------------------------- -function start_simple(model::SemEnsemble; kwargs...) - start_vals = [] - - for sem in model.sems - push!(start_vals, start_simple(sem; kwargs...)) - end - - has_start_val = [.!iszero.(start_val) for start_val in start_vals] +start_simple(model::SemLoss; kwargs...) = + start_simple(observed(model), implied(model); kwargs...) - start_val = similar(start_vals[1]) - start_val .= 0.0 - - for (j, indices) in enumerate(has_start_val) - start_val[indices] .= start_vals[j][indices] - end - - return start_val -end +start_simple(observed::SemObserved, implied::SemImplied; kwargs...) = + start_simple(implied.ram_matrices; kwargs...) function start_simple( ram_matrices::RAMMatrices; @@ -103,3 +80,6 @@ function start_simple( end return start_val end + +# multigroup models +start_simple(model::AbstractSem; kwargs...) = start_values(start_simple, model; kwargs...) diff --git a/src/frontend/finite_diff.jl b/src/frontend/finite_diff.jl new file mode 100644 index 00000000..ee0a9bf9 --- /dev/null +++ b/src/frontend/finite_diff.jl @@ -0,0 +1,35 @@ +_unwrap(wrapper::SemFiniteDiff) = wrapper.model +params(wrapper::SemFiniteDiff) = params(wrapper.model) +loss_terms(wrapper::SemFiniteDiff) = loss_terms(wrapper.model) + +FiniteDiffLossWrappers = Union{LossFiniteDiff, SemLossFiniteDiff} + +_unwrap(term::AbstractLoss) = term +_unwrap(wrapper::FiniteDiffLossWrappers) = wrapper.loss +implied(wrapper::FiniteDiffLossWrappers) = implied(_unwrap(wrapper)) +observed(wrapper::FiniteDiffLossWrappers) = observed(_unwrap(wrapper)) + +FiniteDiffWrapper(model::AbstractSem) = SemFiniteDiff(model) +FiniteDiffWrapper(loss::AbstractLoss) = LossFiniteDiff(loss) +FiniteDiffWrapper(loss::SemLoss) = SemLossFiniteDiff(loss) + +function evaluate!( + objective, + gradient, + hessian, + sem::Union{SemFiniteDiff, FiniteDiffLossWrappers}, + params, +) + wrapped = _unwrap(sem) + obj(p) = _evaluate!( + objective_zero(objective, gradient, hessian), + nothing, + nothing, + wrapped, + p, + ) + isnothing(gradient) || FiniteDiff.finite_difference_gradient!(gradient, obj, params) + isnothing(hessian) || FiniteDiff.finite_difference_hessian!(hessian, obj, params) + # FIXME if objective is not calculated, SemLoss implied states may not correspond to params + return !isnothing(objective) ? obj(params) : nothing +end diff --git a/src/frontend/fit/fitmeasures/RMSEA.jl b/src/frontend/fit/fitmeasures/RMSEA.jl index 8539896f..7406b74c 100644 --- a/src/frontend/fit/fitmeasures/RMSEA.jl +++ b/src/frontend/fit/fitmeasures/RMSEA.jl @@ -19,28 +19,18 @@ for the SEM model. For multigroup models, the correction proposed by J.H. Steiger is applied (see [Steiger, J. H. (1998). *A note on multiple sample extensions of the RMSEA fit index*](https://doi.org/10.1080/10705519809540115)). """ -function RMSEA end - RMSEA(fit::SemFit) = RMSEA(fit, fit.model) -function RMSEA(fit::SemFit, model::AbstractSemSingle) - check_single_lossfun(model; throw_error = true) - return RMSEA(dof(fit), χ²(fit), nsamples(fit)+rmsea_correction(model.loss.functions[1])) -end - -function RMSEA(fit::SemFit, model::SemEnsemble) - check_single_lossfun(model; throw_error = true) - n = nsamples(fit)+model.n*rmsea_correction(model.sems[1].loss.functions[1]) - return sqrt(length(model.sems)) * RMSEA(dof(fit), χ²(fit), n) -end - -function RMSEA(dof, chi2, N⁻) - rmsea = (chi2 - dof) / (N⁻ * dof) - rmsea = rmsea > 0 ? rmsea : 0 - return sqrt(rmsea) +# scaling corrections +RMSEA_corr_scale(::Type{<:SemFIML}) = 0 +RMSEA_corr_scale(::Type{<:SemML}) = -1 +RMSEA_corr_scale(::Type{<:SemWLS}) = -1 + +function RMSEA(fit::SemFit, model::AbstractSem) + term_type = check_single_lossfun(model; throw_error = true) + n = nsamples(fit) + nsem_terms(model) * RMSEA_corr_scale(term_type) + sqrt(nsem_terms(model)) * RMSEA(dof(fit), χ²(fit), n) end -# scaling corrections -rmsea_correction(::SemFIML) = 0 -rmsea_correction(::SemML) = -1 -rmsea_correction(::SemWLS) = -1 +RMSEA(dof::Number, chi2::Number, nsamples::Number) = + sqrt(max((chi2 - dof) / (nsamples * dof), 0.0)) diff --git a/src/frontend/fit/fitmeasures/chi2.jl b/src/frontend/fit/fitmeasures/chi2.jl index 8ce5f079..22d6c2e2 100644 --- a/src/frontend/fit/fitmeasures/chi2.jl +++ b/src/frontend/fit/fitmeasures/chi2.jl @@ -12,57 +12,42 @@ with the *observed* covariance matrix. """ χ²(fit::SemFit) = χ²(fit, fit.model) -############################################################################################ -# Single Models -############################################################################################ +function χ²(fit::SemFit, model::AbstractSem) + terms = sem_terms(model) + isempty(terms) && return 0.0 + + term1 = _unwrap(loss(terms[1])) + L = typeof(term1).name + + # check that all SemLoss terms are of the same class (ML, FIML, WLS etc), ignore typeparams + for (i, term) in enumerate(terms) + lossterm = _unwrap(loss(term)) + @assert lossterm isa SemLoss + if typeof(_unwrap(lossterm)).name != L + @error "SemLoss term #$i is $(typeof(_unwrap(lossterm)).name), expected $L. Heterogeneous loss functions are not supported" + end + end -function χ²(fit::SemFit, model::AbstractSemSingle) - check_single_lossfun(model; throw_error = true) - return χ²(model.loss.functions[1], fit::SemFit, model::AbstractSemSingle) + return χ²(typeof(term1), fit, model) end -χ²(::SemML, fit::SemFit, model::AbstractSemSingle) = - (nsamples(fit) - 1) * - (fit.minimum - logdet(obs_cov(observed(model))) - nobserved_vars(model)) - # bollen, p. 115, only correct for GLS weight matrix -χ²(::SemWLS, fit::SemFit, model::AbstractSemSingle) = - (nsamples(fit) - 1) * fit.minimum - -# FIML -function χ²(::SemFIML, fit::SemFit, model::AbstractSemSingle) - ll_H0 = minus2ll(fit) - ll_H1 = minus2ll(observed(model)) - return ll_H0 - ll_H1 -end - -############################################################################################ -# Collections -############################################################################################ - -function χ²(fit::SemFit, model::SemEnsemble) - check_single_lossfun(model; throw_error = true) - lossfun = model.sems[1].loss.functions[1] - return χ²(lossfun, fit, model) -end - -function χ²(::SemWLS, fit::SemFit, models::SemEnsemble) - return (nsamples(models) - models.n) * fit.minimum -end - -function χ²(::SemML, fit::SemFit, models::SemEnsemble) - F = 0 - for model in models.sems - Fᵢ = objective(model, fit.solution) - Fᵢ -= logdet(obs_cov(observed(model))) + nobserved_vars(model) - Fᵢ *= nsamples(model) - 1 - F += Fᵢ +χ²(::Type{<:SemWLS}, fit::SemFit, model::AbstractSem) = (nsamples(model) - 1) * fit.minimum + +function χ²(::Type{<:SemML}, fit::SemFit, model::AbstractSem) + G = sum(loss_terms(model)) do term + if issemloss(term) + data = observed(term) + something(weight(term), 1.0) * (logdet(obs_cov(data)) + nobserved_vars(data)) + else + return 0.0 + end end - return F + return (nsamples(model) - 1) * (fit.minimum - G) end -function χ²(::SemFIML, fit::SemFit, models::SemEnsemble) +function χ²(::Type{<:SemFIML}, fit::SemFit, model::AbstractSem) ll_H0 = minus2ll(fit) - ll_H1 = sum(minus2ll ∘ observed, models.sems) + ll_H1 = sum(minus2ll ∘ observed, sem_terms(model)) return ll_H0 - ll_H1 end diff --git a/src/frontend/fit/fitmeasures/dof.jl b/src/frontend/fit/fitmeasures/dof.jl index 0e051d02..49b7febf 100644 --- a/src/frontend/fit/fitmeasures/dof.jl +++ b/src/frontend/fit/fitmeasures/dof.jl @@ -18,13 +18,16 @@ dof(fit::SemFit) = dof(fit.model) dof(model::AbstractSem) = n_dp(model) - nparams(model) -function n_dp(model::AbstractSemSingle) - nvars = nobserved_vars(model) +# length of Σ and μ (if present) +function n_dp(implied::SemImplied) + nvars = nobserved_vars(implied) ndp = 0.5(nvars^2 + nvars) - if !isnothing(model.implied.μ) + if !isnothing(implied.μ) ndp += nvars end return ndp end -n_dp(model::SemEnsemble) = sum(n_dp.(model.sems)) +n_dp(term::SemLoss) = n_dp(implied(term)) + +n_dp(model::AbstractSem) = sum(n_dp ∘ loss, sem_terms(model)) diff --git a/src/frontend/fit/fitmeasures/minus2ll.jl b/src/frontend/fit/fitmeasures/minus2ll.jl index c6a954ef..3b353f5c 100644 --- a/src/frontend/fit/fitmeasures/minus2ll.jl +++ b/src/frontend/fit/fitmeasures/minus2ll.jl @@ -6,31 +6,27 @@ Calculate the *-2⋅log(likelihood(fit))*. # See also [`fit_measures`](@ref) """ -minus2ll(fit::SemFit) = minus2ll(fit, fit.model) +minus2ll(fit::SemFit) = minus2ll(fit.model, fit) ############################################################################################ -# Single Models +# Single SEM Terms Models ############################################################################################ -function minus2ll(fit::SemFit, model::AbstractSemSingle) - check_single_lossfun(model; throw_error = true) - F = objective(model, fit.solution) - return minus2ll(model.loss.functions[1], F, model) +function minus2ll(term::SemLoss, fit::SemFit) + minimum = objective(term, fit.solution) + return minus2ll(term, minimum) end -# SemML ------------------------------------------------------------------------------------ -function minus2ll(::SemML, F, model::AbstractSemSingle) - return nsamples(model) * (F + log(2π) * nobserved_vars(model)) -end +minus2ll(term::SemML, minimum::Number) = + nsamples(term) * (minimum + log(2π) * nobserved_vars(term)) -# WLS -------------------------------------------------------------------------------------- -minus2ll(::SemWLS, F, ::AbstractSemSingle) = missing +minus2ll(term::SemWLS, minimum::Number) = missing # compute likelihood for missing data - H0 ------------------------------------------------- -# -2ll = (∑ log(2π)*(nᵢ*mᵢ)) + F*n -function minus2ll(::SemFIML, F, model::AbstractSemSingle) - obs = observed(model)::SemObservedMissing - F *= nsamples(obs) +# -2ll = (∑ log(2π)*(nᵢ + mᵢ)) + F*n +function minus2ll(term::SemFIML, minimum::Number) + obs = observed(term)::SemObservedMissing + F = minimum * nsamples(obs) F += log(2π) * sum(pat -> nsamples(pat) * nmeasured_vars(pat), obs.patterns) return F end @@ -62,10 +58,10 @@ function minus2ll(observed::SemObservedMissing) end ############################################################################################ -# Collection +# Multi-group ############################################################################################ -function minus2ll(fit::SemFit, model::SemEnsemble) +function minus2ll(model::AbstractSem, fit::SemFit) check_single_lossfun(model; throw_error = true) - return sum(Base.Fix1(minus2ll, fit), model.sems) + sum(Base.Fix2(minus2ll, fit) ∘ _unwrap ∘ loss, sem_terms(model)) end diff --git a/src/frontend/fit/standard_errors/hessian.jl b/src/frontend/fit/standard_errors/hessian.jl index 6ae53407..80b96d33 100644 --- a/src/frontend/fit/standard_errors/hessian.jl +++ b/src/frontend/fit/standard_errors/hessian.jl @@ -35,20 +35,21 @@ function se_hessian(fit::SemFit; method = :finitediff) end # Addition functions ------------------------------------------------------------- -function H_scaling(model::AbstractSemSingle) - if length(model.loss.functions) > 1 - @warn "Hessian scaling for multiple loss functions is not implemented yet" - end - return H_scaling(model.loss.functions[1], model) -end - -H_scaling(lossfun::SemML, model::AbstractSemSingle) = 2 / (nsamples(model) - 1) +H_scaling(loss::SemML) = 2 / (nsamples(loss) - 1) -function H_scaling(lossfun::SemWLS, model::AbstractSemSingle) +function H_scaling(loss::SemWLS) @warn "Standard errors for WLS are only correct if a GLS weight matrix (the default) is used." - return 2 / (nsamples(model) - 1) + return 2 / (nsamples(loss) - 1) end -H_scaling(lossfun::SemFIML, model::AbstractSemSingle) = 2 / nsamples(model) +H_scaling(loss::SemFIML) = 2 / nsamples(loss) -H_scaling(model::SemEnsemble) = 2 / nsamples(model) +function H_scaling(model::AbstractSem) + semterms = SEM.sem_terms(model) + if length(semterms) > 1 + #@warn "Hessian scaling for multiple loss functions is not implemented yet" + return 2 / nsamples(model) + else + return length(semterms) >= 1 ? H_scaling(loss(semterms[1])) : 1.0 + end +end diff --git a/src/frontend/pretty_printing.jl b/src/frontend/pretty_printing.jl index 2fa970f2..7b6975f6 100644 --- a/src/frontend/pretty_printing.jl +++ b/src/frontend/pretty_printing.jl @@ -32,9 +32,11 @@ end # Loss Function, Implied, Observed, Optimizer ############################################################## -function Base.show(io::IO, struct_inst::SemLossFunction) - print_type_name(io, struct_inst) - print_field_types(io, struct_inst) +function Base.show(io::IO, sem::SemLoss) + println(io, "Structural Equation Model Loss ($(nameof(typeof(sem))))") + println(io, "- Observed: $(nameof(typeof(observed(sem)))) ($(nsamples(sem)) samples)") + println(io, "- Implied: $(nameof(typeof(implied(sem)))) ($(nparams(sem)) parameters)") + println(io, "- Variables: $(nobserved_vars(sem)) observed, $(nlatent_vars(sem)) latent") end function Base.show(io::IO, struct_inst::SemImplied) diff --git a/src/frontend/specification/Sem.jl b/src/frontend/specification/Sem.jl index a47bad4b..684cfa62 100644 --- a/src/frontend/specification/Sem.jl +++ b/src/frontend/specification/Sem.jl @@ -1,7 +1,167 @@ +losstype(::Type{<:LossTerm{L, W, I}}) where {L, W, I} = L +losstype(term::LossTerm) = losstype(typeof(term)) +loss(term::LossTerm) = term.loss +weight(term::LossTerm) = term.weight +id(term::LossTerm) = term.id + +""" + issemloss(term::LossTerm) -> Bool + +Check if a SEM model term is a SEM loss function ([`SemLoss`](@ref)). +""" +issemloss(term::LossTerm) = isa(loss(term), SemLoss) + +for f in ( + :implied, + :observed, + :nsamples, + :observed_vars, + :nobserved_vars, + :vars, + :nvars, + :latent_vars, + :nlatent_vars, + :params, + :nparams, +) + @eval $f(term::LossTerm) = $f(loss(term)) +end + +function Base.show(io::IO, term::LossTerm) + if !isnothing(id(term)) + print(io, ":$(id(term)): ") + end + print(io, nameof(losstype(term))) + if issemloss(term) + print( + io, + " ($(nsamples(term)) samples, $(nobserved_vars(term)) observed, $(nlatent_vars(term)) latent variables)", + ) + end + if !isnothing(weight(term)) + print(io, " w=$(round(weight(term), digits=3))") + else + print(io, " w=1") + end +end + ############################################################################################ # constructor for Sem types ############################################################################################ +function multigroup_weights(models, n) + nsamples_total = sum(nsamples, models) + uniform_lossfun = check_single_lossfun(models...; throw_error = false) + if !uniform_lossfun + @info """ + Your ensemble model contains heterogeneous loss functions. + Default weights of (#samples per group/#total samples) will be used + """ + return [(nsamples(model)) / (nsamples_total) for model in models] + end + lossfun = models[1].loss.functions[1] + if !applicable(mg_correction, lossfun) + @info """ + We don't know how to choose group weights for the specified loss function. + Default weights of (#samples per group/#total samples) will be used + """ + return [(nsamples(model)) / (nsamples_total) for model in models] + end + c = mg_correction(lossfun) + return [(nsamples(model)+c) / (nsamples_total+n*c) for model in models] +end + +function Sem( + loss_terms...; + params::Union{Vector{Symbol}, Nothing} = nothing, + default_sem_weights = :nsamples, +) + default_sem_weights ∈ [:nsamples, :uniform, :one] || + throw(ArgumentError("Unsupported default_sem_weights=:$default_sem_weights")) + # assemble a list of weighted losses and check params equality + terms = Vector{LossTerm}() + params = !isnothing(params) ? copy(params) : params + has_sem_weights = false + nsems = 0 + for inp_term in loss_terms + if inp_term isa AbstractLoss # term + term = inp_term + term_w = nothing + term_id = nothing + elseif inp_term isa Pair + if inp_term[1] isa AbstractLoss # term => weight + term, term_w = inp_term + term_id = nothing + elseif inp_term[2] isa AbstractLoss # id => term + term_id, term = inp_term + term_w = nothing + elseif inp_term[2] isa Pair # id => term => weight + term_id, (term, term_w) = inp_term + isa(term, AbstractLoss) || throw( + ArgumentError( + "AbstractLoss expected as a second argument of a loss term double pair (id => loss => weight), $(nameof(typeof(term))) found", + ), + ) + end + elseif inp_term isa LossTerm + term_id = id(inp_term) + term = loss(inp_term) + term_w = weight(inp_term) + else + "[id =>] AbstractLoss [=> weight] expected as a loss term, $(nameof(typeof(inp_term))) found" |> + ArgumentError |> + throw + end + + if term isa SemLoss + nsems += 1 + has_sem_weights |= !isnothing(term_w) + # check integrity + if isnothing(params) + params = SEM.params(term) + elseif params != SEM.params(term) + # FIXME the suggestion might no longer be relevant, since ParTable also stores params order + """ + The parameters of your SEM models do not match. + Maybe you tried to specify models of an ensemble via ParameterTables. + In that case, you may use RAMMatrices instead. + """ |> error + end + check_observed_vars(term) + elseif !(term isa AbstractLoss) + "AbstractLoss term expected at $(length(terms)+1) position, $(nameof(typeof(term))) found" |> + ArgumentError |> + throw + end + push!(terms, LossTerm(term, term_id, term_w)) + end + isnothing(params) && throw(ErrorException("No SEM models provided.")) + + if !has_sem_weights && nsems > 1 + # set the weights of SEMs in the ensemble + if default_sem_weights == :nsamples + # weight SEM by the number of samples + nsamples_total = sum(nsamples(term) for term in terms if issemloss(term)) + for (i, term) in enumerate(terms) + if issemloss(term) + terms[i] = + LossTerm(loss(term), id(term), nsamples(term) / nsamples_total) + end + end + elseif default_sem_weights == :uniform # uniform weights + for (i, term) in enumerate(terms) + if issemloss(term) + terms[i] = LossTerm(loss(term), id(term), 1 / nsems) + end + end + elseif default_sem_weights == :one # do nothing + end + end + + terms_tuple = Tuple(terms) + return Sem{typeof(terms_tuple)}(terms_tuple, params) +end + function Sem(; specification = ParameterTable, observed::O = SemObservedData, @@ -13,99 +173,127 @@ function Sem(; set_field_type_kwargs!(kwdict, observed, implied, loss, O, I) - observed, implied, loss = get_fields!(kwdict, specification, observed, implied, loss) - - sem = Sem(observed, implied, loss) + loss = get_fields!(kwdict, specification, observed, implied, loss) - return sem + return Sem(loss...) end +############################################################################################ +# functions +############################################################################################ + +params(model::AbstractSem) = model.params + """ - implied(model::AbstractSemSingle) -> SemImplied + loss_terms(model::AbstractSem) -Returns the [*implied*](@ref SemImplied) part of a model. +Returns a tuple of all [`LossTerm`](@ref) weighted terms in the SEM model. + +See also [`sem_terms`](@ref), [`loss_term`](@ref). """ -implied(model::AbstractSemSingle) = model.implied +loss_terms(model::AbstractSem) = model.loss_terms +nloss_terms(model::AbstractSem) = length(loss_terms(model)) -nvars(model::AbstractSemSingle) = nvars(implied(model)) -nobserved_vars(model::AbstractSemSingle) = nobserved_vars(implied(model)) -nlatent_vars(model::AbstractSemSingle) = nlatent_vars(implied(model)) +""" + sem_terms(model::AbstractSem) -vars(model::AbstractSemSingle) = vars(implied(model)) -observed_vars(model::AbstractSemSingle) = observed_vars(implied(model)) -latent_vars(model::AbstractSemSingle) = latent_vars(implied(model)) +Returns a tuple of all weighted SEM terms in the SEM model. -param_labels(model::AbstractSemSingle) = param_labels(implied(model)) -nparams(model::AbstractSemSingle) = nparams(implied(model)) +In comparison to [`loss_terms`](@ref) that returns all model terms, including e.g. +regularization terms, this function returns only the [`SemLoss`] terms. +See also [`loss_terms`](@ref), [`sem_term`](@ref). """ - observed(model::AbstractSemSingle) -> SemObserved +sem_terms(model::AbstractSem) = Tuple(term for term in loss_terms(model) if issemloss(term)) +nsem_terms(model::AbstractSem) = sum(issemloss, loss_terms(model)) + +nsamples(model::AbstractSem) = + sum(term -> issemloss(term) ? nsamples(term) : 0, loss_terms(model)) -Returns the [*observed*](@ref SemObserved) part of a model. """ -observed(model::AbstractSemSingle) = model.observed + loss_term(model::AbstractSem, id::Any) -> AbstractLoss -nsamples(model::AbstractSemSingle) = nsamples(observed(model)) +Returns the loss term with the specified `id` from the `model`. +Throws an error if the model has no term with the specified `id`. +See also [`loss_terms`](@ref). """ - loss(model::AbstractSemSingle) -> SemLoss +function loss_term(model::AbstractSem, id::Any) + for term in loss_terms(model) + if SEM.id(term) == id + return loss(term) + end + end + error("No loss term with id=$id found") +end -Returns the [*loss*](@ref SemLoss) function of a model. """ -loss(model::AbstractSemSingle) = model.loss - -# sum of samples in all sub-models -nsamples(ensemble::SemEnsemble) = sum(nsamples, ensemble.sems) - -function SemFiniteDiff(; - specification = ParameterTable, - observed::O = SemObservedData, - implied::I = RAM, - loss::L = SemML, - kwargs..., -) where {O, I, L} - kwdict = Dict{Symbol, Any}(kwargs...) + sem_term(model::AbstractSem, [id]) -> SemLoss - set_field_type_kwargs!(kwdict, observed, implied, loss, O, I) +Returns the SEM loss term with the specified `id` from the `model`. +Throws an error if the model has no term with the specified `id` or +if it is not of a [`SemLoss`](@ref) type. - observed, implied, loss = get_fields!(kwdict, specification, observed, implied, loss) +If no `id` is specified and the model contains only one SEM term, the term is returned. +Throws an error if the model contains multiple SEM terms. - sem = SemFiniteDiff(observed, implied, loss) +See also [`loss_term`](@ref), [`sem_terms`](@ref). +""" +function sem_term(model::AbstractSem, id::Any) + term = loss_term(model, id) + issemloss(term) || error("Loss term with id=$id ($(typeof(term))) is not a SEM term") + return term +end - return sem +function sem_term(model::AbstractSem, _::Nothing = nothing) + if nsem_terms(model) != 1 + error( + "Model contains $(nsem_terms(model)) SEM terms, you have to specify a specific term", + ) + end + for term in loss_terms(model) + issemloss(term) && return loss(term) + end + error("Unreachable reached") end -############################################################################################ -# functions -############################################################################################ +# wrappers arounds a single SemLoss term +observed(model::AbstractSem, id::Nothing = nothing) = observed(sem_term(model, id)) +implied(model::AbstractSem, id::Nothing = nothing) = implied(sem_term(model, id)) +vars(model::AbstractSem, id::Nothing = nothing) = vars(implied(model, id)) +observed_vars(model::AbstractSem, id::Nothing = nothing) = observed_vars(implied(model, id)) +latent_vars(model::AbstractSem, id::Nothing = nothing) = latent_vars(implied(model, id)) function set_field_type_kwargs!(kwargs, observed, implied, loss, O, I) kwargs[:observed_type] = O <: Type ? observed : typeof(observed) kwargs[:implied_type] = I <: Type ? implied : typeof(implied) if loss isa SemLoss - kwargs[:loss_types] = [ - lossfun isa SemLossFunction ? typeof(lossfun) : lossfun for - lossfun in loss.functions - ] - elseif applicable(iterate, loss) kwargs[:loss_types] = - [lossfun isa SemLossFunction ? typeof(lossfun) : lossfun for lossfun in loss] + [aloss isa SemLoss ? typeof(aloss) : aloss for aloss in loss.functions] + elseif applicable(iterate, loss) + kwargs[:loss_types] = [aloss isa SemLoss ? typeof(aloss) : aloss for aloss in loss] else - kwargs[:loss_types] = [loss isa SemLossFunction ? typeof(loss) : loss] + kwargs[:loss_types] = [loss isa SemLoss ? typeof(loss) : loss] end end # construct Sem fields -function get_fields!(kwargs, specification, observed, implied, loss) - if !isa(specification, SemSpecification) - specification = specification(; kwargs...) +function get_fields!(kwargs, spec, observed, implied, loss) + if !isa(spec, SemSpecification) + spec = spec(; kwargs...) end # observed if !isa(observed, SemObserved) - observed = observed(; specification, kwargs...) + observed = if spec isa EnsembleParameterTable + Dict( + term_id => observed(; specification = term_spec, kwargs...) for + (term_id, term_spec) in pairs(spec.tables) + ) + else + observed(; specification = spec, kwargs...) + end end - kwargs[:observed] = observed # implied if !isa(implied, SemImplied) @@ -117,95 +305,98 @@ function get_fields!(kwargs, specification, observed, implied, loss) else implied_kwargs = kwargs end - implied = implied(specification; implied_kwargs...) + implied = if spec isa EnsembleParameterTable + Dict( + term_id => implied(term_spec; implied_kwargs...) for + (term_id, term_spec) in pairs(spec.tables) + ) + else + implied(spec; implied_kwargs...) + end end - kwargs[:implied] = implied - kwargs[:nparams] = nparams(implied) - # loss - loss = get_SemLoss(loss; specification, kwargs...) - kwargs[:loss] = loss + loss_kwargs = copy(kwargs) + loss_kwargs[:nparams] = nparams(spec) + loss = build_SemTerms(loss, observed, implied; loss_kwargs...) - return observed, implied, loss + return loss end # construct loss field -function get_SemLoss(loss; kwargs...) +function build_SemTerms(loss, observed, implied; kwargs...) + function build_SemLoss(aloss, observed, implied) + if loss isa AbstractLoss + return loss + elseif aloss <: SemLoss{O, I} where {O, I} + return aloss(observed, implied; kwargs...) + else + return aloss(; kwargs...) + end + end + if loss isa SemLoss - nothing + return loss elseif applicable(iterate, loss) - loss_out = [] - for lossfun in loss - if isa(lossfun, SemLossFunction) - push!(loss_out, lossfun) - else - lossfun = lossfun(; kwargs...) - push!(loss_out, lossfun) - end - end - loss = SemLoss(loss_out...; kwargs...) + return [build_SemLoss(aloss, observed, implied) for aloss in loss] else - if !isa(loss, SemLossFunction) - loss = SemLoss(loss(; kwargs...); kwargs...) + if isa(observed, AbstractDict) && isa(implied, AbstractDict) + observed_ids = Set(keys(observed)) + implied_ids = Set(keys(implied)) + if observed_ids != implied_ids + """" + The term ids of the observed and the implied data do not match. + Observed term ids: $(observed_ids), implied term ids: $(implied_ids) + """ |> + ArgumentError |> + throw + end + loss_out = [ + begin + term_implied = implied[term_id] + if observed_vars(term_observed) != observed_vars(term_implied) + "observed_vars differ between the observed and the implied for the term $term_id" |> + ArgumentError |> + throw + end + LossTerm( + build_SemLoss(loss, term_observed, term_implied), + term_id, + nothing, + ) + end for (term_id, term_observed) in pairs(observed) + ] + return loss_out else - loss = SemLoss(loss; kwargs...) + if observed_vars(observed) != observed_vars(implied) + "observed_vars differ between the observed and the implied" |> + ArgumentError |> + throw + end + return (build_SemLoss(loss, observed, implied),) end end - return loss +end + +function update_observed(sem::Sem, new_observed; kwargs...) + new_terms = Tuple( + update_observed(lossterm.loss, new_observed; kwargs...) for + lossterm in loss_terms(sem) + ) + return Sem(new_terms...) end ############################################################## # pretty printing ############################################################## -function Base.show(io::IO, sem::Sem{O, I, L}) where {O, I, L} - lossfuntypes = @. string(nameof(typeof(sem.loss.functions))) - lossfuntypes = " " .* lossfuntypes .* ("\n") - print(io, "Structural Equation Model \n") - print(io, "- Loss Functions \n") - print(io, lossfuntypes...) - print(io, "- Fields \n") - print(io, " observed: $(nameof(O)) \n") - print(io, " implied: $(nameof(I)) \n") -end - -function Base.show(io::IO, sem::SemFiniteDiff{O, I, L}) where {O, I, L} - lossfuntypes = @. string(nameof(typeof(sem.loss.functions))) - lossfuntypes = " " .* lossfuntypes .* ("\n") - print(io, "Structural Equation Model : Finite Diff Approximation\n") - print(io, "- Loss Functions \n") - print(io, lossfuntypes...) - print(io, "- Fields \n") - print(io, " observed: $(nameof(O)) \n") - print(io, " implied: $(nameof(I)) \n") -end - -function Base.show(io::IO, loss::SemLoss) - lossfuntypes = @. string(nameof(typeof(loss.functions))) - lossfuntypes = " " .* lossfuntypes .* ("\n") - print(io, "SemLoss \n") - print(io, "- Loss Functions \n") - print(io, lossfuntypes...) - print(io, "- Weights \n") - for weight in loss.weights - if isnothing(weight.w) - print(io, " one \n") - else - print(io, "$(round.(weight.w, digits = 2)) \n") - end - end -end - -function Base.show(io::IO, models::SemEnsemble) - print(io, "SemEnsemble \n") - print(io, "- Number of Models: $(models.n) \n") - print(io, "- Weights: $(round.(models.weights, digits = 2)) \n") - - print(io, "\n", "Models: \n") - print(io, "===============================================", "\n") - for (model, i) in zip(models.sems, 1:models.n) - print(io, "---------------------- ", i, " ----------------------", "\n") - print(io, model) +function Base.show(io::IO, sem::AbstractSem) + println(io, "Structural Equation Model ($(nameof(typeof(sem))))") + println(io, "- $(nparams(sem)) parameters") + println(io, "- Loss terms:") + for term in loss_terms(sem) + print(io, " - ") + print(io, term) + println(io) end end diff --git a/src/implied/RAM/generic.jl b/src/implied/RAM/generic.jl index d57500a3..f1c1e08d 100644 --- a/src/implied/RAM/generic.jl +++ b/src/implied/RAM/generic.jl @@ -154,16 +154,11 @@ end ### methods ############################################################################################ -function update!( - targets::EvaluationTargets, - implied::RAM, - model::AbstractSemSingle, - param_labels, -) - materialize!(implied.A, implied.ram_matrices.A, param_labels) - materialize!(implied.S, implied.ram_matrices.S, param_labels) +function update!(targets::EvaluationTargets, implied::RAM, params) + materialize!(implied.A, implied.ram_matrices.A, params) + materialize!(implied.S, implied.ram_matrices.S, params) if !isnothing(implied.M) - materialize!(implied.M, implied.ram_matrices.M, param_labels) + materialize!(implied.M, implied.ram_matrices.M, params) end parent(implied.I_A) .= .-implied.A diff --git a/src/implied/RAM/symbolic.jl b/src/implied/RAM/symbolic.jl index 0f7868da..4c9bda91 100644 --- a/src/implied/RAM/symbolic.jl +++ b/src/implied/RAM/symbolic.jl @@ -176,12 +176,7 @@ end ### objective, gradient, hessian ############################################################################################ -function update!( - targets::EvaluationTargets, - implied::RAMSymbolic, - model::AbstractSemSingle, - par, -) +function update!(targets::EvaluationTargets, implied::RAMSymbolic, par) implied.Σ_eval!(implied.Σ, par) if MeanStruct(implied) === HasMeanStruct implied.μ_eval!(implied.μ, par) @@ -223,7 +218,10 @@ end function eval_Σ_symbolic(S, I_A⁻¹, F; vech::Bool = false, simplify::Bool = false) Σ = F * I_A⁻¹ * S * permutedims(I_A⁻¹) * permutedims(F) Σ = Array(Σ) - vech && (Σ = SEM.vech(Σ)) + if vech + n = size(Σ, 1) + Σ = [Σ[i, j] for j in 1:n for i in j:n] + end if simplify Threads.@threads for i in eachindex(Σ) Σ[i] = Symbolics.simplify(Σ[i]) diff --git a/src/implied/empty.jl b/src/implied/empty.jl index 82a6c946..a327ee13 100644 --- a/src/implied/empty.jl +++ b/src/implied/empty.jl @@ -13,8 +13,8 @@ Empty placeholder for models that don't need an implied part. - `specification`: either a `RAMMatrices` or `ParameterTable` object # Examples -A multigroup model with ridge regularization could be specified as a `SemEnsemble` with one -model per group and an additional model with `ImpliedEmpty` and `SemRidge` for the regularization part. +A multigroup model with ridge regularization could be specified as a `Sem` with one +SEM term (`SemLoss`) per group and an additional `SemRidge` regularization term. # Extended help @@ -45,7 +45,7 @@ end ### methods ############################################################################################ -update!(targets::EvaluationTargets, implied::ImpliedEmpty, par, model) = nothing +update!(targets::EvaluationTargets, implied::ImpliedEmpty, par) = nothing ############################################################################################ ### Recommended methods diff --git a/src/loss/ML/FIML.jl b/src/loss/ML/FIML.jl index 8572b15a..fdedf398 100644 --- a/src/loss/ML/FIML.jl +++ b/src/loss/ML/FIML.jl @@ -75,23 +75,27 @@ Can handle observed data with missing values. # Constructor - SemFIML(; observed::SemObservedMissing, implied::SemImplied, kwargs...) + SemFIML(observed::SemObservedMissing, implied::SemImplied) # Arguments - `observed::SemObservedMissing`: the observed part of the model (see [`SemObservedMissing`](@ref)) - `implied::SemImplied`: the implied part of the model + (see [`SemImplied`](@ref)) # Examples ```julia -my_fiml = SemFIML(observed = my_observed, implied = my_implied) +my_fiml = SemFIML(my_observed, my_implied) ``` # Interfaces Analytic gradients are available. """ -struct SemFIML{T, W} <: SemLossFunction +struct SemFIML{O, I, T, W} <: SemLoss{O, I} hessianeval::ExactHessian + + observed::O + implied::I patterns::Vector{SemFIMLPattern{T}} imp_inv::Matrix{T} # implied inverse @@ -105,7 +109,7 @@ end ### Constructors ############################################################################################ -function SemFIML(; observed::SemObservedMissing, implied, specification, kwargs...) +function SemFIML(observed::SemObservedMissing, implied::SemImplied; kwargs...) if MeanStruct(implied) === NoMeanStruct """ Full information maximum likelihood (FIML) can only be used with a meanstructure. @@ -117,6 +121,8 @@ function SemFIML(; observed::SemObservedMissing, implied, specification, kwargs. return SemFIML( ExactHessian(), + observed, + implied, [SemFIMLPattern(pat) for pat in observed.patterns], zeros(nobserved_vars(observed), nobserved_vars(observed)), CommutationMatrix(nvars(implied)), @@ -128,30 +134,28 @@ end ### methods ############################################################################################ -function evaluate!( - objective, - gradient, - hessian, - loss::SemFIML, - implied::SemImplied, - model::AbstractSemSingle, - params, -) +function evaluate!(objective, gradient, hessian, loss::SemFIML, params) isnothing(hessian) || error("Hessian not implemented for FIML") - if !check(loss, model) + implied = SEM.implied(loss) + observed = SEM.observed(loss) + + copyto!(loss.imp_inv, implied.Σ) + Σ_chol = cholesky!(Symmetric(loss.imp_inv); check = false) + + if !isposdef(Σ_chol) isnothing(objective) || (objective = non_posdef_return(params)) isnothing(gradient) || fill!(gradient, 1) return objective end - prepare!(loss, model) + @inbounds for (patloss, pat) in zip(loss.patterns, observed.patterns) + prepare!(patloss, pat, implied) + end - scale = inv(nsamples(observed(model))) - isnothing(objective) || - (objective = scale * F_FIML(eltype(params), loss, observed(model), model)) - isnothing(gradient) || - (∇F_FIML!(gradient, loss, observed(model), model); gradient .*= scale) + scale = inv(nsamples(observed)) + isnothing(objective) || (objective = scale * F_FIML(eltype(params), loss)) + isnothing(gradient) || (∇F_FIML!(gradient, loss); gradient .*= scale) return objective end @@ -167,27 +171,14 @@ update_observed(loss::SemFIML, observed::SemObserved; kwargs...) = ### additional functions ############################################################################################ -function prepare!(loss::SemFIML, observed::SemObservedMissing, implied::SemImplied) - @inbounds for (patloss, pat) in zip(loss.patterns, observed.patterns) - prepare!(patloss, pat, implied.Σ, implied.μ) - end -end - -prepare!(loss::SemFIML, model::AbstractSemSingle) = - prepare!(loss, observed(model), implied(model)) - -function check(loss::SemFIML, model::AbstractSemSingle) - copyto!(loss.imp_inv, implied(model).Σ) - a = cholesky!(Symmetric(loss.imp_inv); check = false) - return isposdef(a) +function ∇F_fiml_outer!(G, JΣ, Jμ, loss::SemFIML{O, I}) where {O, I <: SemImpliedSymbolic} + mul!(G, loss.implied.∇Σ', JΣ) # should be transposed + mul!(G, loss.implied.∇μ', Jμ, -1, 1) end -function ∇F_fiml_outer!(G, JΣ, Jμ, loss::SemFIML, implied::SemImpliedSymbolic, model) - mul!(G, implied.∇Σ', JΣ) # should be transposed - mul!(G, implied.∇μ', Jμ, -1, 1) -end +function ∇F_fiml_outer!(G, JΣ, Jμ, loss::SemFIML) + implied = loss.implied -function ∇F_fiml_outer!(G, JΣ, Jμ, loss::SemFIML, implied, model) Iₙ = sparse(1.0I, size(implied.A)...) P = kron(implied.F⨉I_A⁻¹, implied.F⨉I_A⁻¹) Q = kron(implied.S * implied.I_A⁻¹', Iₙ) @@ -203,25 +194,20 @@ function ∇F_fiml_outer!(G, JΣ, Jμ, loss::SemFIML, implied, model) mul!(G, ∇μ', Jμ, -1, 1) end -function F_FIML( - ::Type{T}, - loss::SemFIML, - observed::SemObservedMissing, - model::AbstractSemSingle, -) where {T} +function F_FIML(::Type{T}, loss::SemFIML) where {T} F = zero(T) - for (patloss, pat) in zip(loss.patterns, observed.patterns) + for (patloss, pat) in zip(loss.patterns, loss.observed.patterns) F += objective(patloss, pat) end return F end -function ∇F_FIML!(G, loss::SemFIML, observed::SemObservedMissing, model::AbstractSemSingle) - Jμ = zeros(nobserved_vars(model)) - JΣ = zeros(nobserved_vars(model)^2) +function ∇F_FIML!(G, loss::SemFIML) + Jμ = zeros(nobserved_vars(loss)) + JΣ = zeros(nobserved_vars(loss)^2) - for (patloss, pat) in zip(loss.patterns, observed.patterns) + for (patloss, pat) in zip(loss.patterns, loss.observed.patterns) gradient!(JΣ, Jμ, patloss, pat) end - ∇F_fiml_outer!(G, JΣ, Jμ, loss, implied(model), model) + ∇F_fiml_outer!(G, JΣ, Jμ, loss) end diff --git a/src/loss/ML/ML.jl b/src/loss/ML/ML.jl index aae1dada..2d449d73 100644 --- a/src/loss/ML/ML.jl +++ b/src/loss/ML/ML.jl @@ -8,36 +8,41 @@ Maximum likelihood estimation. # Constructor - SemML(; observed, approximate_hessian = false, kwargs...) + SemML(observed, implied; approximate_hessian = false) # Arguments - `observed::SemObserved`: the observed part of the model +- `implied::SemImplied`: [`SemImplied`](@ref) instance - `approximate_hessian::Bool`: if hessian-based optimization is used, should the hessian be swapped for an approximation # Examples ```julia -my_ml = SemML(observed = my_observed) +my_ml = SemML(my_observed, my_implied) ``` # Interfaces Analytic gradients are available, and for models without a meanstructure -and RAMSymbolic implied type, also analytic hessians. +and `RAMSymbolic` implied type, also analytic hessians. """ -struct SemML{HE <: HessianEval, INV, M, M2} <: SemLossFunction +struct SemML{O, I, HE <: HessianEval, INV, M, M2} <: SemLoss{O, I} + observed::O + implied::I hessianeval::HE Σ⁻¹::INV Σ⁻¹Σₒ::M meandiff::M2 - - SemML{HE}(args...) where {HE <: HessianEval} = - new{HE, map(typeof, args)...}(HE(), args...) end ############################################################################################ ### Constructors ############################################################################################ -function SemML(; observed::SemObserved, approximate_hessian::Bool = false, kwargs...) +function SemML( + observed::SemObserved, + implied::SemImplied; + approximate_hessian::Bool = false, + kwargs..., +) if observed isa SemObservedMissing @warn """ ML estimation with `SemObservedMissing` will use an approximate covariance and mean estimated with EM algorithm. @@ -51,12 +56,25 @@ function SemML(; observed::SemObserved, approximate_hessian::Bool = false, kwarg ) """ end + # check integrity + check_observed_vars(observed, implied) + he = approximate_hessian ? ApproxHessian() : ExactHessian() obsmean = obs_mean(observed) obscov = obs_cov(observed) meandiff = isnothing(obsmean) ? nothing : copy(obsmean) - return SemML{approximate_hessian ? ApproxHessian : ExactHessian}( + return SemML{ + typeof(observed), + typeof(implied), + typeof(he), + typeof(obscov), + typeof(obscov), + typeof(meandiff), + }( + observed, + implied, + he, similar(obscov), similar(obscov), meandiff, @@ -74,20 +92,20 @@ function evaluate!( objective, gradient, hessian, - semml::SemML, - implied::SemImpliedSymbolic, - model::AbstractSemSingle, + loss::SemML{<:Any, <:SemImpliedSymbolic}, par, ) + implied = SEM.implied(loss) + if !isnothing(hessian) (MeanStruct(implied) === HasMeanStruct) && throw(DomainError(H, "hessian of ML + meanstructure is not available")) end Σ = implied.Σ - Σₒ = obs_cov(observed(model)) - Σ⁻¹Σₒ = semml.Σ⁻¹Σₒ - Σ⁻¹ = semml.Σ⁻¹ + Σₒ = obs_cov(observed(loss)) + Σ⁻¹Σₒ = loss.Σ⁻¹Σₒ + Σ⁻¹ = loss.Σ⁻¹ copyto!(Σ⁻¹, Σ) Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false) @@ -105,7 +123,7 @@ function evaluate!( if MeanStruct(implied) === HasMeanStruct μ = implied.μ - μₒ = obs_mean(observed(model)) + μₒ = obs_mean(observed(loss)) μ₋ = μₒ - μ isnothing(objective) || (objective += dot(μ₋, Σ⁻¹, μ₋)) @@ -124,7 +142,7 @@ function evaluate!( mul!(gradient, ∇Σ', J') end if !isnothing(hessian) - if HessianEval(semml) === ApproxHessian + if HessianEval(loss) === ApproxHessian mul!(hessian, ∇Σ' * kron(Σ⁻¹, Σ⁻¹), ∇Σ, 2, 0) else ∇²Σ = implied.∇²Σ @@ -143,24 +161,17 @@ end ############################################################################################ ### Non-Symbolic Implied Types -function evaluate!( - objective, - gradient, - hessian, - semml::SemML, - implied::RAM, - model::AbstractSemSingle, - par, -) +function evaluate!(objective, gradient, hessian, loss::SemML, par) if !isnothing(hessian) error("hessian of ML + non-symbolic implied type is not available") end - Σ = implied.Σ - Σₒ = obs_cov(observed(model)) - Σ⁻¹Σₒ = semml.Σ⁻¹Σₒ - Σ⁻¹ = semml.Σ⁻¹ + implied = SEM.implied(loss) + Σ = implied.Σ + Σₒ = obs_cov(observed(loss)) + Σ⁻¹Σₒ = loss.Σ⁻¹Σₒ + Σ⁻¹ = loss.Σ⁻¹ copyto!(Σ⁻¹, Σ) Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false) if !isposdef(Σ_chol) @@ -179,7 +190,7 @@ function evaluate!( if MeanStruct(implied) === HasMeanStruct μ = implied.μ - μₒ = obs_mean(observed(model)) + μₒ = obs_mean(observed(loss)) μ₋ = μₒ - μ objective += dot(μ₋, Σ⁻¹, μ₋) end @@ -198,7 +209,7 @@ function evaluate!( if MeanStruct(implied) === HasMeanStruct μ = implied.μ - μₒ = obs_mean(observed(model)) + μₒ = obs_mean(observed(loss)) ∇M = implied.∇M M = implied.M μ₋ = μₒ - μ @@ -229,16 +240,17 @@ end ### recommended methods ############################################################################################ -update_observed(lossfun::SemML, observed::SemObservedMissing; kwargs...) = +update_observed(loss::SemML, observed::SemObservedMissing; kwargs...) = error("ML estimation does not work with missing data - use FIML instead") -function update_observed(lossfun::SemML, observed::SemObserved; kwargs...) - if size(lossfun.Σ⁻¹) == size(obs_cov(observed)) - return lossfun +function update_observed(loss::SemML, observed::SemObserved; kwargs...) + if (obs_cov(loss) == obs_cov(observed)) && (obs_mean(loss) == obs_mean(observed)) + return loss # no change else - return SemML(; - observed = observed, - approximate_hessian = HessianEval(lossfun) == ApproxHessian, + return SemML( + observed, + loss.implied; + approximate_hessian = HessianEval(loss) == ApproxHessian, kwargs..., ) end diff --git a/src/loss/WLS/WLS.jl b/src/loss/WLS/WLS.jl index 9de011f6..5c4cb252 100644 --- a/src/loss/WLS/WLS.jl +++ b/src/loss/WLS/WLS.jl @@ -9,8 +9,8 @@ At the moment only available with the `RAMSymbolic` implied type. # Constructor - SemWLS(; - observed, implied, + SemWLS( + observed::SemObserved, implied::SemImplied; wls_weight_matrix = nothing, wls_weight_matrix_mean = nothing, approximate_hessian = false, @@ -18,7 +18,7 @@ At the moment only available with the `RAMSymbolic` implied type. # Arguments - `observed`: the `SemObserved` part of the model -- `implied::SemImplied`: the implied part of the model +- `implied`: the `SemImplied` part of the model - `approximate_hessian::Bool`: should the hessian be swapped for an approximation - `wls_weight_matrix`: the weight matrix for weighted least squares. Defaults to GLS estimation (``0.5*(D^T*kron(S,S)*D)`` where D is the duplication matrix @@ -28,29 +28,37 @@ At the moment only available with the `RAMSymbolic` implied type. # Examples ```julia -my_wls = SemWLS(observed = my_observed, implied = my_implied) +my_wls = SemWLS(my_observed, my_implied) ``` # Interfaces Analytic gradients are available, and for models without a meanstructure also analytic hessians. """ -struct SemWLS{HE <: HessianEval, Vt, St, C} <: SemLossFunction +struct SemWLS{O, I, HE <: HessianEval, Vt, St, C} <: SemLoss{O, I} + observed::O + implied::I + hessianeval::HE V::Vt σₒ::St V_μ::C + + SemWLS(observed, implied, ::Type{HE}, args...) where {HE <: HessianEval} = + new{typeof(observed), typeof(implied), HE, map(typeof, args)...}( + observed, + implied, + HE(), + args..., + ) end ############################################################################################ ### Constructors ############################################################################################ -SemWLS{HE}(args...) where {HE <: HessianEval} = - SemWLS{HE, map(typeof, args)...}(HE(), args...) - -function SemWLS(; +function SemWLS( observed::SemObserved, - implied::SemImplied, + implied::SemImplied; wls_weight_matrix::Union{AbstractMatrix, Nothing} = nothing, wls_weight_matrix_mean::Union{AbstractMatrix, Nothing} = nothing, approximate_hessian::Bool = false, @@ -75,14 +83,19 @@ function SemWLS(; ArgumentError |> throw end + # check integrity + check_observed_vars(observed, implied) nobs_vars = nobserved_vars(observed) tril_ind = filter(x -> (x[1] >= x[2]), CartesianIndices(obs_cov(observed))) s = obs_cov(observed)[tril_ind] - size(s) == size(implied.Σ) || - throw(DimensionMismatch("SemWLS requires implied covariance to be in vech-ed form " * - "(vectorized lower triangular part of Σ matrix): $(size(s)) expected, $(size(implied.Σ)) found.\n" * - "$(nameof(typeof(implied))) must be constructed with vech=true.")) + size(s) == size(implied.Σ) || throw( + DimensionMismatch( + "SemWLS requires implied covariance to be in vech-ed form " * + "(vectorized lower triangular part of Σ matrix): $(size(s)) expected, $(size(implied.Σ)) found.\n" * + "$(nameof(typeof(implied))) must be constructed with vech=true.", + ), + ) # compute V here if isnothing(wls_weight_matrix) @@ -101,13 +114,12 @@ function SemWLS(; if MeanStruct(implied) == HasMeanStruct if isnothing(wls_weight_matrix_mean) - @warn "Computing WLS weight matrix for the meanstructure using obs_cov()" + @info "Computing WLS weight matrix for the meanstructure using obs_cov()" wls_weight_matrix_mean = inv(obs_cov(observed)) - else - size(wls_weight_matrix_mean) == (nobs_vars, nobs_vars) || DimensionMismatch( - "wls_weight_matrix_mean has to be of size $(nobs_vars)×$(nobs_vars)", - ) end + size(wls_weight_matrix_mean) == (nobs_vars, nobs_vars) || DimensionMismatch( + "wls_weight_matrix_mean has to be of size $(nobs_vars)×$(nobs_vars)", + ) else isnothing(wls_weight_matrix_mean) || @warn "Ignoring wls_weight_matrix_mean since meanstructure is disabled" @@ -115,31 +127,25 @@ function SemWLS(; end HE = approximate_hessian ? ApproxHessian : ExactHessian - return SemWLS{HE}(wls_weight_matrix, s, wls_weight_matrix_mean) + return SemWLS(observed, implied, HE, wls_weight_matrix, s, wls_weight_matrix_mean) end ############################################################################ ### methods ############################################################################ -function evaluate!( - objective, - gradient, - hessian, - semwls::SemWLS, - implied::SemImpliedSymbolic, - model::AbstractSemSingle, - par, -) +function evaluate!(objective, gradient, hessian, loss::SemWLS, par) + implied = SEM.implied(loss) + if !isnothing(hessian) && (MeanStruct(implied) === HasMeanStruct) error("hessian of WLS with meanstructure is not available") end - V = semwls.V + V = loss.V ∇σ = implied.∇Σ σ = implied.Σ - σₒ = semwls.σₒ + σₒ = loss.σₒ σ₋ = σₒ - σ isnothing(objective) || (objective = dot(σ₋, V, σ₋)) @@ -152,17 +158,17 @@ function evaluate!( gradient .*= -2 end isnothing(hessian) || (mul!(hessian, ∇σ' * V, ∇σ, 2, 0)) - if !isnothing(hessian) && (HessianEval(semwls) === ExactHessian) + if !isnothing(hessian) && (HessianEval(loss) === ExactHessian) ∇²Σ = implied.∇²Σ - J = -2 * (σ₋' * semwls.V)' + J = -2 * (σ₋' * loss.V)' implied.∇²Σ_eval!(∇²Σ, J, par) hessian .+= ∇²Σ end if MeanStruct(implied) === HasMeanStruct μ = implied.μ - μₒ = obs_mean(observed(model)) + μₒ = obs_mean(observed(loss)) μ₋ = μₒ - μ - V_μ = semwls.V_μ + V_μ = loss.V_μ if !isnothing(objective) objective += dot(μ₋, V_μ, μ₋) end @@ -179,23 +185,19 @@ end ############################################################################################ function update_observed( - lossfun::SemWLS, + loss::SemWLS, observed::SemObserved; recompute_V = true, kwargs..., ) if recompute_V - return SemWLS(; - observed = observed, - meanstructure = MeanStruct(kwargs[:implied]) == HasMeanStruct, - kwargs..., - ) + return SemWLS(observed, loss.implied; kwargs...) else - return SemWLS(; - observed = observed, - wls_weight_matrix = lossfun.V, - wls_weight_matrix_mean = lossfun.V_μ, - meanstructure = MeanStruct(kwargs[:implied]) == HasMeanStruct, + return SemWLS( + observed, + loss.implied; + wls_weight_matrix = loss.V, + wls_weight_matrix_mean = loss.V_μ, kwargs..., ) end diff --git a/src/loss/abstract.jl b/src/loss/abstract.jl new file mode 100644 index 00000000..bf8585d6 --- /dev/null +++ b/src/loss/abstract.jl @@ -0,0 +1,42 @@ +""" + observed(loss::SemLoss) -> SemObserved + +Returns the [*observed*](@ref SemObserved) part of a model. +""" +observed(loss::SemLoss) = loss.observed + +""" + implied(loss::SemLoss) -> SemImplied + +Returns the [*implied*](@ref SemImplied) part of a model. +""" +implied(loss::SemLoss) = loss.implied + +for f in (:nsamples, :obs_cov, :obs_mean) + @eval $f(loss::SemLoss) = $f(observed(loss)) +end + +for f in ( + :vars, + :nvars, + :latent_vars, + :nlatent_vars, + :observed_vars, + :nobserved_vars, + :params, + :nparams, +) + @eval $f(loss::SemLoss) = $f(implied(loss)) +end + +function check_observed_vars(observed::SemObserved, implied::SemImplied) + isnothing(observed_vars(implied)) || + observed_vars(observed) == observed_vars(implied) || + throw( + ArgumentError( + "Observed variables defined for \"observed\" and \"implied\" do not match.", + ), + ) +end + +check_observed_vars(sem::SemLoss) = check_observed_vars(observed(sem), implied(sem)) diff --git a/src/loss/constant/constant.jl b/src/loss/constant/constant.jl index 3aed5e27..023076cc 100644 --- a/src/loss/constant/constant.jl +++ b/src/loss/constant/constant.jl @@ -4,6 +4,8 @@ ### Types ############################################################################################ """ + SemConstant{C <: Number} <: AbstractLoss + Constant loss term. Can be used for comparability to other packages. # Constructor @@ -15,37 +17,27 @@ Constant loss term. Can be used for comparability to other packages. # Examples ```julia - my_constant = SemConstant(constant_loss = 42.0) + my_constant = SemConstant(42.0) ``` # Interfaces Analytic gradients and hessians are available. """ -struct SemConstant{C} <: SemLossFunction +struct SemConstant{C <: Number} <: AbstractLoss hessianeval::ExactHessian c::C -end -############################################################################################ -### Constructors -############################################################################################ - -function SemConstant(; constant_loss, kwargs...) - return SemConstant(ExactHessian(), constant_loss) + SemConstant(c::Number) = new{typeof(c)}(ExactHessian(), c) end -############################################################################################ -### methods -############################################################################################ +SemConstant(; constant_loss::Number, kwargs...) = SemConstant(constant_loss) -objective(constant::SemConstant, model::AbstractSem, par) = constant.c -gradient(constant::SemConstant, model::AbstractSem, par) = zero(par) -hessian(constant::SemConstant, model::AbstractSem, par) = - zeros(eltype(par), length(par), length(par)) +objective(loss::SemConstant, par) = convert(eltype(par), loss.c) +gradient(loss::SemConstant, par) = zero(par) +hessian(loss::SemConstant, par) = zeros(eltype(par), length(par), length(par)) ############################################################################################ ### Recommended methods ############################################################################################ -update_observed(loss_function::SemConstant, observed::SemObserved; kwargs...) = - loss_function +update_observed(loss::SemConstant, observed::SemObserved; kwargs...) = loss diff --git a/src/loss/regularization/ridge.jl b/src/loss/regularization/ridge.jl index 90cbcc23..3e2cfbff 100644 --- a/src/loss/regularization/ridge.jl +++ b/src/loss/regularization/ridge.jl @@ -25,7 +25,7 @@ my_ridge = SemRidge(;α_ridge = 0.02, which_ridge = [:λ₁, :λ₂, :ω₂₃], # Interfaces Analytic gradients and hessians are available. """ -struct SemRidge{P, W1, W2, GT, HT} <: SemLossFunction +struct SemRidge{P, W1, W2, GT, HT} <: AbstractLoss hessianeval::ExactHessian α::P which::W1 @@ -74,15 +74,14 @@ end ### methods ############################################################################################ -objective(ridge::SemRidge, model::AbstractSem, par) = - @views ridge.α * sum(abs2, par[ridge.which]) +objective(ridge::SemRidge, par) = @views ridge.α * sum(abs2, par[ridge.which]) -function gradient(ridge::SemRidge, model::AbstractSem, par) +function gradient(ridge::SemRidge, par) @views ridge.gradient[ridge.which] .= (2 * ridge.α) * par[ridge.which] return ridge.gradient end -function hessian(ridge::SemRidge, model::AbstractSem, par) +function hessian(ridge::SemRidge, par) @views @. ridge.hessian[ridge.which_H] .= 2 * ridge.α return ridge.hessian end diff --git a/src/objective_gradient_hessian.jl b/src/objective_gradient_hessian.jl index 69915ffa..23cef4e6 100644 --- a/src/objective_gradient_hessian.jl +++ b/src/objective_gradient_hessian.jl @@ -24,68 +24,61 @@ is_hessian_required(::EvaluationTargets{<:Any, <:Any, H}) where {H} = H (targets::EvaluationTargets)(arg_tuple::Tuple) = targets(arg_tuple...) """ - evaluate!(objective, gradient, hessian [, lossfun], model, params) + evaluate!(objective, gradient, hessian, loss::AbstractLoss, params) + evaluate!(objective, gradient, hessian, model::AbstractSem, params) Evaluates the objective, gradient, and/or Hessian at the given parameter vector. -If a loss function is passed, only this specific loss function is evaluated, otherwise, -the sum of all loss functions in the model is evaluated. + +If a single loss term (`loss`) is passed, only this specific term is evaluated, +otherwise, if the entire SEM `model` is passed, the weighted sum of all loss terms +in the model is evaluated. If objective, gradient or hessian are `nothing`, they are not evaluated. For example, since many numerical optimization algorithms don't require a Hessian, -the computation will be turned off by setting `hessian` to `nothing`. +its computation will be turned off by setting `hessian` to `nothing`. + +During the evaluation, the internal state of the loss term or of the model +could be modified. # Arguments - `objective`: a Number if the objective should be evaluated, otherwise `nothing` - `gradient`: a pre-allocated vector the gradient should be written to, otherwise `nothing` - `hessian`: a pre-allocated matrix the Hessian should be written to, otherwise `nothing` -- `lossfun::SemLossFunction`: loss function to evaluate +- `loss::AbstractLoss`: loss function to evaluate - `model::AbstractSem`: model to evaluate - `params`: vector of parameters # Implementing a new loss function -To implement a new loss function, a new method for `evaluate!` has to be defined. +To implement a new loss (subtype of `SemLoss` for SEM terms, or of `AbstractLoss` for +regularization terms), a new method for `evaluate!` has to be defined. This is explained in the online documentation on [Custom loss functions](@ref). """ function evaluate! end -# dispatch on SemImplied -evaluate!(objective, gradient, hessian, loss::SemLossFunction, model::AbstractSem, params) = - evaluate!(objective, gradient, hessian, loss, implied(model), model, params) - # fallback method -function evaluate!( - obj, - grad, - hess, - loss::SemLossFunction, - implied::SemImplied, - model, - params, -) - isnothing(obj) || (obj = objective(loss, implied, model, params)) - isnothing(grad) || copyto!(grad, gradient(loss, implied, model, params)) - isnothing(hess) || copyto!(hess, hessian(loss, implied, model, params)) +function evaluate!(obj, grad, hess, loss::AbstractLoss, params) + isnothing(obj) || (obj = objective(loss, params)) + isnothing(grad) || copyto!(grad, gradient(loss, params)) + isnothing(hess) || copyto!(hess, hessian(loss, params)) return obj end -# fallback methods -objective(f::SemLossFunction, implied::SemImplied, model, params) = - objective(f, model, params) -gradient(f::SemLossFunction, implied::SemImplied, model, params) = - gradient(f, model, params) -hessian(f::SemLossFunction, implied::SemImplied, model, params) = hessian(f, model, params) +evaluate!(obj, grad, hess, term::LossTerm, params) = + evaluate!(obj, grad, hess, loss(term), params) # fallback method for SemImplied that calls update_xxx!() methods -function update!(targets::EvaluationTargets, implied::SemImplied, model, params) - is_objective_required(targets) && update_objective!(implied, model, params) - is_gradient_required(targets) && update_gradient!(implied, model, params) - is_hessian_required(targets) && update_hessian!(implied, model, params) +function update!(targets::EvaluationTargets, implied::SemImplied, params) + is_objective_required(targets) && update_objective!(implied, params) + is_gradient_required(targets) && update_gradient!(implied, params) + is_hessian_required(targets) && update_hessian!(implied, params) end +const AbstractSemOrLoss = Union{AbstractSem, AbstractLoss} + # guess objective type -objective_type(model::AbstractSem, params::Any) = Float64 -objective_type(model::AbstractSem, params::AbstractVector{T}) where {T <: Number} = T -objective_zero(model::AbstractSem, params::Any) = zero(objective_type(model, params)) +objective_type(model::AbstractSemOrLoss, params::Any) = Float64 +objective_type(model::AbstractSemOrLoss, params::AbstractVector{T}) where {T <: Number} = T +objective_zero(model::AbstractSemOrLoss, params::Any) = zero(objective_type(model, params)) objective_type(objective::T, gradient, hessian) where {T <: Number} = T objective_type( @@ -101,145 +94,151 @@ objective_type( objective_zero(objective, gradient, hessian) = zero(objective_type(objective, gradient, hessian)) +evaluate!(objective, gradient, hessian, model::AbstractSem, params) = + error("evaluate!() for $(typeof(model)) is not implemented") + ############################################################################################ -# methods for AbstractSem +# methods for Sem ############################################################################################ -function evaluate!(objective, gradient, hessian, model::AbstractSemSingle, params) - targets = EvaluationTargets(objective, gradient, hessian) - # update implied state, its gradient and hessian (if required) - update!(targets, implied(model), model, params) - return evaluate!( - !isnothing(objective) ? zero(objective) : nothing, - gradient, - hessian, - loss(model), - model, - params, - ) -end +function evaluate!(objective, gradient, hessian, model::Sem, params) + # reset output + isnothing(objective) || (objective = objective_zero(objective, gradient, hessian)) + isnothing(gradient) || fill!(gradient, zero(eltype(gradient))) + isnothing(hessian) || fill!(hessian, zero(eltype(hessian))) -############################################################################################ -# methods for SemFiniteDiff -# (approximate gradient and hessian with finite differences of objective) -############################################################################################ + # gradient and hessian for individual terms + t_grad = isnothing(gradient) ? nothing : similar(gradient) + t_hess = isnothing(hessian) ? nothing : similar(hessian) + + # update implied states of all SemLoss terms before term calculation loop + # to make sure all terms use updated implied states + targets = EvaluationTargets(objective, gradient, hessian) + for term in loss_terms(model) + issemloss(term) && update!(targets, implied(term), params) + end -function evaluate!(objective, gradient, hessian, model::SemFiniteDiff, params) - function obj(p) - # recalculate implied state for p - update!(EvaluationTargets{true, false, false}(), implied(model), model, p) - evaluate!( - objective_zero(objective, gradient, hessian), - nothing, - nothing, - loss(model), - model, - p, + for term in loss_terms(model) + t_obj = evaluate!(objective, t_grad, t_hess, term, params) + #@show nameof(typeof(term)) t_obj + objective = accumulate_loss!( + objective, + gradient, + hessian, + weight(term), + t_obj, + t_grad, + t_hess, ) end - isnothing(gradient) || FiniteDiff.finite_difference_gradient!(gradient, obj, params) - isnothing(hessian) || FiniteDiff.finite_difference_hessian!(hessian, obj, params) - return !isnothing(objective) ? obj(params) : nothing + return objective end -objective(model::AbstractSem, params) = - evaluate!(objective_zero(model, params), nothing, nothing, model, params) - -############################################################################################ -# methods for SemLoss (weighted sum of individual SemLossFunctions) -############################################################################################ +# internal function to accumulate loss objective, gradient and hessian +function accumulate_loss!( + total_objective, + total_gradient, + total_hessian, + weight::Nothing, + objective, + gradient, + hessian, +) + isnothing(total_gradient) || (total_gradient .+= gradient) + isnothing(total_hessian) || (total_hessian .+= hessian) + return isnothing(total_objective) ? total_objective : (total_objective + objective) +end -function evaluate!(objective, gradient, hessian, loss::SemLoss, model::AbstractSem, params) - isnothing(objective) || (objective = zero(objective)) - isnothing(gradient) || fill!(gradient, zero(eltype(gradient))) - isnothing(hessian) || fill!(hessian, zero(eltype(hessian))) - f_grad = isnothing(gradient) ? nothing : similar(gradient) - f_hess = isnothing(hessian) ? nothing : similar(hessian) - for (f, weight) in zip(loss.functions, loss.weights) - f_obj = evaluate!(objective, f_grad, f_hess, f, model, params) - isnothing(objective) || (objective += weight * f_obj) - isnothing(gradient) || (gradient .+= weight * f_grad) - isnothing(hessian) || (hessian .+= weight * f_hess) - end - return objective +function accumulate_loss!( + total_objective, + total_gradient, + total_hessian, + weight::Number, + objective, + gradient, + hessian, +) + isnothing(total_gradient) || axpy!(weight, gradient, total_gradient) + isnothing(total_hessian) || axpy!(weight, hessian, total_hessian) + return isnothing(total_objective) ? total_objective : + (total_objective + weight * objective) end ############################################################################################ -# methods for SemEnsemble (weighted sum of individual AbstractSemSingle models) +# methods for SemFiniteDiff +# (approximate gradient and hessian with finite differences of objective) ############################################################################################ -function evaluate!(objective, gradient, hessian, ensemble::SemEnsemble, params) - isnothing(objective) || (objective = zero(objective)) - isnothing(gradient) || fill!(gradient, zero(eltype(gradient))) - isnothing(hessian) || fill!(hessian, zero(eltype(hessian))) - sem_grad = isnothing(gradient) ? nothing : similar(gradient) - sem_hess = isnothing(hessian) ? nothing : similar(hessian) - for (sem, weight) in zip(ensemble.sems, ensemble.weights) - sem_obj = evaluate!(objective, sem_grad, sem_hess, sem, params) - isnothing(objective) || (objective += weight * sem_obj) - isnothing(gradient) || (gradient .+= weight * sem_grad) - isnothing(hessian) || (hessian .+= weight * sem_hess) - end - return objective +# evaluate!() wrapper that does some housekeeping, if necessary +_evaluate!(args...) = evaluate!(args...) + +# update implied state, its gradient and hessian +function _evaluate!(objective, gradient, hessian, loss::SemLoss, params) + # note that any other Sem loss terms that are dependent on implied + # should be enumerated after the SemLoss term + # otherwise they would be using outdated implied state + update!(EvaluationTargets(objective, gradient, hessian), implied(loss), params) + return evaluate!(objective, gradient, hessian, loss, params) end +objective(model::AbstractSemOrLoss, params) = + _evaluate!(objective_zero(model, params), nothing, nothing, model, params) + +# throw an error by default if gradient! and hessian! are not implemented + +#= gradient!(model::AbstractSemOrLoss, par, model) = + throw(ArgumentError("gradient for $(nameof(typeof(model))) is not available")) + +hessian!(model::AbstractSemOrLoss, par, model) = + throw(ArgumentError("hessian for $(nameof(typeof(model))) is not available")) =# + ############################################################################################ # Documentation ############################################################################################ """ objective!(model::AbstractSem, params) -Returns the objective value at `params`. -The model object can be modified. +Calculates the objective value at `params`. -# Implementation -To implement a new `SemImplied` or `SemLossFunction` subtype, you need to add a method for - objective!(newtype::MyNewType, params, model::AbstractSemSingle) +The model object can be modified during calculation. -To implement a new `AbstractSem` subtype, you need to add a method for - objective!(model::MyNewType, params) +See also [`evaluate!`](@ref). """ function objective! end """ gradient!(gradient, model::AbstractSem, params) -Writes the gradient value at `params` to `gradient`. +Calculates the model's gradient at `params` and writes it to `gradient`. -# Implementation -To implement a new `SemImplied` or `SemLossFunction` type, you can add a method for - gradient!(newtype::MyNewType, params, model::AbstractSemSingle) +The model object can be modified during calculation. -To implement a new `AbstractSem` subtype, you can add a method for - gradient!(gradient, model::MyNewType, params) +See also [`evaluate!`](@ref). """ function gradient! end """ hessian!(hessian, model::AbstractSem, params) -Writes the hessian value at `params` to `hessian`. +Calculates the model's hessian at `params` and writes it to `hessian`. -# Implementation -To implement a new `SemImplied` or `SemLossFunction` type, you can add a method for - hessian!(newtype::MyNewType, params, model::AbstractSemSingle) +The model object can be modified during calculation. -To implement a new `AbstractSem` subtype, you can add a method for - hessian!(hessian, model::MyNewType, params) +See also [`evaluate!`](@ref). """ function hessian! end objective!(model::AbstractSem, params) = - evaluate!(objective_zero(model, params), nothing, nothing, model, params) + _evaluate!(objective_zero(model, params), nothing, nothing, model, params) gradient!(gradient, model::AbstractSem, params) = - evaluate!(nothing, gradient, nothing, model, params) + _evaluate!(nothing, gradient, nothing, model, params) hessian!(hessian, model::AbstractSem, params) = - evaluate!(nothing, nothing, hessian, model, params) + _evaluate!(nothing, nothing, hessian, model, params) objective_gradient!(gradient, model::AbstractSem, params) = - evaluate!(objective_zero(model, params), gradient, nothing, model, params) + _evaluate!(objective_zero(model, params), gradient, nothing, model, params) objective_hessian!(hessian, model::AbstractSem, params) = - evaluate!(objective_zero(model, params), nothing, hessian, model, params) + _evaluate!(objective_zero(model, params), nothing, hessian, model, params) gradient_hessian!(gradient, hessian, model::AbstractSem, params) = - evaluate!(nothing, gradient, hessian, model, params) + _evaluate!(nothing, gradient, hessian, model, params) objective_gradient_hessian!(gradient, hessian, model::AbstractSem, params) = - evaluate!(objective_zero(model, params), gradient, hessian, model, params) + _evaluate!(objective_zero(model, params), gradient, hessian, model, params) diff --git a/src/optimizer/abstract.jl b/src/optimizer/abstract.jl index 0c7913c4..6774e549 100644 --- a/src/optimizer/abstract.jl +++ b/src/optimizer/abstract.jl @@ -137,13 +137,13 @@ fit(model::AbstractSem; engine::Symbol = :Optim, start_val = nothing, kwargs...) fit(optim::SemOptimizer, model::AbstractSem, start_params; kwargs...) = error("Optimizer $(optim) support not implemented.") -# FABIN3 is the default method for single models -prepare_start_params(start_val::Nothing, model::AbstractSemSingle; kwargs...) = - start_fabin3(model; kwargs...) - -# simple algorithm is the default method for ensembles -prepare_start_params(start_val::Nothing, model::AbstractSem; kwargs...) = - start_simple(model; kwargs...) +# defaults when no starting parameters are specified +function prepare_start_params(start_val::Nothing, model::AbstractSem; kwargs...) + sems = sem_terms(model) + # FABIN3 for single models, simple algorithm for ensembles + return length(sems) == 1 ? start_fabin3(loss(sems[1]); kwargs...) : + start_simple(model; kwargs...) +end # first argument is a function prepare_start_params(start_val, model::AbstractSem; kwargs...) = start_val(model; kwargs...) diff --git a/src/types.jl b/src/types.jl index 3a6b5fdf..87b733cf 100644 --- a/src/types.jl +++ b/src/types.jl @@ -1,14 +1,6 @@ ############################################################################################ # Define the basic type system ############################################################################################ -"Most abstract supertype for all SEMs" -abstract type AbstractSem end - -"Supertype for all single SEMs, e.g. SEMs that have at least the fields `observed`, `implied`, `loss`" -abstract type AbstractSemSingle{O, I, L} <: AbstractSem end - -"Supertype for all collections of multiple SEMs" -abstract type AbstractSemCollection <: AbstractSem end "Meanstructure trait for `SemImplied` subtypes" abstract type MeanStruct end @@ -36,48 +28,8 @@ HessianEval(::Type{T}) where {T} = HessianEval(semobj) = HessianEval(typeof(semobj)) -"Supertype for all loss functions of SEMs. If you want to implement a custom loss function, it should be a subtype of `SemLossFunction`." -abstract type SemLossFunction end - -""" - SemLoss(args...; loss_weights = nothing, ...) - -Constructs the loss field of a SEM. Can contain multiple `SemLossFunction`s, the model is optimized over their sum. -See also [`SemLossFunction`](@ref). - -# Arguments -- `args...`: Multiple `SemLossFunction`s. -- `loss_weights::Vector`: Weights for each loss function. Defaults to unweighted optimization. - -# Examples -```julia -my_ml_loss = SemML(...) -my_ridge_loss = SemRidge(...) -my_loss = SemLoss(SemML, SemRidge; loss_weights = [1.0, 2.0]) -``` -""" -mutable struct SemLoss{F <: Tuple, T} - functions::F - weights::T -end - -function SemLoss(functions...; loss_weights = nothing, kwargs...) - if !isnothing(loss_weights) - loss_weights = SemWeight.(loss_weights) - else - loss_weights = Tuple(SemWeight(nothing) for _ in 1:length(functions)) - end - - return SemLoss(functions, loss_weights) -end - -# weights for loss functions or models. If the weight is nothing, multiplication returns the second argument -struct SemWeight{T} - w::T -end - -Base.:*(x::SemWeight{Nothing}, y) = y -Base.:*(x::SemWeight, y) = x.w * y +"Supertype for all loss functions of SEMs. If you want to implement a custom loss function, it should be a subtype of `AbstractLoss`." +abstract type AbstractLoss end abstract type SemOptimizer{E} end @@ -85,6 +37,8 @@ abstract type SemOptimizer{E} end abstract type SemOptimizerResult{O <: SemOptimizer} end """ + abstract type SemObserved + Supertype of all objects that can serve as the observed field of a SEM. Pre-processes data and computes sufficient statistics for example. If you have a special kind of data, e.g. ordinal data, you should implement a subtype of SemObserved. @@ -103,169 +57,90 @@ abstract type SemImplied end abstract type SemImpliedSymbolic <: SemImplied end """ - Sem(;observed = SemObservedData, implied = RAM, loss = SemML, kwargs...) + abstract type SemLoss{O <: SemObserved, I <: SemImplied} <: AbstractLoss -Constructor for the basic `Sem` type. -All additional kwargs are passed down to the constructors for the observed, implied, and loss fields. +The base type for calculating the loss of the implied SEM model when explaining the observed data. -# Arguments -- `observed`: object of subtype `SemObserved` or a constructor. -- `implied`: object of subtype `SemImplied` or a constructor. -- `loss`: object of subtype `SemLossFunction`s or constructor; or a tuple of such. - -Returns a Sem with fields -- `observed::SemObserved`: Stores observed data, sample statistics, etc. See also [`SemObserved`](@ref). -- `implied::SemImplied`: Computes model implied statistics, like Σ, μ, etc. See also [`SemImplied`](@ref). -- `loss::SemLoss`: Computes the objective and gradient of a sum of loss functions. See also [`SemLoss`](@ref). +All subtypes of `SemLoss` should have the following fields: +- `observed::O`: object of subtype [`SemObserved`](@ref). +- `implied::I`: object of subtype [`SemImplied`](@ref). """ -mutable struct Sem{O <: SemObserved, I <: SemImplied, L <: SemLoss} <: - AbstractSemSingle{O, I, L} - observed::O - implied::I - loss::L -end +abstract type SemLoss{O <: SemObserved, I <: SemImplied} <: AbstractLoss end -############################################################################################ -# automatic differentiation -############################################################################################ """ - SemFiniteDiff(;observed = SemObservedData, implied = RAM, loss = SemML, kwargs...) + abstract type AbstractSem -A wrapper around [`Sem`](@ref) that substitutes dedicated evaluation of gradient and hessian with -finite difference approximation. +The base type for all SEMs. +""" +abstract type AbstractSem end -# Arguments -- `observed`: object of subtype `SemObserved` or a constructor. -- `implied`: object of subtype `SemImplied` or a constructor. -- `loss`: object of subtype `SemLossFunction`s or constructor; or a tuple of such. - -Returns a Sem with fields -- `observed::SemObserved`: Stores observed data, sample statistics, etc. See also [`SemObserved`](@ref). -- `implied::SemImplied`: Computes model implied statistics, like Σ, μ, etc. See also [`SemImplied`](@ref). -- `loss::SemLoss`: Computes the objective and gradient of a sum of loss functions. See also [`SemLoss`](@ref). """ -struct SemFiniteDiff{O <: SemObserved, I <: SemImplied, L <: SemLoss} <: - AbstractSemSingle{O, I, L} - observed::O - implied::I + struct LossTerm{L, I, W} + +A term of a [`Sem`](@ref) model that wraps [`AbstractLoss`](@ref) loss function of type `L`. +Loss term can have an optional *id* of type `I` and *weight* of numeric type `W`. +""" +struct LossTerm{L <: AbstractLoss, I <: Union{Symbol, Nothing}, W <: Union{Number, Nothing}} loss::L + id::I + weight::W end -############################################################################################ -# ensemble models -############################################################################################ """ - (1) SemEnsemble(models...; weights = nothing, groups = nothing, kwargs...) + Sem(loss_terms...; [params], kwargs...) - (2) SemEnsemble(;specification, data, groups, column = :group, kwargs...) +SEM model (including multi-group SEMs) that combines all the data, implied SEM structure +and regularization terms. -Constructor for ensemble models. (2) can be used to conveniently specify multigroup models. +All terms of the `Sem` object share the same set of parameters. +`Sem` implements the calculation of the weighted sum of its terms (the *objective* +function), as well as the gradient and Hessian of this sum. # Arguments -- `models...`: `AbstractSem`s. -- `weights::Vector`: Weights for each model. Defaults to the number of observed data points. -- `specification::EnsembleParameterTable`: Model specification. -- `data::DataFrame`: Observed data. Must contain a `column` of type `Vector{Symbol}` that contains the group. -- `groups::Vector{Symbol}`: Group names. -- `column::Symbol`: Name of the column in `data` that contains the group. - -All additional kwargs are passed down to the model parts. - -Returns a SemEnsemble with fields -- `n::Int`: Number of models. -- `sems::Tuple`: `AbstractSem`s. -- `weights::Vector`: Weights for each model. -- `param_labels::Vector`: Stores parameter labels and their position. - -For instructions on multigroup models, see the online documentation. +- `loss_terms...`: [`AbstractLoss`](@ref) objects, including SEM losses ([`SemLoss`](@ref)), + optionally can be a pair of a loss object and its numeric weight + +# Fields +- `loss_terms::Tuple`: a tuple of all loss functions and their weights +- `params::Vector{Symbol}`: the vector of parameter ids shared by all loss functions. """ -struct SemEnsemble{N, T <: Tuple, V <: AbstractVector, I, G <: Vector{Symbol}} <: - AbstractSemCollection - n::N - sems::T - weights::V - param_labels::I - groups::G +struct Sem{L <: Tuple} <: AbstractSem + loss_terms::L + params::Vector{Symbol} end -# constructor from multiple models -function SemEnsemble(models...; weights = nothing, groups = nothing, kwargs...) - n = length(models) - # default weights - weights = isnothing(weights) ? multigroup_weights(models, n) : weights - # default group labels - groups = isnothing(groups) ? Symbol.(:g, 1:n) : groups - # check parameters equality - param_labels = SEM.param_labels(models[1]) - for model in models - if param_labels != SEM.param_labels(model) - throw(ErrorException("The parameters of your models do not match. \n - Maybe you tried to specify models of an ensemble via ParameterTables. \n - In that case, you may use RAMMatrices instead.")) - end - end - - return SemEnsemble(n, models, weights, param_labels, groups) -end +############################################################################################ +# automatic differentiation +############################################################################################ -# constructor from EnsembleParameterTable and data set -function SemEnsemble(; specification, data, groups, column = :group, kwargs...) - if specification isa EnsembleParameterTable - specification = convert(Dict{Symbol, RAMMatrices}, specification) - end - models = [] - for group in groups - ram_matrices = specification[group] - data_group = select(filter(r -> r[column] == group, data), Not(column)) - if iszero(nrow(data_group)) - error("Your data does not contain any observations from group `$(group)`.") - end - model = Sem(; specification = ram_matrices, data = data_group, kwargs...) - push!(models, model) - end - return SemEnsemble(models...; groups = groups, kwargs...) -end +""" + SemFiniteDiff(model::AbstractSem) -function multigroup_weights(models, n) - nsamples_total = sum(nsamples, models) - uniform_lossfun = check_single_lossfun(models...; throw_error = false) - if !uniform_lossfun - @info "Your ensemble model contains heterogeneous loss functions. - Default weights of (#samples per group/#total samples) will be used." - return [(nsamples(model)) / (nsamples_total) for model in models] - end - lossfun = models[1].loss.functions[1] - if !applicable(mg_correction, lossfun) - @info "We don't know how to choose group weights for the specified loss function. - Default weights of (#samples per group/#total samples) will be used." - return [(nsamples(model)) / (nsamples_total) for model in models] - end - c = mg_correction(lossfun) - return [(nsamples(model)+c) / (nsamples_total+n*c) for model in models] -end +A wrapper around [`AbstractSem`](@ref) that substitutes dedicated evaluation of gradient and +hessian with finite difference approximation. -param_labels(ensemble::SemEnsemble) = ensemble.param_labels +`SemFiniteDiff` could be used to enable gradient-based optimization of the SEM models +when the dedicated calculation of gradient and hessian are not available. +For approximation, it uses the *FiniteDiff.jl* package. +# Arguments +- `model::Sem`: the SEM model to wrap """ - n_models(ensemble::SemEnsemble) -> Integer +struct SemFiniteDiff{S <: AbstractSem} <: AbstractSem + model::S +end -Returns the number of models in an ensemble model. -""" -n_models(ensemble::SemEnsemble) = ensemble.n -""" - models(ensemble::SemEnsemble) -> Tuple{AbstractSem} +struct LossFiniteDiff{L <: AbstractLoss} <: AbstractLoss + loss::L +end -Returns the models in an ensemble model. -""" -models(ensemble::SemEnsemble) = ensemble.sems -""" - weights(ensemble::SemEnsemble) -> Vector +struct SemLossFiniteDiff{O, I, L <: SemLoss{O, I}} <: SemLoss{O, I} + loss::L +end -Returns the weights of an ensemble model. """ -weights(ensemble::SemEnsemble) = ensemble.weights + abstract type SemSpecification end -""" Base type for all SEM specifications. """ abstract type SemSpecification end diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl index bb7db3b5..48723fbe 100644 --- a/test/examples/multigroup/build_models.jl +++ b/test/examples/multigroup/build_models.jl @@ -4,34 +4,36 @@ const SEM = StructuralEquationModels # ML estimation ############################################################################################ -model_g1 = Sem(specification = specification_g1, data = dat_g1, implied = RAMSymbolic) +obs_g1 = SemObservedData(data = dat_g1, observed_vars = SEM.observed_vars(specification_g1)) +obs_g2 = SemObservedData(data = dat_g2, observed_vars = SEM.observed_vars(specification_g2)) -model_g2 = Sem(specification = specification_g2, data = dat_g2, implied = RAM) +model_ml_multigroup = Sem( + :Pasteur => SemML(obs_g1, RAMSymbolic(specification_g1)), + :Grant_White => SemML(obs_g2, RAM(specification_g2)), +) -@test SEM.param_labels(model_g1.implied.ram_matrices) == - SEM.param_labels(model_g2.implied.ram_matrices) +@testset "Sem API" begin + @test SEM.nsamples(model_ml_multigroup) == nsamples(obs_g1) + nsamples(obs_g2) + @test SEM.nsem_terms(model_ml_multigroup) == 2 + @test length(SEM.sem_terms(model_ml_multigroup)) == 2 +end -# test the different constructors -model_ml_multigroup = SemEnsemble(model_g1, model_g2; groups = [:Pasteur, :Grant_White]) -model_ml_multigroup2 = SemEnsemble( - specification = partable, - data = dat, - column = :school, - groups = [:Pasteur, :Grant_White], - loss = SemML, +# replace observed using Dict of data matrices +model_ml_multigroup3 = replace_observed( + model_ml_multigroup, + Dict(:Pasteur => dat_g1, :Grant_White => dat_g2), ) -model_ml_multigroup3 = replace_observed( - model_ml_multigroup2, - column = :school, - specification = partable, - data = dat, +# replace observed using DataFrame with group column +model_ml_multigroup4 = replace_observed( + model_ml_multigroup, + dat; + semterm_column = :school, ) # gradients @testset "ml_gradients_multigroup" begin test_gradient(model_ml_multigroup, start_test; atol = 1e-9) - test_gradient(model_ml_multigroup2, start_test; atol = 1e-9) end # fit @@ -44,50 +46,18 @@ end atol = 1e-4, lav_groups = Dict(:Pasteur => 1, :Grant_White => 2), ) - solution = fit(semoptimizer, model_ml_multigroup2) - update_estimate!(partable, solution) - test_estimates( - partable, - solution_lav[:parameter_estimates_ml]; - atol = 1e-4, - lav_groups = Dict(:Pasteur => 1, :Grant_White => 2), - ) end @testset "replace_observed_multigroup" begin sem_fit_1 = fit(semoptimizer, model_ml_multigroup) - sem_fit_2 = fit(semoptimizer, model_ml_multigroup3) - @test sem_fit_1.solution ≈ sem_fit_2.solution + sem_fit_3 = fit(semoptimizer, model_ml_multigroup3) + @test sem_fit_1.solution ≈ sem_fit_3.solution + sem_fit_4 = fit(semoptimizer, model_ml_multigroup4) + @test sem_fit_1.solution ≈ sem_fit_4.solution end @testset "fitmeasures/se_ml" begin - solution_ml = fit(model_ml_multigroup) - test_fitmeasures( - fit_measures(solution_ml), - solution_lav[:fitmeasures_ml]; - rtol = 1e-2, - atol = 1e-7, - ) - test_fitmeasures( - Dict(:CFI => CFI(solution_ml)), - solution_lav[:fitmeasures_ml]; - fitmeasure_names = Dict(:CFI => "cfi"), - ) - - update_se_hessian!(partable, solution_ml) - test_estimates( - partable, - solution_lav[:parameter_estimates_ml]; - atol = 1e-3, - col = :se, - lav_col = :se, - lav_groups = Dict(:Pasteur => 1, :Grant_White => 2), - ) - - test_bootstrap(solution_ml, partable; rtol_hessian = 0.3) - smoketest_CI_z(solution_ml, partable) - - solution_ml = fit(model_ml_multigroup2) + solution_ml = fit(semoptimizer, model_ml_multigroup) test_fitmeasures( fit_measures(solution_ml), solution_lav[:fitmeasures_ml]; @@ -118,15 +88,19 @@ end partable_s = sort_vars(partable) specification_s = convert(Dict{Symbol, RAMMatrices}, partable_s) +obs_g1_s = SemObservedData( + data = dat_g1, + observed_vars = SEM.observed_vars(specification_s[:Pasteur]), +) +obs_g2_s = SemObservedData( + data = dat_g2, + observed_vars = SEM.observed_vars(specification_s[:Grant_White]), +) -specification_g1_s = specification_s[:Pasteur] -specification_g2_s = specification_s[:Grant_White] - -model_g1 = Sem(specification = specification_g1_s, data = dat_g1, implied = RAMSymbolic) - -model_g2 = Sem(specification = specification_g2_s, data = dat_g2, implied = RAM) - -model_ml_multigroup = SemEnsemble(model_g1, model_g2; optimizer = semoptimizer) +model_ml_multigroup = Sem( + SemML(obs_g1_s, RAMSymbolic(specification_s[:Pasteur])), + SemML(obs_g2_s, RAM(specification_s[:Grant_White])), +) # gradients @testset "ml_gradients_multigroup | sorted" begin @@ -142,7 +116,7 @@ grad_fd = FiniteDiff.finite_difference_gradient( # fit @testset "ml_solution_multigroup | sorted" begin - solution = fit(model_ml_multigroup) + solution = fit(semoptimizer, model_ml_multigroup) update_estimate!(partable_s, solution) test_estimates( partable_s, @@ -153,7 +127,7 @@ grad_fd = FiniteDiff.finite_difference_gradient( end @testset "fitmeasures/se_ml | sorted" begin - solution_ml = fit(model_ml_multigroup) + solution_ml = fit(semoptimizer, model_ml_multigroup) test_fitmeasures( fit_measures(solution_ml), solution_lav[:fitmeasures_ml]; @@ -178,28 +152,26 @@ end end @testset "sorted | LowerTriangular A" begin - @test implied(model_ml_multigroup.sems[2]).A isa LowerTriangular + @test implied(SEM.sem_terms(model_ml_multigroup)[2]).A isa LowerTriangular end ############################################################################################ # ML estimation - user defined loss function ############################################################################################ -struct UserSemML <: SemLossFunction +struct UserSemML{O, I} <: SemLoss{O, I} hessianeval::ExactHessian - UserSemML() = new(ExactHessian()) -end - -############################################################################################ -### functors -############################################################################################ + observed::O + implied::I -using LinearAlgebra: isposdef, logdet, tr, inv + UserSemML(observed::SemObserved, implied::SemImplied) = + new{typeof(observed), typeof(implied)}(ExactHessian(), observed, implied) +end -function SEM.objective(ml::UserSemML, model::AbstractSem, params) - Σ = implied(model).Σ - Σₒ = SEM.obs_cov(observed(model)) +function SEM.objective(ml::UserSemML, params) + Σ = implied(ml).Σ + Σₒ = SEM.obs_cov(observed(ml)) if !isposdef(Σ) return Inf else @@ -208,24 +180,18 @@ function SEM.objective(ml::UserSemML, model::AbstractSem, params) end # models -model_g1 = Sem(specification = specification_g1, data = dat_g1, implied = RAMSymbolic) - -model_g2 = SemFiniteDiff( - specification = specification_g2, - data = dat_g2, - implied = RAMSymbolic, - loss = UserSemML(), +model_ml_multigroup = Sem( + SemML(obs_g1, RAMSymbolic(specification_g1)), + SEM.FiniteDiffWrapper(UserSemML(obs_g2, RAMSymbolic(specification_g2))), ) -model_ml_multigroup = SemEnsemble(model_g1, model_g2; optimizer = semoptimizer) - @testset "gradients_user_defined_loss" begin test_gradient(model_ml_multigroup, start_test; atol = 1e-9) end # fit @testset "solution_user_defined_loss" begin - solution = fit(model_ml_multigroup) + solution = fit(semoptimizer, model_ml_multigroup) update_estimate!(partable, solution) test_estimates( partable, @@ -239,25 +205,9 @@ end # GLS estimation ############################################################################################ -model_ls_g1 = Sem( - specification = specification_g1, - data = dat_g1, - implied = RAMSymbolic, - loss = SemWLS, -) - -model_ls_g2 = Sem( - specification = specification_g2, - data = dat_g2, - implied = RAMSymbolic, - loss = SemWLS, -) - -model_ls_multigroup = SemEnsemble( - model_ls_g1, - model_ls_g2; - groups = [:Pasteur, :Grant_White], - optimizer = semoptimizer, +model_ls_multigroup = Sem( + SemWLS(obs_g1, RAMSymbolic(specification_g1, vech = true)), + SemWLS(obs_g2, RAMSymbolic(specification_g2, vech = true)), ) @testset "ls_gradients_multigroup" begin @@ -265,7 +215,7 @@ model_ls_multigroup = SemEnsemble( end @testset "ls_solution_multigroup" begin - solution = fit(model_ls_multigroup) + solution = fit(semoptimizer, model_ls_multigroup) update_estimate!(partable, solution) test_estimates( partable, @@ -276,7 +226,7 @@ end end @testset "fitmeasures/se_ls" begin - solution_ls = fit(model_ls_multigroup) + solution_ls = fit(semoptimizer, model_ls_multigroup) test_fitmeasures( fit_measures(solution_ls), solution_lav[:fitmeasures_ls]; @@ -308,40 +258,27 @@ end ############################################################################################ if !isnothing(specification_miss_g1) - model_g1 = Sem( - specification = specification_miss_g1, - observed = SemObservedMissing, - loss = SemFIML, - data = dat_miss_g1, - implied = RAM, - meanstructure = true, - ) - - model_g2 = Sem( - specification = specification_miss_g2, - observed = SemObservedMissing, - loss = SemFIML, - data = dat_miss_g2, - implied = RAM, - meanstructure = true, - ) - - model_ml_multigroup = SemEnsemble(model_g1, model_g2) - model_ml_multigroup2 = SemEnsemble( - specification = partable_miss, - data = dat_missing, - column = :school, - groups = [:Pasteur, :Grant_White], - loss = SemFIML, - observed = SemObservedMissing, - meanstructure = true, + model_ml_multigroup = Sem( + SemFIML( + SemObservedMissing( + data = dat_miss_g1, + observed_vars = SEM.observed_vars(specification_miss_g1), + ), + RAM(specification_miss_g1), + ), + SemFIML( + SemObservedMissing( + data = dat_miss_g2, + observed_vars = SEM.observed_vars(specification_miss_g2), + ), + RAM(specification_miss_g2), + ), ) - model_ml_varonly = SemEnsemble( + model_ml_varonly = Sem( specification = partable_varonly, data = dat_missing, - column = :school, - groups = [:Pasteur, :Grant_White], + semterm_column = :school, loss = SemFIML, observed = SemObservedMissing, meanstructure = true, @@ -373,7 +310,6 @@ if !isnothing(specification_miss_g1) @testset "fiml_gradients_multigroup" begin test_gradient(model_ml_multigroup, start_test; atol = 1e-7) - test_gradient(model_ml_multigroup2, start_test; atol = 1e-7) end @testset "fiml_solution_multigroup" begin @@ -385,14 +321,6 @@ if !isnothing(specification_miss_g1) atol = 1e-4, lav_groups = Dict(:Pasteur => 1, :Grant_White => 2), ) - solution = fit(semoptimizer, model_ml_multigroup2) - update_estimate!(partable_miss, solution) - test_estimates( - partable_miss, - solution_lav[:parameter_estimates_fiml]; - atol = 1e-4, - lav_groups = Dict(:Pasteur => 1, :Grant_White => 2), - ) end @testset "fitmeasures/se_fiml" begin diff --git a/test/examples/political_democracy/by_parts.jl b/test/examples/political_democracy/by_parts.jl index d2d468a9..6866eead 100644 --- a/test/examples/political_democracy/by_parts.jl +++ b/test/examples/political_democracy/by_parts.jl @@ -2,103 +2,83 @@ ### models w.o. meanstructure ############################################################################################ -# observed --------------------------------------------------------------------------------- -observed = SemObservedData(specification = spec, data = dat) +semoptimizer = SemOptimizer(engine = opt_engine) -# implied -implied_ram = RAM(specification = spec) +model_ml = Sem(specification = spec, data = dat) +@test SEM.params(model_ml) == SEM.params(spec) -implied_ram_sym = RAMSymbolic(specification = spec) +model_ls_sym = + Sem(specification = spec, data = dat, implied = RAMSymbolic, vech = true, loss = SemWLS) -# loss functions --------------------------------------------------------------------------- -ml = SemML(observed = observed) +model_ml_sym = Sem(specification = spec, data = dat, implied = RAMSymbolic) -wls = SemWLS(observed = observed) - -ridge = SemRidge(α_ridge = 0.001, which_ridge = 16:20, nparams = 31) - -constant = SemConstant(constant_loss = 3.465) - -# loss ------------------------------------------------------------------------------------- -loss_ml = SemLoss(ml) - -loss_wls = SemLoss(wls) - -# optimizer ------------------------------------------------------------------------------------- -optimizer_obj = SemOptimizer(engine = opt_engine) - -# models ----------------------------------------------------------------------------------- - -model_ml = Sem(observed, implied_ram, loss_ml) - -model_ls_sym = Sem(observed, RAMSymbolic(specification = spec, vech = true), loss_wls) - -model_ml_sym = Sem(observed, implied_ram_sym, loss_ml) - -model_ridge = Sem(observed, implied_ram, SemLoss(ml, ridge)) +model_ml_ridge = Sem( + specification = spec, + data = dat, + loss = (SemML, SemRidge), + α_ridge = 0.001, + which_ridge = 16:20, +) -model_constant = Sem(observed, implied_ram, SemLoss(ml, constant)) +model_ml_const = Sem( + specification = spec, + data = dat, + loss = (SemML, SemConstant), + constant_loss = 3.465, +) -model_ml_weighted = - Sem(observed, implied_ram, SemLoss(ml; loss_weights = [nsamples(model_ml)])) +model_ml_weighted = Sem(SemML(SemObservedData(data = dat), RAM(spec)) => nsamples(model_ml)) ############################################################################################ ### test gradients ############################################################################################ -models = - [model_ml, model_ls_sym, model_ridge, model_constant, model_ml_sym, model_ml_weighted] -model_names = ["ml", "ls_sym", "ridge", "constant", "ml_sym", "ml_weighted"] +models = Dict( + "ml" => model_ml, + "ls_sym" => model_ls_sym, + "ml_ridge" => model_ml_ridge, + "ml_const" => model_ml_const, + "ml_sym" => model_ml_sym, + "ml_weighted" => model_ml_weighted, +) -for (model, name) in zip(models, model_names) - try - @testset "$(name)_gradient" begin - test_gradient(model, start_test; rtol = 1e-9) - end - catch - end +@testset "$(id)_gradient" for (id, model) in pairs(models) + test_gradient(model, start_test; rtol = 1e-9) end ############################################################################################ ### test solution ############################################################################################ -models = [model_ml, model_ls_sym, model_ml_sym, model_constant] -model_names = ["ml", "ls_sym", "ml_sym", "constant"] -solution_names = Symbol.("parameter_estimates_" .* ["ml", "ls", "ml", "ml"]) - -for (model, name, solution_name) in zip(models, model_names, solution_names) - try - @testset "$(name)_solution" begin - solution = fit(optimizer_obj, model) - update_estimate!(partable, solution) - test_estimates(partable, solution_lav[solution_name]; atol = 1e-2) - end - catch - end +@testset "$(id)_solution" for id in ["ml", "ls_sym", "ml_sym", "ml_const"] + model = models[id] + solution = fit(semoptimizer, model) + sol_name = Symbol("parameter_estimates_", replace(id, r"_.+$" => "")) + update_estimate!(partable, solution) + test_estimates(partable, solution_lav[sol_name]; atol = 1e-2) end @testset "ridge_solution" begin - solution_ridge = fit(optimizer_obj, model_ridge) - solution_ml = fit(optimizer_obj, model_ml) - # solution_ridge_id = fit(optimizer_obj, model_ridge_id) - @test solution_ridge.minimum < solution_ml.minimum + 1 + solution_ridge = fit(semoptimizer, model_ml_ridge) + solution_ml = fit(semoptimizer, model_ml) + # solution_ridge_id = fit(model_ridge_id) + @test abs(solution_ridge.minimum - solution_ml.minimum) < 1 end # test constant objective value @testset "constant_objective_and_gradient" begin - @test (objective!(model_constant, start_test) - 3.465) ≈ + @test (objective!(model_ml_const, start_test) - 3.465) ≈ objective!(model_ml, start_test) grad = similar(start_test) grad2 = similar(start_test) - gradient!(grad, model_constant, start_test) + gradient!(grad, model_ml_const, start_test) gradient!(grad2, model_ml, start_test) @test grad ≈ grad2 end @testset "ml_solution_weighted" begin - solution_ml = fit(optimizer_obj, model_ml) - solution_ml_weighted = fit(optimizer_obj, model_ml_weighted) + solution_ml = fit(semoptimizer, model_ml) + solution_ml_weighted = fit(semoptimizer, model_ml_weighted) @test solution(solution_ml) ≈ solution(solution_ml_weighted) rtol = 1e-3 @test nsamples(model_ml) * StructuralEquationModels.minimum(solution_ml) ≈ StructuralEquationModels.minimum(solution_ml_weighted) rtol = 1e-6 @@ -109,7 +89,7 @@ end ############################################################################################ @testset "fitmeasures/se_ml" begin - solution_ml = fit(optimizer_obj, model_ml) + solution_ml = fit(semoptimizer, model_ml) test_fitmeasures(fit_measures(solution_ml), solution_lav[:fitmeasures_ml]; atol = 1e-3) test_fitmeasures( Dict(:CFI => CFI(solution_ml)), @@ -128,20 +108,15 @@ end end @testset "fitmeasures/se_ls" begin - solution_ls = fit(optimizer_obj, model_ls_sym) + solution_ls = fit(semoptimizer, model_ls_sym) fm = fit_measures(solution_ls) test_fitmeasures( - fm, + merge(fm, Dict(:CFI => CFI(solution_ls))), solution_lav[:fitmeasures_ls]; atol = 1e-3, - fitmeasure_names = fitmeasure_names_ls, - ) - test_fitmeasures( - Dict(:CFI => CFI(solution_ls)), - solution_lav[:fitmeasures_ls]; - fitmeasure_names = Dict(:CFI => "cfi"), + fitmeasure_names = merge(fitmeasure_names_ls, Dict(:CFI => "cfi")) ) - @test (fm[:AIC] === missing) & (fm[:BIC] === missing) & (fm[:minus2ll] === missing) + @test ismissing(fm[:AIC]) && ismissing(fm[:BIC]) && ismissing(fm[:minus2ll]) @suppress update_se_hessian!(partable, solution_ls) test_estimates( @@ -160,22 +135,22 @@ end if opt_engine == :Optim using Optim, LineSearches - optimizer_obj = SemOptimizer( - engine = opt_engine, - algorithm = Newton(; - linesearch = BackTracking(order = 3), - alphaguess = InitialHagerZhang(), - ), + model_ls = Sem( + data = dat, + specification = spec, + implied = RAMSymbolic, + loss = SemWLS, + vech = true, + hessian = true, ) - implied_sym_hessian_vech = - RAMSymbolic(specification = spec, vech = true, hessian = true) - - implied_sym_hessian = RAMSymbolic(specification = spec, hessian = true) - - model_ls = Sem(observed, implied_sym_hessian_vech, loss_wls) - - model_ml = Sem(observed, implied_sym_hessian, loss_ml) + model_ml = Sem( + data = dat, + specification = spec, + implied = RAMSymbolic, + loss = SemML, + hessian = true, + ) @testset "ml_hessians" begin test_hessian(model_ml, start_test; atol = 1e-4) @@ -186,13 +161,23 @@ if opt_engine == :Optim end @testset "ml_solution_hessian" begin - solution = fit(optimizer_obj, model_ml) + solution = fit(SemOptimizer(engine = :Optim, algorithm = Newton()), model_ml) + update_estimate!(partable, solution) test_estimates(partable, solution_lav[:parameter_estimates_ml]; atol = 1e-2) end @testset "ls_solution_hessian" begin - solution = fit(optimizer_obj, model_ls) + solution = fit( + SemOptimizer( + engine = :Optim, + algorithm = Newton( + linesearch = BackTracking(order = 3), + alphaguess = InitialHagerZhang(), + ), + ), + model_ls, + ) update_estimate!(partable, solution) test_estimates( partable, @@ -207,69 +192,47 @@ end ### meanstructure ############################################################################################ -# observed --------------------------------------------------------------------------------- -observed = SemObservedData(specification = spec_mean, data = dat, meanstructure = true) - -# implied -implied_ram = RAM(specification = spec_mean, meanstructure = true) - -implied_ram_sym = RAMSymbolic(specification = spec_mean, meanstructure = true) - -# loss functions --------------------------------------------------------------------------- -ml = SemML(observed = observed, meanstructure = true) - -wls = SemWLS(observed = observed, meanstructure = true) - -# loss ------------------------------------------------------------------------------------- -loss_ml = SemLoss(ml) - -loss_wls = SemLoss(wls) - -# optimizer ------------------------------------------------------------------------------------- -optimizer_obj = SemOptimizer(engine = opt_engine) +# models +model_ls = Sem( + data = dat, + specification = spec_mean, + implied = RAMSymbolic, + loss = SemWLS, + vech = true, +) -# models ----------------------------------------------------------------------------------- -model_ml = Sem(observed, implied_ram, loss_ml) +model_ml = Sem(data = dat, specification = spec_mean, implied = RAM, loss = SemML) -model_ls = Sem( - observed, - RAMSymbolic(specification = spec_mean, meanstructure = true, vech = true), - loss_wls, +model_ml_cov = Sem( + specification = spec, + observed = SemObservedCovariance, + obs_cov = cov(Matrix(dat)), + observed_vars = Symbol.(names(dat)), + nsamples = 75, ) -model_ml_sym = Sem(observed, implied_ram_sym, loss_ml) +model_ml_sym = + Sem(data = dat, specification = spec_mean, implied = RAMSymbolic, loss = SemML) ############################################################################################ ### test gradients ############################################################################################ -models = [model_ml, model_ls, model_ml_sym] -model_names = ["ml", "ls_sym", "ml_sym"] +models = Dict("ml" => model_ml, "ls_sym" => model_ls, "ml_sym" => model_ml_sym) -for (model, name) in zip(models, model_names) - try - @testset "$(name)_gradient_mean" begin - test_gradient(model, start_test_mean; rtol = 1e-9) - end - catch - end +@testset "$(id)_gradient_mean" for (id, model) in pairs(models) + test_gradient(model, start_test_mean; rtol = 1e-9) end ############################################################################################ ### test solution ############################################################################################ -solution_names = Symbol.("parameter_estimates_" .* ["ml", "ls", "ml"] .* "_mean") - -for (model, name, solution_name) in zip(models, model_names, solution_names) - try - @testset "$(name)_solution_mean" begin - solution = fit(optimizer_obj, model) - update_estimate!(partable_mean, solution) - test_estimates(partable_mean, solution_lav[solution_name]; atol = 1e-2) - end - catch - end +@testset "$(id)_solution_mean" for (id, model) in pairs(models) + solution = fit(semoptimizer, model, start_val = start_test_mean) + update_estimate!(partable_mean, solution) + sol_name = Symbol("parameter_estimates_", replace(id, r"_.+$" => ""), "_mean") + test_estimates(partable_mean, solution_lav[sol_name]; atol = 1e-2) end ############################################################################################ @@ -277,7 +240,7 @@ end ############################################################################################ @testset "fitmeasures/se_ml_mean" begin - solution_ml = fit(optimizer_obj, model_ml) + solution_ml = fit(semoptimizer, model_ml) test_fitmeasures( fit_measures(solution_ml), solution_lav[:fitmeasures_ml_mean]; @@ -300,20 +263,15 @@ end end @testset "fitmeasures/se_ls_mean" begin - solution_ls = fit(optimizer_obj, model_ls) + solution_ls = fit(semoptimizer, model_ls) fm = fit_measures(solution_ls) test_fitmeasures( - fm, + merge(fm, Dict(:CFI => CFI(solution_ls))), solution_lav[:fitmeasures_ls_mean]; atol = 1e-3, - fitmeasure_names = fitmeasure_names_ls, - ) - test_fitmeasures( - Dict(:CFI => CFI(solution_ls)), - solution_lav[:fitmeasures_ls_mean]; - fitmeasure_names = Dict(:CFI => "cfi"), + fitmeasure_names = merge(fitmeasure_names_ls, Dict(:CFI => "cfi")), ) - @test (fm[:AIC] === missing) & (fm[:BIC] === missing) & (fm[:minus2ll] === missing) + @test ismissing(fm[:AIC]) && ismissing(fm[:BIC]) && ismissing(fm[:minus2ll]) @suppress update_se_hessian!(partable_mean, solution_ls) test_estimates( @@ -329,16 +287,22 @@ end ### fiml ############################################################################################ -observed = - SemObservedMissing(specification = spec_mean, data = dat_missing, rtol_em = 1e-10) - -fiml = SemFIML(observed = observed, specification = spec_mean) - -loss_fiml = SemLoss(fiml) - -model_ml = Sem(observed, implied_ram, loss_fiml) +# models +model_ml = Sem( + data = dat_missing, + observed = SemObservedMissing, + specification = spec_mean, + implied = RAM, + loss = SemFIML, +) -model_ml_sym = Sem(observed, implied_ram_sym, loss_fiml) +model_ml_sym = Sem( + data = dat_missing, + observed = SemObservedMissing, + specification = spec_mean, + implied = RAMSymbolic, + loss = SemFIML, +) ############################################################################################ ### test gradients @@ -357,13 +321,13 @@ end ############################################################################################ @testset "fiml_solution" begin - solution = fit(optimizer_obj, model_ml) + solution = fit(semoptimizer, model_ml) update_estimate!(partable_mean, solution) test_estimates(partable_mean, solution_lav[:parameter_estimates_fiml]; atol = 1e-2) end @testset "fiml_solution_symbolic" begin - solution = fit(optimizer_obj, model_ml_sym) + solution = fit(semoptimizer, model_ml_sym, start_val = start_test_mean) update_estimate!(partable_mean, solution) test_estimates(partable_mean, solution_lav[:parameter_estimates_fiml]; atol = 1e-2) end @@ -373,7 +337,7 @@ end ############################################################################################ @testset "fitmeasures/se_fiml" begin - solution_ml = fit(optimizer_obj, model_ml) + solution_ml = fit(semoptimizer, model_ml) test_fitmeasures( fit_measures(solution_ml), solution_lav[:fitmeasures_fiml]; @@ -384,7 +348,7 @@ end test_estimates( partable_mean, solution_lav[:parameter_estimates_fiml]; - atol = 1e-3, + atol = 0.002, col = :se, lav_col = :se, ) diff --git a/test/examples/political_democracy/constraints.jl b/test/examples/political_democracy/constraints.jl index 7a6670fa..0291e7ea 100644 --- a/test/examples/political_democracy/constraints.jl +++ b/test/examples/political_democracy/constraints.jl @@ -50,7 +50,7 @@ end @test solution_constrained.solution[31] * solution_constrained.solution[30] >= (0.6 - 1e-8) - @test all(abs.(solution_constrained.solution) .< 10) - @test solution_constrained.optimization_result.result[3] == :FTOL_REACHED + @test all(p -> abs(p) < 10, solution_constrained.solution) + @test solution_constrained.optimization_result.result[3] == :FTOL_REACHED skip = true @test solution_constrained.minimum <= 21.21 + 0.01 end diff --git a/test/examples/political_democracy/constructor.jl b/test/examples/political_democracy/constructor.jl index 1c1c42e5..25a6da91 100644 --- a/test/examples/political_democracy/constructor.jl +++ b/test/examples/political_democracy/constructor.jl @@ -5,21 +5,22 @@ semoptimizer = SemOptimizer(engine = opt_engine) model_ml = Sem(specification = spec, data = dat) -@test SEM.param_labels(model_ml.implied.ram_matrices) == SEM.param_labels(spec) +@test SEM.param_labels(model_ml) == SEM.param_labels(spec) model_ml_cov = Sem( specification = spec, observed = SemObservedCovariance, obs_cov = cov(Matrix(dat)), - obs_colnames = Symbol.(names(dat)), + observed_vars = Symbol.(names(dat)), nsamples = 75, ) -model_ls_sym = Sem(specification = spec, data = dat, implied = RAMSymbolic, loss = SemWLS) +model_ls_sym = + Sem(specification = spec, data = dat, implied = RAMSymbolic, vech = true, loss = SemWLS) model_ml_sym = Sem(specification = spec, data = dat, implied = RAMSymbolic) -model_ridge = Sem( +model_ml_ridge = Sem( specification = spec, data = dat, loss = (SemML, SemRidge), @@ -27,7 +28,7 @@ model_ridge = Sem( which_ridge = 16:20, ) -model_constant = Sem( +model_ml_const = Sem( specification = spec, data = dat, loss = (SemML, SemConstant), @@ -35,65 +36,52 @@ model_constant = Sem( ) model_ml_weighted = - Sem(specification = partable, data = dat, loss_weights = (nsamples(model_ml),)) + Sem(SemML(SemObservedData(data = dat), RAMSymbolic(spec)) => nsamples(model_ml)) ############################################################################################ ### test gradients ############################################################################################ -models = [ - model_ml, - model_ml_cov, - model_ls_sym, - model_ridge, - model_constant, - model_ml_sym, - model_ml_weighted, -] -model_names = ["ml", "ml_cov", "ls_sym", "ridge", "constant", "ml_sym", "ml_weighted"] - -for (model, name) in zip(models, model_names) - try - @testset "$(name)_gradient" begin - test_gradient(model, start_test; rtol = 1e-9) - end - catch - end +models = Dict( + "ml" => model_ml, + "ml_cov" => model_ml_cov, + "ls_sym" => model_ls_sym, + "ridge" => model_ml_ridge, + "ml_const" => model_ml_const, + "ml_sym" => model_ml_sym, + "ml_weighted" => model_ml_weighted, +) + +@testset "$(id)_gradient" for (id, model) in pairs(models) + test_gradient(model, start_test; rtol = 1e-9) end ############################################################################################ ### test solution ############################################################################################ -models = [model_ml, model_ml_cov, model_ls_sym, model_ml_sym, model_constant] -model_names = ["ml", "ml_cov", "ls_sym", "ml_sym", "constant"] -solution_names = Symbol.("parameter_estimates_" .* ["ml", "ml", "ls", "ml", "ml"]) - -for (model, name, solution_name) in zip(models, model_names, solution_names) - try - @testset "$(name)_solution" begin - solution = fit(semoptimizer, model) - update_estimate!(partable, solution) - test_estimates(partable, solution_lav[solution_name]; atol = 1e-2) - end - catch - end +@testset "$(id)_solution" for id in ["ml", "ml_cov", "ls_sym", "ml_sym", "ml_const"] + model = models[id] + solution = fit(semoptimizer, model) + sol_name = Symbol("parameter_estimates_", replace(id, r"_.+$" => "")) + update_estimate!(partable, solution) + test_estimates(partable, solution_lav[sol_name]; atol = 1e-2) end @testset "ridge_solution" begin - solution_ridge = fit(semoptimizer, model_ridge) + solution_ridge = fit(semoptimizer, model_ml_ridge) solution_ml = fit(semoptimizer, model_ml) - # solution_ridge_id = fit(semoptimizer, model_ridge_id) + # solution_ridge_id = fit(model_ridge_id) @test abs(solution_ridge.minimum - solution_ml.minimum) < 1 end # test constant objective value @testset "constant_objective_and_gradient" begin - @test (objective!(model_constant, start_test) - 3.465) ≈ + @test (objective!(model_ml_const, start_test) - 3.465) ≈ objective!(model_ml, start_test) grad = similar(start_test) grad2 = similar(start_test) - gradient!(grad, model_constant, start_test) + gradient!(grad, model_ml_const, start_test) gradient!(grad2, model_ml, start_test) @test grad ≈ grad2 end @@ -101,12 +89,9 @@ end @testset "ml_solution_weighted" begin solution_ml = fit(semoptimizer, model_ml) solution_ml_weighted = fit(semoptimizer, model_ml_weighted) - @test isapprox(solution(solution_ml), solution(solution_ml_weighted), rtol = 1e-3) - @test isapprox( - nsamples(model_ml) * StructuralEquationModels.minimum(solution_ml), - StructuralEquationModels.minimum(solution_ml_weighted), - rtol = 1e-6, - ) + @test solution(solution_ml) ≈ solution(solution_ml_weighted) rtol = 1e-3 + @test nsamples(model_ml) * StructuralEquationModels.minimum(solution_ml) ≈ + StructuralEquationModels.minimum(solution_ml_weighted) rtol = 1e-6 end ############################################################################################ @@ -181,19 +166,14 @@ end ) # set seed for simulation Random.seed!(83472834) - colnames = Symbol.(names(example_data("political_democracy"))) # simulate data model_ml_new = replace_observed( model_ml, - data = rand(model_ml, params, 1_000_000), - specification = spec, - obs_colnames = colnames, + rand(model_ml, params, 1_000_000), ) model_ml_sym_new = replace_observed( model_ml_sym, - data = rand(model_ml_sym, params, 1_000_000), - specification = spec, - obs_colnames = colnames, + rand(model_ml_sym, params, 1_000_000), ) # fit models sol_ml = solution(fit(semoptimizer, model_ml_new)) @@ -211,23 +191,19 @@ if opt_engine == :Optim using Optim, LineSearches model_ls = Sem( - specification = spec, data = dat, - implied = RAMSymbolic, + specification = spec, + observed = SemObservedData, + implied = RAMSymbolic(spec, vech = true, hessian = true), loss = SemWLS, - hessian = true, - algorithm = Newton(; - linesearch = BackTracking(order = 3), - alphaguess = InitialHagerZhang(), - ), ) model_ml = Sem( - specification = spec, data = dat, - implied = RAMSymbolic, - hessian = true, - algorithm = Newton(), + specification = spec, + observed = SemObservedData, + implied = RAMSymbolic(spec, hessian = true), + loss = SemML, ) @testset "ml_hessians" begin @@ -239,13 +215,23 @@ if opt_engine == :Optim end @testset "ml_solution_hessian" begin - solution = fit(semoptimizer, model_ml) + solution = fit(SemOptimizer(engine = :Optim, algorithm = Newton()), model_ml) + update_estimate!(partable, solution) test_estimates(partable, solution_lav[:parameter_estimates_ml]; atol = 1e-2) end @testset "ls_solution_hessian" begin - solution = fit(semoptimizer, model_ls) + solution = fit( + SemOptimizer( + engine = :Optim, + algorithm = Newton( + linesearch = BackTracking(order = 3), + alphaguess = InitialHagerZhang(), + ), + ), + model_ls, + ) update_estimate!(partable, solution) test_estimates( partable, @@ -266,6 +252,7 @@ model_ls = Sem( specification = spec_mean, data = dat, implied = RAMSymbolic, + vech = true, loss = SemWLS, meanstructure = true, ) @@ -277,7 +264,7 @@ model_ml_cov = Sem( observed = SemObservedCovariance, obs_cov = cov(Matrix(dat)), obs_mean = vcat(mean(Matrix(dat), dims = 1)...), - obs_colnames = Symbol.(names(dat)), + observed_vars = Symbol.(names(dat)), meanstructure = true, nsamples = 75, ) @@ -289,33 +276,26 @@ model_ml_sym = ### test gradients ############################################################################################ -models = [model_ml, model_ml_cov, model_ls, model_ml_sym] -model_names = ["ml", "ml_cov", "ls_sym", "ml_sym"] +models = Dict( + "ml" => model_ml, + "ml_cov" => model_ml_cov, + "ls_sym" => model_ls, + "ml_sym" => model_ml_sym, +) -for (model, name) in zip(models, model_names) - try - @testset "$(name)_gradient_mean" begin - test_gradient(model, start_test_mean; rtol = 1e-9) - end - catch - end +@testset "$(id)_gradient_mean" for (id, model) in pairs(models) + test_gradient(model, start_test_mean; rtol = 1e-9) end ############################################################################################ ### test solution ############################################################################################ -solution_names = Symbol.("parameter_estimates_" .* ["ml", "ml", "ls", "ml"] .* "_mean") - -for (model, name, solution_name) in zip(models, model_names, solution_names) - try - @testset "$(name)_solution_mean" begin - solution = fit(semoptimizer, model) - update_estimate!(partable_mean, solution) - test_estimates(partable_mean, solution_lav[solution_name]; atol = 1e-2) - end - catch - end +@testset "$(id)_solution_mean" for (id, model) in pairs(models) + solution = fit(semoptimizer, model, start_val = start_test_mean) + update_estimate!(partable_mean, solution) + sol_name = Symbol("parameter_estimates_", replace(id, r"_.+$" => ""), "_mean") + test_estimates(partable_mean, solution_lav[sol_name]; atol = 1e-2) end ############################################################################################ @@ -395,21 +375,14 @@ end ) # set seed for simulation Random.seed!(83472834) - colnames = Symbol.(names(example_data("political_democracy"))) # simulate data model_ml_new = replace_observed( model_ml, - data = rand(model_ml, params, 1_000_000), - specification = spec, - obs_colnames = colnames, - meanstructure = true, + rand(model_ml, params, 1_000_000), ) model_ml_sym_new = replace_observed( model_ml_sym, - data = rand(model_ml_sym, params, 1_000_000), - specification = spec, - obs_colnames = colnames, - meanstructure = true, + rand(model_ml_sym, params, 1_000_000), ) # fit models sol_ml = solution(fit(semoptimizer, model_ml_new)) @@ -474,7 +447,7 @@ end end @testset "fiml_solution_symbolic" begin - solution = fit(semoptimizer, model_ml_sym) + solution = fit(semoptimizer, model_ml_sym, start_val = start_test_mean) update_estimate!(partable_mean, solution) test_estimates(partable_mean, solution_lav[:parameter_estimates_fiml]; atol = 1e-2) end diff --git a/test/examples/recover_parameters/recover_parameters_twofact.jl b/test/examples/recover_parameters/recover_parameters_twofact.jl index a4bd7d5f..ebaaae83 100644 --- a/test/examples/recover_parameters/recover_parameters_twofact.jl +++ b/test/examples/recover_parameters/recover_parameters_twofact.jl @@ -1,5 +1,7 @@ using StructuralEquationModels, Distributions, Random, Optim, LineSearches +SEM = StructuralEquationModels + include( joinpath( chop(dirname(pathof(StructuralEquationModels)), tail = 3), @@ -7,7 +9,7 @@ include( ), ) -x = Symbol.("x", 1:13) +pars = Symbol.("x", 1:13) S = [ :x1 0 0 0 0 0 0 0 @@ -40,7 +42,7 @@ A = [ 0 0 0 0 0 0 0 0 ] -ram_matrices = RAMMatrices(; A = A, S = S, F = F, param_labels = x, vars = nothing) +ram_matrices = RAMMatrices(; A = A, S = S, F = F, param_labels = pars, vars = nothing) true_val = [ repeat([1], 8) @@ -53,19 +55,19 @@ start = [ repeat([0.5], 4) ] -implied_ml = RAMSymbolic(ram_matrices; start_val = start) +implied_sym = RAMSymbolic(ram_matrices) -implied_ml.Σ_eval!(implied_ml.Σ, true_val) +implied_sym.Σ_eval!(implied_sym.Σ, true_val) -true_dist = MultivariateNormal(implied_ml.Σ) +true_dist = MultivariateNormal(implied_sym.Σ) Random.seed!(1234) -x = transpose(rand(true_dist, 100_000)) -semobserved = SemObservedData(data = x, specification = nothing) +x = permutedims(rand(true_dist, 10^5), (2, 1)) + +observed = SemObservedData(data = x, specification = ram_matrices) -loss_ml = SemLoss(SemML(; observed = semobserved, nparams = length(start))) +model_ml = Sem(SemML(observed, implied_sym)) -model_ml = Sem(semobserved, implied_ml, loss_ml) objective!(model_ml, true_val) optimizer = SemOptimizer( @@ -73,6 +75,6 @@ optimizer = SemOptimizer( Optim.Options(; f_reltol = 1e-10, x_abstol = 1.5e-8), ) -solution_ml = fit(optimizer, model_ml) +solution_ml = fit(optimizer, model_ml, start_val = start) -@test true_val ≈ solution(solution_ml) atol = 0.05 +@test solution(solution_ml) ≈ true_val atol = 0.05 diff --git a/test/unit_tests/model.jl b/test/unit_tests/model.jl index fbe2a937..87812fba 100644 --- a/test/unit_tests/model.jl +++ b/test/unit_tests/model.jl @@ -68,9 +68,8 @@ end test_vars_api(implied(model), ram_matrices) test_params_api(implied(model), ram_matrices) - @test @inferred(loss(model)) isa SemLoss - semloss = loss(model).functions[1] - @test semloss isa SemML + @test @inferred(sem_term(model)) isa SemLoss + @test sem_term(model) isa losstype @test @inferred(nsamples(model)) == nsamples(obs) end From bab1317c939a3d33bfc44b5e74f7f3fb54a9c697 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 21 Mar 2026 11:17:31 -0700 Subject: [PATCH 05/18] params/param_labels(): use both as synonyms for now --- src/frontend/specification/Sem.jl | 1 + src/implied/abstract.jl | 1 + 2 files changed, 2 insertions(+) diff --git a/src/frontend/specification/Sem.jl b/src/frontend/specification/Sem.jl index 684cfa62..d89606b6 100644 --- a/src/frontend/specification/Sem.jl +++ b/src/frontend/specification/Sem.jl @@ -183,6 +183,7 @@ end ############################################################################################ params(model::AbstractSem) = model.params +param_labels(model::AbstractSem) = params(model) # alias """ loss_terms(model::AbstractSem) diff --git a/src/implied/abstract.jl b/src/implied/abstract.jl index d4868d74..e41e79f6 100644 --- a/src/implied/abstract.jl +++ b/src/implied/abstract.jl @@ -8,6 +8,7 @@ nobserved_vars(implied::SemImplied) = nobserved_vars(implied.ram_matrices) nlatent_vars(implied::SemImplied) = nlatent_vars(implied.ram_matrices) param_labels(implied::SemImplied) = param_labels(implied.ram_matrices) +params(implied::SemImplied) = param_labels(implied) nparams(implied::SemImplied) = nparams(implied.ram_matrices) # checks if the A matrix is acyclic From f7f74520ea765ebea72d5b54098e1ab31f7bda7d Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 21 Mar 2026 11:17:31 -0700 Subject: [PATCH 06/18] check_same_semterm_type(): refactor check_single_lossfun() --- src/additional_functions/helper.jl | 29 --------------- src/frontend/fit/fitmeasures/RMSEA.jl | 2 +- src/frontend/fit/fitmeasures/chi2.jl | 17 ++------- src/frontend/fit/fitmeasures/minus2ll.jl | 2 +- src/frontend/specification/Sem.jl | 47 +++++++++++++++++++++++- 5 files changed, 50 insertions(+), 47 deletions(-) diff --git a/src/additional_functions/helper.jl b/src/additional_functions/helper.jl index b3f5212b..8f2342c3 100644 --- a/src/additional_functions/helper.jl +++ b/src/additional_functions/helper.jl @@ -116,35 +116,6 @@ function nonunique(values::AbstractVector) return res end -# check that a model only has a single lossfun -function check_single_lossfun(model::AbstractSemSingle; throw_error) - if (length(model.loss.functions) > 1) & throw_error - @error "The model has $(length(sem.loss.functions)) loss functions. - Only a single loss function is supported." - end - return isone(length(model.loss.functions)) -end - -# check that all models use the same single loss function -function check_single_lossfun(models::AbstractSemSingle...; throw_error) - uniform = true - lossfun = models[1].loss.functions[1] - L = typeof(lossfun) - for (i, model) in enumerate(models) - uniform &= check_single_lossfun(model; throw_error = throw_error) - cur_lossfun = model.loss.functions[1] - if !isa(cur_lossfun, L) & throw_error - @error "Loss function for group #$i model is $(typeof(cur_lossfun)), expected $L. - Heterogeneous loss functions are not supported." - end - uniform &= isa(cur_lossfun, L) - end - return uniform -end - -check_single_lossfun(model::SemEnsemble; throw_error) = - check_single_lossfun(model.sems...; throw_error) - # scaling corrections for multigroup models mg_correction(::SemFIML) = 0 mg_correction(::SemML) = 0 diff --git a/src/frontend/fit/fitmeasures/RMSEA.jl b/src/frontend/fit/fitmeasures/RMSEA.jl index 7406b74c..ac2d890d 100644 --- a/src/frontend/fit/fitmeasures/RMSEA.jl +++ b/src/frontend/fit/fitmeasures/RMSEA.jl @@ -27,7 +27,7 @@ RMSEA_corr_scale(::Type{<:SemML}) = -1 RMSEA_corr_scale(::Type{<:SemWLS}) = -1 function RMSEA(fit::SemFit, model::AbstractSem) - term_type = check_single_lossfun(model; throw_error = true) + term_type = check_same_semterm_type(model; throw_error = true) n = nsamples(fit) + nsem_terms(model) * RMSEA_corr_scale(term_type) sqrt(nsem_terms(model)) * RMSEA(dof(fit), χ²(fit), n) end diff --git a/src/frontend/fit/fitmeasures/chi2.jl b/src/frontend/fit/fitmeasures/chi2.jl index 22d6c2e2..c56b9a2a 100644 --- a/src/frontend/fit/fitmeasures/chi2.jl +++ b/src/frontend/fit/fitmeasures/chi2.jl @@ -14,21 +14,10 @@ with the *observed* covariance matrix. function χ²(fit::SemFit, model::AbstractSem) terms = sem_terms(model) - isempty(terms) && return 0.0 + @assert !isempty(terms) - term1 = _unwrap(loss(terms[1])) - L = typeof(term1).name - - # check that all SemLoss terms are of the same class (ML, FIML, WLS etc), ignore typeparams - for (i, term) in enumerate(terms) - lossterm = _unwrap(loss(term)) - @assert lossterm isa SemLoss - if typeof(_unwrap(lossterm)).name != L - @error "SemLoss term #$i is $(typeof(_unwrap(lossterm)).name), expected $L. Heterogeneous loss functions are not supported" - end - end - - return χ²(typeof(term1), fit, model) + L = check_same_semterm_type(model; throw_error = true) + return χ²(L, fit, model) end # bollen, p. 115, only correct for GLS weight matrix diff --git a/src/frontend/fit/fitmeasures/minus2ll.jl b/src/frontend/fit/fitmeasures/minus2ll.jl index 3b353f5c..1cdf5c07 100644 --- a/src/frontend/fit/fitmeasures/minus2ll.jl +++ b/src/frontend/fit/fitmeasures/minus2ll.jl @@ -62,6 +62,6 @@ end ############################################################################################ function minus2ll(model::AbstractSem, fit::SemFit) - check_single_lossfun(model; throw_error = true) + check_same_semterm_type(model; throw_error = true) sum(Base.Fix2(minus2ll, fit) ∘ _unwrap ∘ loss, sem_terms(model)) end diff --git a/src/frontend/specification/Sem.jl b/src/frontend/specification/Sem.jl index d89606b6..c3e4bd2b 100644 --- a/src/frontend/specification/Sem.jl +++ b/src/frontend/specification/Sem.jl @@ -51,8 +51,8 @@ end function multigroup_weights(models, n) nsamples_total = sum(nsamples, models) - uniform_lossfun = check_single_lossfun(models...; throw_error = false) - if !uniform_lossfun + semloss_type = check_same_semterm_type(semterms; throw_error = false) + if isnothing(semloss_type) @info """ Your ensemble model contains heterogeneous loss functions. Default weights of (#samples per group/#total samples) will be used @@ -258,6 +258,49 @@ function sem_term(model::AbstractSem, _::Nothing = nothing) error("Unreachable reached") end +# check that all models use the same single loss function +# returns the type of the single SEM loss function, SemLoss if there are multiple different SEM losses, +# nothing if there are no SEM terms. +# If throw_error=true, throws an error if there are multiple different SEM loss functions +check_same_semterm_type(model::AbstractSem; throw_error::Bool = true) = + check_same_semterm_type(sem_terms(model); throw_error = throw_error) + +# check that all models use the same single loss function +# returns the type of the single SEM loss function, +# nothing if there are multiple different SEM losses or no SEM terms. +# If throw_error=true, throws an error if there are multiple different SEM loss functions +function check_same_semterm_type(terms::Tuple; throw_error::Bool = true) + isempty(terms) && return nothing + + _semloss(term::SemLoss) = _unwrap(term) + _semloss(term::LossTerm) = _semloss(loss(term)) + _semloss(term) = throw(ArgumentError("SemLoss term expected, $(typeof(term)) found")) + _semloss_label(i::Integer, _::Union{SemLoss, LossTerm{<:SemLoss, Nothing}}) = "#$i" + _semloss_label(i::Integer, term::LossTerm{<:SemLoss, Symbol}) = "#$i ($(SEM.id(term)))" + + term1 = _semloss(terms[1]) + L = typeof(term1).name + + # check that all SemLoss terms are of the same class (ML, FIML, WLS etc), ignore typeparams + for (i, term) in enumerate(terms) + lossterm = _semloss(term) + @assert lossterm isa SemLoss + if typeof(lossterm).name != L + if throw_error + error( + "SemLoss term $(_semloss_label(i, term)) is $(typeof(lossterm).name), expected $L. Heterogeneous loss functions are not supported", + ) + else + return nothing + end + end + end + + # return the type of the first SEM term + # note that type params of the SEM terms might be different + return typeof(term1) +end + # wrappers arounds a single SemLoss term observed(model::AbstractSem, id::Nothing = nothing) = observed(sem_term(model, id)) implied(model::AbstractSem, id::Nothing = nothing) = implied(sem_term(model, id)) From 961a3c8931941c7c5e0c150b95a74d05119da537 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 21 Mar 2026 11:17:31 -0700 Subject: [PATCH 07/18] update multi-group correction deduplicate the correction scale methods and move to Sem.jl --- src/additional_functions/helper.jl | 5 --- src/frontend/fit/fitmeasures/RMSEA.jl | 7 +--- src/frontend/specification/Sem.jl | 46 +++++++++++++++++---------- 3 files changed, 31 insertions(+), 27 deletions(-) diff --git a/src/additional_functions/helper.jl b/src/additional_functions/helper.jl index 8f2342c3..5442357f 100644 --- a/src/additional_functions/helper.jl +++ b/src/additional_functions/helper.jl @@ -115,8 +115,3 @@ function nonunique(values::AbstractVector) end return res end - -# scaling corrections for multigroup models -mg_correction(::SemFIML) = 0 -mg_correction(::SemML) = 0 -mg_correction(::SemWLS) = -1 diff --git a/src/frontend/fit/fitmeasures/RMSEA.jl b/src/frontend/fit/fitmeasures/RMSEA.jl index ac2d890d..9d33e47e 100644 --- a/src/frontend/fit/fitmeasures/RMSEA.jl +++ b/src/frontend/fit/fitmeasures/RMSEA.jl @@ -21,14 +21,9 @@ For multigroup models, the correction proposed by J.H. Steiger is applied """ RMSEA(fit::SemFit) = RMSEA(fit, fit.model) -# scaling corrections -RMSEA_corr_scale(::Type{<:SemFIML}) = 0 -RMSEA_corr_scale(::Type{<:SemML}) = -1 -RMSEA_corr_scale(::Type{<:SemWLS}) = -1 - function RMSEA(fit::SemFit, model::AbstractSem) term_type = check_same_semterm_type(model; throw_error = true) - n = nsamples(fit) + nsem_terms(model) * RMSEA_corr_scale(term_type) + n = nsamples(fit) + nsem_terms(model) * multigroup_correction_scale(term_type) sqrt(nsem_terms(model)) * RMSEA(dof(fit), χ²(fit), n) end diff --git a/src/frontend/specification/Sem.jl b/src/frontend/specification/Sem.jl index c3e4bd2b..d8696e82 100644 --- a/src/frontend/specification/Sem.jl +++ b/src/frontend/specification/Sem.jl @@ -45,32 +45,46 @@ function Base.show(io::IO, term::LossTerm) end end -############################################################################################ -# constructor for Sem types -############################################################################################ +# scaling corrections for multigroup models + +# fallback method for non-standard SemLoss type +multigroup_correction_scale(::Type{<:SemLoss}) = nothing -function multigroup_weights(models, n) - nsamples_total = sum(nsamples, models) +multigroup_correction_scale(::Type{<:SemFIML}) = 0 +multigroup_correction_scale(::Type{<:SemML}) = 0 +multigroup_correction_scale(::Type{<:SemWLS}) = -1 + +multigroup_correction_scale(loss::SemLoss) = multigroup_correction_scale(typeof(loss)) + +# calculate sem term weights for multigroup models +# correcting for the number of samples and the loss type +function multigroup_weights(semterms...) + n = length(semterms) + nsamples_total = sum(nsamples, semterms) semloss_type = check_same_semterm_type(semterms; throw_error = false) if isnothing(semloss_type) @info """ Your ensemble model contains heterogeneous loss functions. Default weights of (#samples per group/#total samples) will be used """ - return [(nsamples(model)) / (nsamples_total) for model in models] - end - lossfun = models[1].loss.functions[1] - if !applicable(mg_correction, lossfun) - @info """ - We don't know how to choose group weights for the specified loss function. - Default weights of (#samples per group/#total samples) will be used - """ - return [(nsamples(model)) / (nsamples_total) for model in models] + c = 0 + else + c = multigroup_correction_scale(semloss_type) + if isnothing(c) + @info """ + We don't know how to choose group weights for the specified loss function. + Default weights of (#samples per group/#total samples) will be used + """ + c = 0 + end end - c = mg_correction(lossfun) - return [(nsamples(model)+c) / (nsamples_total+n*c) for model in models] + return [(nsamples(term)+c) / (nsamples_total+n*c) for term in semterms] end +############################################################################################ +# constructor for Sem types +############################################################################################ + function Sem( loss_terms...; params::Union{Vector{Symbol}, Nothing} = nothing, From a9ee00bdafcbd63df536eeaa929ebb4a218fe6f3 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 21 Mar 2026 11:18:19 -0700 Subject: [PATCH 08/18] replace_observed(): simplify & refactor remove update_observed!() --- ext/SEMNLOptExt/NLopt.jl | 7 -- ext/SEMProximalOptExt/ProximalAlgorithms.jl | 7 -- src/StructuralEquationModels.jl | 1 - src/additional_functions/simulation.jl | 115 ------------------- src/frontend/specification/Sem.jl | 116 +++++++++++++++++++- src/implied/RAM/generic.jl | 17 --- src/implied/RAM/symbolic.jl | 20 ---- src/implied/empty.jl | 6 - src/loss/ML/FIML.jl | 8 -- src/loss/ML/ML.jl | 20 ---- src/loss/WLS/WLS.jl | 23 ---- src/loss/abstract.jl | 29 +++++ src/loss/constant/constant.jl | 6 - src/loss/regularization/ridge.jl | 6 - src/optimizer/Empty.jl | 6 - src/optimizer/optim.jl | 6 - 16 files changed, 140 insertions(+), 253 deletions(-) diff --git a/ext/SEMNLOptExt/NLopt.jl b/ext/SEMNLOptExt/NLopt.jl index 90004b90..87601030 100644 --- a/ext/SEMNLOptExt/NLopt.jl +++ b/ext/SEMNLOptExt/NLopt.jl @@ -107,13 +107,6 @@ function SemOptimizerNLopt(; ) end -############################################################################################ -### Recommended methods -############################################################################################ - -SEM.update_observed(optimizer::SemOptimizerNLopt, observed::SemObserved; kwargs...) = - optimizer - ############################################################################################ ### additional methods ############################################################################################ diff --git a/ext/SEMProximalOptExt/ProximalAlgorithms.jl b/ext/SEMProximalOptExt/ProximalAlgorithms.jl index 3ec32453..0937ee04 100644 --- a/ext/SEMProximalOptExt/ProximalAlgorithms.jl +++ b/ext/SEMProximalOptExt/ProximalAlgorithms.jl @@ -34,13 +34,6 @@ SemOptimizerProximal(; SEM.sem_optimizer_subtype(::Val{:Proximal}) = SemOptimizerProximal -############################################################################################ -### Recommended methods -############################################################################################ - -SEM.update_observed(optimizer::SemOptimizerProximal, observed::SemObserved; kwargs...) = - optimizer - ############################################################################ ### Model fitting ############################################################################ diff --git a/src/StructuralEquationModels.jl b/src/StructuralEquationModels.jl index d98e7925..0dbcd16a 100644 --- a/src/StructuralEquationModels.jl +++ b/src/StructuralEquationModels.jl @@ -205,7 +205,6 @@ export AbstractSem, z_test!, example_data, replace_observed, - update_observed, @StenoGraph, →, ←, diff --git a/src/additional_functions/simulation.jl b/src/additional_functions/simulation.jl index 6d694c97..e85e9d5c 100644 --- a/src/additional_functions/simulation.jl +++ b/src/additional_functions/simulation.jl @@ -1,118 +1,3 @@ -""" - (1) replace_observed(model::AbstractSemSingle; kwargs...) - - (2) replace_observed(model::AbstractSemSingle, observed; kwargs...) - - (3) replace_observed(model::SemEnsemble; column = :group, weights = nothing, kwargs...) - -Return a new model with swaped observed part. - -# Arguments -- `model::AbstractSemSingle`: model to swap the observed part of. -- `kwargs`: additional keyword arguments; typically includes `data` and `specification` -- `observed`: Either an object of subtype of `SemObserved` or a subtype of `SemObserved` - -# For SemEnsemble models: -- `column`: if a DataFrame is passed as `data = ...`, which column signifies the group? -- `weights`: how to weight the different sub-models, - defaults to number of samples per group in the new data -- `kwargs`: has to be a dict with keys equal to the group names. - For `data` can also be a DataFrame with `column` containing the group information, - and for `specification` can also be an `EnsembleParameterTable`. - -# Examples -See the online documentation on [Replace observed data](@ref). -""" -function replace_observed end - -""" - update_observed(to_update, observed::SemObserved; kwargs...) - -Update a `SemImplied`, `SemLossFunction` or `SemOptimizer` object to use a `SemObserved` object. - -# Examples -See the online documentation on [Replace observed data](@ref). - -# Implementation -You can provide a method for this function when defining a new type, for more information -on this see the online developer documentation on [Update observed data](@ref). -""" -function update_observed end - -############################################################################################ -# change observed (data) without reconstructing the whole model -############################################################################################ - -# don't change non-SEM terms -replace_observed(loss::AbstractLoss; kwargs...) = loss - -# use the same observed type as before -replace_observed(loss::SemLoss; kwargs...) = - replace_observed(loss, typeof(SEM.observed(loss)).name.wrapper; kwargs...) - -# construct a new observed type -replace_observed(loss::SemLoss, observed_type; kwargs...) = - replace_observed(loss, observed_type(; kwargs...); kwargs...) - -function replace_observed(loss::SemLoss, new_observed::SemObserved; kwargs...) - kwargs = Dict{Symbol, Any}(kwargs...) - old_observed = SEM.observed(loss) - implied = SEM.implied(loss) - - # get field types - kwargs[:observed_type] = typeof(new_observed) - kwargs[:old_observed_type] = typeof(old_observed) - - # update implied - new_implied = update_observed(implied, new_observed; kwargs...) - kwargs[:implied] = new_implied - kwargs[:implied_type] = typeof(new_implied) - kwargs[:nparams] = nparams(new_implied) - - # update loss - return update_observed(loss, new_observed; kwargs...) -end - -replace_observed(loss::LossTerm; kwargs...) = - LossTerm(replace_observed(loss.loss; kwargs...), loss.id, loss.weight) - -function replace_observed(sem::Sem; kwargs...) - updated_terms = Tuple(replace_observed(term; kwargs...) for term in loss_terms(sem)) - return Sem(updated_terms...) -end - -function replace_observed( - emodel::SemEnsemble; - column = :group, - weights = nothing, - kwargs..., -) - kwargs = Dict{Symbol, Any}(kwargs...) - # allow for EnsembleParameterTable to be passed as specification - if haskey(kwargs, :specification) && isa(kwargs[:specification], EnsembleParameterTable) - kwargs[:specification] = convert(Dict{Symbol, RAMMatrices}, kwargs[:specification]) - end - # allow for DataFrame with group variable "column" to be passed as new data - if haskey(kwargs, :data) && isa(kwargs[:data], DataFrame) - kwargs[:data] = Dict( - group => - select(filter(r -> r[column] == group, kwargs[:data]), Not(column)) for - group in emodel.groups - ) - end - # update each model for new data - models = emodel.sems - new_models = Tuple( - replace_observed(m; group_kwargs(g, kwargs)...) for - (m, g) in zip(models, emodel.groups) - ) - return SemEnsemble(new_models...; weights = weights, groups = emodel.groups) -end - -function group_kwargs(g, kwargs) - return Dict(k => kwargs[k][g] for k in keys(kwargs)) -end - ############################################################################################ # simulate data ############################################################################################ diff --git a/src/frontend/specification/Sem.jl b/src/frontend/specification/Sem.jl index d8696e82..42ff2d3e 100644 --- a/src/frontend/specification/Sem.jl +++ b/src/frontend/specification/Sem.jl @@ -436,12 +436,118 @@ function build_SemTerms(loss, observed, implied; kwargs...) end end -function update_observed(sem::Sem, new_observed; kwargs...) - new_terms = Tuple( - update_observed(lossterm.loss, new_observed; kwargs...) for - lossterm in loss_terms(sem) +############################################################## +# replace_observed: Sem level +############################################################## + +""" + replace_observed(model::Sem, observed::SemObserved) + replace_observed(model::Sem, data::AbstractDict{Symbol}) + replace_observed(model::Sem, data::AbstractDataFrame; [semterm_column]) + replace_observed(loss::SemLoss, observed::SemObserved) + replace_observed(loss::SemLoss, data::Union{AbstractMatrix, DataFrame}) + +Construct a new SEM model or SEM loss with replaced observed data. + +The SEM structure (implied covariance, loss type) is preserved; +only the observed data is swapped. + +# Single-term models + +Pass a `SemObserved` object, a data matrix, or a `DataFrame`: +```julia +replace_observed(model, new_data_matrix) +replace_observed(model, new_sem_observed) +replace_observed(model, new_df) +``` + +# Multi-term models + +Pass a `Dict{Symbol}` mapping term ids to data or `SemObserved` objects: +```julia +replace_observed(model, Dict(:g1 => data1, :g2 => data2)) +``` + +Or pass a `DataFrame` with a `semterm_column` identifying the group: +```julia +replace_observed(model, new_df; semterm_column = :group) +``` +""" +function replace_observed end + +function replace_observed(sem::Sem, data::Union{SemObserved, AbstractMatrix}) + nsem_terms(sem) > 1 && throw( + ArgumentError( + "Model contains $(nsem_terms(sem)) SEM terms. " * + "Use a Dict{Symbol} or a DataFrame with `semterm_column` to provide per-term data.", + ), + ) + updated_terms = Tuple(replace_observed(term, data) for term in loss_terms(sem)) + return Sem(updated_terms...) +end + +function replace_observed(sem::Sem, data::AbstractDict{Symbol}) + term_ids = Set( + if !isnothing(id(term)) + id(term) + else + "Multigroup replace_observed(sem, data::Dict) requires all SEM terms to have ids." |> + ArgumentError |> + throw + end for term in loss_terms(sem) if issemloss(term) + ) + # check for extra ids + extra_term_ids = setdiff(keys(data), term_ids) + isempty(extra_term_ids) || + @warn "Ignoring data with ids=$(collect(extra_term_ids)): no such SEM terms exist in the model" + + updated_terms = map(loss_terms(sem)) do term + issemloss(term) || return term + tid = id(term) + term_data = get(data, tid, nothing) + isnothing(term_data) && + throw(ArgumentError("No data provided for SEM term :$tid")) + return replace_observed(term, term_data) + end + return Sem(Tuple(updated_terms)...) +end + +function replace_observed(sem::Sem, data::AbstractVector) + nsem = nsem_terms(sem) + nsem == length(data) || throw( + ArgumentError( + "Length of data ($(length(data))) does not match number of SEM terms ($nsem)", + ), + ) + updated_terms = map(enumerate(loss_terms(sem))) do (i, term) + issemloss(term) ? replace_observed(term, data[i]) : term + end + return Sem(Tuple(updated_terms)...) +end + +function replace_observed( + sem::Sem, + data::AbstractDataFrame; + semterm_column::Union{Symbol, Nothing} = nothing, +) + if isnothing(semterm_column) + # single-term shortcut + nsem_terms(sem) > 1 && throw( + ArgumentError( + "Model contains $(nsem_terms(sem)) SEM terms. " * + "Provide `semterm_column` to specify which DataFrame column identifies the groups.", + ), + ) + updated_terms = Tuple(replace_observed(term, data) for term in loss_terms(sem)) + return Sem(updated_terms...) + end + + # multi-term: split DataFrame by semterm_column + terms_data = Dict( + g[semterm_column] => group_data for + (g, group_data) in pairs(groupby(data, semterm_column)) ) - return Sem(new_terms...) + return replace_observed(sem, terms_data) end ############################################################## diff --git a/src/implied/RAM/generic.jl b/src/implied/RAM/generic.jl index f1c1e08d..1569b341 100644 --- a/src/implied/RAM/generic.jl +++ b/src/implied/RAM/generic.jl @@ -179,20 +179,3 @@ function update!(targets::EvaluationTargets, implied::RAM, params) mul!(implied.μ, implied.F⨉I_A⁻¹, implied.M) end end - -############################################################################################ -### Recommended methods -############################################################################################ - -function update_observed(implied::RAM, observed::SemObserved; kwargs...) - if nobserved_vars(observed) == nobserved_vars(implied) - return implied - else - return RAM(; - observed = observed, - gradient_required = !isnothing(implied.∇A), - meanstructure = MeanStruct(implied) == HasMeanStruct, - kwargs..., - ) - end -end diff --git a/src/implied/RAM/symbolic.jl b/src/implied/RAM/symbolic.jl index 4c9bda91..52a192e6 100644 --- a/src/implied/RAM/symbolic.jl +++ b/src/implied/RAM/symbolic.jl @@ -190,26 +190,6 @@ function update!(targets::EvaluationTargets, implied::RAMSymbolic, par) end end -############################################################################################ -### Recommended methods -############################################################################################ - -function update_observed(implied::RAMSymbolic, observed::SemObserved; kwargs...) - if nobserved_vars(observed) == nobserved_vars(implied) - return implied - else - return RAMSymbolic(; - observed = observed, - vech = implied.Σ isa Vector, - gradient = !isnothing(implied.∇Σ), - hessian = !isnothing(implied.∇²Σ), - meanstructure = MeanStruct(implied) == HasMeanStruct, - approximate_hessian = isnothing(implied.∇²Σ), - kwargs..., - ) - end -end - ############################################################################################ ### additional functions ############################################################################################ diff --git a/src/implied/empty.jl b/src/implied/empty.jl index a327ee13..a650a07a 100644 --- a/src/implied/empty.jl +++ b/src/implied/empty.jl @@ -46,9 +46,3 @@ end ############################################################################################ update!(targets::EvaluationTargets, implied::ImpliedEmpty, par) = nothing - -############################################################################################ -### Recommended methods -############################################################################################ - -update_observed(implied::ImpliedEmpty, observed::SemObserved; kwargs...) = implied diff --git a/src/loss/ML/FIML.jl b/src/loss/ML/FIML.jl index fdedf398..15081e20 100644 --- a/src/loss/ML/FIML.jl +++ b/src/loss/ML/FIML.jl @@ -159,14 +159,6 @@ function evaluate!(objective, gradient, hessian, loss::SemFIML, params) return objective end - -############################################################################################ -### Recommended methods -############################################################################################ - -update_observed(loss::SemFIML, observed::SemObserved; kwargs...) = - SemFIML(; observed = observed, kwargs...) - ############################################################################################ ### additional functions ############################################################################################ diff --git a/src/loss/ML/ML.jl b/src/loss/ML/ML.jl index 2d449d73..9f327544 100644 --- a/src/loss/ML/ML.jl +++ b/src/loss/ML/ML.jl @@ -235,23 +235,3 @@ function non_posdef_return(par) return typemax(eltype(par)) end end - -############################################################################################ -### recommended methods -############################################################################################ - -update_observed(loss::SemML, observed::SemObservedMissing; kwargs...) = - error("ML estimation does not work with missing data - use FIML instead") - -function update_observed(loss::SemML, observed::SemObserved; kwargs...) - if (obs_cov(loss) == obs_cov(observed)) && (obs_mean(loss) == obs_mean(observed)) - return loss # no change - else - return SemML( - observed, - loss.implied; - approximate_hessian = HessianEval(loss) == ApproxHessian, - kwargs..., - ) - end -end diff --git a/src/loss/WLS/WLS.jl b/src/loss/WLS/WLS.jl index 5c4cb252..8f4a109c 100644 --- a/src/loss/WLS/WLS.jl +++ b/src/loss/WLS/WLS.jl @@ -179,26 +179,3 @@ function evaluate!(objective, gradient, hessian, loss::SemWLS, par) return objective end - -############################################################################################ -### Recommended methods -############################################################################################ - -function update_observed( - loss::SemWLS, - observed::SemObserved; - recompute_V = true, - kwargs..., -) - if recompute_V - return SemWLS(observed, loss.implied; kwargs...) - else - return SemWLS( - observed, - loss.implied; - wls_weight_matrix = loss.V, - wls_weight_matrix_mean = loss.V_μ, - kwargs..., - ) - end -end diff --git a/src/loss/abstract.jl b/src/loss/abstract.jl index bf8585d6..bcd6d62b 100644 --- a/src/loss/abstract.jl +++ b/src/loss/abstract.jl @@ -40,3 +40,32 @@ function check_observed_vars(observed::SemObserved, implied::SemImplied) end check_observed_vars(sem::SemLoss) = check_observed_vars(observed(sem), implied(sem)) + +############################################################################################ +# replace_observed: SemLoss, AbstractLoss, LossTerm +############################################################################################ + +function replace_observed(loss::SemLoss, new_observed::SemObserved) + old_obs = SEM.observed(loss) + observed_vars(old_obs) == observed_vars(new_observed) || throw( + ArgumentError( + "observed_vars of the new data do not match the model: " * + "expected $(observed_vars(old_obs)), got $(observed_vars(new_observed))", + ), + ) + return typeof(loss).name.wrapper(new_observed, SEM.implied(loss)) +end + +function replace_observed(loss::SemLoss, data::Union{AbstractMatrix, DataFrame}) + old_obs = SEM.observed(loss) + new_observed = + typeof(old_obs).name.wrapper(data = data, observed_vars = observed_vars(old_obs)) + return replace_observed(loss, new_observed) +end + +# non-SEM loss terms are unchanged +replace_observed(loss::AbstractLoss, ::Any) = loss + +# LossTerm: delegate to inner loss +replace_observed(term::LossTerm, data) = + LossTerm(replace_observed(loss(term), data), id(term), weight(term)) diff --git a/src/loss/constant/constant.jl b/src/loss/constant/constant.jl index 023076cc..2aff0156 100644 --- a/src/loss/constant/constant.jl +++ b/src/loss/constant/constant.jl @@ -35,9 +35,3 @@ SemConstant(; constant_loss::Number, kwargs...) = SemConstant(constant_loss) objective(loss::SemConstant, par) = convert(eltype(par), loss.c) gradient(loss::SemConstant, par) = zero(par) hessian(loss::SemConstant, par) = zeros(eltype(par), length(par), length(par)) - -############################################################################################ -### Recommended methods -############################################################################################ - -update_observed(loss::SemConstant, observed::SemObserved; kwargs...) = loss diff --git a/src/loss/regularization/ridge.jl b/src/loss/regularization/ridge.jl index 3e2cfbff..813aff11 100644 --- a/src/loss/regularization/ridge.jl +++ b/src/loss/regularization/ridge.jl @@ -85,9 +85,3 @@ function hessian(ridge::SemRidge, par) @views @. ridge.hessian[ridge.which_H] .= 2 * ridge.α return ridge.hessian end - -############################################################################################ -### Recommended methods -############################################################################################ - -update_observed(loss::SemRidge, observed::SemObserved; kwargs...) = loss diff --git a/src/optimizer/Empty.jl b/src/optimizer/Empty.jl index f95c067c..fd36acb5 100644 --- a/src/optimizer/Empty.jl +++ b/src/optimizer/Empty.jl @@ -11,12 +11,6 @@ struct SemOptimizerEmpty <: SemOptimizer{:Empty} end sem_optimizer_subtype(::Val{:Empty}) = SemOptimizerEmpty -############################################################################################ -### Recommended methods -############################################################################################ - -update_observed(optimizer::SemOptimizerEmpty, observed::SemObserved; kwargs...) = optimizer - ############################################################################################ ### Pretty Printing ############################################################################################ diff --git a/src/optimizer/optim.jl b/src/optimizer/optim.jl index 70413193..a0aae22a 100644 --- a/src/optimizer/optim.jl +++ b/src/optimizer/optim.jl @@ -57,12 +57,6 @@ SemOptimizerOptim(; sem_optimizer_subtype(::Val{:Optim}) = SemOptimizerOptim -############################################################################################ -### Recommended methods -############################################################################################ - -update_observed(optimizer::SemOptimizerOptim, observed::SemObserved; kwargs...) = optimizer - ############################################################################################ ### additional methods ############################################################################################ From 84c66530d59da0c645b9fed8aa1578f1747f0795 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 21 Mar 2026 17:56:25 -0700 Subject: [PATCH 09/18] bootstrap: sync with Sem updates --- src/frontend/fit/standard_errors/bootstrap.jl | 348 +++++++++--------- test/examples/helper.jl | 27 +- test/examples/multigroup/build_models.jl | 4 +- .../political_democracy/constructor.jl | 14 +- 4 files changed, 199 insertions(+), 194 deletions(-) diff --git a/src/frontend/fit/standard_errors/bootstrap.jl b/src/frontend/fit/standard_errors/bootstrap.jl index 845a209e..ce84e923 100644 --- a/src/frontend/fit/standard_errors/bootstrap.jl +++ b/src/frontend/fit/standard_errors/bootstrap.jl @@ -1,35 +1,118 @@ +# base type for accumulators of intermediate bootstrap results +abstract type BootstrapAccumulator end + +# internal function to run bootstrap +function bootstrap!( + acc::BootstrapAccumulator, + fitted::SemFit; + data = nothing, + engine = :Optim, + parallel = false, + fit_kwargs = Dict(), +) + sem = model(fitted) + data = isnothing(data) ? _bootstrap_data(sem) : data + start = solution(fitted) + + n_boot = n_bootstrap(acc) + + # fit to bootstrap samples + if !parallel + for i in 1:n_boot + new_fit = _fit_bootstrap_sample(sem, data, start; engine, fit_kwargs) + update!(acc, i, new_fit, nothing) + end + else + n_threads = Threads.nthreads() + # Pre-create one independent model copy per thread via deepcopy. + model_pool = Channel(n_threads) + for _ in 1:n_threads + put!(model_pool, deepcopy(sem)) + end + lk = ReentrantLock() + Threads.@threads for i in 1:n_boot + thread_model = take!(model_pool) + new_fit = _fit_bootstrap_sample(thread_model, data, start; engine, fit_kwargs) + update!(acc, i, new_fit, lk) + put!(model_pool, thread_model) + end + end + + return acc +end + +# a simple accumulator that just stores the statistic for each sample and whether it converged +struct SimpleBootstrapAccumulator{F} <: BootstrapAccumulator + statistic::F + samples::Vector{Any} + converged_mask::Vector{Bool} +end + +SimpleBootstrapAccumulator(statistic, n_boot::Integer) = + SimpleBootstrapAccumulator(statistic, Vector{Any}(undef, n_boot), fill(false, n_boot)) + +n_bootstrap(acc::SimpleBootstrapAccumulator) = length(acc.samples) + +function update!(acc::SimpleBootstrapAccumulator, i::Integer, fit::SemFit, _) + acc.samples[i] = acc.statistic(fit) + acc.converged_mask[i] = converged(fit) +end + +""" + struct BootstrapResult{T} + +Stores the output of a [`bootstrap`](@ref) call. +""" +struct BootstrapResult{T} + samples::Vector{T} + converged_mask::BitVector + n_boot::Int + n_converged::Int +end + +function Base.show(io::IO, result::BootstrapResult{T}) where {T} + println( + io, + "BootstrapResult{$(T)} with $(result.n_converged)/$(result.n_boot) converged samples", + ) +end + """ bootstrap( - fitted::SemFit, - specification::SemSpecification; + fitted::SemFit; statistic = solution, n_boot = 3000, data = nothing, engine = :Optim, parallel = false, - fit_kwargs = Dict(), - replace_kwargs = Dict()) + fit_kwargs = Dict() + ) -> BootstrapResult + +Bootstrap the samples and apply `statistic` function to each. -Return bootstrap samples for `statistic`. +Returns a [`BootstrapResult`](@ref) object containing the results of `statistic` +applied to each bootstrapped sample. + +Supports both single-group and multi-group models. +For multi-group models, each group is resampled independently. # Arguments - `fitted`: a fitted SEM. -- `specification`: a `ParameterTable` or `RAMMatrices` object passed to `replace_observed`. - `statistic`: any function that can be called on a `SemFit` object. The output will be returned as the bootstrap sample. - `n_boot`: number of boostrap samples -- `data`: data to sample from. Only needed if different than the data from `sem_fit` +- `data`: data to sample from. Only needed if different than the fitted model. + For multi-group models, pass a `Dict{Symbol}` mapping term ids to data matrices. - `engine`: optimizer engine, passed to `fit`. - `parallel`: if `true`, run bootstrap samples in parallel on all available threads. The number of threads is controlled by the `JULIA_NUM_THREADS` environment variable or the `--threads` flag when starting Julia. - `fit_kwargs` : a `Dict` controlling model fitting for each bootstrap sample, - passed to `fit` -- `replace_kwargs`: a `Dict` passed to `replace_observed` + passed to [`fit`](@ref) # Example ```julia -# 1000 boostrap samples of the minimum, fitted with :Optim +# 1000 bootstrap samples of the minimum, fitted with :Optim bootstrap( fitted; statistic = StructuralEquationModels.minimum, @@ -40,95 +123,74 @@ bootstrap( """ function bootstrap( fitted::SemFit, - specification::SemSpecification; - statistic = solution, + statistic = solution; n_boot = 3000, data = nothing, engine = :Optim, parallel = false, fit_kwargs = Dict(), - replace_kwargs = Dict(), ) - # access data and convert to matrix - data = prepare_data_bootstrap(data, fitted.model) - start = solution(fitted) - # pre-allocations - out = Vector{Any}(nothing, n_boot) - conv = fill(false, n_boot) - # fit to bootstrap samples - if !parallel - for i in 1:n_boot - sample_data = bootstrap_sample(data) - new_model = replace_observed( - fitted.model; - data = sample_data, - specification = specification, - replace_kwargs..., - ) - new_fit = fit(new_model; start_val = start, engine = engine, fit_kwargs...) - sample = statistic(new_fit) - c = converged(new_fit) - out[i] = sample - conv[i] = c - end - else - n_threads = Threads.nthreads() - # Pre-create one independent model copy per thread via deepcopy. - model_pool = Channel(n_threads) - for _ in 1:n_threads - put!(model_pool, deepcopy(fitted.model)) - end - # fit models in parallel - lk = ReentrantLock() - Threads.@threads for i in 1:n_boot - thread_model = take!(model_pool) - sample_data = bootstrap_sample(data) - new_model = replace_observed( - thread_model; - data = sample_data, - specification = specification, - replace_kwargs..., - ) - new_fit = fit(new_model; start_val = start, engine = engine, fit_kwargs...) - sample = statistic(new_fit) - c = converged(new_fit) - out[i] = sample - conv[i] = c - put!(model_pool, thread_model) - end - end - return Dict( - :samples => collect(a for a in out), - :n_boot => n_boot, - :n_converged => sum(conv), - :converged => conv, + acc = SimpleBootstrapAccumulator(statistic, n_boot) + bootstrap!(acc, fitted; data, engine, parallel, fit_kwargs) + return BootstrapResult( + [s for s in acc.samples], + convert(BitVector, acc.converged_mask), + n_bootstrap(acc), + sum(acc.converged_mask), ) end +# bootstrap accumulator for se_bootstrap() +# accumulates per-parameter sum and sum of squares across bootstrap samples +struct StdErrBootstrapAccumulator <: BootstrapAccumulator + n_boot::Int + sum::Vector{Float64} + squared_sum::Vector{Float64} + n_converged::Ref{Int} +end + +n_bootstrap(acc::StdErrBootstrapAccumulator) = acc.n_boot + +StdErrBootstrapAccumulator(n_params::Integer, n_boot::Integer) = + StdErrBootstrapAccumulator(n_boot, zeros(n_params), zeros(n_params), Ref(0)) + +function update!( + acc::StdErrBootstrapAccumulator, + i::Integer, + fit::SemFit, + lk::Union{Base.AbstractLock, Nothing}, +) + conv = converged(fit) + if conv + sol = solution(fit) + isnothing(lk) || lock(lk) + acc.n_converged[] += 1 + @. acc.sum += sol + @. acc.squared_sum += abs2(sol) + isnothing(lk) || unlock(lk) + end +end + """ - se_bootstrap( - fitted::SemFit, - specification::SemSpecification; - n_boot = 3000, - data = nothing, - parallel = false, - fit_kwargs = Dict(), - replace_kwargs = Dict()) + se_bootstrap(fitted::SemFit; n_boot = 3000, kwargs...) -Return bootstrap standard errors. +Calculate standard errors using bootstrap approach. + +Supports both single-group and multi-group models. +For multi-group models, each group is resampled independently. # Arguments - `fitted`: a fitted SEM. -- `specification`: a `ParameterTable` or `RAMMatrices` object passed to `replace_observed`. - `n_boot`: number of boostrap samples -- `data`: data to sample from. Only needed if different than the data from `sem_fit` +- `data`: data to sample from. Only needed if different than the fitted model. + For multi-group models, pass a `Dict{Symbol}` mapping term ids to data matrices. - `engine`: optimizer engine, passed to `fit`. - `parallel`: if `true`, run bootstrap samples in parallel on all available threads. The number of threads is controlled by the `JULIA_NUM_THREADS` environment variable or the `--threads` flag when starting Julia. -- `fit_kwargs` : a `Dict` controlling model fitting for each bootstrap sample, - passed to `sem_fit` -- `replace_kwargs`: a `Dict` passed to `replace_observed` +- `fit_kwargs` : a `Dict` controlling model fitting for each bootstrap sample, + passed to [`fit`](@ref) + # Example ```julia @@ -142,109 +204,53 @@ se_bootstrap( ) ``` """ -function se_bootstrap( - fitted::SemFit, - specification::SemSpecification; - n_boot = 3000, - data = nothing, - engine = :Optim, - parallel = false, - fit_kwargs = Dict(), - replace_kwargs = Dict(), -) - # access data and convert to matrix - data = prepare_data_bootstrap(data, fitted.model) - start = solution(fitted) - # pre-allocations - total_sum = zero(start) - total_squared_sum = zero(start) - n_conv = Ref(0) - # fit to bootstrap samples - if !parallel - for _ in 1:n_boot - sample_data = bootstrap_sample(data) - new_model = replace_observed( - fitted.model; - data = sample_data, - specification = specification, - replace_kwargs..., - ) - new_fit = fit(new_model; start_val = start, engine = engine, fit_kwargs...) - sol = solution(new_fit) - conv = converged(new_fit) - if conv - n_conv[] += 1 - @. total_sum += sol - @. total_squared_sum += sol^2 - end - end +function se_bootstrap(fitted::SemFit; n_boot = 3000, kwargs...) + acc = StdErrBootstrapAccumulator(nparams(fitted), n_boot) + bootstrap!(acc, fitted; kwargs...) + n_conv = acc.n_converged[] + @info "$n_conv models converged" + + if n_conv == 0 + @warn "No bootstrap samples converged. Returning NaN." + return fill(NaN, length(acc.sum)) else - n_threads = Threads.nthreads() - # Pre-create one independent model copy per thread via deepcopy. - model_pool = Channel(n_threads) - for _ in 1:n_threads - put!(model_pool, deepcopy(fitted.model)) - end - # fit models in parallel - lk = ReentrantLock() - Threads.@threads for _ in 1:n_boot - thread_model = take!(model_pool) - sample_data = bootstrap_sample(data) - new_model = replace_observed( - thread_model; - data = sample_data, - specification = specification, - replace_kwargs..., - ) - new_fit = fit(new_model; start_val = start, engine = engine, fit_kwargs...) - sol = solution(new_fit) - conv = converged(new_fit) - if conv - lock(lk) do - n_conv[] += 1 - @. total_sum += sol - @. total_squared_sum += sol^2 - end - end - put!(model_pool, thread_model) - end + return sqrt.(acc.squared_sum ./ n_conv - abs2.(acc.sum / n_conv)) end - # compute parameters - n_conv = n_conv[] - sd = sqrt.(total_squared_sum / n_conv - (total_sum / n_conv) .^ 2) - @info string(n_conv)*" models converged" - return sd end ############################################################################################ ### Helper Functions ############################################################################################ -function bootstrap_sample(data::Matrix) - nobs = size(data, 1) - index_new = rand(1:nobs, nobs) - data_new = data[index_new, :] - return data_new -end +""" + resample_with_replacement(data::AbstractMatrix) + resample_with_replacement(data::AbstractVector{<:AbstractMatrix}) -bootstrap_sample(data::Dict) = Dict(k => bootstrap_sample(data[k]) for k in keys(data)) +Resample rows of a data matrix with replacement (bootstrap sample). +For a vector of matrices (multi-group models), independently resamples each matrix. +""" +function resample_with_replacement(data::AbstractMatrix) + n = size(data, 1) + return data[rand(1:n, n), :] +end -function prepare_data_bootstrap(data, model::AbstractSemSingle) - if isnothing(data) - data = samples(observed(model)) - end - data = Matrix(data) - return data +function resample_with_replacement(data::AbstractVector{<:AbstractMatrix}) + return [resample_with_replacement(term_data) for term_data in data] end -function prepare_data_bootstrap(data, model::SemEnsemble) - sems = model.sems - groups = model.groups - if isnothing(data) - data = Dict(g => samples(observed(m)) for (g, m) in zip(groups, sems)) +# Extract data from a model for bootstrap resampling. +function _bootstrap_data(sem::AbstractSem) + terms = sem_terms(sem) + if length(terms) == 1 + return samples(observed(loss(terms[1]))) + else + return [samples(observed(loss(term))) for term in terms] end - data = Dict(k => Matrix(data[k]) for k in keys(data)) - return data end - +# Fit one bootstrap replicate: resample, replace observed data, fit. +function _fit_bootstrap_sample(sem_model, data, start; engine, fit_kwargs) + boot_data = resample_with_replacement(data) + boot_model = replace_observed(sem_model, boot_data) + return fit(boot_model; start_val = start, engine = engine, fit_kwargs...) +end diff --git a/test/examples/helper.jl b/test/examples/helper.jl index c4191fdb..f14fec62 100644 --- a/test/examples/helper.jl +++ b/test/examples/helper.jl @@ -138,17 +138,17 @@ function test_estimates( end function test_bootstrap( - model_fit, - spec; + model_fit::SemFit; compare_hessian = true, rtol_hessian = 0.2, compare_bs = true, rtol_bs = 0.1, n_boot = 500, + seed = 32432, ) - @testset rng = Random.seed!(32432) "bootstrap" begin - se_bs = @suppress se_bootstrap(model_fit, spec; n_boot = n_boot) - # hessian and bootstrap se are close + @testset rng = Random.seed!(seed) "bootstrap" begin + se_bs = @suppress se_bootstrap(model_fit; n_boot = n_boot) + # hessian-based and bootstrap-based std.errors are close if compare_hessian se_he = @suppress se_hessian(model_fit) #println(maximum(abs.(se_he - se_bs))) @@ -156,10 +156,9 @@ function test_bootstrap( end # se_bootstrap and bootstrap |> se are close if compare_bs - bs_samples = bootstrap(model_fit, spec; n_boot = n_boot) - @test bs_samples[:n_converged] >= 0.95*n_boot - bs_samples = - cat(bs_samples[:samples][BitVector(bs_samples[:converged])]..., dims = 2) + bs_samples = bootstrap(model_fit; n_boot = n_boot) + @test bs_samples.n_converged >= 0.95*n_boot + bs_samples = reduce(hcat, bs_samples.samples[bs_samples.converged_mask]) se_bs_2 = sqrt.(var(bs_samples, corrected = false, dims = 2)) #println(maximum(abs.(se_bs_2 - se_bs))) @test isapprox(se_bs_2, se_bs, rtol = rtol_bs) @@ -167,14 +166,14 @@ function test_bootstrap( end end -function smoketest_bootstrap(model_fit, spec; n_boot = 5) - # hessian and bootstrap se are close - se_bs = se_bootstrap(model_fit, spec; n_boot = n_boot) - bs_samples = bootstrap(model_fit, spec; n_boot = n_boot) +function smoketest_bootstrap(model_fit::SemFit; n_boot = 5) + # just test that both methods succeed + se_bs = se_bootstrap(model_fit; n_boot = n_boot) + bs_samples = bootstrap(model_fit; n_boot = n_boot) return se_bs, bs_samples end -function smoketest_CI_z(model_fit, partable) +function smoketest_CI_z(model_fit::SemFit, partable) se_he = @suppress se_hessian(model_fit) normal_CI!(partable, model_fit, se_he) z_test!(partable, model_fit, se_he) diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl index 48723fbe..462deab6 100644 --- a/test/examples/multigroup/build_models.jl +++ b/test/examples/multigroup/build_models.jl @@ -249,7 +249,7 @@ end lav_col = :se, lav_groups = Dict(:Pasteur => 1, :Grant_White => 2), ) - # test_bootstrap(solution_ls, partable; compare_bs = false, rtol_hessian = 0.3) + # test_bootstrap(solution_ls; compare_bs = false, rtol_hessian = 0.3) smoketest_CI_z(solution_ls, partable) end @@ -360,7 +360,7 @@ if !isnothing(specification_miss_g1) fitmeasure_names = Dict(:CFI => "cfi"), ) - test_bootstrap(solution, partable_miss; compare_bs = false, rtol_hessian = 0.5) + test_bootstrap(solution; compare_bs = false, rtol_hessian = 0.5) smoketest_CI_z(solution, partable_miss) update_se_hessian!(partable_miss, solution) diff --git a/test/examples/political_democracy/constructor.jl b/test/examples/political_democracy/constructor.jl index 25a6da91..759875b2 100644 --- a/test/examples/political_democracy/constructor.jl +++ b/test/examples/political_democracy/constructor.jl @@ -116,7 +116,7 @@ end lav_col = :se, ) - test_bootstrap(solution_ml, partable) + test_bootstrap(solution_ml) smoketest_CI_z(solution_ml, partable) end @@ -146,7 +146,7 @@ end lav_col = :se, ) - test_bootstrap(solution_ls, partable; compare_bs = false) + test_bootstrap(solution_ls; compare_bs = false) smoketest_CI_z(solution_ls, partable) end @@ -324,7 +324,7 @@ end lav_col = :se, ) - test_bootstrap(solution_ml, partable_mean) + test_bootstrap(solution_ml) smoketest_CI_z(solution_ml, partable_mean) end @@ -353,8 +353,8 @@ end lav_col = :se, ) - test_bootstrap(solution_ls, partable_mean, compare_bs = false) - # smoketest_bootstrap(solution_ls, partable_mean) + test_bootstrap(solution_ls, compare_bs = false) + # smoketest_bootstrap(solution_ls) smoketest_CI_z(solution_ls, partable_mean) end @@ -481,7 +481,7 @@ end lav_col = :se, ) - # test_bootstrap(solution_ml, partable_mean) # too much compute - smoketest_bootstrap(solution_ml, partable_mean) + # test_bootstrap(solution_ml) # too much compute + smoketest_bootstrap(solution_ml) smoketest_CI_z(solution_ml, partable_mean) end From 24261d52055b8400400b8271774a7db483bd6906 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 21 Mar 2026 17:56:25 -0700 Subject: [PATCH 10/18] CFI: sync with Sem refactor --- src/frontend/fit/fitmeasures/CFI.jl | 47 ++++++++++++++--------------- 1 file changed, 22 insertions(+), 25 deletions(-) diff --git a/src/frontend/fit/fitmeasures/CFI.jl b/src/frontend/fit/fitmeasures/CFI.jl index 9f3c5a2d..e2bcb8c7 100644 --- a/src/frontend/fit/fitmeasures/CFI.jl +++ b/src/frontend/fit/fitmeasures/CFI.jl @@ -5,11 +5,11 @@ Calculate the Comparative Fit Index (CFI). -The CFI ranges from 0-1 and measures how much better the model +The CFI ranges from 0-1 and measures how much better the model fits the data compared to a baseline model. If no baseline model is provided, a model with unconstrained variances (and means) is compaired against. -For multigroup models, variances (and means) per group are free +For multigroup models, variances (and means) per group are free without any equality constraints between groups. """ function CFI end @@ -35,34 +35,31 @@ function CFI(χ², dof, χ²₀, dof₀) end ### -function χ²_varonly(model::AbstractSemSingle) - check_single_lossfun(model; throw_error = true) - return χ²_varonly(model.loss.functions[1], model) -end - -function χ²_varonly(model::SemEnsemble) - check_single_lossfun(model; throw_error = true) - return sum(χ²_varonly, model.sems) +function χ²_varonly(model::AbstractSem) + check_same_semterm_type(model; throw_error = true) + return sum(sem_terms(model)) do semterm + χ²_varonly(_unwrap(loss(semterm))) + end end -function χ²_varonly(::SemML, model::AbstractSemSingle) - N⁻ = (nsamples(model) - 1) - S = obs_cov(observed(model)) +function χ²_varonly(loss::SemML) + N⁻ = (nsamples(loss) - 1) + S = obs_cov(observed(loss)) Σ₀ = Diagonal(S) - p = nobserved_vars(model) + p = nobserved_vars(loss) return N⁻*(logdet(Σ₀) + tr(inv(Σ₀)*S) - logdet(S) - p) end # for the optimal variance only model, we have to solve 1/2 tr((I-XS⁻¹)^2) with X diagonal -function χ²_varonly(::SemWLS, model) - N⁻ = (nsamples(model) - 1) - S⁻¹ = inv((obs_cov(observed(model)))) +function χ²_varonly(loss::SemWLS) + N⁻ = (nsamples(loss) - 1) + S⁻¹ = inv((obs_cov(observed(loss)))) Σ₀ = Diagonal(inv(S⁻¹ .* S⁻¹)*diag(S⁻¹)) return N⁻*0.5*tr((I - Σ₀*S⁻¹)^2) end # For FIML, an explicit bl model has to be passed -function χ²_varonly(::SemFIML, model) +function χ²_varonly(loss::SemFIML) """ Computing the CFI with FIML requires explicitely passing a fitted baseline model as CFI(fit::SemFit, fit_baseline::SemFit) @@ -71,12 +68,12 @@ function χ²_varonly(::SemFIML, model) throw end -function dof_varonly(model::AbstractSemSingle) - nparams_varonly = nobserved_vars(model) - if MeanStruct(model.implied) === HasMeanStruct - nparams_varonly *= 2 +function dof_varonly(model::AbstractSem) + return sum(sem_terms(model)) do semterm + nparams_varonly = nobserved_vars(semterm) + if MeanStruct(implied(semterm)) === HasMeanStruct + nparams_varonly *= 2 + end + return n_dp(loss(semterm)) - nparams_varonly end - return n_dp(model) - nparams_varonly end - -dof_varonly(model::SemEnsemble) = sum(dof_varonly, model.sems) From e4d38e582e19d0ee0b6a8d13e407928a99305f24 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 21 Mar 2026 17:57:46 -0700 Subject: [PATCH 11/18] test/build_models: remove redundant model --- test/examples/multigroup/build_models.jl | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl index 462deab6..329c5502 100644 --- a/test/examples/multigroup/build_models.jl +++ b/test/examples/multigroup/build_models.jl @@ -347,19 +347,6 @@ if !isnothing(specification_miss_g1) lav_groups = Dict(:Pasteur => 1, :Grant_White => 2), ) - solution = fit(semoptimizer, model_ml_multigroup2) - test_fitmeasures( - fit_measures(solution), - solution_lav[:fitmeasures_fiml]; - rtol = 1e-3, - atol = 0, - ) - test_fitmeasures( - Dict(:CFI => CFI(solution, solution_varonly)), - solution_lav[:fitmeasures_fiml]; - fitmeasure_names = Dict(:CFI => "cfi"), - ) - test_bootstrap(solution; compare_bs = false, rtol_hessian = 0.5) smoketest_CI_z(solution, partable_miss) From cb9b1e772392ed6697f4c641c8bff1830fae9ce3 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 21 Mar 2026 17:57:46 -0700 Subject: [PATCH 12/18] revert using --- test/examples/multigroup/multigroup.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/examples/multigroup/multigroup.jl b/test/examples/multigroup/multigroup.jl index dd654731..c8ae8c1f 100644 --- a/test/examples/multigroup/multigroup.jl +++ b/test/examples/multigroup/multigroup.jl @@ -1,5 +1,5 @@ using StructuralEquationModels, Test, FiniteDiff, Suppressor -using LinearAlgebra: diagind, LowerTriangular +using LinearAlgebra: diagind, isposdef, logdet, tr, LowerTriangular using Statistics: var using Random From afac0b4b56dcbfa3599b6d6ba43f1490729eb6bb Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 21 Mar 2026 17:57:46 -0700 Subject: [PATCH 13/18] WLS: verbose option to suppress info about inv(obs_cov) --- src/loss/WLS/WLS.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/loss/WLS/WLS.jl b/src/loss/WLS/WLS.jl index 8f4a109c..d04bc346 100644 --- a/src/loss/WLS/WLS.jl +++ b/src/loss/WLS/WLS.jl @@ -62,6 +62,7 @@ function SemWLS( wls_weight_matrix::Union{AbstractMatrix, Nothing} = nothing, wls_weight_matrix_mean::Union{AbstractMatrix, Nothing} = nothing, approximate_hessian::Bool = false, + verbose::Bool = false, kwargs..., ) if observed isa SemObservedMissing @@ -114,7 +115,8 @@ function SemWLS( if MeanStruct(implied) == HasMeanStruct if isnothing(wls_weight_matrix_mean) - @info "Computing WLS weight matrix for the meanstructure using obs_cov()" + verbose && + @info "Computing WLS weight matrix for the meanstructure using obs_cov()" wls_weight_matrix_mean = inv(obs_cov(observed)) end size(wls_weight_matrix_mean) == (nobs_vars, nobs_vars) || DimensionMismatch( From 53a615a3b602e4485a300987c24f9f3035249ff2 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sat, 21 Mar 2026 17:57:46 -0700 Subject: [PATCH 14/18] docs: sync with Sem refactor --- docs/src/developer/implied.md | 30 ++++------- docs/src/developer/loss.md | 21 ++------ docs/src/developer/optimizer.md | 8 --- docs/src/developer/sem.md | 17 +++---- docs/src/internals/types.md | 20 +++++--- docs/src/performance/mixed_differentiation.md | 18 +++---- docs/src/performance/simulation.md | 17 +------ docs/src/tutorials/collection/collection.md | 51 +++++++++++++------ docs/src/tutorials/collection/multigroup.md | 17 ++++--- 9 files changed, 91 insertions(+), 108 deletions(-) diff --git a/docs/src/developer/implied.md b/docs/src/developer/implied.md index 056cd663..6321decb 100644 --- a/docs/src/developer/implied.md +++ b/docs/src/developer/implied.md @@ -13,9 +13,9 @@ end and a method to update!: ```julia -import StructuralEquationModels: objective! +import StructuralEquationModels: update! -function update!(targets::EvaluationTargets, implied::MyImplied, model::AbstractSemSingle, params) +function update!(targets::EvaluationTargets, implied::MyImplied, params) if is_objective_required(targets) ... @@ -31,11 +31,9 @@ function update!(targets::EvaluationTargets, implied::MyImplied, model::Abstract end ``` -As you can see, `update` gets passed as a first argument `targets`, which is telling us whether the objective value, gradient, and/or hessian are needed. +As you can see, `update!` gets passed as a first argument `targets`, which is telling us whether the objective value, gradient, and/or hessian are needed. We can then use the functions `is_..._required` and conditional on what the optimizer needs, we can compute and store things we want to make available to the loss functions. For example, as we have seen in [Second example - maximum likelihood](@ref), the `RAM` implied type computes the model-implied covariance matrix and makes it available via `implied.Σ`. - - Just as described in [Custom loss functions](@ref), you may define a constructor. Typically, this will depend on the `specification = ...` argument that can be a `ParameterTable` or a `RAMMatrices` object. We implement an `ImpliedEmpty` type in our package that does nothing but serving as an `implied` field in case you are using a loss function that does not need any implied type at all. You may use it as a template for defining your own implied type, as it also shows how to handle the specification objects: @@ -56,7 +54,7 @@ Empty placeholder for models that don't need an implied part. - `specification`: either a `RAMMatrices` or `ParameterTable` object # Examples -A multigroup model with ridge regularization could be specified as a `SemEnsemble` with one +A multigroup model with ridge regularization could be specified as a `Sem` with one model per group and an additional model with `ImpliedEmpty` and `SemRidge` for the regularization part. # Extended help @@ -75,26 +73,20 @@ end ### Constructors ############################################################################################ -function ImpliedEmpty(; - specification, - meanstruct = NoMeanStruct(), - hessianeval = ExactHessian(), +function ImpliedEmpty( + spec::SemSpecification; + hessianeval::HessianApprox = ExactHessian(), kwargs..., ) - return ImpliedEmpty(hessianeval, meanstruct, convert(RAMMatrices, specification)) + ram_matrices = convert(RAMMatrices, spec) + return ImpliedEmpty(hessianeval, MeanStruct(ram_matrices), ram_matrices) end ############################################################################################ ### methods ############################################################################################ -update!(targets::EvaluationTargets, implied::ImpliedEmpty, par, model) = nothing - -############################################################################################ -### Recommended methods -############################################################################################ - -update_observed(implied::ImpliedEmpty, observed::SemObserved; kwargs...) = implied +update!(targets::EvaluationTargets, implied::ImpliedEmpty, par) = nothing ``` -As you see, similar to [Custom loss functions](@ref) we implement a method for `update_observed`. \ No newline at end of file +As you see, similar to [Custom loss functions](@ref) we implement a constructor. \ No newline at end of file diff --git a/docs/src/developer/loss.md b/docs/src/developer/loss.md index d6949842..aa6a1e17 100644 --- a/docs/src/developer/loss.md +++ b/docs/src/developer/loss.md @@ -11,9 +11,9 @@ Since we allow for the optimization of sums of loss functions, and the maximum l using StructuralEquationModels ``` -To define a new loss function, you have to define a new type that is a subtype of `SemLossFunction`: +To define a new loss function, you have to define a new type that is a subtype of `AbstractLoss`: ```@example loss -struct Ridge <: SemLossFunction +struct MyRidge <: AbstractLoss α I end @@ -25,8 +25,8 @@ Additionaly, we need to define a *method* of the function `evaluate!` to compute ```@example loss import StructuralEquationModels: evaluate! -evaluate!(objective::Number, gradient::Nothing, hessian::Nothing, ridge::Ridge, model::AbstractSem, par) = - ridge.α * sum(i -> par[i]^2, ridge.I) +evaluate!(objective::Number, gradient::Nothing, hessian::Nothing, ridge::MyRidge, par) = + ridge.α * sum(i -> abs2(par[i]), ridge.I) ``` The function `evaluate!` recognizes by the types of the arguments `objective`, `gradient` and `hessian` whether it should compute the objective value, gradient or hessian of the model w.r.t. the parameters. @@ -98,7 +98,7 @@ function evaluate!(objective, gradient, hessian::Nothing, ridge::Ridge, model::A gradient[ridge.I] .= 2 * ridge.α * par[ridge.I] end # compute objective - if !isnothing(objective) + if !isnothing(objective) return ridge.α * sum(i -> par[i]^2, ridge.I) end end @@ -166,17 +166,6 @@ end ## Additional functionality -### Update observed data - -If you are planing a simulation study where you have to fit the **same model** to many **different datasets**, it is computationally beneficial to not build the whole model completely new everytime you change your data. -Therefore, we provide a function to update the data of your model, `replace_observed(model(semfit); data = new_data)`. However, we can not know beforehand in what way your loss function depends on the specific datasets. The solution is to provide a method for `update_observed`. Since `Ridge` does not depend on the data at all, this is quite easy: - -```julia -import StructuralEquationModels: update_observed - -update_observed(ridge::Ridge, observed::SemObserved; kwargs...) = ridge -``` - ### Access additional information If you want to provide a way to query information about loss functions of your type, you can provide functions for that: diff --git a/docs/src/developer/optimizer.md b/docs/src/developer/optimizer.md index b5c9a6e0..4659ba5d 100644 --- a/docs/src/developer/optimizer.md +++ b/docs/src/developer/optimizer.md @@ -25,12 +25,6 @@ struct MyoptResult{O <: SemOptimizerMyopt} <: SEM.SemOptimizerResult{O} ... end -############################################################################################ -### Recommended methods -############################################################################################ - -update_observed(optimizer::SemOptimizerMyopt, observed::SemObserved; kwargs...) = optimizer - ############################################################################################ ### additional methods ############################################################################################ @@ -43,8 +37,6 @@ and `SEM.sem_optimizer_subtype(::Val{:Myopt})` returns `SemOptimizerMyopt`. This instructs *SEM.jl* to use `SemOptimizerMyopt` when `:Myopt` is specified as the engine for model fitting: `fit(..., engine = :Myopt)`. -A method for `update_observed` and additional methods might be usefull, but are not necessary. - Now comes the essential part: we need to provide the [`fit`](@ref) method with `SemOptimizerMyopt` as the first positional argument. diff --git a/docs/src/developer/sem.md b/docs/src/developer/sem.md index c54ff26a..bb077f43 100644 --- a/docs/src/developer/sem.md +++ b/docs/src/developer/sem.md @@ -1,13 +1,14 @@ # Custom model types -The abstract supertype for all models is `AbstractSem`, which has two subtypes, `AbstractSemSingle{O, I, L}` and `AbstractSemCollection`. Currently, there are 2 subtypes of `AbstractSemSingle`: `Sem`, `SemFiniteDiff`. All subtypes of `AbstractSemSingle` should have at least observed, implied, loss and optimizer fields, and share their types (`{O, I, L}`) with the parametric abstract supertype. For example, the `SemFiniteDiff` type is implemented as +The abstract supertype for all models is [`AbstractSem`](@ref). Currently, there are 2 concrete subtypes: +`Sem{L <: Tuple}` and `SemFiniteDiff{S <: AbstractSem}`. +A `Sem` model holds a tuple of `LossTerm`s (each wrapping an `AbstractLoss`) and a vector of parameter labels. Both single-group and multigroup models are represented as `Sem`. + +`SemFiniteDiff` wraps any `AbstractSem` and substitutes dedicated gradient/hessian evaluation with finite difference approximation: ```julia -struct SemFiniteDiff{O <: SemObserved, I <: SemImplied, L <: SemLoss} <: - AbstractSemSingle{O, I, L} - observed::O - implied::I - loss::L +struct SemFiniteDiff{S <: AbstractSem} <: AbstractSem + model::S end ``` @@ -17,6 +18,4 @@ Additionally, you can change how objective/gradient/hessian values are computed evaluate!(objective, gradient, hessian, model::SemFiniteDiff, params) = ... ``` -Additionally, we can define constructors like the one in `"src/frontend/specification/Sem.jl"`. - -It is also possible to add new subtypes for `AbstractSemCollection`. \ No newline at end of file +Additionally, we can define constructors like the one in `"src/frontend/specification/Sem.jl"`. \ No newline at end of file diff --git a/docs/src/internals/types.md b/docs/src/internals/types.md index e70a52ca..4b4cd4fa 100644 --- a/docs/src/internals/types.md +++ b/docs/src/internals/types.md @@ -2,12 +2,16 @@ The type hierarchy is implemented in `"src/types.jl"`. -`AbstractSem`: the most abstract type in our package -- `AbstractSemSingle{O, I, L} <: AbstractSem` is an abstract parametric type that is a supertype of all single models - - `Sem`: models that do not need automatic differentiation or finite difference approximation - - `SemFiniteDiff`: models whose gradients and/or hessians should be computed via finite difference approximation -- `AbstractSemCollection <: AbstractSem` is an abstract supertype of all models that contain multiple `AbstractSem` submodels +[`AbstractLoss`](@ref): is the base abstract type for all loss functions: +- `SemLoss{O <: SemObserved, I <: SemImplied}`: is the subtype of `AbstractLoss`, which is the + base for all SEM-specific loss functions ([`SemML`](@ref), [`SemWLS`](@ref) etc) that + evaluate how closely the implied covariation structure (represented by the object of type `I`) + matches the observed one (contained in the object of type `O`); +- regularizing terms (e.g. [`SemRidge`](@ref)) are implemented as subtypes of `AbstractLoss`. -Every `AbstractSemSingle` has to have `SemObserved`, `SemImplied`, and `SemLoss` fields (and can have additional fields). - -`SemLoss` is a container for multiple `SemLossFunctions`. \ No newline at end of file +[`AbstractSem`](@ref) is the base abstract type for all SEM models. It has two concrete subtypes: +- `Sem{L <: Tuple} <: AbstractSem`: the main SEM model type that implements a list of weighted +loss terms (using [`LossTerm`](@ref) wrapper around `AbstractLoss`) and allows modeling both single +and multi-group SEMs and combining them with regularization terms. +- `SemFiniteDiff{S <: AbstractSem} <: AbstractSem`: a wrapper around any `AbstractSem` that + substitutes dedicated gradient/hessian evaluation with finite difference approximation. diff --git a/docs/src/performance/mixed_differentiation.md b/docs/src/performance/mixed_differentiation.md index b7ae333b..f33fa6ab 100644 --- a/docs/src/performance/mixed_differentiation.md +++ b/docs/src/performance/mixed_differentiation.md @@ -2,22 +2,20 @@ This way of specifying our model is not ideal, however, because now also the maximum likelihood loss function lives inside a `SemFiniteDiff` model, and this means even though we have defined analytical gradients for it, we do not make use of them. -A more efficient way is therefore to specify our model as an ensemble model: +A more efficient way is therefore to specify our model as a combined model with multiple loss terms: ```julia -model_ml = Sem( - specification = partable, - data = data, - loss = SemML +ml_term = SemML( + SemObservedData(data = data, specification = partable), + RAMSymbolic(partable) ) -model_ridge = SemFiniteDiff( - specification = partable, - data = data, - loss = myridge +ridge_term = SemRidge( + α_ridge = 0.01, + which_ridge = params(ml_term) ) -model_ml_ridge = SemEnsemble(model_ml, model_ridge) +model_ml_ridge = Sem(ml_term, ridge_term) model_ml_ridge_fit = fit(model_ml_ridge) ``` diff --git a/docs/src/performance/simulation.md b/docs/src/performance/simulation.md index 85a0c0a0..61a9d5ad 100644 --- a/docs/src/performance/simulation.md +++ b/docs/src/performance/simulation.md @@ -57,19 +57,7 @@ model = Sem( data = data_1 ) -model_updated = replace_observed(model; data = data_2, specification = partable) -``` - -If you are building your models by parts, you can also update each part seperately with the function `update_observed`. -For example, - -```@example replace_observed - -new_observed = SemObservedData(;data = data_2, specification = partable) - -my_optimizer = SemOptimizer() - -new_optimizer = update_observed(my_optimizer, new_observed) +model_updated = replace_observed(model, data_2) ``` ## Multithreading @@ -90,7 +78,7 @@ model1 = Sem( data = data_1 ) -model2 = deepcopy(replace_observed(model; data = data_2, specification = partable)) +model2 = deepcopy(replace_observed(model, data_2)) models = [model1, model2] fits = Vector{SemFit}(undef, 2) @@ -104,5 +92,4 @@ end ```@docs replace_observed -update_observed ``` \ No newline at end of file diff --git a/docs/src/tutorials/collection/collection.md b/docs/src/tutorials/collection/collection.md index f60b7312..2a8ea92c 100644 --- a/docs/src/tutorials/collection/collection.md +++ b/docs/src/tutorials/collection/collection.md @@ -1,31 +1,52 @@ # Collections -With StructuralEquationModels.jl, you can fit weighted sums of structural equation models. -The most common use case for this are [Multigroup models](@ref). -Another use case may be optimizing the sum of loss functions for some of which you do know the analytic gradient, but not for others. -In this case, you can optimize the sum of a `Sem` and a `SemFiniteDiff` (or any other differentiation method). +With *StructuralEquationModels.jl*, you can fit weighted sums of structural equation models. +The most common use case for this are [Multigroup models](@ref). +Another use case may be optimizing the sum of loss functions for some of which you do know the analytic gradient, but not for others. +In this case, [`FiniteDiffWrapper`](@ref) can generate a wrapper around the specific `SemLoss` term. The wrapper loss term will +only use the objective of the original term to calculate its gradient using finite difference approximation. -To use this feature, you have to construct a `SemEnsemble` model, which is actually quite easy: +```julia +loss_1 = SemML(observed_1, implied_1) +loss_2 = SemML(observed_2, implied_2) +loss_2_findiff = FiniteDiffWrapper(loss_2) +``` + +To construct `Sem` from the the individual `SemLoss` (or other `AbstractLoss`) terms, they are +just passed to the `Sem` constructor: ```julia -# models -model_1 = Sem(...) +model = Sem(loss_1, loss_2) +model_findiff = Sem(loss_1, loss_2_findiff) +``` + +It is also possible to use finite difference for the entire `Sem` model: -model_2 = SemFiniteDiff(...) +```julia +model_findiff2 = FiniteDiffWrapper(model) +``` -model_3 = Sem(...) +The weighting scheme of the SEM loss terms is specified using `default_set_weights` argument of the `Sem` constructor. +The `:nsamples` scheme (the default) weights SEM terms by ``N_{term}/N_{total}``, i.e. each term is weighted by the number +of observations in its data (which matches the formula for multigroup models). +The weights for the loss terms (both SEM and regularization) can be explicitly specified the pair syntax `loss => weight`: -model_ensemble = SemEnsemble(model_1, model_2, model_3) +```julia +model_weighted = Sem(loss_1 => 0.5, loss_2 => 1.0) ``` -So you just construct the individual models (however you like) and pass them to `SemEnsemble`. -You may also pass a vector of weigths to `SemEnsemble`. By default, those are set to ``N_{model}/N_{total}``, i.e. each model is weighted by the number of observations in it's data (which matches the formula for multigroup models). +`Sem` support assigning unique identifier to each loss term, which is essential for complex multi-term model. +The syntax is `id => loss`, or `id => loss => weight`: -Multigroup models can also be specified via the graph interface; for an example, see [Multigroup models](@ref). +```julia +model2 = Sem(:main => loss_1, :alt => loss_2) +model2_weighted = Sem(:main => loss_1 => 0.5, :alt => loss_2 => 1.0) +``` # API - collections ```@docs -SemEnsemble -AbstractSemCollection +Sem +LossTerm +FiniteDiffWrapper ``` \ No newline at end of file diff --git a/docs/src/tutorials/collection/multigroup.md b/docs/src/tutorials/collection/multigroup.md index 16d3dcd7..04f1893d 100644 --- a/docs/src/tutorials/collection/multigroup.md +++ b/docs/src/tutorials/collection/multigroup.md @@ -4,19 +4,20 @@ using StructuralEquationModels ``` -As an example, we will fit the model from [the `lavaan` tutorial](https://lavaan.ugent.be/tutorial/groups.html) with loadings constrained to equality across groups. +As an example, we will fit the model from [the `lavaan` tutorial](https://lavaan.ugent.be/tutorial/groups.html) +with loadings constrained to equality across groups. -We first load the example data. +We first load the example data. We have to make sure that the column indicating the group (here called `school`) is a vector of `Symbol`s, not strings - so we convert it. ```@setup mg dat = example_data("holzinger_swineford") -dat.school = ifelse.(dat.school .== "Pasteur", :Pasteur, :Grant_White) +dat.school = Symbol.(replace.(dat.school, "-" => "_")) ``` ```julia dat = example_data("holzinger_swineford") -dat.school = ifelse.(dat.school .== "Pasteur", :Pasteur, :Grant_White) +dat.school = Symbol.(replace.(dat.school, "-" => "_")) ``` We then specify our model via the graph interface: @@ -59,19 +60,19 @@ You can then use the resulting graph to specify an `EnsembleParameterTable` groups = [:Pasteur, :Grant_White] partable = EnsembleParameterTable( - graph, + graph, observed_vars = observed_vars, latent_vars = latent_vars, groups = groups) ``` -The parameter table can be used to create a `SemEnsemble` model: +The parameter table can be used to create a multigroup `Sem` model: ```@example mg; ansicolor = true -model_ml_multigroup = SemEnsemble( +model_ml_multigroup = Sem( specification = partable, data = dat, - column = :school, + semterm_column = :school, groups = groups) ``` From 240e3cd42058aa3ccecb2c44dcc521de0bfe116c Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 22 Mar 2026 13:19:46 -0700 Subject: [PATCH 15/18] test: fix formatting --- test/examples/multigroup/build_models.jl | 12 +++-------- .../political_democracy/constructor.jl | 20 ++++--------------- 2 files changed, 7 insertions(+), 25 deletions(-) diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl index 329c5502..6811cb40 100644 --- a/test/examples/multigroup/build_models.jl +++ b/test/examples/multigroup/build_models.jl @@ -19,17 +19,11 @@ model_ml_multigroup = Sem( end # replace observed using Dict of data matrices -model_ml_multigroup3 = replace_observed( - model_ml_multigroup, - Dict(:Pasteur => dat_g1, :Grant_White => dat_g2), -) +model_ml_multigroup3 = + replace_observed(model_ml_multigroup, Dict(:Pasteur => dat_g1, :Grant_White => dat_g2)) # replace observed using DataFrame with group column -model_ml_multigroup4 = replace_observed( - model_ml_multigroup, - dat; - semterm_column = :school, -) +model_ml_multigroup4 = replace_observed(model_ml_multigroup, dat; semterm_column = :school) # gradients @testset "ml_gradients_multigroup" begin diff --git a/test/examples/political_democracy/constructor.jl b/test/examples/political_democracy/constructor.jl index 759875b2..48ba1b96 100644 --- a/test/examples/political_democracy/constructor.jl +++ b/test/examples/political_democracy/constructor.jl @@ -167,14 +167,8 @@ end # set seed for simulation Random.seed!(83472834) # simulate data - model_ml_new = replace_observed( - model_ml, - rand(model_ml, params, 1_000_000), - ) - model_ml_sym_new = replace_observed( - model_ml_sym, - rand(model_ml_sym, params, 1_000_000), - ) + model_ml_new = replace_observed(model_ml, rand(model_ml, params, 1_000_000)) + model_ml_sym_new = replace_observed(model_ml_sym, rand(model_ml_sym, params, 1_000_000)) # fit models sol_ml = solution(fit(semoptimizer, model_ml_new)) sol_ml_sym = solution(fit(semoptimizer, model_ml_sym_new)) @@ -376,14 +370,8 @@ end # set seed for simulation Random.seed!(83472834) # simulate data - model_ml_new = replace_observed( - model_ml, - rand(model_ml, params, 1_000_000), - ) - model_ml_sym_new = replace_observed( - model_ml_sym, - rand(model_ml_sym, params, 1_000_000), - ) + model_ml_new = replace_observed(model_ml, rand(model_ml, params, 1_000_000)) + model_ml_sym_new = replace_observed(model_ml_sym, rand(model_ml_sym, params, 1_000_000)) # fit models sol_ml = solution(fit(semoptimizer, model_ml_new)) sol_ml_sym = solution(fit(semoptimizer, model_ml_sym_new)) From a277cb0e87a77e9b6772996aaac10dc0b588e78a Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 22 Mar 2026 20:54:23 -0700 Subject: [PATCH 16/18] fit_measures(): support vectors of funcs also add CFI to the list --- src/frontend/fit/fitmeasures/fit_measures.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/frontend/fit/fitmeasures/fit_measures.jl b/src/frontend/fit/fitmeasures/fit_measures.jl index 185b348c..7fabc950 100644 --- a/src/frontend/fit/fitmeasures/fit_measures.jl +++ b/src/frontend/fit/fitmeasures/fit_measures.jl @@ -1,6 +1,8 @@ -fit_measures(fit) = fit_measures(fit, nparams, dof, AIC, BIC, RMSEA, χ², p_value, minus2ll) +const DEFAULT_FIT_MEASURES = [AIC, BIC, dof, χ², p_value, nparams, RMSEA, CFI] -fit_measures(fit, measures...) = Dict(Symbol(fn) => fn(fit) for fn in measures) +fit_measures(fit, measures::AbstractVector) = Dict(Symbol(fn) => fn(fit) for fn in measures) +fit_measures(fit, measures...) = fit_measures(fit, measures) +fit_measures(fit) = fit_measures(fit, DEFAULT_FIT_MEASURES) """ fit_measures(fit::SemFit, measures...) -> Dict{Symbol} @@ -20,6 +22,7 @@ fit_measures(semfit, nparams, dof, p_value) ``` # See also -[`AIC`](@ref), [`BIC`](@ref), [`RMSEA`](@ref), [`χ²`](@ref), [`p_value`](@ref), [`minus2ll`](@ref) +[`AIC`](@ref), [`BIC`](@ref), [`RMSEA`](@ref), [`χ²`](@ref), [`p_value`](@ref), +[`minus2ll`](@ref), [`CFI`](@ref) """ fit_measures From 60dbdc7a250a9e289d928e2847823f8eede82834 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 22 Mar 2026 20:58:00 -0700 Subject: [PATCH 17/18] test_fitmeasures(): refactor/simplify --- test/examples/helper.jl | 34 ++++++---- test/examples/multigroup/build_models.jl | 47 ++------------ test/examples/political_democracy/by_parts.jl | 36 ++--------- .../political_democracy/constructor.jl | 64 ++++--------------- 4 files changed, 42 insertions(+), 139 deletions(-) diff --git a/test/examples/helper.jl b/test/examples/helper.jl index f14fec62..fed95f3c 100644 --- a/test/examples/helper.jl +++ b/test/examples/helper.jl @@ -49,7 +49,8 @@ function test_hessian(model, params; rtol = 1e-4, atol = 0) @test hessian ≈ true_hessian rtol = rtol atol = atol end -fitmeasure_names_ml = Dict( +# map from the SEM.jl name of the fit measure to the lavaan's one +fitmeasure_semjl_to_lavaan = Dict( :AIC => "aic", :BIC => "bic", :dof => "df", @@ -57,26 +58,31 @@ fitmeasure_names_ml = Dict( :p_value => "pvalue", :nparams => "npar", :RMSEA => "rmsea", -) - -fitmeasure_names_ls = Dict( - :dof => "df", - :χ² => "chisq", - :p_value => "pvalue", - :nparams => "npar", - :RMSEA => "rmsea", + :CFI => "cfi", ) function test_fitmeasures( - measures, + fitted::SemFit, measures_lav; + fitmeasures::AbstractVector = SEM.DEFAULT_FIT_MEASURES, + fitted_baseline::Union{SemFit, Nothing} = nothing, rtol = 1e-4, atol = 0, - fitmeasure_names = fitmeasure_names_ml, ) - @testset "$name" for (key, name) in pairs(fitmeasure_names) - measure_lav = measures_lav.x[findfirst(==(name), measures_lav[!, 1])] - @test measures[key] ≈ measure_lav rtol = rtol atol = atol + @testset "$fn" for fn in fitmeasures + name = Symbol(fn) + # FIML CFI requires the baseline model + measure = + fn != CFI || isnothing(fitted_baseline) ? fn(fitted) : + fn(fitted, fitted_baseline) + lav_name = fitmeasure_semjl_to_lavaan[name] + lav_ix = findfirst(==(lav_name), measures_lav[!, 1]) + if isnothing(lav_ix) + @test ismissing(measure) + else + measure_lav = measures_lav.x[lav_ix] + @test measure ≈ measure_lav rtol = rtol atol = atol + end end end diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl index 6811cb40..71128760 100644 --- a/test/examples/multigroup/build_models.jl +++ b/test/examples/multigroup/build_models.jl @@ -52,17 +52,7 @@ end @testset "fitmeasures/se_ml" begin solution_ml = fit(semoptimizer, model_ml_multigroup) - test_fitmeasures( - fit_measures(solution_ml), - solution_lav[:fitmeasures_ml]; - rtol = 1e-2, - atol = 1e-7, - ) - test_fitmeasures( - Dict(:CFI => CFI(solution_ml)), - solution_lav[:fitmeasures_ml]; - fitmeasure_names = Dict(:CFI => "cfi"), - ) + test_fitmeasures(solution_ml, solution_lav[:fitmeasures_ml]; rtol = 1e-2, atol = 1e-7) update_se_hessian!(partable, solution_ml) test_estimates( @@ -122,17 +112,7 @@ end @testset "fitmeasures/se_ml | sorted" begin solution_ml = fit(semoptimizer, model_ml_multigroup) - test_fitmeasures( - fit_measures(solution_ml), - solution_lav[:fitmeasures_ml]; - rtol = 1e-2, - atol = 1e-7, - ) - test_fitmeasures( - Dict(:CFI => CFI(solution_ml)), - solution_lav[:fitmeasures_ml]; - fitmeasure_names = Dict(:CFI => "cfi"), - ) + test_fitmeasures(solution_ml, solution_lav[:fitmeasures_ml]; rtol = 1e-2, atol = 1e-7) update_se_hessian!(partable_s, solution_ml) test_estimates( @@ -221,18 +201,7 @@ end @testset "fitmeasures/se_ls" begin solution_ls = fit(semoptimizer, model_ls_multigroup) - test_fitmeasures( - fit_measures(solution_ls), - solution_lav[:fitmeasures_ls]; - fitmeasure_names = fitmeasure_names_ls, - rtol = 1e-2, - atol = 1e-5, - ) - test_fitmeasures( - Dict(:CFI => CFI(solution_ls)), - solution_lav[:fitmeasures_ls]; - fitmeasure_names = Dict(:CFI => "cfi"), - ) + test_fitmeasures(solution_ls, solution_lav[:fitmeasures_ls]; rtol = 1e-2, atol = 1e-5) @suppress update_se_hessian!(partable, solution_ls) test_estimates( @@ -319,18 +288,14 @@ if !isnothing(specification_miss_g1) @testset "fitmeasures/se_fiml" begin solution = fit(semoptimizer, model_ml_multigroup) + solution_varonly = fit(semoptimizer, model_ml_varonly) test_fitmeasures( - fit_measures(solution), + solution, solution_lav[:fitmeasures_fiml]; + fitted_baseline = solution_varonly, rtol = 1e-3, atol = 0, ) - solution_varonly = fit(semoptimizer, model_ml_varonly) - test_fitmeasures( - Dict(:CFI => CFI(solution, solution_varonly)), - solution_lav[:fitmeasures_fiml]; - fitmeasure_names = Dict(:CFI => "cfi"), - ) update_se_hessian!(partable_miss, solution) test_estimates( partable_miss, diff --git a/test/examples/political_democracy/by_parts.jl b/test/examples/political_democracy/by_parts.jl index 6866eead..ef634a59 100644 --- a/test/examples/political_democracy/by_parts.jl +++ b/test/examples/political_democracy/by_parts.jl @@ -90,12 +90,7 @@ end @testset "fitmeasures/se_ml" begin solution_ml = fit(semoptimizer, model_ml) - test_fitmeasures(fit_measures(solution_ml), solution_lav[:fitmeasures_ml]; atol = 1e-3) - test_fitmeasures( - Dict(:CFI => CFI(solution_ml)), - solution_lav[:fitmeasures_ml]; - fitmeasure_names = Dict(:CFI => "cfi"), - ) + test_fitmeasures(solution_ml, solution_lav[:fitmeasures_ml]; atol = 1e-3) update_se_hessian!(partable, solution_ml) test_estimates( @@ -109,14 +104,7 @@ end @testset "fitmeasures/se_ls" begin solution_ls = fit(semoptimizer, model_ls_sym) - fm = fit_measures(solution_ls) - test_fitmeasures( - merge(fm, Dict(:CFI => CFI(solution_ls))), - solution_lav[:fitmeasures_ls]; - atol = 1e-3, - fitmeasure_names = merge(fitmeasure_names_ls, Dict(:CFI => "cfi")) - ) - @test ismissing(fm[:AIC]) && ismissing(fm[:BIC]) && ismissing(fm[:minus2ll]) + test_fitmeasures(solution_ls, solution_lav[:fitmeasures_ls]; atol = 1e-3) @suppress update_se_hessian!(partable, solution_ls) test_estimates( @@ -241,16 +229,7 @@ end @testset "fitmeasures/se_ml_mean" begin solution_ml = fit(semoptimizer, model_ml) - test_fitmeasures( - fit_measures(solution_ml), - solution_lav[:fitmeasures_ml_mean]; - atol = 1e-3, - ) - test_fitmeasures( - Dict(:CFI => CFI(solution_ml)), - solution_lav[:fitmeasures_ml_mean]; - fitmeasure_names = Dict(:CFI => "cfi"), - ) + test_fitmeasures(solution_ml, solution_lav[:fitmeasures_ml_mean]; atol = 1e-3) update_se_hessian!(partable_mean, solution_ml) test_estimates( @@ -264,14 +243,7 @@ end @testset "fitmeasures/se_ls_mean" begin solution_ls = fit(semoptimizer, model_ls) - fm = fit_measures(solution_ls) - test_fitmeasures( - merge(fm, Dict(:CFI => CFI(solution_ls))), - solution_lav[:fitmeasures_ls_mean]; - atol = 1e-3, - fitmeasure_names = merge(fitmeasure_names_ls, Dict(:CFI => "cfi")), - ) - @test ismissing(fm[:AIC]) && ismissing(fm[:BIC]) && ismissing(fm[:minus2ll]) + test_fitmeasures(solution_ls, solution_lav[:fitmeasures_ls_mean]; atol = 1e-3) @suppress update_se_hessian!(partable_mean, solution_ls) test_estimates( diff --git a/test/examples/political_democracy/constructor.jl b/test/examples/political_democracy/constructor.jl index 48ba1b96..2efa5abe 100644 --- a/test/examples/political_democracy/constructor.jl +++ b/test/examples/political_democracy/constructor.jl @@ -100,12 +100,7 @@ end @testset "fitmeasures/se_ml" begin solution_ml = fit(semoptimizer, model_ml) - test_fitmeasures(fit_measures(solution_ml), solution_lav[:fitmeasures_ml]; atol = 1e-3) - test_fitmeasures( - Dict(:CFI => CFI(solution_ml)), - solution_lav[:fitmeasures_ml]; - fitmeasure_names = Dict(:CFI => "cfi"), - ) + test_fitmeasures(solution_ml, solution_lav[:fitmeasures_ml]; atol = 1e-3) update_se_hessian!(partable, solution_ml) test_estimates( @@ -122,20 +117,7 @@ end @testset "fitmeasures/se_ls" begin solution_ls = fit(semoptimizer, model_ls_sym) - fm = fit_measures(solution_ls) - test_fitmeasures( - fm, - solution_lav[:fitmeasures_ls]; - atol = 1e-3, - fitmeasure_names = fitmeasure_names_ls, - ) - test_fitmeasures( - Dict(:CFI => CFI(solution_ls)), - solution_lav[:fitmeasures_ls]; - fitmeasure_names = Dict(:CFI => "cfi"), - ) - - @test ismissing(fm[:AIC]) && ismissing(fm[:BIC]) && ismissing(fm[:minus2ll]) + test_fitmeasures(solution_ls, solution_lav[:fitmeasures_ls]; atol = 1e-3) @suppress update_se_hessian!(partable, solution_ls) test_estimates( @@ -298,16 +280,7 @@ end @testset "fitmeasures/se_ml_mean" begin solution_ml = fit(semoptimizer, model_ml) - test_fitmeasures( - fit_measures(solution_ml), - solution_lav[:fitmeasures_ml_mean]; - atol = 0.002, - ) - test_fitmeasures( - Dict(:CFI => CFI(solution_ml)), - solution_lav[:fitmeasures_ml_mean]; - fitmeasure_names = Dict(:CFI => "cfi"), - ) + test_fitmeasures(solution_ml, solution_lav[:fitmeasures_ml_mean]; atol = 0.002) update_se_hessian!(partable_mean, solution_ml) test_estimates( @@ -324,19 +297,7 @@ end @testset "fitmeasures/se_ls_mean" begin solution_ls = fit(semoptimizer, model_ls) - fm = fit_measures(solution_ls) - test_fitmeasures( - fm, - solution_lav[:fitmeasures_ls_mean]; - atol = 1e-3, - fitmeasure_names = fitmeasure_names_ls, - ) - @test ismissing(fm[:AIC]) && ismissing(fm[:BIC]) && ismissing(fm[:minus2ll]) - test_fitmeasures( - Dict(:CFI => CFI(solution_ls)), - solution_lav[:fitmeasures_ls_mean]; - fitmeasure_names = Dict(:CFI => "cfi"), - ) + test_fitmeasures(solution_ls, solution_lav[:fitmeasures_ls_mean]; atol = 1e-3) @suppress update_se_hessian!(partable_mean, solution_ls) test_estimates( @@ -410,6 +371,8 @@ if !ismissing(spec_varonly) loss = SemFIML, meanstructure = true, ) +else + model_varonly = nothing end ############################################################################################ @@ -446,19 +409,16 @@ end @testset "fitmeasures/se_fiml" begin solution_ml = fit(semoptimizer, model_ml) + solution_varonly = + !isnothing(model_varonly) ? fit(semoptimizer, model_varonly) : nothing test_fitmeasures( - fit_measures(solution_ml), + solution_ml, solution_lav[:fitmeasures_fiml]; + fitted_baseline = solution_varonly, + fitmeasures = !isnothing(solution_varonly) ? SEM.DEFAULT_FIT_MEASURES : + filter(!=(CFI), SEM.DEFAULT_FIT_MEASURES), atol = 1e-3, ) - if !ismissing(spec_varonly) - solution_varonly = fit(semoptimizer, model_varonly) - test_fitmeasures( - Dict(:CFI => CFI(solution_ml, solution_varonly)), - solution_lav[:fitmeasures_fiml]; - fitmeasure_names = Dict(:CFI => "cfi"), - ) - end update_se_hessian!(partable_mean, solution_ml) test_estimates( From 05abcd9c7f2170d8fe89cb3ffd1edc3d2f5035da Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Sun, 22 Mar 2026 23:51:27 -0700 Subject: [PATCH 18/18] test/multigroup: small tweaks --- test/examples/multigroup/build_models.jl | 2 +- test/examples/multigroup/multigroup.jl | 17 +++++++++-------- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl index 71128760..6c22a453 100644 --- a/test/examples/multigroup/build_models.jl +++ b/test/examples/multigroup/build_models.jl @@ -240,7 +240,7 @@ if !isnothing(specification_miss_g1) model_ml_varonly = Sem( specification = partable_varonly, - data = dat_missing, + data = dat_miss, semterm_column = :school, loss = SemFIML, observed = SemObservedMissing, diff --git a/test/examples/multigroup/multigroup.jl b/test/examples/multigroup/multigroup.jl index c8ae8c1f..35fe20e6 100644 --- a/test/examples/multigroup/multigroup.jl +++ b/test/examples/multigroup/multigroup.jl @@ -10,17 +10,18 @@ const SEM = StructuralEquationModels include(joinpath(chop(dirname(pathof(SEM)), tail = 3), "test/examples/helper.jl")) dat = example_data("holzinger_swineford") -dat_missing = example_data("holzinger_swineford_missing") -solution_lav = example_data("holzinger_swineford_solution") +dat.school = Symbol.(replace.(dat.school, "-" => "_")) + +dat_miss = example_data("holzinger_swineford_missing") +dat_miss.school = Symbol.(replace.(dat_miss.school, "-" => "_")) -dat_g1 = dat[dat.school .== "Pasteur", :] -dat_g2 = dat[dat.school .== "Grant-White", :] +solution_lav = example_data("holzinger_swineford_solution") -dat_miss_g1 = dat_missing[dat_missing.school .== "Pasteur", :] -dat_miss_g2 = dat_missing[dat_missing.school .== "Grant-White", :] +dat_g1 = dat[dat.school .== :Pasteur, :] +dat_g2 = dat[dat.school .== :Grant_White, :] -dat.school = ifelse.(dat.school .== "Pasteur", :Pasteur, :Grant_White) -dat_missing.school = ifelse.(dat_missing.school .== "Pasteur", :Pasteur, :Grant_White) +dat_miss_g1 = dat_miss[dat_miss.school .== :Pasteur, :] +dat_miss_g2 = dat_miss[dat_miss.school .== :Grant_White, :] ############################################################################################ ### specification - RAMMatrices