Convergence, bootstrap, p-values

ext/SEMNLOptExt/NLopt.jl

@@ -134,6 +134,10 @@ end
 SEM.algorithm_name(res::NLoptResult) = res.problem.algorithm
 SEM.n_iterations(res::NLoptResult) = res.problem.numevals
 SEM.convergence(res::NLoptResult) = res.result[3]
+function SEM.converged(res::NLoptResult)
+    flag = res.result[3]
+    return flag ∈ [:SUCCESS, :STOPVAL_REACHED, :FTOL_REACHED, :XTOL_REACHED]
+end
 
 # construct NLopt.jl problem
 function NLopt_problem(algorithm, options, npar)

ext/SEMProximalOptExt/ProximalAlgorithms.jl

@@ -90,10 +90,10 @@ SEM.algorithm_name(res::ProximalResult) = SEM.algorithm_name(res.optimizer.algor
 SEM.algorithm_name(
     ::ProximalAlgorithms.IterativeAlgorithm{I, H, S, D, K},
 ) where {I, H, S, D, K} = nameof(I)
-
 SEM.convergence(
     ::ProximalResult,
 ) = "No standard convergence criteria for proximal \n algorithms available."
+SEM.converged(::ProximalResult) = missing
 SEM.n_iterations(res::ProximalResult) = res.n_iterations
 
 ############################################################################################

src/StructuralEquationModels.jl

@@ -90,6 +90,8 @@ include("frontend/fit/fitmeasures/fit_measures.jl")
 # standard errors
 include("frontend/fit/standard_errors/hessian.jl")
 include("frontend/fit/standard_errors/bootstrap.jl")
+include("frontend/fit/standard_errors/z_test.jl")
+include("frontend/fit/standard_errors/confidence_intervals.jl")
 
 export AbstractSem,
     AbstractSemSingle,
@@ -129,6 +131,7 @@ export AbstractSem,
     optimizer_engine_doc,
     optimizer_engines,
     n_iterations,
+    converged,
     convergence,
     SemObserved,
     SemObservedData,
@@ -191,7 +194,12 @@ export AbstractSem,
     CFI,
     EmMVNModel,
     se_hessian,
+    bootstrap,
     se_bootstrap,
+    normal_CI,
+    normal_CI!,
+    z_test,
+    z_test!,
     example_data,
     replace_observed,
     update_observed,

src/frontend/fit/SemFit.jl

@@ -15,7 +15,8 @@ Fitted structural equation model.
 
 - `algorithm_name(::SemFit)` -> optimization algorithm
 - `n_iterations(::SemFit)` -> number of iterations
-- `convergence(::SemFit)` -> convergence properties
+- `convergence(::SemFit)` -> convergence flags
+- `converged(::SemFit)` -> convergence success
 """
 mutable struct SemFit{Mi, So, St, Mo, O}
     minimum::Mi
@@ -71,3 +72,4 @@ optimizer_engine(sem_fit::SemFit) = optimizer_engine(optimization_result(sem_fit
 algorithm_name(sem_fit::SemFit) = algorithm_name(optimization_result(sem_fit))
 n_iterations(sem_fit::SemFit) = n_iterations(optimization_result(sem_fit))
 convergence(sem_fit::SemFit) = convergence(optimization_result(sem_fit))
+converged(sem_fit::SemFit) = converged(optimization_result(sem_fit))

src/frontend/fit/standard_errors/bootstrap.jl

@@ -1,39 +1,166 @@
+"""
+    bootstrap(
+        fitted::SemFit,
+        specification::SemSpecification;
+        statistic = solution,
+        n_boot = 3000,
+        data = nothing,
+        engine = :Optim,
+        parallel = false,
+        fit_kwargs = Dict(),
+        replace_kwargs = Dict())
+
+Return bootstrap samples for `statistic`.
+
+# 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`
+- `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`
+
+# Example
+```julia
+# 1000 boostrap samples of the minimum, fitted with :Optim
+bootstrap(
+    fitted;
+    statistic = StructuralEquationModels.minimum,
+    n_boot = 1000,
+    engine = :Optim,
+)
+```
+"""
+function bootstrap(
+    fitted::SemFit,
+    specification::SemSpecification;
+    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 = []
+    conv = []
+    # 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...)
+            sample = statistic(new_fit)
+            c = converged(new_fit)
+            push!(out, sample)
+            push!(conv, 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 _ 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)
+            lock(lk) do
+                push!(out, sample)
+                push!(conv, c)
+            end
+            put!(model_pool, thread_model)
+        end
+    end
+    return Dict(
+        :samples => out,
+        :n_boot => n_boot,
+        :n_converged => isempty(conv) ? 0 : sum(conv),
+        :converged => conv,
+    )
+end
+
 """
     se_bootstrap(
-        sem_fit::SemFit;
+        fitted::SemFit,
+        specification::SemSpecification;
         n_boot = 3000,
         data = nothing,
-        specification = nothing,
         parallel = false,
-        kwargs...)
+        fit_kwargs = Dict(),
+        replace_kwargs = Dict())
 
 Return bootstrap standard errors.
 
 # 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`
-- `specification`: a `ParameterTable` or `RAMMatrices` object passed down to `replace_observed`.
-   Necessary for FIML models.
+- `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.
-- `kwargs...`: passed down to `replace_observed`
+- `fit_kwargs` : a `Dict` controlling model fitting for each bootstrap sample, 
+  passed to `sem_fit`
+- `replace_kwargs`: a `Dict` passed to `replace_observed`
+
+# Example
+```julia
+# 1000 boostrap samples, fitted with :NLopt
+using NLopt
+
+se_bootstrap(
+    fitted;
+    n_boot = 1000,
+    engine = :NLopt,
+)
+```
 """
 function se_bootstrap(
-    fitted::SemFit{Mi, So, St, Mo, O};
+    fitted::SemFit,
+    specification::SemSpecification;
     n_boot = 3000,
     data = nothing,
-    specification = nothing,
+    engine = :Optim,
     parallel = false,
-    kwargs...,
-) where {Mi, So, St, Mo, O}
+    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_failed          = Ref(0)
+    n_conv            = Ref(0)
     # fit to bootstrap samples
     if !parallel
         for _ in 1:n_boot
@@ -42,14 +169,15 @@ function se_bootstrap(
                 fitted.model;
                 data = sample_data,
                 specification = specification,
-                kwargs...,
+                replace_kwargs...,
             )
-            try
-                sol = solution(fit(new_model; start_val = start))
+            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
-            catch
-                n_failed[] += 1
             end
         end
     else
@@ -63,37 +191,37 @@ function se_bootstrap(
         lk = ReentrantLock()
         Threads.@threads for _ in 1:n_boot
             thread_model = take!(model_pool)
-            try
-                sample_data = bootstrap_sample(data)
-                new_model = replace_observed(
-                    thread_model;
-                    data = sample_data,
-                    specification = specification,
-                    kwargs...,
-                )
-                sol = solution(fit(new_model; start_val = start))
+            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
-            catch
-                lock(lk) do
-                    n_failed[] += 1
-                end
-            finally
-                put!(model_pool, thread_model)
             end
+            put!(model_pool, thread_model)
         end
     end
     # compute parameters
-    n_conv = n_boot - n_failed[]
+    n_conv = n_conv[]
     sd = sqrt.(total_squared_sum / n_conv - (total_sum / n_conv) .^ 2)
-    if !iszero(n_failed[])
-        @warn "During bootstrap sampling, "*string(n_failed[])*" models did not converge"
-    end
+    @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)

src/frontend/fit/standard_errors/confidence_intervals.jl

@@ -0,0 +1,43 @@
+_doc_normal_CI = """
+    (1) normal_CI(fitted, se; α = 0.05, name_lower = :ci_lower, name_upper = :ci_upper)
+
+    (2) normal_CI!(partable, fitted, se; α = 0.05, name_lower = :ci_lower, name_upper = :ci_upper)
+
+Return normal-theory confidence intervals for all model parameters.
+`normal_CI!` additionally writes the result into `partable`.
+
+# Arguments
+- `fitted`: a fitted SEM.
+- `se`: standard errors for each parameter, e.g. from [`se_hessian`](@ref) or
+  [`se_bootstrap`](@ref).
+- `partable`: a [`ParameterTable`](@ref) to write confidence intervals to.
+- `α`: significance level. Defaults to `0.05` (95% intervals).
+- `name_lower`: column name for the lower bound in `partable`. Defaults to `:ci_lower`.
+- `name_upper`: column name for the upper bound in `partable`. Defaults to `:ci_upper`.
+
+# Returns
+- a `Dict` with keys `name_lower` and `name_upper`, each mapping to a vector of bounds 
+  over all parameters.
+"""
+
+@doc "$(_doc_normal_CI)"
+function normal_CI(fitted, se; α = 0.05, name_lower = :ci_lower, name_upper = :ci_upper)
+    qnt = quantile(Normal(0, 1), 1-α/2);
+    sol = solution(fitted)
+    return Dict(name_lower => sol - qnt*se, name_upper => sol + qnt*se)
+end
+
+@doc "$(_doc_normal_CI)"
+function normal_CI!(
+    partable,
+    fitted,
+    se;
+    α = 0.05,
+    name_lower = :ci_lower,
+    name_upper = :ci_upper,
+)
+    cis = normal_CI(fitted, se; α, name_lower, name_upper)
+    update_partable!(partable, name_lower, fitted, cis[name_lower])
+    update_partable!(partable, name_upper, fitted, cis[name_upper])
+    return cis
+end