diff --git a/ext/SEMNLOptExt/NLopt.jl b/ext/SEMNLOptExt/NLopt.jl
index 909dbbfc1..90004b907 100644
--- a/ext/SEMNLOptExt/NLopt.jl
+++ b/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)
diff --git a/ext/SEMProximalOptExt/ProximalAlgorithms.jl b/ext/SEMProximalOptExt/ProximalAlgorithms.jl
index 1d7f83632..3ec324530 100644
--- a/ext/SEMProximalOptExt/ProximalAlgorithms.jl
+++ b/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
 
 ############################################################################################
diff --git a/src/StructuralEquationModels.jl b/src/StructuralEquationModels.jl
index f144c98bc..19dd6f43a 100644
--- a/src/StructuralEquationModels.jl
+++ b/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,
diff --git a/src/frontend/fit/SemFit.jl b/src/frontend/fit/SemFit.jl
index 9c2d114e7..1d2e82a60 100644
--- a/src/frontend/fit/SemFit.jl
+++ b/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))
diff --git a/src/frontend/fit/standard_errors/bootstrap.jl b/src/frontend/fit/standard_errors/bootstrap.jl
index eb7aefa7b..b2a909e9b 100644
--- a/src/frontend/fit/standard_errors/bootstrap.jl
+++ b/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)
diff --git a/src/frontend/fit/standard_errors/confidence_intervals.jl b/src/frontend/fit/standard_errors/confidence_intervals.jl
new file mode 100644
index 000000000..20bf58a73
--- /dev/null
+++ b/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
diff --git a/src/frontend/fit/standard_errors/z_test.jl b/src/frontend/fit/standard_errors/z_test.jl
new file mode 100644
index 000000000..27bebf147
--- /dev/null
+++ b/src/frontend/fit/standard_errors/z_test.jl
@@ -0,0 +1,36 @@
+_doc_z_test = """
+    (1) z_test(fitted, se)
+
+    (2) z_test!(partable, fitted, se, name = :p_value)
+
+Return two-sided p-values from a z-test for each model parameter.
+
+Tests the null hypothesis that each parameter is zero using the test statistic
+`z = estimate / se`, which is compared against a standard normal distribution.
+`z_test!` 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 p-values to.
+- `name`: column name for the p-values in `partable`. Defaults to `:p_value`.
+
+# Returns
+- a vector of p-values.
+"""
+
+@doc "$(_doc_z_test)"
+function z_test(fitted, se)
+    dev = solution(fitted) ./ se
+    dist = Normal(0, 1)
+    p = 2*ccdf.(dist, abs.(dev))
+    return p
+end
+
+@doc "$(_doc_z_test)"
+function z_test!(partable, fitted, se, name = :p_value)
+    p = z_test(fitted, se)
+    update_partable!(partable, name, fitted, p)
+    return p
+end
diff --git a/src/implied/RAM/generic.jl b/src/implied/RAM/generic.jl
index 37591c232..3b8596874 100644
--- a/src/implied/RAM/generic.jl
+++ b/src/implied/RAM/generic.jl
@@ -196,9 +196,13 @@ end
 ############################################################################################
 
 function update_observed(implied::RAM, observed::SemObserved; kwargs...)
-    if nobserved_vars(observed) == size(implied.Σ, 1)
+    if nobserved_vars(observed) == nobserved_vars(implied)
         return implied
     else
-        return RAM(; observed = observed, kwargs...)
+        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 436f339b7..0d3ba5e11 100644
--- a/src/implied/RAM/symbolic.jl
+++ b/src/implied/RAM/symbolic.jl
@@ -210,10 +210,17 @@ end
 ############################################################################################
 
 function update_observed(implied::RAMSymbolic, observed::SemObserved; kwargs...)
-    if nobserved_vars(observed) == size(implied.Σ, 1)
+    if nobserved_vars(observed) == nobserved_vars(implied)
         return implied
     else
-        return RAMSymbolic(; observed = observed, kwargs...)
+        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
 
diff --git a/src/loss/ML/ML.jl b/src/loss/ML/ML.jl
index 6461ba087..67d4fe524 100644
--- a/src/loss/ML/ML.jl
+++ b/src/loss/ML/ML.jl
@@ -237,6 +237,9 @@ function update_observed(lossfun::SemML, observed::SemObserved; kwargs...)
     if size(lossfun.Σ⁻¹) == size(obs_cov(observed))
         return lossfun
     else
-        return SemML(; observed = observed, kwargs...)
+        return SemML(;
+        observed = observed,
+        approximate_hessian = HessianEval(lossfun) == ApproxHessian,
+        kwargs...)
     end
 end
diff --git a/src/loss/WLS/WLS.jl b/src/loss/WLS/WLS.jl
index b2aed17c0..2b10d7b47 100644
--- a/src/loss/WLS/WLS.jl
+++ b/src/loss/WLS/WLS.jl
@@ -173,5 +173,7 @@ end
 ### Recommended methods
 ############################################################################################
 
-update_observed(lossfun::SemWLS, observed::SemObserved; kwargs...) =
-    SemWLS(; observed = observed, kwargs...)
+update_observed(lossfun::SemWLS, observed::SemObserved; kwargs...) = SemWLS(;
+    observed = observed,
+    meanstructure = MeanStruct(kwargs[:implied]) == HasMeanStruct,
+    kwargs...)
diff --git a/src/optimizer/optim.jl b/src/optimizer/optim.jl
index 83ebbe5e1..704131938 100644
--- a/src/optimizer/optim.jl
+++ b/src/optimizer/optim.jl
@@ -77,7 +77,17 @@ end
 
 algorithm_name(res::SemOptimResult) = Optim.summary(res.result)
 n_iterations(res::SemOptimResult) = Optim.iterations(res.result)
-convergence(res::SemOptimResult) = Optim.converged(res.result)
+function convergence(res::SemOptimResult)
+    flags = res.result.stopped_by
+    active_flags = Symbol[]
+    for key in keys(flags)
+        if flags[key]
+            push!(active_flags, key)
+        end
+    end
+    return active_flags
+end
+converged(res::SemOptimResult) = Optim.converged(res.result)
 
 function fit(
     optim::SemOptimizerOptim,
diff --git a/test/examples/helper.jl b/test/examples/helper.jl
index acc3ccd08..31d65679a 100644
--- a/test/examples/helper.jl
+++ b/test/examples/helper.jl
@@ -135,3 +135,40 @@ function test_estimates(
         @test actual ≈ expected rtol = rtol atol = atol norm = Base.Fix2(norm, Inf)
     end
 end
+
+function test_bootstrap(
+        model_fit,
+        spec;
+        compare_hessian = true,
+        rtol_hessian = 0.2,
+        compare_bs = true,
+        rtol_bs = 0.1,
+        n_boot = 500)
+    se_bs = se_bootstrap(model_fit, spec; n_boot = n_boot)
+    # hessian and bootstrap se are close
+    if compare_hessian
+        se_he = se_hessian(model_fit)
+        @test isapprox(se_bs, se_he, rtol = rtol_hessian)
+    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)
+        se_bs_2 = sqrt.(var(bs_samples, corrected = false, dims = 2))
+        @test isapprox(se_bs_2, se_bs, rtol = rtol_bs)
+    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)
+    return se_bs, bs_samples
+end
+
+function smoketest_CI_z(model_fit, partable)
+    se_he = se_hessian(model_fit)
+    normal_CI!(partable, model_fit, se_he)
+    z_test!(partable, model_fit, se_he)
+end
diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl
index e10a8a058..5cbb345c7 100644
--- a/test/examples/multigroup/build_models.jl
+++ b/test/examples/multigroup/build_models.jl
@@ -12,7 +12,7 @@ model_g2 = Sem(specification = specification_g2, data = dat_g2, implied = RAM)
       SEM.param_labels(model_g2.implied.ram_matrices)
 
 # test the different constructors
-model_ml_multigroup = SemEnsemble(model_g1, model_g2)
+model_ml_multigroup = SemEnsemble(model_g1, model_g2; groups = [:Pasteur, :Grant_White])
 model_ml_multigroup2 = SemEnsemble(
     specification = partable,
     data = dat,
@@ -83,6 +83,8 @@ end
         lav_col = :se,
         lav_groups = Dict(:Pasteur => 1, :Grant_White => 2),
     )
+    test_bootstrap(solution_ml, partable; rtol_hessian = 0.3, rtol_bs = 0.2, n_boot = 1_000)
+    smoketest_CI_z(solution_ml, partable)
 
     solution_ml = fit(model_ml_multigroup2)
     test_fitmeasures(
@@ -250,7 +252,7 @@ model_ls_g2 = Sem(
     loss = SemWLS,
 )
 
-model_ls_multigroup = SemEnsemble(model_ls_g1, model_ls_g2; optimizer = semoptimizer)
+model_ls_multigroup = SemEnsemble(model_ls_g1, model_ls_g2; groups = [:Pasteur, :Grant_White], optimizer = semoptimizer)
 
 @testset "ls_gradients_multigroup" begin
     test_gradient(model_ls_multigroup, start_test; atol = 1e-9)
@@ -291,6 +293,8 @@ end
         lav_col = :se,
         lav_groups = Dict(:Pasteur => 1, :Grant_White => 2),
     )
+    test_bootstrap(solution_ls, partable; compare_bs = false, rtol_hessian = 0.3)
+    smoketest_CI_z(solution_ls, partable)
 end
 
 ############################################################################################
@@ -409,6 +413,7 @@ if !isnothing(specification_miss_g1)
             lav_groups = Dict(:Pasteur => 1, :Grant_White => 2),
         )
 
+
         solution = fit(semoptimizer, model_ml_multigroup2)
         test_fitmeasures(
             fit_measures(solution),
@@ -422,6 +427,9 @@ if !isnothing(specification_miss_g1)
             fitmeasure_names = Dict(:CFI => "cfi"),
         )
 
+        test_bootstrap(solution, partable_miss; compare_bs = false, rtol_hessian = 0.3)
+        smoketest_CI_z(solution, partable_miss)
+
         update_se_hessian!(partable_miss, solution)
         test_estimates(
             partable_miss,
diff --git a/test/examples/multigroup/multigroup.jl b/test/examples/multigroup/multigroup.jl
index 1f89714d8..dd654731d 100644
--- a/test/examples/multigroup/multigroup.jl
+++ b/test/examples/multigroup/multigroup.jl
@@ -1,5 +1,9 @@
 using StructuralEquationModels, Test, FiniteDiff, Suppressor
 using LinearAlgebra: diagind, LowerTriangular
+using Statistics: var
+using Random
+
+Random.seed!(948723)
 
 const SEM = StructuralEquationModels
 
diff --git a/test/examples/political_democracy/constructor.jl b/test/examples/political_democracy/constructor.jl
index 45de75d13..1c1c42e54 100644
--- a/test/examples/political_democracy/constructor.jl
+++ b/test/examples/political_democracy/constructor.jl
@@ -130,6 +130,9 @@ end
         col = :se,
         lav_col = :se,
     )
+
+    test_bootstrap(solution_ml, partable)
+    smoketest_CI_z(solution_ml, partable)
 end
 
 @testset "fitmeasures/se_ls" begin
@@ -157,6 +160,9 @@ end
         col = :se,
         lav_col = :se,
     )
+
+    test_bootstrap(solution_ls, partable; compare_bs = false)
+    smoketest_CI_z(solution_ls, partable)
 end
 
 ############################################################################################
@@ -337,6 +343,9 @@ end
         col = :se,
         lav_col = :se,
     )
+
+    test_bootstrap(solution_ml, partable_mean)
+    smoketest_CI_z(solution_ml, partable_mean)
 end
 
 @testset "fitmeasures/se_ls_mean" begin
@@ -363,6 +372,10 @@ end
         col = :se,
         lav_col = :se,
     )
+
+    test_bootstrap(solution_ls, partable_mean, compare_bs = false)
+    # smoketest_bootstrap(solution_ls, partable_mean)
+    smoketest_CI_z(solution_ls, partable_mean)
 end
 
 ############################################################################################
@@ -494,4 +507,8 @@ end
         col = :se,
         lav_col = :se,
     )
+
+    # test_bootstrap(solution_ml, partable_mean) # too much compute
+    smoketest_bootstrap(solution_ml, partable_mean)
+    smoketest_CI_z(solution_ml, partable_mean)
 end
diff --git a/test/examples/political_democracy/political_democracy.jl b/test/examples/political_democracy/political_democracy.jl
index 9c8cc2a7b..8c8d6c36e 100644
--- a/test/examples/political_democracy/political_democracy.jl
+++ b/test/examples/political_democracy/political_democracy.jl
@@ -1,7 +1,9 @@
 using StructuralEquationModels, Test, Suppressor, FiniteDiff
-using Statistics: cov, mean
+using Statistics: cov, mean, var
 using Random, NLopt
 
+Random.seed!(464577)
+
 SEM = StructuralEquationModels
 
 include(