Documentation and bugfixes for 0.5.0

Project.toml

@@ -27,22 +27,23 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 
 [compat]
-julia = "1.9, 1.10, 1.11"
+julia = "1.10, 1.11, 1.12"
 StenoGraphs = "0.5"
 DataFrames = "1"
 Distributions = "0.25"
 FiniteDiff = "2"
-InteractiveUtils = "1.11.0"
 LineSearches = "7"
-NLSolversBase = "7"
+NLSolversBase = "8"
 NLopt = "0.6, 1"
-Optim = "1"
+Optim = "2"
 PrettyTables = "3"
 ProximalAlgorithms = "0.7"
 StatsBase = "0.33, 0.34"
 Symbolics = "4, 5, 6, 7"
 SymbolicUtils = "1.4 - 1.5, 1.7, 2, 3, 4"
 StatsAPI = "1"
+DelimitedFiles = "1"
+Statistics = "1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

docs/make.jl

@@ -52,11 +52,6 @@ makedocs(
             "Starting values" => "performance/starting_values.md",
             "Parametric Types" => "performance/parametric.md",
         ],
-        "Internals and design" => [
-            "Internals and design" => "internals/internals.md",
-            "files" => "internals/files.md",
-            "types" => "internals/types.md",
-        ],
     ],
     format = Documenter.HTML(
         prettyurls = get(ENV, "CI", nothing) == "true",

docs/src/developer/extending.md

@@ -1,10 +1,12 @@
 # Extending the package
 
-As discussed in the section on [Model Construction](@ref), every Structural Equation Model (`Sem`) consists of three (four with the optimizer) parts:
+As discussed in [Our Concept of a Structural Equation Model](@ref), a `Sem` is built from one or more loss
+terms, and each SEM loss function bundles an *observed* and an *implied* part:
 
 ![SEM concept typed](../assets/concept_typed.svg)
 
-On the following pages, we will explain how you can define your own custom parts and "plug them in". There are certain things you **have to do** to define custom parts and some things you **can do** to have a more pleasent experience. In general, these requirements fall into the categories
+On the following pages, we will explain how you can define your own custom parts (a loss function, an observed
+type, an implied type, or an optimizer) and "plug them in". There are certain things you **have to do** to define custom parts and some things you **can do** to have a more pleasent experience. In general, these requirements fall into the categories
 - minimal (to use your custom part and fit a `Sem` with it)
 - use the outer constructor to build a model in a more convenient way
 - use additional functionality like standard errors, fit measures, etc.

docs/src/developer/implied.md

@@ -43,8 +43,7 @@ We implement an `ImpliedEmpty` type in our package that does nothing but serving
 ### Types
 ############################################################################################
 """
-Empty placeholder for models that don't need an implied part.
-(For example, models that only regularize parameters.)
+Empty placeholder for loss functions that don't need an implied part.
 
 # Constructor
 
@@ -55,7 +54,7 @@ Empty placeholder for models that don't need an implied part.
 
 # Examples
 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.
+SEM term (`SemLoss`) per group and an additional `SemRidge` regularization term.
 
 # Extended help
 
@@ -73,13 +72,13 @@ end
 ### Constructors
 ############################################################################################
 
-function ImpliedEmpty(
-    spec::SemSpecification;
-    hessianeval::HessianApprox = ExactHessian(),
+function ImpliedEmpty(;
+    specification,
+    meanstruct = NoMeanStruct(),
+    hessianeval = ExactHessian(),
     kwargs...,
 )
-    ram_matrices = convert(RAMMatrices, spec)
-    return ImpliedEmpty(hessianeval, MeanStruct(ram_matrices), ram_matrices)
+    return ImpliedEmpty(hessianeval, meanstruct, convert(RAMMatrices, specification))
 end
 
 ############################################################################################

docs/src/developer/loss.md

@@ -71,7 +71,7 @@ partable = ParameterTable(
 )
 
 parameter_indices = getindex.([param_indices(partable)], [:a, :b, :c])
-myridge = Ridge(0.01, parameter_indices)
+myridge = MyRidge(0.01, parameter_indices)
 
 model = SemFiniteDiff(
     specification = partable,
@@ -82,16 +82,16 @@ model = SemFiniteDiff(
 model_fit = fit(model)
 ```
 
-This is one way of specifying the model - we now have **one model** with **multiple loss functions**. Because we did not provide a gradient for `Ridge`, we have to specify a `SemFiniteDiff` model that computes numerical gradients with finite difference approximation.
+This is one way of specifying the model - we now have **one model** with **multiple loss functions**. Because we did not provide a gradient for `MyRidge`, we have to specify a `SemFiniteDiff` model that computes numerical gradients with finite difference approximation.
 
-Note that the last argument to the `objective!` method is the whole model. Therefore, we can access everything that is stored inside our model everytime we compute the objective value for our loss function. Since ridge regularization is a very easy case, we do not need to do this. But maximum likelihood estimation for example depends on both the observed and the model implied covariance matrix. See [Second example - maximum likelihood](@ref) for information on how to do that.
+Ridge regularization only depends on the parameters, so the `evaluate!` method above does not need anything else. Other loss functions, however, depend on the observed data and on what the model implies about it. Loss functions that compare the implied and the observed structure are subtypes of [`SemLoss`](@ref) and store their own `observed` and `implied` parts, which can be accessed inside `evaluate!` via `observed(loss)` and `implied(loss)`. See [Second example - maximum likelihood](@ref) for information on how to do that.
 
 ### Improve performance
 
 By far the biggest improvements in performance will result from specifying analytical gradients. We can do this for our example:
 
 ```@example loss
-function evaluate!(objective, gradient, hessian::Nothing, ridge::Ridge, model::AbstractSem, par)
+function evaluate!(objective, gradient, hessian::Nothing, ridge::MyRidge, par)
     # compute gradient
     if !isnothing(gradient)
         fill!(gradient, 0)
@@ -136,37 +136,43 @@ Additionally, you may provide analytic hessians by writing a respective method f
 
 ## Convenient
 
-To be able to build the loss term, it needs a constructor.
-Every `SemLoss` subtype should provide a constructor with 3 positional arguments:
-  * `observed::SemObserved`: the observed part of the model
-  * `implied::SemImplied`: the implied part of the model
-  * `refloss::Union{MyLoss, Nothing} = nothing`: optional loss term of the same type
-    to use as a reference for any loss-specific configuration.
+In the minimal example above we built `myridge` ourselves and passed the ready-made instance to
+the model via `loss = (SemML, myridge)`. Alternatively, you can let the outer [`Sem`](@ref)
+constructor build the loss term for you: pass the loss *type* instead of an instance and provide
+a keyword constructor.
 
-Any additional loss configuration details should be passed as optional keyword arguments.
-If both `refloss` and the keyword arguments are provided, the keyword arguments take
-precedence. This constructor is used internally by the functions like [`replace_observed`](@ref)
-to rebuild the loss term with new observed data while preserving the implied state.
+```julia
+MyRidge(; α_ridge, which_ridge, kwargs...) = MyRidge(α_ridge, which_ridge)
+```
+
+Any keyword arguments passed to `Sem(...)` are forwarded to this constructor (along with some that
+the model supplies automatically, such as `nparams`), so the loss can be configured directly from
+the model call:
 
 ```julia
-function MyLoss(
-    observed::SemObserved, implied::SemImplied, refloss::Union{MyLoss, Nothing} = nothing;
-    kwarg1 = ..., kwarg2 = ..., kwargs...
+model = SemFiniteDiff(
+    specification = partable,
+    data = example_data("political_democracy"),
+    loss = (SemML, MyRidge),
+    α_ridge = 0.01,
+    which_ridge = parameter_indices,
 )
-    ...
-    return MyLoss(...) # internal MyLoss constructor
-end
 ```
 
+Note that, being a plain `AbstractLoss`, `MyRidge` neither stores nor receives an `observed` or
+`implied` part — it depends only on the parameters. 
+SEM-specific loss functions are constructed differently; see
+[Second example - maximum likelihood](@ref).
+
 ## Additional functionality
 
 ### 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:
 
 ```julia
-hyperparameter(ridge::Ridge) = ridge.α
-regularization_indices(ridge::Ridge) = ridge.I
+hyperparameter(ridge::MyRidge) = ridge.α
+regularization_indices(ridge::MyRidge) = ridge.I
 ```
 
 # Second example - maximum likelihood
@@ -178,7 +184,9 @@ To keep it simple, we only cover models without a meanstructure. The maximum lik
 F_{ML} = \log \det \Sigma_i + \mathrm{tr}\left(\Sigma_{i}^{-1} \Sigma_o \right)
 ```
 
-where ``\Sigma_i`` is the model implied covariance matrix and ``\Sigma_o`` is the observed covariance matrix. We can query the model implied covariance matrix from the `implied` par of our model, and the observed covariance matrix from the `observed` path of our model.
+where ``\Sigma_i`` is the model implied covariance matrix and ``\Sigma_o`` is the observed covariance matrix. We can query the model implied covariance matrix from the `implied` part of our loss term, and the observed covariance matrix from the `observed` part of our loss term.
+
+Since this loss function compares the implied and the observed structure, it is a subtype of [`SemLoss`](@ref) rather than a plain `AbstractLoss`. Every `SemLoss` stores its own `observed` and `implied` parts, which can be accessed inside `evaluate!` via `observed(loss)` and `implied(loss)`.
 
 To get information on what we can access from a certain `implied` or `observed` type, we can check it`s documentation an the pages [API - model parts](@ref) or via the help mode of the REPL:
 
@@ -190,20 +198,33 @@ help?> RAM
 help?> SemObservedData
 ```
 
-We see that the model implied covariance matrix can be assessed as `Σ(implied)` and the observed covariance matrix as `obs_cov(observed)`.
+We see that the model implied covariance matrix can be assessed as `implied(loss).Σ` and the observed covariance matrix as `obs_cov(observed(loss))`.
+
+Unlike a plain `AbstractLoss`, a `SemLoss` subtype stores its `observed` and `implied` parts (in its first two fields), and the [`Sem`](@ref) constructor builds it for you. To support this, every `SemLoss` subtype should provide a constructor with three positional arguments:
+  * `observed::SemObserved`: the observed part of the loss term
+  * `implied::SemImplied`: the implied part of the loss term
+  * `refloss::Union{MaximumLikelihood, Nothing} = nothing`: an optional existing loss term of the
+    same type, used as a reference for any loss-specific configuration.
 
-With this information, we write can implement maximum likelihood optimization as
+Any additional configuration is passed as optional keyword arguments; if both `refloss` and keyword arguments are given, the keyword arguments take precedence. This constructor is also used by [`replace_observed`](@ref) to rebuild the loss term with new observed data while sharing the implied state. With this, we can implement maximum likelihood optimization as
 
 ```@example loss
-struct MaximumLikelihood <: SemLossFunction end
+struct MaximumLikelihood{O <: SemObserved, I <: SemImplied} <: SemLoss{O, I}
+    observed::O
+    implied::I
+end
+
+# constructor used by the `Sem` constructor to build the loss term
+MaximumLikelihood(observed::SemObserved, implied::SemImplied, refloss = nothing; kwargs...) =
+    MaximumLikelihood{typeof(observed), typeof(implied)}(observed, implied)
 
 using LinearAlgebra
-import StructuralEquationModels: obs_cov, evaluate!
+import StructuralEquationModels: evaluate!
 
-function evaluate!(objective::Number, gradient::Nothing, hessian::Nothing, semml::MaximumLikelihood, model::AbstractSem, par)
+function evaluate!(objective::Number, gradient::Nothing, hessian::Nothing, semml::MaximumLikelihood, par)
     # access the model implied and observed covariance matrices
-    Σᵢ = implied(model).Σ
-    Σₒ = obs_cov(observed(model))
+    Σᵢ = implied(semml).Σ
+    Σₒ = obs_cov(observed(semml))
     # compute the objective
     if isposdef(Symmetric(Σᵢ)) # is the model implied covariance matrix positive definite?
         return logdet(Σᵢ) + tr(inv(Σᵢ)*Σₒ)
@@ -221,7 +242,7 @@ Let's specify and fit a model:
 model_ml = SemFiniteDiff(
     specification = partable,
     data = example_data("political_democracy"),
-    loss = MaximumLikelihood()
+    loss = MaximumLikelihood
 )
 
 model_fit = fit(model_ml)

docs/src/performance/parametric.md

@@ -25,7 +25,7 @@ giving each field a type and adding them as parameters to our type declaration.
 Recall our example from [Custom loss functions](@ref):
 
 ```julia
-struct Ridge <: SemLossFunction
+struct Ridge <: AbstractLoss
     α
     I
 end
@@ -34,7 +34,7 @@ end
 We could also declare it as a parametric type:
 
 ```julia
-struct ParametricRidge{X, Y} <: SemLossFunction
+struct ParametricRidge{X, Y} <: AbstractLoss
     α::X
     I::Y
 end

docs/src/tutorials/collection/collection.md

@@ -3,7 +3,7 @@
 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
+In this case, [`FiniteDiffWrapper`](@ref StructuralEquationModels.FiniteDiffWrapper) 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.
 
 ```julia
@@ -26,16 +26,22 @@ It is also possible to use finite difference for the entire `Sem` model:
 model_findiff2 = FiniteDiffWrapper(model)
 ```
 
-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`:
+The weighting scheme of the SEM loss terms is specified using the `default_sem_weights` argument of the `Sem` constructor.
+The available schemes are:
+- `:nsamples_corrected` (the default): like `:nsamples`, but applies a loss-type-specific correction
+  to the sample counts (e.g. ``N_{term} - 1`` for maximum likelihood and weighted least squares). For FIML the correction is zero, so it coincides with `:nsamples`,
+- `:nsamples`: weights each SEM term by ``N_{term}/N_{total}``, i.e. by the (uncorrected) number of
+  observations in its data,
+- `:uniform`: weights each of the ``k`` SEM terms by ``1/k``,
+- `:one`: leaves all SEM terms unweighted.
+
+The weights for the loss terms (both SEM and regularization) can also be explicitly specified using the pair syntax `loss => weight`:
 
 ```julia
 model_weighted = Sem(loss_1 => 0.5, loss_2 => 1.0)
 ```
 
-`Sem` support assigning unique identifier to each loss term, which is essential for complex multi-term model.
+`Sem` supports assigning a unique identifier to each loss term, which is useful for complex multi-term models.
 The syntax is `id => loss`, or `id => loss => weight`:
 
 ```julia
@@ -47,6 +53,8 @@ model2_weighted = Sem(:main => loss_1 => 0.5, :alt => loss_2 => 1.0)
 
 ```@docs
 Sem
-LossTerm
-FiniteDiffWrapper
+SemFiniteDiff
+AbstractSem
+StructuralEquationModels.LossTerm
+StructuralEquationModels.FiniteDiffWrapper
 ```

docs/src/tutorials/collection/multigroup.md

@@ -52,7 +52,7 @@ graph = @StenoGraph begin
 end
 ```
 
-You can pass multiple arguments to `fix()` and `label()` for each group. Parameters with the same label (within and across groups) are constrained to be equal. To fix a parameter in one group, but estimate it freely in the other, you may write `fix(NaN, 4.3)`.
+You can pass multiple arguments to `fixed()` and `label()` for each group. Parameters with the same label (within and across groups) are constrained to be equal. To fix a parameter in one group, but estimate it freely in the other, you may write `fixed(NaN, 4.3)`.
 
 You can then use the resulting graph to specify an `EnsembleParameterTable`
 
@@ -72,8 +72,7 @@ The parameter table can be used to create a multigroup `Sem` model:
 model_ml_multigroup = Sem(
     specification = partable,
     data = dat,
-    semterm_column = :school,
-    groups = groups)
+    semterm_column = :school)
 ```
 
 !!! note "A different way to specify"