Merge pull request #304 from StructuralEquationModels/fiml_meanstructure_argument

Fiml meanstructure argument

Author: Maximilian-Stefan-Ernst

src/implied/RAM/generic.jl

@@ -98,6 +98,8 @@ function RAM(;
 )
     ram_matrices = convert(RAMMatrices, specification)
 
+    check_meanstructure_specification(meanstructure, ram_matrices)
+
     # get dimensions of the model
     n_par = nparams(ram_matrices)
     n_obs = nobserved_vars(ram_matrices)
@@ -126,11 +128,6 @@ function RAM(;
     # μ
     if meanstructure
         MS = HasMeanStruct
-        !isnothing(ram_matrices.M) || throw(
-            ArgumentError(
-                "You set `meanstructure = true`, but your model specification contains no mean parameters.",
-            ),
-        )
         M_pre = materialize(ram_matrices.M, rand_params)
         ∇M = gradient_required ? sparse_gradient(ram_matrices.M) : nothing
         μ = zeros(n_obs)

src/implied/RAM/symbolic.jl

@@ -92,6 +92,8 @@ function RAMSymbolic(;
 )
     ram_matrices = convert(RAMMatrices, specification)
 
+    check_meanstructure_specification(meanstructure, ram_matrices)
+
     n_par = nparams(ram_matrices)
     par = (Symbolics.@variables θ[1:n_par])[1]
 

src/implied/abstract.jl

@@ -31,3 +31,17 @@ 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

src/loss/ML/FIML.jl

@@ -43,7 +43,14 @@ end
 ### Constructors
 ############################################################################################
 
-function SemFIML(; observed::SemObservedMissing, specification, kwargs...)
+function SemFIML(; observed::SemObservedMissing, implied, specification, kwargs...)
+
+    if implied.meanstruct isa NoMeanStruct
+        throw(ArgumentError(
+            "Full information maximum likelihood (FIML) can only be used with a meanstructure.
+            Did you forget to set `Sem(..., meanstructure = true)`?"))
+    end
+
     inverses =
         [zeros(nmeasured_vars(pat), nmeasured_vars(pat)) for pat in observed.patterns]
     choleskys = Array{Cholesky{Float64, Array{Float64, 2}}, 1}(undef, length(inverses))

src/loss/ML/ML.jl

@@ -39,6 +39,21 @@ end
 ############################################################################################
 
 function SemML(; observed::SemObserved, approximate_hessian::Bool = false, kwargs...)
+
+    if observed isa SemObservedMissing
+        throw(ArgumentError(
+            "Normal maximum likelihood estimation can't be used with `SemObservedMissing`.
+            Use full information maximum likelihood (FIML) estimation or remove missing 
+            values in your data.
+            A FIML model can be constructed with 
+            Sem(
+                ...,
+                observed = SemObservedMissing,
+                loss = SemFIML,
+                meanstructure = true
+            )"))
+    end
+
     obsmean = obs_mean(observed)
     obscov = obs_cov(observed)
     meandiff = isnothing(obsmean) ? nothing : copy(obsmean)

src/loss/WLS/WLS.jl

@@ -51,12 +51,33 @@ SemWLS{HE}(args...) where {HE <: HessianEval} =
 
 function SemWLS(;
     observed,
+    implied,
     wls_weight_matrix = nothing,
     wls_weight_matrix_mean = nothing,
     approximate_hessian = false,
     meanstructure = false,
     kwargs...,
 )
+
+    if observed isa SemObservedMissing
+        throw(ArgumentError(
+            "WLS estimation can't be used with `SemObservedMissing`.
+            Use full information maximum likelihood (FIML) estimation or remove missing 
+            values in your data.
+            A FIML model can be constructed with 
+            Sem(
+                ...,
+                observed = SemObservedMissing,
+                loss = SemFIML,
+                meanstructure = true
+            )"))
+    end
+
+    if !(implied isa RAMSymbolic) 
+        throw(ArgumentError(
+            "WLS estimation is only available with the implied type RAMSymbolic at the moment."))
+    end
+
     nobs_vars = nobserved_vars(observed)
     tril_ind = filter(x -> (x[1] >= x[2]), CartesianIndices(obs_cov(observed)))
     s = obs_cov(observed)[tril_ind]

src/observed/data.jl

@@ -38,10 +38,24 @@ function SemObservedData(;
     observed_var_prefix::Union{Symbol, AbstractString} = :obs,
     kwargs...,
 )
+
     data, obs_vars, _ =
         prepare_data(data, observed_vars, specification; observed_var_prefix)
     obs_mean, obs_cov = mean_and_cov(data, 1)
 
+    if any(ismissing.(data))
+        throw(ArgumentError(
+            "Your dataset contains missing values.
+            Remove missing values or use full information maximum likelihood (FIML) estimation.
+            A FIML model can be constructed with 
+            Sem(
+                ...,
+                observed = SemObservedMissing,
+                loss = SemFIML,
+                meanstructure = true
+            )"))
+    end
+
     return SemObservedData(data, obs_vars, obs_cov, vec(obs_mean), size(data, 1))
 end
 

test/unit_tests/model.jl

@@ -46,14 +46,13 @@ function test_params_api(semobj, spec::SemSpecification)
     @test @inferred(param_labels(semobj)) == param_labels(spec)
 end
 
-@testset "Sem(implied=$impliedtype, loss=$losstype)" for impliedtype in (RAM, RAMSymbolic),
-    losstype in (SemML, SemWLS)
+@testset "Sem(implied=$impliedtype, loss=SemML)" for impliedtype in (RAM, RAMSymbolic)
 
     model = Sem(
         specification = ram_matrices,
         observed = obs,
         implied = impliedtype,
-        loss = losstype,
+        loss = SemML,
     )
 
     @test model isa Sem
@@ -68,7 +67,33 @@ end
 
     @test @inferred(loss(model)) isa SemLoss
     semloss = loss(model).functions[1]
-    @test semloss isa losstype
+    @test semloss isa SemML
 
     @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