diff --git a/docs/src/tutorials/concept.md b/docs/src/tutorials/concept.md index 035144d62..2b453925a 100644 --- a/docs/src/tutorials/concept.md +++ b/docs/src/tutorials/concept.md @@ -50,8 +50,8 @@ Available loss functions are - [`SemRidge`](@ref): ridge regularization ## The optimizer part aka `SemOptimizer` -The optimizer part of a model connects to the numerical optimization backend used to fit the model. -It can be used to control options like the optimization algorithm, linesearch, stopping criteria, etc. +The optimizer part of a model connects to the numerical optimization backend used to fit the model. +It can be used to control options like the optimization algorithm, linesearch, stopping criteria, etc. There are currently three available backends, [`SemOptimizerOptim`](@ref) connecting to the [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl) backend, [`SemOptimizerNLopt`](@ref) connecting to the [NLopt.jl](https://github.com/JuliaOpt/NLopt.jl) backend and [`SemOptimizerProximal`](@ref) connecting to [ProximalAlgorithms.jl](https://github.com/JuliaFirstOrder/ProximalAlgorithms.jl). For more information about the available options see also the tutorials about [Using Optim.jl](@ref) and [Using NLopt.jl](@ref), as well as [Constrained optimization](@ref) and [Regularization](@ref) . diff --git a/docs/src/tutorials/constraints/constraints.md b/docs/src/tutorials/constraints/constraints.md index 3d9b3d1a1..c433240a9 100644 --- a/docs/src/tutorials/constraints/constraints.md +++ b/docs/src/tutorials/constraints/constraints.md @@ -38,7 +38,7 @@ end partable = ParameterTable( graph, - latent_vars = latent_vars, + latent_vars = latent_vars, observed_vars = observed_vars) data = example_data("political_democracy") @@ -64,7 +64,7 @@ Let's introduce some constraints: (Of course those constaints only serve an illustratory purpose.) -We first need to get the indices of the respective parameters that are invoved in the constraints. +We first need to get the indices of the respective parameters that are invoved in the constraints. We can look up their labels in the output above, and retrieve their indices as ```@example constraints @@ -112,7 +112,7 @@ end ``` If the algorithm needs gradients at an iteration, it will pass the vector `gradient` that is of the same size as the parameters. -With `if length(gradient) > 0` we check if the algorithm needs gradients, and if it does, we fill the `gradient` vector with the gradients +With `if length(gradient) > 0` we check if the algorithm needs gradients, and if it does, we fill the `gradient` vector with the gradients of the constraint w.r.t. the parameters. In NLopt, vector-valued constraints are also possible, but we refer to the documentation for that. @@ -134,7 +134,7 @@ constrained_optimizer = SemOptimizerNLopt( ``` As you see, the equality constraints and inequality constraints are passed as keyword arguments, and the bounds are passed as options for the (outer) optimization algorithm. -Additionally, for equality and inequality constraints, a feasibility tolerance can be specified that controls if a solution can be accepted, even if it violates the constraints by a small amount. +Additionally, for equality and inequality constraints, a feasibility tolerance can be specified that controls if a solution can be accepted, even if it violates the constraints by a small amount. Especially for equality constraints, it is recommended to allow for a small positive tolerance. In this example, we set both tolerances to `1e-8`. @@ -142,7 +142,7 @@ In this example, we set both tolerances to `1e-8`. We have often observed that the default convergence criteria in NLopt lead to non-convergence flags. Indeed, this example does not convergence with default criteria. As you see above, we used a realively liberal absolute tolerance in the optimization parameters of 1e-4. - This should not be a problem in most cases, as the sampling variance in (almost all) structural equation models + This should not be a problem in most cases, as the sampling variance in (almost all) structural equation models should lead to uncertainty in the parameter estimates that are orders of magnitude larger. We nontheless recommend choosing a convergence criterion with care (i.e. w.r.t. the scale of your parameters), inspecting the solutions for plausibility, and comparing them to unconstrained solutions. @@ -162,14 +162,14 @@ As you can see, the optimizer converged (`:XTOL_REACHED`) and investigating the update_partable!( partable, :estimate_constr, - model_fit_constrained, - solution(model_fit_constrained), - ) + model_fit_constrained, + solution(model_fit_constrained), +) details(partable) ``` -As we can see, the constrained solution is very close to the original solution (compare the columns estimate and estimate_constr), with the difference that the constrained parameters fulfill their constraints. +As we can see, the constrained solution is very close to the original solution (compare the columns estimate and estimate_constr), with the difference that the constrained parameters fulfill their constraints. As all parameters are estimated simultaneously, it is expexted that some unconstrained parameters are also affected (e.g., the constraint on `dem60 → y2` leads to a higher estimate of the residual variance `y2 ↔ y2`). ## Using the Optim.jl backend diff --git a/docs/src/tutorials/construction/build_by_parts.md b/docs/src/tutorials/construction/build_by_parts.md index 6b6b59ac9..680e28804 100644 --- a/docs/src/tutorials/construction/build_by_parts.md +++ b/docs/src/tutorials/construction/build_by_parts.md @@ -40,7 +40,7 @@ end partable = ParameterTable( graph, - latent_vars = lat_vars, + latent_vars = lat_vars, observed_vars = obs_vars) ``` diff --git a/docs/src/tutorials/fitting/fitting.md b/docs/src/tutorials/fitting/fitting.md index fff06abaa..d7353c9f9 100644 --- a/docs/src/tutorials/fitting/fitting.md +++ b/docs/src/tutorials/fitting/fitting.md @@ -7,19 +7,19 @@ model_fit = fit(model) # output -Fitted Structural Equation Model -=============================================== ---------------------- Model ------------------- +Fitted Structural Equation Model +=============================================== +--------------------- Model ------------------- -Structural Equation Model -- Loss Functions +Structural Equation Model +- Loss Functions SemML -- Fields - observed: SemObservedData - implied: RAM - optimizer: SemOptimizerOptim +- Fields + observed: SemObservedData + implied: RAM + optimizer: SemOptimizerOptim -------------- Optimization result ------------- +------------- Optimization result ------------- * Status: success diff --git a/docs/src/tutorials/meanstructure.md b/docs/src/tutorials/meanstructure.md index b2da5029a..4e6d2a36a 100644 --- a/docs/src/tutorials/meanstructure.md +++ b/docs/src/tutorials/meanstructure.md @@ -40,7 +40,7 @@ end partable = ParameterTable( graph, - latent_vars = latent_vars, + latent_vars = latent_vars, observed_vars = observed_vars) ``` @@ -78,7 +78,7 @@ end partable = ParameterTable( graph, - latent_vars = latent_vars, + latent_vars = latent_vars, observed_vars = observed_vars) ``` diff --git a/docs/src/tutorials/regularization/regularization.md b/docs/src/tutorials/regularization/regularization.md index bcb9b7747..2b2c6df30 100644 --- a/docs/src/tutorials/regularization/regularization.md +++ b/docs/src/tutorials/regularization/regularization.md @@ -2,7 +2,7 @@ ## Setup -For ridge regularization, you can simply use `SemRidge` as an additional loss function +For ridge regularization, you can simply use `SemRidge` as an additional loss function (for example, a model with the loss functions `SemML` and `SemRidge` corresponds to ridge-regularized maximum likelihood estimation). For lasso, elastic net and (far) beyond, you can load the `ProximalAlgorithms.jl` and `ProximalOperators.jl` packages alongside `StructuralEquationModels`: @@ -22,7 +22,7 @@ using StructuralEquationModels, ProximalAlgorithms, ProximalOperators ## `SemOptimizerProximal` To estimate regularized models, we provide a "building block" for the optimizer part, called `SemOptimizerProximal`. -It connects our package to the [`ProximalAlgorithms.jl`](https://github.com/JuliaFirstOrder/ProximalAlgorithms.jl) optimization backend, providing so-called proximal optimization algorithms. +It connects our package to the [`ProximalAlgorithms.jl`](https://github.com/JuliaFirstOrder/ProximalAlgorithms.jl) optimization backend, providing so-called proximal optimization algorithms. Those can handle, amongst other things, various forms of regularization. It can be used as @@ -32,7 +32,7 @@ SemOptimizerProximal( algorithm = ProximalAlgorithms.PANOC(), operator_g, operator_h = nothing - ) +) ``` The proximal operator (aka the regularization function) can be passed as `operator_g`. @@ -69,7 +69,7 @@ end partable = ParameterTable( graph, - latent_vars = latent_vars, + latent_vars = latent_vars, observed_vars = observed_vars ) @@ -85,7 +85,7 @@ We labeled the covariances between the items because we want to regularize those ```@example reg ind = getindex.( - [param_indices(model)], + Ref(param_indices(model)), [:cov_15, :cov_24, :cov_26, :cov_37, :cov_48, :cov_68]) ``` @@ -107,7 +107,7 @@ and use `SemOptimizerProximal`. ```@example reg optimizer_lasso = SemOptimizerProximal( operator_g = NormL1(λ) - ) +) model_lasso = Sem( specification = partable, @@ -158,7 +158,7 @@ prox_operator = SlicedSeparableSum((NormL0(20.0), NormL1(0.02), NormL0(0.0)), ([ model_mixed = Sem( specification = partable, - data = data, + data = data, ) fit_mixed = fit(model_mixed; engine = :Proximal, operator_g = prox_operator) diff --git a/docs/src/tutorials/specification/graph_interface.md b/docs/src/tutorials/specification/graph_interface.md index 75e1d1b6d..62eeef00b 100644 --- a/docs/src/tutorials/specification/graph_interface.md +++ b/docs/src/tutorials/specification/graph_interface.md @@ -1,6 +1,6 @@ # Graph interface -## Workflow +## Workflow As discussed before, when using the graph interface, you can specify your model as a graph ```julia @@ -17,7 +17,7 @@ lat_vars = ... partable = ParameterTable( graph, - latent_vars = lat_vars, + latent_vars = lat_vars, observed_vars = obs_vars) model = Sem( @@ -32,7 +32,7 @@ In general, there are two different types of parameters: **directed** and **indi We allow multiple variables on both sides of an arrow, for example `x → [y z]` or `[a b] → [c d]`. The later specifies element wise edges; that is its the same as `a → c; b → d`. If you want edges corresponding to the cross-product, we have the double lined arrow `[a b] ⇒ [c d]`, corresponding to `a → c; a → d; b → c; b → d`. The undirected arrows ↔ (element-wise) and ⇔ (crossproduct) behave the same way. !!! note "Unicode symbols in julia" - The `→` symbol is a unicode symbol allowed in julia (among many others; see this [list](https://docs.julialang.org/en/v1/manual/unicode-input/)). You can enter it in the julia REPL or the vscode IDE by typing `\to` followed by hitting `tab`. Similarly, + The `→` symbol is a unicode symbol allowed in julia (among many others; see this [list](https://docs.julialang.org/en/v1/manual/unicode-input/)). You can enter it in the julia REPL or the vscode IDE by typing `\to` followed by hitting `tab`. Similarly, - `←` = `\leftarrow`, - `↔` = `\leftrightarrow`, - `⇒` = `\Rightarrow`, @@ -54,7 +54,7 @@ graph = @StenoGraph begin ξ₃ ↔ fixed(1.0)*ξ₃ end ``` -would +would - fix the directed effects from `ξ₁` to `x1` and from `ξ₂` to `x2` to `1` - leave the directed effect from `ξ₃` to `x7` free but instead restrict the variance of `ξ₃` to `1` - give the effect from `ξ₁` to `x3` the label `:a` (which can be convenient later if you want to retrieve information from your model about that specific parameter) @@ -66,7 +66,7 @@ As you saw above and in the [A first model](@ref) example, the graph object need ```julia partable = ParameterTable( graph, - latent_vars = lat_vars, + latent_vars = lat_vars, observed_vars = obs_vars) ``` @@ -85,7 +85,7 @@ The variable names (`:x1`) have to be symbols, the syntax `:something` creates a _(lat_vars) ⇔ _(lat_vars) end ``` -creates undirected effects coresponding to +creates undirected effects coresponding to 1. the variances of all observed variables and 2. the variances plus covariances of all latent variables So if you want to work with a subset of variables, simply specify a vector of symbols `somevars = [...]`, and inside the graph specification, refer to them as `_(somevars)`. diff --git a/docs/src/tutorials/specification/parameter_table.md b/docs/src/tutorials/specification/parameter_table.md index c328a3b1a..62c45c9a4 100644 --- a/docs/src/tutorials/specification/parameter_table.md +++ b/docs/src/tutorials/specification/parameter_table.md @@ -5,5 +5,5 @@ As lavaan also uses parameter tables to store model specifications, we are worki ## Convert from and to RAMMatrices -To convert a RAMMatrices object to a ParameterTable, simply use `partable = ParameterTable(rammatrices)`. +To convert a RAMMatrices object to a ParameterTable, simply use `partable = ParameterTable(rammatrices)`. To convert an object of type `ParameterTable` to RAMMatrices, you can use `ram_matrices = RAMMatrices(partable)`. \ No newline at end of file diff --git a/docs/src/tutorials/specification/ram_matrices.md b/docs/src/tutorials/specification/ram_matrices.md index abe76ea6f..2730ff4bf 100644 --- a/docs/src/tutorials/specification/ram_matrices.md +++ b/docs/src/tutorials/specification/ram_matrices.md @@ -1,6 +1,6 @@ # RAMMatrices interface -Models can also be specified by an object of type `RAMMatrices`. +Models can also be specified by an object of type `RAMMatrices`. The RAM (reticular action model) specification corresponds to three matrices; the `A` matrix containing all directed parameters, the `S` matrix containing all undirected parameters, and the `F` matrix filtering out latent variables from the model implied covariance. The model implied covariance matrix for the observed variables of a SEM is then computed as @@ -56,9 +56,9 @@ A =[0 0 0 0 0 0 0 0 0 0 0 1.0 0 0 θ = Symbol.(:θ, 1:31) spec = RAMMatrices(; - A = A, - S = S, - F = F, + A = A, + S = S, + F = F, param_labels = θ, vars = [:x1, :x2, :x3, :y1, :y2, :y3, :y4, :y5, :y6, :y7, :y8, :ind60, :dem60, :dem65] ) @@ -71,9 +71,9 @@ model = Sem( Let's look at this step by step: -First, we specify the `A`, `S` and `F`-Matrices. -For a free parameter, we write a `Symbol` like `:θ1` (or any other symbol we like) to the corresponding place in the respective matrix, to constrain parameters to be equal we just use the same `Symbol` in the respective entries. -To fix a parameter (as in the `A`-Matrix above), we just write down the number we want to fix it to. +First, we specify the `A`, `S` and `F`-Matrices. +For a free parameter, we write a `Symbol` like `:θ1` (or any other symbol we like) to the corresponding place in the respective matrix, to constrain parameters to be equal we just use the same `Symbol` in the respective entries. +To fix a parameter (as in the `A`-Matrix above), we just write down the number we want to fix it to. All other entries are 0. Second, we specify a vector of symbols containing our parameters: @@ -82,14 +82,14 @@ Second, we specify a vector of symbols containing our parameters: θ = Symbol.(:θ, 1:31) ``` -Third, we construct an object of type `RAMMatrices`, passing our matrices and parameters, as well as the column names of our matrices. +Third, we construct an object of type `RAMMatrices`, passing our matrices and parameters, as well as the column names of our matrices. Those are quite important, as they will be used to rearrange your data to match it to your `RAMMatrices` specification. ```julia spec = RAMMatrices(; - A = A, - S = S, - F = F, + A = A, + S = S, + F = F, param_labels = θ, vars = [:x1, :x2, :x3, :y1, :y2, :y3, :y4, :y5, :y6, :y7, :y8, :ind60, :dem60, :dem65] ) @@ -109,7 +109,7 @@ According to the RAM, model implied mean values of the observed variables are co ```math \mu = F(I-A)^{-1}M ``` -where `M` is a vector of mean parameters. To estimate the means of the observed variables in our example (and set the latent means to `0`), we would specify the model just as before but add +where `M` is a vector of mean parameters. To estimate the means of the observed variables in our example (and set the latent means to `0`), we would specify the model just as before but add ```julia ... @@ -128,5 +128,5 @@ spec = RAMMatrices(; ## Convert from and to ParameterTables -To convert a RAMMatrices object to a ParameterTable, simply use `partable = ParameterTable(ram_matrices)`. +To convert a RAMMatrices object to a ParameterTable, simply use `partable = ParameterTable(ram_matrices)`. To convert an object of type `ParameterTable` to RAMMatrices, you can use `ram_matrices = RAMMatrices(partable)`. \ No newline at end of file diff --git a/ext/SEMNLOptExt/NLopt.jl b/ext/SEMNLOptExt/NLopt.jl index c5e0ad6cb..27bc30039 100644 --- a/ext/SEMNLOptExt/NLopt.jl +++ b/ext/SEMNLOptExt/NLopt.jl @@ -19,7 +19,7 @@ function SemOptimizerNLopt(; applicable(iterate, equality_constraints) && !isa(equality_constraints, NamedTuple) || (equality_constraints = [equality_constraints]) applicable(iterate, inequality_constraints) && - !isa(inequality_constraints, NamedTuple) || + !isa(inequality_constraints, NamedTuple) || (inequality_constraints = [inequality_constraints]) return SemOptimizerNLopt( algorithm, diff --git a/src/additional_functions/helper.jl b/src/additional_functions/helper.jl index 6cbcb0573..d6a1fc6c8 100644 --- a/src/additional_functions/helper.jl +++ b/src/additional_functions/helper.jl @@ -21,7 +21,7 @@ function batch_inv!(fun, model) end # computes A*S*B -> C, where ind gives the entries of S that are 1 -function sparse_outer_mul!(C, A, B, ind) +function sparse_outer_mul!(C, A, B, ind) fill!(C, 0.0) for i in 1:length(ind) BLAS.ger!(1.0, A[:, ind[i][1]], B[ind[i][2], :], C) @@ -29,14 +29,14 @@ function sparse_outer_mul!(C, A, B, ind) end # computes A*∇m, where ∇m ind gives the entries of ∇m that are 1 -function sparse_outer_mul!(C, A, ind) +function sparse_outer_mul!(C, A, ind) fill!(C, 0.0) @views C .= sum(A[:, ind], dims = 2) return C end # computes A*S*B -> C, where ind gives the entries of S that are 1 -function sparse_outer_mul!(C, A, B::Vector, ind) +function sparse_outer_mul!(C, A, B::Vector, ind) fill!(C, 0.0) @views @inbounds for i in 1:length(ind) C .+= B[ind[i][2]] .* A[:, ind[i][1]] diff --git a/src/additional_functions/params_array.jl b/src/additional_functions/params_array.jl index 1031e349e..8cd89032e 100644 --- a/src/additional_functions/params_array.jl +++ b/src/additional_functions/params_array.jl @@ -199,11 +199,7 @@ materialize!( kwargs..., ) = materialize!(parent(dest), src, params; kwargs...) -function sparse_materialize( - ::Type{T}, - arr::ParamsMatrix, - params::AbstractVector, -) where {T} +function sparse_materialize(::Type{T}, arr::ParamsMatrix, params::AbstractVector) where {T} nparams(arr) == length(params) || throw( DimensionMismatch( "Number of values ($(length(params))) does not match the number of parameter ($(nparams(arr)))", @@ -251,8 +247,8 @@ sparse_gradient(arr::ParamsArray{T}) where {T} = sparse_gradient(T, arr) # range of parameters that are referenced in the matrix function params_range(arr::ParamsArray; allow_gaps::Bool = false) - first_i = findfirst(i -> arr.param_ptr[i+1] > arr.param_ptr[i], 1:nparams(arr)-1) - last_i = findlast(i -> arr.param_ptr[i+1] > arr.param_ptr[i], 1:nparams(arr)-1) + first_i = findfirst(i -> arr.param_ptr[i+1] > arr.param_ptr[i], 1:(nparams(arr)-1)) + last_i = findlast(i -> arr.param_ptr[i+1] > arr.param_ptr[i], 1:(nparams(arr)-1)) if !allow_gaps && !isnothing(first_i) && !isnothing(last_i) for i in first_i:last_i diff --git a/src/additional_functions/simulation.jl b/src/additional_functions/simulation.jl index a787516b2..4839bc272 100644 --- a/src/additional_functions/simulation.jl +++ b/src/additional_functions/simulation.jl @@ -48,7 +48,7 @@ replace_observed(model::AbstractSemSingle; kwargs...) = replace_observed(model, typeof(observed(model)).name.wrapper; kwargs...) function replace_observed(model::AbstractSemSingle, observed_type; kwargs...) - new_observed = observed_type(;kwargs...) + new_observed = observed_type(; kwargs...) kwargs = Dict{Symbol, Any}(kwargs...) # get field types @@ -65,11 +65,7 @@ function replace_observed(model::AbstractSemSingle, observed_type; kwargs...) # update loss new_loss = update_observed(model.loss, new_observed; kwargs...) - return Sem( - new_observed, - new_implied, - new_loss - ) + return Sem(new_observed, new_implied, new_loss) end function update_observed(loss::SemLoss, new_observed; kwargs...) @@ -79,7 +75,6 @@ function update_observed(loss::SemLoss, new_observed; kwargs...) return SemLoss(new_functions, loss.weights) end - function replace_observed( emodel::SemEnsemble; column = :group, @@ -94,16 +89,16 @@ function replace_observed( # 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) + 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) + 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 diff --git a/src/frontend/StatsAPI.jl b/src/frontend/StatsAPI.jl index edd677e34..b00c451af 100644 --- a/src/frontend/StatsAPI.jl +++ b/src/frontend/StatsAPI.jl @@ -13,11 +13,7 @@ Note that the function combines the duplicate occurences of the same parameter in `partable` and will raise an error if the values do not match. """ -function params!( - out::AbstractVector, - partable::ParameterTable, - col::Symbol = :estimate, -) +function params!(out::AbstractVector, partable::ParameterTable, col::Symbol = :estimate) (length(out) == nparams(partable)) || throw( DimensionMismatch( "The length of parameter values vector ($(length(out))) does not match the number of parameters ($(nparams(partable)))", @@ -75,4 +71,8 @@ Synonymous to [`nsamples`](@ref). """ nobs(model::AbstractSem) = nsamples(model) -coeftable(model::AbstractSem; level::Real=0.95) = throw(ArgumentError("StructuralEquationModels does not support the `CoefTable` interface; see [`ParameterTable`](@ref) instead.")) \ No newline at end of file +coeftable(model::AbstractSem; level::Real = 0.95) = throw( + ArgumentError( + "StructuralEquationModels does not support the `CoefTable` interface; see [`ParameterTable`](@ref) instead.", + ), +) diff --git a/src/frontend/fit/fitmeasures/RMSEA.jl b/src/frontend/fit/fitmeasures/RMSEA.jl index b9fff648e..f9dae84ed 100644 --- a/src/frontend/fit/fitmeasures/RMSEA.jl +++ b/src/frontend/fit/fitmeasures/RMSEA.jl @@ -1,15 +1,16 @@ """ - RMSEA(sem_fit::SemFit) + RMSEA(fit::SemFit) Return the RMSEA. """ function RMSEA end -RMSEA(sem_fit::SemFit{Mi, So, St, Mo, O} where {Mi, So, St, Mo <: AbstractSemSingle, O}) = - RMSEA(dof(sem_fit), χ²(sem_fit), nsamples(sem_fit)) +RMSEA(fit::SemFit) = RMSEA(fit, fit.model) -RMSEA(sem_fit::SemFit{Mi, So, St, Mo, O} where {Mi, So, St, Mo <: SemEnsemble, O}) = - sqrt(length(sem_fit.model.sems)) * RMSEA(dof(sem_fit), χ²(sem_fit), nsamples(sem_fit)) +RMSEA(fit::SemFit, model::AbstractSemSingle) = RMSEA(dof(fit), χ²(fit), nsamples(fit)) + +RMSEA(fit::SemFit, model::SemEnsemble) = + sqrt(length(model.sems)) * RMSEA(dof(fit), χ²(fit), nsamples(fit)) function RMSEA(dof, chi2, nsamples) rmsea = (chi2 - dof) / (nsamples * dof) diff --git a/src/frontend/fit/fitmeasures/chi2.jl b/src/frontend/fit/fitmeasures/chi2.jl index 333783f95..dc19467fc 100644 --- a/src/frontend/fit/fitmeasures/chi2.jl +++ b/src/frontend/fit/fitmeasures/chi2.jl @@ -1,89 +1,70 @@ """ - χ²(sem_fit::SemFit) + χ²(fit::SemFit) Return the χ² value. """ -function χ² end +χ²(fit::SemFit) = χ²(fit, fit.model) ############################################################################################ # Single Models ############################################################################################ -# SemFit splices loss functions ------------------------------------------------------------ -χ²(sem_fit::SemFit{Mi, So, St, Mo, O} where {Mi, So, St, Mo <: AbstractSemSingle, O}) = χ²( - sem_fit, - sem_fit.model.observed, - sem_fit.model.implied, - sem_fit.model.loss.functions..., -) +χ²(fit::SemFit, model::AbstractSemSingle) = + sum(loss -> χ²(loss, fit, model), model.loss.functions) # RAM + SemML -χ²(sem_fit::SemFit, observed, imp::Union{RAM, RAMSymbolic}, loss_ml::SemML) = - (nsamples(sem_fit) - 1) * - (sem_fit.minimum - logdet(observed.obs_cov) - nobserved_vars(observed)) +χ²(lossfun::SemML, fit::SemFit, model::AbstractSemSingle) = + (nsamples(fit) - 1) * + (fit.minimum - logdet(obs_cov(observed(model))) - nobserved_vars(observed(model))) # bollen, p. 115, only correct for GLS weight matrix -χ²(sem_fit::SemFit, observed, imp::Union{RAM, RAMSymbolic}, loss_ml::SemWLS) = - (nsamples(sem_fit) - 1) * sem_fit.minimum +χ²(lossfun::SemWLS, fit::SemFit, model::AbstractSemSingle) = + (nsamples(fit) - 1) * fit.minimum # FIML -function χ²(sem_fit::SemFit, observed::SemObservedMissing, imp, loss_ml::SemFIML) - ll_H0 = minus2ll(sem_fit) - ll_H1 = minus2ll(observed) - chi2 = ll_H0 - ll_H1 - return chi2 +function χ²(lossfun::SemFIML, fit::SemFit, model::AbstractSemSingle) + ll_H0 = minus2ll(fit) + ll_H1 = minus2ll(observed(model)) + return ll_H0 - ll_H1 end ############################################################################################ # Collections ############################################################################################ -# SemFit splices loss functions ------------------------------------------------------------ -χ²(sem_fit::SemFit{Mi, So, St, Mo, O} where {Mi, So, St, Mo <: SemEnsemble, O}) = - χ²(sem_fit, sem_fit.model, sem_fit.model.sems[1].loss.functions[1]) +function χ²(fit::SemFit, models::SemEnsemble) + isempty(models.sems) && return 0.0 -function χ²(sem_fit::SemFit, model::SemEnsemble, lossfun::L) where {L <: SemWLS} - check_ensemble_length(model) - check_lossfun_types(model, L) - return (nsamples(model) - 1) * sem_fit.minimum -end + lossfun = models.sems[1].loss.functions[1] + # check that all models use the same single loss function + L = typeof(lossfun) + for (i, sem) in enumerate(models.sems) + if length(sem.loss.functions) > 1 + @error "Model for group #$i has $(length(sem.loss.functions)) loss functions. Only the single one is supported" + end + cur_lossfun = sem.loss.functions[1] + if !isa(cur_lossfun, L) + @error "Loss function for group #$i model is $(typeof(cur_lossfun)), expected $L. Heterogeneous loss functions are not supported" + end + end -function χ²(sem_fit::SemFit, model::SemEnsemble, lossfun::L) where {L <: SemML} - check_ensemble_length(model) - check_lossfun_types(model, L) - F_G = sem_fit.minimum - F_G -= sum([ - w * (logdet(m.observed.obs_cov) + nobserved_vars(m.observed)) for - (w, m) in zip(model.weights, model.sems) - ]) - return (nsamples(model) - 1) * F_G + return χ²(lossfun, fit, models) end -function χ²(sem_fit::SemFit, model::SemEnsemble, lossfun::L) where {L <: SemFIML} - check_ensemble_length(model) - check_lossfun_types(model, L) - - ll_H0 = minus2ll(sem_fit) - ll_H1 = sum(minus2ll.(observed.(models(model)))) - chi2 = ll_H0 - ll_H1 - - return chi2 +function χ²(lossfun::SemWLS, fit::SemFit, models::SemEnsemble) + return (nsamples(models) - 1) * fit.minimum end -function check_ensemble_length(model) - for sem in model.sems - if length(sem.loss.functions) > 1 - @error "A model for one of the groups contains multiple loss functions." - end +function χ²(lossfun::SemML, fit::SemFit, models::SemEnsemble) + G = sum(zip(models.weights, models.sems)) do (w, model) + data = observed(model) + w * (logdet(obs_cov(data)) + nobserved_vars(data)) end + return (nsamples(models) - 1) * (fit.minimum - G) end -function check_lossfun_types(model, type) - for sem in model.sems - for lossfun in sem.loss.functions - if !isa(lossfun, type) - @error "Your model(s) contain multiple lossfunctions with differing types." - end - end - end +function χ²(lossfun::SemFIML, fit::SemFit, models::SemEnsemble) + ll_H0 = minus2ll(fit) + ll_H1 = sum(minus2ll ∘ observed, models.sems) + return ll_H0 - ll_H1 end diff --git a/src/frontend/fit/fitmeasures/minus2ll.jl b/src/frontend/fit/fitmeasures/minus2ll.jl index ab4d24e53..9b211fb44 100644 --- a/src/frontend/fit/fitmeasures/minus2ll.jl +++ b/src/frontend/fit/fitmeasures/minus2ll.jl @@ -9,30 +9,31 @@ function minus2ll end # Single Models ############################################################################################ -# SemFit splices loss functions ------------------------------------------------------------ -minus2ll( - sem_fit::SemFit{Mi, So, St, Mo, O} where {Mi, So, St, Mo <: AbstractSemSingle, O}, -) = minus2ll( - sem_fit, - sem_fit.model.observed, - sem_fit.model.implied, - sem_fit.model.loss.functions..., -) - -minus2ll(sem_fit::SemFit, obs, imp, args...) = minus2ll(sem_fit.minimum, obs, imp, args...) +minus2ll(fit::SemFit) = minus2ll(fit, fit.model) + +function minus2ll(fit::SemFit, model::AbstractSemSingle) + minimum = objective(model, fit.solution) + return minus2ll(minimum, model) +end + +minus2ll(minimum::Number, model::AbstractSemSingle) = + sum(lossfun -> minus2ll(lossfun, minimum, model), model.loss.functions) # SemML ------------------------------------------------------------------------------------ -minus2ll(minimum::Number, obs, imp::Union{RAM, RAMSymbolic}, loss_ml::SemML) = - nsamples(obs) * (minimum + log(2π) * nobserved_vars(obs)) +function minus2ll(lossfun::SemML, minimum::Number, model::AbstractSemSingle) + obs = observed(model) + return nsamples(obs) * (minimum + log(2π) * nobserved_vars(obs)) +end # WLS -------------------------------------------------------------------------------------- -minus2ll(minimum::Number, obs, imp::Union{RAM, RAMSymbolic}, loss_ml::SemWLS) = missing +minus2ll(lossfun::SemWLS, minimum::Number, model::AbstractSemSingle) = missing # compute likelihood for missing data - H0 ------------------------------------------------- # -2ll = (∑ log(2π)*(nᵢ + mᵢ)) + F*n -function minus2ll(minimum::Number, observed, imp::Union{RAM, RAMSymbolic}, loss_ml::SemFIML) - F = minimum * nsamples(observed) - F += log(2π) * sum(pat -> nsamples(pat) * nmeasured_vars(pat), observed.patterns) +function minus2ll(lossfun::SemFIML, minimum::Number, model::AbstractSemSingle) + obs = observed(model)::SemObservedMissing + F = minimum * nsamples(obs) + F += log(2π) * sum(pat -> nsamples(pat) * nmeasured_vars(pat), obs.patterns) return F end @@ -66,16 +67,4 @@ end # Collection ############################################################################################ -minus2ll(minimum, model::AbstractSemSingle) = - minus2ll(minimum, model.observed, model.implied, model.loss.functions...) - -function minus2ll( - sem_fit::SemFit{Mi, So, St, Mo, O} where {Mi, So, St, Mo <: SemEnsemble, O}, -) - m2ll = 0.0 - for sem in sem_fit.model.sems - minimum = objective!(sem, sem_fit.solution) - m2ll += minus2ll(minimum, sem) - end - return m2ll -end +minus2ll(fit::SemFit, model::SemEnsemble) = sum(Base.Fix1(minus2ll, fit), model.sems) diff --git a/src/frontend/fit/fitmeasures/p.jl b/src/frontend/fit/fitmeasures/p.jl index 8c69d5ec2..da9bedaf6 100644 --- a/src/frontend/fit/fitmeasures/p.jl +++ b/src/frontend/fit/fitmeasures/p.jl @@ -3,4 +3,4 @@ Return the p value computed from the χ² test statistic. """ -p_value(sem_fit::SemFit) = 1 - cdf(Chisq(dof(sem_fit)), χ²(sem_fit)) +p_value(sem_fit::SemFit) = ccdf(Chisq(dof(sem_fit)), χ²(sem_fit)) diff --git a/src/frontend/fit/standard_errors/bootstrap.jl b/src/frontend/fit/standard_errors/bootstrap.jl index 4b3e302bb..bebad2935 100644 --- a/src/frontend/fit/standard_errors/bootstrap.jl +++ b/src/frontend/fit/standard_errors/bootstrap.jl @@ -14,8 +14,7 @@ function se_bootstrap( data = nothing, specification = nothing, kwargs..., - ) where {Mi, So, St, Mo <: AbstractSemSingle, O} - +) where {Mi, So, St, Mo <: AbstractSemSingle, O} if isnothing(data) data = samples(observed(model(semfit))) end @@ -67,11 +66,10 @@ function se_bootstrap( data = nothing, specification = nothing, kwargs..., - ) where {Mi, So, St, Mo <: SemEnsemble, O} - +) where {Mi, So, St, Mo <: SemEnsemble, O} models = semfit.model.sems groups = semfit.model.groups - + if isnothing(data) data = Dict(g => samples(observed(m)) for (g, m) in zip(groups, models)) end diff --git a/src/frontend/specification/EnsembleParameterTable.jl b/src/frontend/specification/EnsembleParameterTable.jl index 14169dd94..227e902ff 100644 --- a/src/frontend/specification/EnsembleParameterTable.jl +++ b/src/frontend/specification/EnsembleParameterTable.jl @@ -30,7 +30,7 @@ function EnsembleParameterTable( param_labels = if isnothing(param_labels) # collect all SEM parameters in ensemble if not specified # and apply the set to all partables - unique(mapreduce(SEM.param_labels, vcat, values(spec_ensemble), init = Vector{Symbol}())) + mapreduce(SEM.param_labels, vcat, values(spec_ensemble), init = Symbol[]) |> unique else copy(param_labels) end diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl index 5f8aff1eb..600c9058c 100644 --- a/src/frontend/specification/ParameterTable.jl +++ b/src/frontend/specification/ParameterTable.jl @@ -34,7 +34,9 @@ function ParameterTable( latent_vars::Union{AbstractVector{Symbol}, Nothing} = nothing, param_labels::Union{AbstractVector{Symbol}, Nothing} = nothing, ) - param_labels = isnothing(param_labels) ? unique!(filter(!=(:const), columns[:label])) : copy(param_labels) + param_labels = + isnothing(param_labels) ? unique!(filter(!=(:const), columns[:label])) : + copy(param_labels) check_param_labels(param_labels, columns[:label]) return ParameterTable( columns, @@ -400,7 +402,6 @@ function update_se_hessian!( return update_partable!(partable, :se, param_labels(fit), se) end - """ lavaan_params!(out::AbstractVector, partable_lav, partable::ParameterTable, @@ -449,8 +450,8 @@ function lavaan_params!( lav_ind = findallrows( r -> r[:lhs] == String(to) && - r[:op] == "~1" && - (isnothing(lav_group) || r[:group] == lav_group), + r[:op] == "~1" && + (isnothing(lav_group) || r[:group] == lav_group), partable_lav, ) else @@ -469,20 +470,20 @@ function lavaan_params!( lav_ind = findallrows( r -> ( - (r[:lhs] == String(from) && r[:rhs] == String(to)) || - (r[:lhs] == String(to) && r[:rhs] == String(from)) - ) && - r[:op] == lav_type && - (isnothing(lav_group) || r[:group] == lav_group), + (r[:lhs] == String(from) && r[:rhs] == String(to)) || + (r[:lhs] == String(to) && r[:rhs] == String(from)) + ) && + r[:op] == lav_type && + (isnothing(lav_group) || r[:group] == lav_group), partable_lav, ) else lav_ind = findallrows( r -> r[:lhs] == String(from) && - r[:rhs] == String(to) && - r[:op] == lav_type && - (isnothing(lav_group) || r[:group] == lav_group), + r[:rhs] == String(to) && + r[:op] == lav_type && + (isnothing(lav_group) || r[:group] == lav_group), partable_lav, ) end @@ -535,10 +536,4 @@ lavaan_params( partable::ParameterTable, lav_col::Symbol = :est, lav_group = nothing, -) = lavaan_params!( - fill(NaN, nparams(partable)), - partable_lav, - partable, - lav_col, - lav_group, -) +) = lavaan_params!(fill(NaN, nparams(partable)), partable_lav, partable, lav_col, lav_group) diff --git a/src/frontend/specification/checks.jl b/src/frontend/specification/checks.jl index 5ef41c59d..2d00be26d 100644 --- a/src/frontend/specification/checks.jl +++ b/src/frontend/specification/checks.jl @@ -4,8 +4,11 @@ function check_param_labels( param_refs::Union{AbstractVector{Symbol}, Nothing}, ) dup_param_labels = nonunique(param_labels) - isempty(dup_param_labels) || - throw(ArgumentError("Duplicate parameter labels detected: $(join(dup_param_labels, ", "))")) + isempty(dup_param_labels) || throw( + ArgumentError( + "Duplicate parameter labels detected: $(join(dup_param_labels, ", "))", + ), + ) any(==(:const), param_labels) && throw(ArgumentError("Parameters constain reserved :const name")) diff --git a/src/frontend/specification/documentation.jl b/src/frontend/specification/documentation.jl index e46620fbc..a3a8d2659 100644 --- a/src/frontend/specification/documentation.jl +++ b/src/frontend/specification/documentation.jl @@ -37,7 +37,6 @@ Return the vector of parameter labels (in the same order as [`params`](@ref)). """ param_labels(spec::SemSpecification) = spec.param_labels - """ `ParameterTable`s contain the specification of a structural equation model. diff --git a/src/implied/RAM/generic.jl b/src/implied/RAM/generic.jl index fd2ef61b5..6d4a18194 100644 --- a/src/implied/RAM/generic.jl +++ b/src/implied/RAM/generic.jl @@ -152,7 +152,7 @@ function RAM(; F⨉I_A⁻¹, F⨉I_A⁻¹S, I_A, - copy(I_A), + similar(I_A), ∇A, ∇S, ∇M, @@ -163,7 +163,12 @@ end ### methods ############################################################################################ -function update!(targets::EvaluationTargets, implied::RAM, model::AbstractSemSingle, param_labels) +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) if !isnothing(implied.M) diff --git a/src/implied/RAM/symbolic.jl b/src/implied/RAM/symbolic.jl index 9634bfa89..39efe453b 100644 --- a/src/implied/RAM/symbolic.jl +++ b/src/implied/RAM/symbolic.jl @@ -37,13 +37,13 @@ Jacobians (only available in gradient! calls) - `ram.∇Σ` -> ``∂vec(Σ)/∂θᵀ`` - `ram.∇μ` -> ``∂μ/∂θᵀ`` -- `ram.∇Σ_function` -> function to overwrite `∇Σ` in place, - i.e. `∇Σ_function(∇Σ, θ)`. Typically, you do not want to use this but simply +- `∇Σ_eval!(::RAMSymbolic)` -> function to evaluate `∇Σ` in place, + i.e. `∇Σ_eval!(∇Σ, θ)`. Typically, you do not want to use this but simply query `ram.∇Σ`. Hessians The computation of hessians is more involved. -Therefore, we desribe it in the online documentation, +Therefore, we desribe it in the online documentation, and the respective interfaces are omitted here. ## RAM notation @@ -56,23 +56,19 @@ and for models with a meanstructure, the model implied means are computed as \mu = F(I-A)^{-1}M ``` """ -struct RAMSymbolic{MS, F1, F2, F3, A1, A2, A3, S1, S2, S3, V2, F4, A4, F5, A5} <: - SemImpliedSymbolic +struct RAMSymbolic{MS, F1, F2, F3, A1, A2, A3, V2, F4, A4, F5, A5} <: SemImpliedSymbolic meanstruct::MS hessianeval::ExactHessian - Σ_function::F1 - ∇Σ_function::F2 - ∇²Σ_function::F3 + Σ_eval!::F1 + ∇Σ_eval!::F2 + ∇²Σ_eval!::F3 Σ::A1 ∇Σ::A2 ∇²Σ::A3 - Σ_symbolic::S1 - ∇Σ_symbolic::S2 - ∇²Σ_symbolic::S3 ram_matrices::V2 - μ_function::F4 + μ_eval!::F4 μ::A4 - ∇μ_function::F5 + ∇μ_eval!::F5 ∇μ::A5 RAMSymbolic{MS}(args...) where {MS <: MeanStruct} = @@ -111,81 +107,75 @@ function RAMSymbolic(; I_A⁻¹ = neumann_series(A) # Σ - Σ_symbolic = eval_Σ_symbolic(S, I_A⁻¹, F; vech = vech, simplify = simplify_symbolics) - #print(Symbolics.build_function(Σ_symbolic)[2]) - Σ_function = Symbolics.build_function(Σ_symbolic, par, expression = Val{false})[2] - Σ = zeros(size(Σ_symbolic)) - precompile(Σ_function, (typeof(Σ), Vector{Float64})) + Σ_sym = eval_Σ_symbolic(S, I_A⁻¹, F; vech, simplify = simplify_symbolics) + #print(Symbolics.build_function(Σ_sym)[2]) + Σ_eval! = Symbolics.build_function(Σ_sym, par, expression = Val{false})[2] + Σ = zeros(size(Σ_sym)) + precompile(Σ_eval!, (typeof(Σ), Vector{Float64})) # ∇Σ if gradient - ∇Σ_symbolic = Symbolics.sparsejacobian(vec(Σ_symbolic), [par...]) - ∇Σ_function = Symbolics.build_function(∇Σ_symbolic, par, expression = Val{false})[2] - constr = findnz(∇Σ_symbolic) - ∇Σ = sparse(constr[1], constr[2], fill(1.0, nnz(∇Σ_symbolic)), size(∇Σ_symbolic)...) - precompile(∇Σ_function, (typeof(∇Σ), Vector{Float64})) + ∇Σ_sym = Symbolics.sparsejacobian(vec(Σ_sym), [par...]) + ∇Σ_eval! = Symbolics.build_function(∇Σ_sym, par, expression = Val{false})[2] + constr = findnz(∇Σ_sym) + ∇Σ = sparse(constr[1], constr[2], fill(1.0, nnz(∇Σ_sym)), size(∇Σ_sym)...) + precompile(∇Σ_eval!, (typeof(∇Σ), Vector{Float64})) else - ∇Σ_symbolic = nothing - ∇Σ_function = nothing + ∇Σ_eval! = nothing ∇Σ = nothing end if hessian && !approximate_hessian - n_sig = length(Σ_symbolic) - ∇²Σ_symbolic_vec = [Symbolics.sparsehessian(σᵢ, [par...]) for σᵢ in vec(Σ_symbolic)] + n_sig = length(Σ_sym) + ∇²Σ_sym_vec = [Symbolics.sparsehessian(σᵢ, [par...]) for σᵢ in vec(Σ_sym)] @variables J[1:n_sig] - ∇²Σ_symbolic = zeros(Num, n_par, n_par) + ∇²Σ_sym = zeros(Num, n_par, n_par) for i in 1:n_sig - ∇²Σ_symbolic += J[i] * ∇²Σ_symbolic_vec[i] + ∇²Σ_sym += J[i] * ∇²Σ_sym_vec[i] end - ∇²Σ_function = - Symbolics.build_function(∇²Σ_symbolic, J, par, expression = Val{false})[2] + ∇²Σ_eval! = Symbolics.build_function(∇²Σ_sym, J, par, expression = Val{false})[2] ∇²Σ = zeros(n_par, n_par) else - ∇²Σ_symbolic = nothing - ∇²Σ_function = nothing + ∇²Σ_sym = nothing + ∇²Σ_eval! = nothing ∇²Σ = nothing end # μ if meanstructure MS = HasMeanStruct - μ_symbolic = eval_μ_symbolic(M, I_A⁻¹, F; simplify = simplify_symbolics) - μ_function = Symbolics.build_function(μ_symbolic, par, expression = Val{false})[2] - μ = zeros(size(μ_symbolic)) + μ_sym = eval_μ_symbolic(M, I_A⁻¹, F; simplify = simplify_symbolics) + μ_eval! = Symbolics.build_function(μ_sym, par, expression = Val{false})[2] + μ = zeros(size(μ_sym)) if gradient - ∇μ_symbolic = Symbolics.jacobian(μ_symbolic, [par...]) - ∇μ_function = - Symbolics.build_function(∇μ_symbolic, par, expression = Val{false})[2] + ∇μ_sym = Symbolics.jacobian(μ_sym, [par...]) + ∇μ_eval! = Symbolics.build_function(∇μ_sym, par, expression = Val{false})[2] ∇μ = zeros(size(F, 1), size(par, 1)) else - ∇μ_function = nothing + ∇μ_eval! = nothing ∇μ = nothing end else MS = NoMeanStruct - μ_function = nothing + μ_eval! = nothing μ = nothing - ∇μ_function = nothing + ∇μ_eval! = nothing ∇μ = nothing end return RAMSymbolic{MS}( - Σ_function, - ∇Σ_function, - ∇²Σ_function, + Σ_eval!, + ∇Σ_eval!, + ∇²Σ_eval!, Σ, ∇Σ, ∇²Σ, - Σ_symbolic, - ∇Σ_symbolic, - ∇²Σ_symbolic, ram_matrices, - μ_function, + μ_eval!, μ, - ∇μ_function, + ∇μ_eval!, ∇μ, ) end @@ -200,15 +190,15 @@ function update!( model::AbstractSemSingle, par, ) - implied.Σ_function(implied.Σ, par) + implied.Σ_eval!(implied.Σ, par) if MeanStruct(implied) === HasMeanStruct - implied.μ_function(implied.μ, par) + implied.μ_eval!(implied.μ, par) end if is_gradient_required(targets) || is_hessian_required(targets) - implied.∇Σ_function(implied.∇Σ, par) + implied.∇Σ_eval!(implied.∇Σ, par) if MeanStruct(implied) === HasMeanStruct - implied.∇μ_function(implied.∇μ, par) + implied.∇μ_eval!(implied.∇μ, par) end end end diff --git a/src/loss/ML/ML.jl b/src/loss/ML/ML.jl index ec5eb997c..f5a42e3f4 100644 --- a/src/loss/ML/ML.jl +++ b/src/loss/ML/ML.jl @@ -114,10 +114,9 @@ function evaluate!( if HessianEval(semml) === ApproxHessian mul!(hessian, ∇Σ' * kron(Σ⁻¹, Σ⁻¹), ∇Σ, 2, 0) else - ∇²Σ_function! = implied.∇²Σ_function ∇²Σ = implied.∇²Σ # inner - ∇²Σ_function!(∇²Σ, J, par) + implied.∇²Σ_eval!(∇²Σ, J, par) # outer H_outer = kron(2Σ⁻¹ΣₒΣ⁻¹ - Σ⁻¹, Σ⁻¹) mul!(hessian, ∇Σ' * H_outer, ∇Σ) diff --git a/src/loss/WLS/WLS.jl b/src/loss/WLS/WLS.jl index dd5be4874..59508ce16 100644 --- a/src/loss/WLS/WLS.jl +++ b/src/loss/WLS/WLS.jl @@ -127,10 +127,9 @@ function evaluate!( end isnothing(hessian) || (mul!(hessian, ∇σ' * V, ∇σ, 2, 0)) if !isnothing(hessian) && (HessianEval(semwls) === ExactHessian) - ∇²Σ_function! = implied.∇²Σ_function ∇²Σ = implied.∇²Σ J = -2 * (σ₋' * semwls.V)' - ∇²Σ_function!(∇²Σ, J, par) + implied.∇²Σ_eval!(∇²Σ, J, par) hessian .+= ∇²Σ end if MeanStruct(implied) === HasMeanStruct diff --git a/src/observed/EM.jl b/src/observed/EM.jl index 46d0622be..288082ccc 100644 --- a/src/observed/EM.jl +++ b/src/observed/EM.jl @@ -26,7 +26,6 @@ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. =# -# outer function --------------------------------------------------------------------------- """ em_mvn!(; observed::SemObservedMissing, @@ -37,6 +36,12 @@ THE SOFTWARE. =# Estimates the covariance matrix and mean vector of the normal distribution via expectation maximization for `observed`. Overwrites the statistics stored in `observed`. + +Uses the EM algorithm for MVN-distributed data with missing values +adapted from the supplementary material to the book *Machine Learning: A Probabilistic Perspective*, +copyright (2010) Kevin Murphy and Matt Dunham: see +[*gaussMissingFitEm.m*](https://github.com/probml/pmtk3/blob/master/toolbox/BasicModels/gauss/sub/gaussMissingFitEm.m) and +[*emAlgo.m*](https://github.com/probml/pmtk3/blob/master/toolbox/Algorithms/optimization/emAlgo.m) scripts. """ function em_mvn!( observed::SemObservedMissing; diff --git a/src/optimizer/abstract.jl b/src/optimizer/abstract.jl index 3b1e98842..c1ad72592 100644 --- a/src/optimizer/abstract.jl +++ b/src/optimizer/abstract.jl @@ -41,7 +41,7 @@ function fit(optim::SemOptimizer, model::AbstractSem; start_val = nothing, kwarg end fit(model::AbstractSem; engine::Symbol = :Optim, start_val = nothing, kwargs...) = -fit(SemOptimizer(; engine, kwargs...), model; start_val, kwargs...) + fit(SemOptimizer(; engine, kwargs...), model; start_val, kwargs...) # fallback method fit(optim::SemOptimizer, model::AbstractSem, start_params; kwargs...) = @@ -56,8 +56,7 @@ prepare_start_params(start_val::Nothing, model::AbstractSem; kwargs...) = start_simple(model; kwargs...) # first argument is a function -prepare_start_params(start_val, model::AbstractSem; kwargs...) = - start_val(model; kwargs...) +prepare_start_params(start_val, model::AbstractSem; kwargs...) = start_val(model; kwargs...) function prepare_start_params(start_val::AbstractVector, model::AbstractSem; kwargs...) (length(start_val) == nparams(model)) || throw( diff --git a/src/optimizer/optim.jl b/src/optimizer/optim.jl index bd57942d8..2d782473a 100644 --- a/src/optimizer/optim.jl +++ b/src/optimizer/optim.jl @@ -67,7 +67,7 @@ SemOptimizer{:Optim}(args...; kwargs...) = SemOptimizerOptim(args...; kwargs...) SemOptimizerOptim(; algorithm = LBFGS(), - options = Optim.Options(;f_reltol = 1e-10, x_abstol = 1.5e-8), + options = Optim.Options(; f_reltol = 1e-10, x_abstol = 1.5e-8), kwargs..., ) = SemOptimizerOptim(algorithm, options) diff --git a/src/types.jl b/src/types.jl index 660c1c43f..0e279e5b2 100644 --- a/src/types.jl +++ b/src/types.jl @@ -192,7 +192,8 @@ Returns a SemEnsemble with fields For instructions on multigroup models, see the online documentation. """ -struct SemEnsemble{N, T <: Tuple, V <: AbstractVector, I, G <: Vector{Symbol}} <: AbstractSemCollection +struct SemEnsemble{N, T <: Tuple, V <: AbstractVector, I, G <: Vector{Symbol}} <: + AbstractSemCollection n::N sems::T weights::V diff --git a/test/examples/helper.jl b/test/examples/helper.jl index 4ff9bd507..acc3ccd08 100644 --- a/test/examples/helper.jl +++ b/test/examples/helper.jl @@ -90,12 +90,8 @@ function test_estimates( skip::Bool = false, ) actual = StructuralEquationModels.params(partable, col) - expected = StructuralEquationModels.lavaan_params( - partable_lav, - partable, - lav_col, - lav_group, - ) + expected = + StructuralEquationModels.lavaan_params(partable_lav, partable, lav_col, lav_group) @test !any(isnan, actual) @test !any(isnan, expected) diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl index f5ea0b5d7..e12d1c7c0 100644 --- a/test/examples/multigroup/build_models.jl +++ b/test/examples/multigroup/build_models.jl @@ -8,7 +8,8 @@ model_g1 = Sem(specification = specification_g1, data = dat_g1, implied = RAMSym model_g2 = Sem(specification = specification_g2, data = dat_g2, implied = RAM) -@test SEM.param_labels(model_g1.implied.ram_matrices) == SEM.param_labels(model_g2.implied.ram_matrices) +@test SEM.param_labels(model_g1.implied.ram_matrices) == + SEM.param_labels(model_g2.implied.ram_matrices) # test the different constructors model_ml_multigroup = SemEnsemble(model_g1, model_g2) diff --git a/test/examples/multigroup/multigroup.jl b/test/examples/multigroup/multigroup.jl index 239bf713c..43de554ce 100644 --- a/test/examples/multigroup/multigroup.jl +++ b/test/examples/multigroup/multigroup.jl @@ -9,11 +9,11 @@ dat = example_data("holzinger_swineford") dat_missing = example_data("holzinger_swineford_missing") solution_lav = example_data("holzinger_swineford_solution") -dat_g1 = dat[dat.school.=="Pasteur", :] -dat_g2 = dat[dat.school.=="Grant-White", :] +dat_g1 = dat[dat.school .== "Pasteur", :] +dat_g2 = dat[dat.school .== "Grant-White", :] -dat_miss_g1 = dat_missing[dat_missing.school.=="Pasteur", :] -dat_miss_g2 = dat_missing[dat_missing.school.=="Grant-White", :] +dat_miss_g1 = dat_missing[dat_missing.school .== "Pasteur", :] +dat_miss_g2 = dat_missing[dat_missing.school .== "Grant-White", :] dat.school = ifelse.(dat.school .== "Pasteur", :Pasteur, :Grant_White) dat_missing.school = ifelse.(dat_missing.school .== "Pasteur", :Pasteur, :Grant_White) diff --git a/test/examples/political_democracy/constructor.jl b/test/examples/political_democracy/constructor.jl index 7a8adc72e..4045141ce 100644 --- a/test/examples/political_democracy/constructor.jl +++ b/test/examples/political_democracy/constructor.jl @@ -1,6 +1,3 @@ -using Statistics: cov, mean -using Random, NLopt - ############################################################################################ ### models w.o. meanstructure ############################################################################################ diff --git a/test/examples/political_democracy/political_democracy.jl b/test/examples/political_democracy/political_democracy.jl index ad06e0fcd..cbdf7ea74 100644 --- a/test/examples/political_democracy/political_democracy.jl +++ b/test/examples/political_democracy/political_democracy.jl @@ -1,4 +1,6 @@ using StructuralEquationModels, Test, Suppressor, FiniteDiff +using Statistics: cov, mean +using Random, NLopt SEM = StructuralEquationModels diff --git a/test/examples/proximal/l0.jl b/test/examples/proximal/l0.jl index 374f8e58a..f74dfb2d1 100644 --- a/test/examples/proximal/l0.jl +++ b/test/examples/proximal/l0.jl @@ -1,5 +1,3 @@ -using StructuralEquationModels, Test, ProximalAlgorithms, ProximalOperators - # load data dat = example_data("political_democracy") @@ -62,6 +60,7 @@ fit_prox = fit(model_prox, engine = :Proximal, operator_g = prox_operator) @test fit_prox.optimization_result.result[:iterations] < 1000 @test solution(fit_prox)[31] == 0.0 @test abs( - StructuralEquationModels.minimum(fit_prox) - StructuralEquationModels.minimum(sem_fit), + StructuralEquationModels.minimum(fit_prox) - + StructuralEquationModels.minimum(sem_fit), ) < 1.0 end diff --git a/test/examples/proximal/lasso.jl b/test/examples/proximal/lasso.jl index beb5cf529..356ac6188 100644 --- a/test/examples/proximal/lasso.jl +++ b/test/examples/proximal/lasso.jl @@ -1,5 +1,3 @@ -using StructuralEquationModels, Test, ProximalAlgorithms, ProximalOperators - # load data dat = example_data("political_democracy") diff --git a/test/examples/proximal/proximal.jl b/test/examples/proximal/proximal.jl index 40e72a1ef..84a9162cb 100644 --- a/test/examples/proximal/proximal.jl +++ b/test/examples/proximal/proximal.jl @@ -1,3 +1,5 @@ +using StructuralEquationModels, Test, ProximalAlgorithms, ProximalOperators, Suppressor + @testset "Ridge" begin include("ridge.jl") end diff --git a/test/examples/proximal/ridge.jl b/test/examples/proximal/ridge.jl index fd7ae113d..61b7fa12a 100644 --- a/test/examples/proximal/ridge.jl +++ b/test/examples/proximal/ridge.jl @@ -1,5 +1,3 @@ -using StructuralEquationModels, Test, ProximalAlgorithms, ProximalOperators, Suppressor - # load data dat = example_data("political_democracy") diff --git a/test/examples/recover_parameters/recover_parameters_twofact.jl b/test/examples/recover_parameters/recover_parameters_twofact.jl index a7b4cec9a..ce7dc61ff 100644 --- a/test/examples/recover_parameters/recover_parameters_twofact.jl +++ b/test/examples/recover_parameters/recover_parameters_twofact.jl @@ -55,7 +55,7 @@ start = [ implied_ml = RAMSymbolic(; specification = ram_matrices, start_val = start) -implied_ml.Σ_function(implied_ml.Σ, true_val) +implied_ml.Σ_eval!(implied_ml.Σ, true_val) true_dist = MultivariateNormal(implied_ml.Σ) diff --git a/test/unit_tests/StatsAPI.jl b/test/unit_tests/StatsAPI.jl index 8648fc363..5907ee7b5 100644 --- a/test/unit_tests/StatsAPI.jl +++ b/test/unit_tests/StatsAPI.jl @@ -5,10 +5,7 @@ end partable = ParameterTable(graph, observed_vars = [:a, :b], latent_vars = Symbol[]) update_partable!(partable, :estimate, param_labels(partable), [3.1415]) data = randn(100, 2) -model = Sem( - specification = partable, - data = data -) +model = Sem(specification = partable, data = data) model_fit = fit(model) @testset "params" begin @@ -25,5 +22,5 @@ end end @testset "coeftable" begin - @test_throws "StructuralEquationModels does not support" coeftable(model) -end \ No newline at end of file + @test_throws "StructuralEquationModels does not support" coeftable(model) +end diff --git a/test/unit_tests/model.jl b/test/unit_tests/model.jl index 2bf5dedaf..d2c1b7a02 100644 --- a/test/unit_tests/model.jl +++ b/test/unit_tests/model.jl @@ -25,7 +25,6 @@ graph = @StenoGraph begin y8 ↔ y4 + y6 end - ram_matrices = RAMMatrices(ParameterTable(graph, observed_vars = obs_vars, latent_vars = lat_vars)) diff --git a/test/unit_tests/unit_tests.jl b/test/unit_tests/unit_tests.jl index 7189addd4..4d7dad7cf 100644 --- a/test/unit_tests/unit_tests.jl +++ b/test/unit_tests/unit_tests.jl @@ -7,7 +7,7 @@ available_tests = Dict( "data_input_formats" => "SemObserved", "specification" => "SemSpecification", "model" => "Sem model", - "StatsAPI" => "StatsAPI" + "StatsAPI" => "StatsAPI", ) # Determine which tests to run based on command-line arguments diff --git a/test/unit_tests/unit_tests_interactive.jl b/test/unit_tests/unit_tests_interactive.jl index cf082fa60..10a384403 100644 --- a/test/unit_tests/unit_tests_interactive.jl +++ b/test/unit_tests/unit_tests_interactive.jl @@ -7,4 +7,4 @@ try catch e @warn "Error initializing Test Env" exception=(e, catch_backtrace()) end -include("unit_tests.jl") \ No newline at end of file +include("unit_tests.jl")