diff --git a/README.md b/README.md
index 79c11da21..9754a8c20 100644
--- a/README.md
+++ b/README.md
@@ -15,7 +15,7 @@ Models you can fit include
 - Multigroup SEM
 - Sums of arbitrary loss functions (everything the optimizer can handle).
 
-# What are the merrits?
+# What are the merits?
 
 We provide fast objective functions, gradients, and for some cases hessians as well as approximations thereof.
 As a user, you can easily define custom loss functions.
diff --git a/docs/src/developer/loss.md b/docs/src/developer/loss.md
index 57a7b485d..931c2d0e5 100644
--- a/docs/src/developer/loss.md
+++ b/docs/src/developer/loss.md
@@ -79,7 +79,7 @@ model = SemFiniteDiff(
     loss = (SemML, myridge)
 )
 
-model_fit = sem_fit(model)
+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.
@@ -117,7 +117,7 @@ model_new = Sem(
     loss = (SemML, myridge)
 )
 
-model_fit = sem_fit(model_new)
+model_fit = fit(model_new)
 ```
 
 The results are the same, but we can verify that the computational costs are way lower (for this, the julia package `BenchmarkTools` has to be installed):
@@ -125,9 +125,9 @@ The results are the same, but we can verify that the computational costs are way
 ```julia
 using BenchmarkTools
 
-@benchmark sem_fit(model)
+@benchmark fit(model)
 
-@benchmark sem_fit(model_new)
+@benchmark fit(model_new)
 ```
 
 The exact results of those benchmarks are of course highly depended an your system (processor, RAM, etc.), but you should see that the median computation time with analytical gradients drops to about 5% of the computation without analytical gradients.
@@ -241,7 +241,7 @@ model_ml = SemFiniteDiff(
     loss = MaximumLikelihood()
 )
 
-model_fit = sem_fit(model_ml)
+model_fit = fit(model_ml)
 ```
 
 If you want to differentiate your own loss functions via automatic differentiation, check out the [AutoDiffSEM](https://github.com/StructuralEquationModels/AutoDiffSEM) package.
diff --git a/docs/src/developer/optimizer.md b/docs/src/developer/optimizer.md
index 82ec594d8..a651ec636 100644
--- a/docs/src/developer/optimizer.md
+++ b/docs/src/developer/optimizer.md
@@ -34,7 +34,7 @@ algorithm(optimizer::SemOptimizerName) = optimizer.algorithm
 options(optimizer::SemOptimizerName) = optimizer.options
 ```
 
-Note that your optimizer is a subtype of `SemOptimizer{:Name}`, where you can choose a `:Name` that can later be used as a keyword argument to `sem_fit(engine = :Name)`.
+Note that your optimizer is a subtype of `SemOptimizer{:Name}`, where you can choose a `:Name` that can later be used as a keyword argument to `fit(engine = :Name)`.
 Similarly, `SemOptimizer{:Name}(args...; kwargs...) = SemOptimizerName(args...; kwargs...)` should be defined as well as a constructor that uses only keyword arguments:
 
 ´´´julia
@@ -46,10 +46,10 @@ SemOptimizerName(;
 ´´´
 A method for `update_observed` and additional methods might be usefull, but are not necessary.
 
-Now comes the substantive part: We need to provide a method for `sem_fit`:
+Now comes the substantive part: We need to provide a method for `fit`:
 
 ```julia
-function sem_fit(
+function fit(
     optim::SemOptimizerName,
     model::AbstractSem,
     start_params::AbstractVector;
diff --git a/docs/src/internals/files.md b/docs/src/internals/files.md
index 0872c2b02..90ceceaaf 100644
--- a/docs/src/internals/files.md
+++ b/docs/src/internals/files.md
@@ -11,7 +11,7 @@ Source code is in the `"src"` folder:
 - `"types.jl"` defines all abstract types and the basic type hierarchy
 - `"objective_gradient_hessian.jl"` contains methods for computing objective, gradient and hessian values for different model types as well as generic fallback methods
 - The four folders `"observed"`, `"implied"`, `"loss"` and `"diff"` contain implementations of specific subtypes (for example, the `"loss"` folder contains a file `"ML.jl"` that implements the `SemML` loss function).
-- `"optimizer"` contains connections to different optimization backends (aka methods for `sem_fit`)
+- `"optimizer"` contains connections to different optimization backends (aka methods for `fit`)
     - `"optim.jl"`: connection to the `Optim.jl` package
 - `"frontend"` contains user-facing functions
     - `"specification"` contains functionality for model specification
diff --git a/docs/src/performance/mixed_differentiation.md b/docs/src/performance/mixed_differentiation.md
index 2ac937077..b7ae333b5 100644
--- a/docs/src/performance/mixed_differentiation.md
+++ b/docs/src/performance/mixed_differentiation.md
@@ -19,7 +19,7 @@ model_ridge = SemFiniteDiff(
 
 model_ml_ridge = SemEnsemble(model_ml, model_ridge)
 
-model_ml_ridge_fit = sem_fit(model_ml_ridge)
+model_ml_ridge_fit = fit(model_ml_ridge)
 ```
 
 The results of both methods will be the same, but we can verify that the computation costs differ (the package `BenchmarkTools` has to be installed for this):
@@ -27,7 +27,7 @@ The results of both methods will be the same, but we can verify that the computa
 ```julia
 using BenchmarkTools
 
-@benchmark sem_fit(model)
+@benchmark fit(model)
 
-@benchmark sem_fit(model_ml_ridge)
+@benchmark fit(model_ml_ridge)
 ```
\ No newline at end of file
diff --git a/docs/src/performance/mkl.md b/docs/src/performance/mkl.md
index 0d5467658..4361ab445 100644
--- a/docs/src/performance/mkl.md
+++ b/docs/src/performance/mkl.md
@@ -27,9 +27,9 @@ To check the performance implications for fitting a SEM, you can use the [`Bench
 ```julia
 using BenchmarkTools
 
-@benchmark sem_fit($your_model)
+@benchmark fit($your_model)
 
 using MKL
 
-@benchmark sem_fit($your_model)
+@benchmark fit($your_model)
 ```
\ No newline at end of file
diff --git a/docs/src/performance/simulation.md b/docs/src/performance/simulation.md
index 881da6222..0cb2ea25d 100644
--- a/docs/src/performance/simulation.md
+++ b/docs/src/performance/simulation.md
@@ -100,7 +100,7 @@ models = [model1, model2]
 fits = Vector{SemFit}(undef, 2)
 
 Threads.@threads for i in 1:2
-    fits[i] = sem_fit(models[i])
+    fits[i] = fit(models[i])
 end
 ```
 
diff --git a/docs/src/performance/starting_values.md b/docs/src/performance/starting_values.md
index ba7b4f41d..2df8d94d4 100644
--- a/docs/src/performance/starting_values.md
+++ b/docs/src/performance/starting_values.md
@@ -1,9 +1,9 @@
 # Starting values
 
-The `sem_fit` function has a keyword argument that takes either a vector of starting values or a function that takes a model as input to compute starting values. Current options are `start_fabin3` for fabin 3 starting values [^Hägglund82] or `start_simple` for simple starting values. Additional keyword arguments to `sem_fit` are passed to the starting value function. For example,
+The `fit` function has a keyword argument that takes either a vector of starting values or a function that takes a model as input to compute starting values. Current options are `start_fabin3` for fabin 3 starting values [^Hägglund82] or `start_simple` for simple starting values. Additional keyword arguments to `fit` are passed to the starting value function. For example,
 
 ```julia
-    sem_fit(
+    fit(
         model; 
         start_val = start_simple,
         start_covariances_latent = 0.5
diff --git a/docs/src/tutorials/collection/multigroup.md b/docs/src/tutorials/collection/multigroup.md
index 23c13b950..1007f4563 100644
--- a/docs/src/tutorials/collection/multigroup.md
+++ b/docs/src/tutorials/collection/multigroup.md
@@ -81,7 +81,7 @@ model_ml_multigroup = SemEnsemble(
 We now fit the model and inspect the parameter estimates:
 
 ```@example mg; ansicolor = true
-fit = sem_fit(model_ml_multigroup)
+fit = fit(model_ml_multigroup)
 update_estimate!(partable, fit)
 details(partable)
 ```
diff --git a/docs/src/tutorials/constraints/constraints.md b/docs/src/tutorials/constraints/constraints.md
index 66f3def19..338803cb3 100644
--- a/docs/src/tutorials/constraints/constraints.md
+++ b/docs/src/tutorials/constraints/constraints.md
@@ -48,7 +48,7 @@ model = Sem(
     data = data
 )
 
-model_fit = sem_fit(model)
+model_fit = fit(model)
 
 update_estimate!(partable, model_fit)
 
@@ -153,7 +153,7 @@ model_constrained = Sem(
     data = data
 )
 
-model_fit_constrained = sem_fit(constrained_optimizer, model_constrained)
+model_fit_constrained = fit(constrained_optimizer, model_constrained)
 ```
 
 As you can see, the optimizer converged (`:XTOL_REACHED`) and investigating the solution yields
diff --git a/docs/src/tutorials/construction/build_by_parts.md b/docs/src/tutorials/construction/build_by_parts.md
index 606a6576e..45d2a2ea1 100644
--- a/docs/src/tutorials/construction/build_by_parts.md
+++ b/docs/src/tutorials/construction/build_by_parts.md
@@ -65,5 +65,5 @@ optimizer = SemOptimizerOptim()
 
 model_ml = Sem(observed, implied_ram, loss_ml)
 
-sem_fit(optimizer, model_ml)
+fit(optimizer, model_ml)
 ```
\ No newline at end of file
diff --git a/docs/src/tutorials/construction/outer_constructor.md b/docs/src/tutorials/construction/outer_constructor.md
index 6a3cd2cef..a1c0b8ad3 100644
--- a/docs/src/tutorials/construction/outer_constructor.md
+++ b/docs/src/tutorials/construction/outer_constructor.md
@@ -131,4 +131,4 @@ model = SemFiniteDiff(
 )
 ```
 
-constructs a model that will use finite difference approximation if you estimate the parameters via `sem_fit(model)`.
\ No newline at end of file
+constructs a model that will use finite difference approximation if you estimate the parameters via `fit(model)`.
\ No newline at end of file
diff --git a/docs/src/tutorials/first_model.md b/docs/src/tutorials/first_model.md
index 5b7284649..e8048966c 100644
--- a/docs/src/tutorials/first_model.md
+++ b/docs/src/tutorials/first_model.md
@@ -110,7 +110,7 @@ model = Sem(
 We can now fit the model via
 
 ```@example high_level; ansicolor = true
-model_fit = sem_fit(model)
+model_fit = fit(model)
 ```
 
 and compute fit measures as
diff --git a/docs/src/tutorials/fitting/fitting.md b/docs/src/tutorials/fitting/fitting.md
index a3e4b9b91..fff06abaa 100644
--- a/docs/src/tutorials/fitting/fitting.md
+++ b/docs/src/tutorials/fitting/fitting.md
@@ -3,7 +3,7 @@
 As we saw in [A first model](@ref), after you have build a model, you can fit it via
 
 ```julia
-model_fit = sem_fit(model)
+model_fit = fit(model)
 
 # output
 
@@ -45,24 +45,24 @@ Structural Equation Model
 
 ## Choosing an optimizer
 
-To choose a different optimizer, you can call `sem_fit` with the keyword argument `engine = ...`, and pass additional keyword arguments:
+To choose a different optimizer, you can call `fit` with the keyword argument `engine = ...`, and pass additional keyword arguments:
 
 ```julia
 using Optim
 
-model_fit = sem_fit(model; engine = :Optim, algorithm = BFGS())
+model_fit = fit(model; engine = :Optim, algorithm = BFGS())
 ```
 
 Available options for engine are `:Optim`, `:NLopt` and `:Proximal`, where `:NLopt` and `:Proximal` are only available if the `NLopt.jl` and `ProximalAlgorithms.jl` packages are loaded respectively.
 
 The available keyword arguments are listed in the sections [Using Optim.jl](@ref), [Using NLopt.jl](@ref) and [Regularization](@ref).
 
-Alternative, you can also explicitely define a `SemOptimizer` and pass it as the first argument to `sem_fit`:
+Alternative, you can also explicitely define a `SemOptimizer` and pass it as the first argument to `fit`:
 
 ```julia
 my_optimizer = SemOptimizerOptim(algorithm = BFGS())
 
-sem_fit(my_optimizer, model)
+fit(my_optimizer, model)
 ```
 
 You may also optionally specify [Starting values](@ref).
@@ -70,5 +70,5 @@ You may also optionally specify [Starting values](@ref).
 # API - model fitting
 
 ```@docs
-sem_fit
+fit
 ```
\ No newline at end of file
diff --git a/docs/src/tutorials/inspection/inspection.md b/docs/src/tutorials/inspection/inspection.md
index 628df62d0..abd416c1c 100644
--- a/docs/src/tutorials/inspection/inspection.md
+++ b/docs/src/tutorials/inspection/inspection.md
@@ -42,13 +42,13 @@ model = Sem(
     data = data
 )
 
-model_fit = sem_fit(model)
+model_fit = fit(model)
 ```
 
 After you fitted a model,
 
 ```julia
-model_fit = sem_fit(model)
+model_fit = fit(model)
 ```
 
 you end up with an object of type [`SemFit`](@ref).
@@ -126,7 +126,7 @@ fit_measures
 AIC
 BIC
 χ²
-df
+dof
 minus2ll
 nobserved_vars
 nsamples
diff --git a/docs/src/tutorials/meanstructure.md b/docs/src/tutorials/meanstructure.md
index 60578224a..b2da5029a 100644
--- a/docs/src/tutorials/meanstructure.md
+++ b/docs/src/tutorials/meanstructure.md
@@ -96,7 +96,7 @@ model = Sem(
     meanstructure = true
 )
 
-sem_fit(model)
+fit(model)
 ```
 
 If we build the model by parts, we have to pass the `meanstructure = true` argument to every part that requires it (when in doubt, simply consult the documentation for the respective part).
@@ -112,5 +112,5 @@ ml = SemML(observed = observed, meanstructure = true)
 
 model = Sem(observed, implied_ram, SemLoss(ml))
 
-sem_fit(model)
+fit(model)
 ```
\ No newline at end of file
diff --git a/docs/src/tutorials/regularization/regularization.md b/docs/src/tutorials/regularization/regularization.md
index 34d7c509a..3d82fcfba 100644
--- a/docs/src/tutorials/regularization/regularization.md
+++ b/docs/src/tutorials/regularization/regularization.md
@@ -120,13 +120,13 @@ Let's fit the regularized model
 
 ```@example reg
 
-fit_lasso = sem_fit(optimizer_lasso, model_lasso)
+fit_lasso = fit(optimizer_lasso, model_lasso)
 ```
 
 and compare the solution to unregularizted estimates:
 
 ```@example reg
-fit = sem_fit(model)
+fit = fit(model)
 
 update_estimate!(partable, fit)
 
@@ -135,10 +135,10 @@ update_partable!(partable, :estimate_lasso, param_labels(fit_lasso), solution(fi
 details(partable)
 ```
 
-Instead of explicitely defining a `SemOptimizerProximal` object, you can also pass `engine = :Proximal` and additional keyword arguments to `sem_fit`:
+Instead of explicitely defining a `SemOptimizerProximal` object, you can also pass `engine = :Proximal` and additional keyword arguments to `fit`:
 
 ```@example reg
-fit = sem_fit(model; engine = :Proximal, operator_g = NormL1(λ))
+fit = fit(model; engine = :Proximal, operator_g = NormL1(λ))
 ```
 
 ## Second example - mixed l1 and l0 regularization
@@ -162,7 +162,7 @@ model_mixed = Sem(
     data = data,    
 )
 
-fit_mixed = sem_fit(model_mixed; engine = :Proximal, operator_g = prox_operator)
+fit_mixed = fit(model_mixed; engine = :Proximal, operator_g = prox_operator)
 ```
 
 Let's again compare the different results:
diff --git a/ext/SEMNLOptExt/NLopt.jl b/ext/SEMNLOptExt/NLopt.jl
index a614c501b..c5e0ad6cb 100644
--- a/ext/SEMNLOptExt/NLopt.jl
+++ b/ext/SEMNLOptExt/NLopt.jl
@@ -71,8 +71,8 @@ function SemFit_NLopt(optimization_result, model::AbstractSem, start_val, opt)
     )
 end
 
-# sem_fit method
-function SEM.sem_fit(
+# fit method
+function SEM.fit(
     optim::SemOptimizerNLopt,
     model::AbstractSem,
     start_params::AbstractVector;
diff --git a/ext/SEMProximalOptExt/ProximalAlgorithms.jl b/ext/SEMProximalOptExt/ProximalAlgorithms.jl
index 2f1775e85..0d4748e3a 100644
--- a/ext/SEMProximalOptExt/ProximalAlgorithms.jl
+++ b/ext/SEMProximalOptExt/ProximalAlgorithms.jl
@@ -40,7 +40,7 @@ mutable struct ProximalResult
     result::Any
 end
 
-function SEM.sem_fit(
+function SEM.fit(
     optim::SemOptimizerProximal,
     model::AbstractSem,
     start_params::AbstractVector;
diff --git a/src/StructuralEquationModels.jl b/src/StructuralEquationModels.jl
index 13591d72c..f6068dc50 100644
--- a/src/StructuralEquationModels.jl
+++ b/src/StructuralEquationModels.jl
@@ -16,7 +16,7 @@ using LinearAlgebra,
     DelimitedFiles,
     DataFrames
 
-import StatsAPI: params, coef, coefnames
+import StatsAPI: params, coef, coefnames, dof, fit, nobs, coeftable
 
 export StenoGraphs, @StenoGraph, meld
 
@@ -78,7 +78,7 @@ include("additional_functions/simulation.jl")
 include("frontend/fit/fitmeasures/AIC.jl")
 include("frontend/fit/fitmeasures/BIC.jl")
 include("frontend/fit/fitmeasures/chi2.jl")
-include("frontend/fit/fitmeasures/df.jl")
+include("frontend/fit/fitmeasures/dof.jl")
 include("frontend/fit/fitmeasures/minus2ll.jl")
 include("frontend/fit/fitmeasures/p.jl")
 include("frontend/fit/fitmeasures/RMSEA.jl")
@@ -95,6 +95,7 @@ export AbstractSem,
     AbstractSemCollection,
     coef,
     coefnames,
+    coeftable,
     Sem,
     SemFiniteDiff,
     SemEnsemble,
@@ -135,8 +136,9 @@ export AbstractSem,
     obs_cov,
     obs_mean,
     nsamples,
+    nobs,
     samples,
-    sem_fit,
+    fit,
     SemFit,
     minimum,
     solution,
@@ -179,7 +181,7 @@ export AbstractSem,
     AIC,
     BIC,
     χ²,
-    df,
+    dof,
     fit_measures,
     minus2ll,
     p_value,
diff --git a/src/frontend/StatsAPI.jl b/src/frontend/StatsAPI.jl
index 156ef4941..edd677e34 100644
--- a/src/frontend/StatsAPI.jl
+++ b/src/frontend/StatsAPI.jl
@@ -57,13 +57,22 @@ params(partable::ParameterTable, col::Symbol = :estimate) =
 """
     coef(x::ParameterTable)
 
-For a `SEM`, this function is equivalent to `params(x)`.
-It returns the parameters for the given model.
+For a `ParameterTable`, this function is synonymous to [`params`](@ref).
 """
 coef(x::ParameterTable) = params(x)
 
 """
     coefnames(x::ParameterTable)
-To maintain compatibility with the `lavaan` package, this function is a synonym for `param_labels(x)`.
+
+Synonymous to [`param_labels`](@ref param_labels).
 """
 coefnames(x::ParameterTable) = param_labels(x)
+
+"""
+    nobs(model::AbstractSem) -> Int
+
+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
diff --git a/src/frontend/fit/fitmeasures/RMSEA.jl b/src/frontend/fit/fitmeasures/RMSEA.jl
index b91e81d3e..b9fff648e 100644
--- a/src/frontend/fit/fitmeasures/RMSEA.jl
+++ b/src/frontend/fit/fitmeasures/RMSEA.jl
@@ -6,13 +6,13 @@ Return the RMSEA.
 function RMSEA end
 
 RMSEA(sem_fit::SemFit{Mi, So, St, Mo, O} where {Mi, So, St, Mo <: AbstractSemSingle, O}) =
-    RMSEA(df(sem_fit), χ²(sem_fit), nsamples(sem_fit))
+    RMSEA(dof(sem_fit), χ²(sem_fit), nsamples(sem_fit))
 
 RMSEA(sem_fit::SemFit{Mi, So, St, Mo, O} where {Mi, So, St, Mo <: SemEnsemble, O}) =
-    sqrt(length(sem_fit.model.sems)) * RMSEA(df(sem_fit), χ²(sem_fit), nsamples(sem_fit))
+    sqrt(length(sem_fit.model.sems)) * RMSEA(dof(sem_fit), χ²(sem_fit), nsamples(sem_fit))
 
-function RMSEA(df, chi2, nsamples)
-    rmsea = (chi2 - df) / (nsamples * df)
+function RMSEA(dof, chi2, nsamples)
+    rmsea = (chi2 - dof) / (nsamples * dof)
     rmsea > 0 ? nothing : rmsea = 0
     return sqrt(rmsea)
 end
diff --git a/src/frontend/fit/fitmeasures/df.jl b/src/frontend/fit/fitmeasures/dof.jl
similarity index 62%
rename from src/frontend/fit/fitmeasures/df.jl
rename to src/frontend/fit/fitmeasures/dof.jl
index 4d9025601..3df49d89d 100644
--- a/src/frontend/fit/fitmeasures/df.jl
+++ b/src/frontend/fit/fitmeasures/dof.jl
@@ -1,14 +1,14 @@
 """
-    df(sem_fit::SemFit)
-    df(model::AbstractSem)
+    dof(sem_fit::SemFit)
+    dof(model::AbstractSem)
 
 Return the degrees of freedom.
 """
-function df end
+function dof end
 
-df(sem_fit::SemFit) = df(sem_fit.model)
+dof(sem_fit::SemFit) = dof(sem_fit.model)
 
-df(model::AbstractSem) = n_dp(model) - nparams(model)
+dof(model::AbstractSem) = n_dp(model) - nparams(model)
 
 function n_dp(model::AbstractSemSingle)
     nvars = nobserved_vars(model)
diff --git a/src/frontend/fit/fitmeasures/fit_measures.jl b/src/frontend/fit/fitmeasures/fit_measures.jl
index 40e3caae0..2fc4dfba0 100644
--- a/src/frontend/fit/fitmeasures/fit_measures.jl
+++ b/src/frontend/fit/fitmeasures/fit_measures.jl
@@ -1,5 +1,5 @@
 fit_measures(sem_fit) =
-    fit_measures(sem_fit, nparams, df, AIC, BIC, RMSEA, χ², p_value, minus2ll)
+    fit_measures(sem_fit, nparams, dof, AIC, BIC, RMSEA, χ², p_value, minus2ll)
 
 function fit_measures(sem_fit, args...)
     measures = Dict{Symbol, Union{Float64, Missing}}()
diff --git a/src/frontend/fit/fitmeasures/p.jl b/src/frontend/fit/fitmeasures/p.jl
index 3d4275f95..8c69d5ec2 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(df(sem_fit)), χ²(sem_fit))
+p_value(sem_fit::SemFit) = 1 - cdf(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 e8d840d0c..4589dc020 100644
--- a/src/frontend/fit/standard_errors/bootstrap.jl
+++ b/src/frontend/fit/standard_errors/bootstrap.jl
@@ -1,5 +1,5 @@
 """
-    se_bootstrap(semfit::SemFit; n_boot = 3000, data = nothing, kwargs...)
+    se_bootstrap(sem_fit::SemFit; n_boot = 3000, data = nothing, kwargs...)
 
 Return boorstrap standard errors.
 Only works for single models.
@@ -52,7 +52,7 @@ function se_bootstrap(
         new_solution .= 0.0
 
         try
-            new_solution = solution(sem_fit(new_model; start_val = start))
+            new_solution = solution(fit(new_model; start_val = start))
         catch
             n_failed += 1
         end
diff --git a/src/frontend/specification/Sem.jl b/src/frontend/specification/Sem.jl
index e7c545ca7..7ba8f7fb7 100644
--- a/src/frontend/specification/Sem.jl
+++ b/src/frontend/specification/Sem.jl
@@ -45,6 +45,11 @@ Returns the [*observed*](@ref SemObserved) part of a model.
 """
 observed(model::AbstractSemSingle) = model.observed
 
+"""
+    nsamples(model::AbstractSem) -> Int
+
+Returns the number of samples from the [*observed*](@ref SemObserved) part of a model.
+"""
 nsamples(model::AbstractSemSingle) = nsamples(observed(model))
 
 """
diff --git a/src/optimizer/abstract.jl b/src/optimizer/abstract.jl
index 68bcc04ad..2487b7c52 100644
--- a/src/optimizer/abstract.jl
+++ b/src/optimizer/abstract.jl
@@ -1,5 +1,5 @@
 """
-    sem_fit([optim::SemOptimizer], model::AbstractSem;
+    fit([optim::SemOptimizer], model::AbstractSem;
             [engine::Symbol], start_val = start_val, kwargs...)
 
 Return the fitted `model`.
@@ -20,25 +20,25 @@ the online documentation on [Starting values](@ref).
 
 # Examples
 ```julia
-sem_fit(
+fit(
     my_model;
     start_val = start_simple,
     start_covariances_latent = 0.5)
 ```
 """
-function sem_fit(optim::SemOptimizer, model::AbstractSem; start_val = nothing, kwargs...)
+function fit(optim::SemOptimizer, model::AbstractSem; start_val = nothing, kwargs...)
     start_params = prepare_start_params(start_val, model; kwargs...)
     @assert start_params isa AbstractVector
     @assert length(start_params) == nparams(model)
 
-    sem_fit(optim, model, start_params; kwargs...)
+    fit(optim, model, start_params; kwargs...)
 end
 
-sem_fit(model::AbstractSem; engine::Symbol = :Optim, start_val = nothing, kwargs...) =
-    sem_fit(SemOptimizer(; engine, kwargs...), model; start_val, kwargs...)
+fit(model::AbstractSem; engine::Symbol = :Optim, start_val = nothing, kwargs...) =
+fit(SemOptimizer(; engine, kwargs...), model; start_val, kwargs...)
 
 # fallback method
-sem_fit(optim::SemOptimizer, model::AbstractSem, start_params; kwargs...) =
+fit(optim::SemOptimizer, model::AbstractSem, start_params; kwargs...) =
     error("Optimizer $(optim) support not implemented.")
 
 # FABIN3 is the default method for single models
diff --git a/src/optimizer/optim.jl b/src/optimizer/optim.jl
index cec37a77a..8f5404bc2 100644
--- a/src/optimizer/optim.jl
+++ b/src/optimizer/optim.jl
@@ -102,7 +102,7 @@ optimizer(res::Optim.MultivariateOptimizationResults) = Optim.summary(res)
 n_iterations(res::Optim.MultivariateOptimizationResults) = Optim.iterations(res)
 convergence(res::Optim.MultivariateOptimizationResults) = Optim.converged(res)
 
-function sem_fit(
+function fit(
     optim::SemOptimizerOptim,
     model::AbstractSem,
     start_params::AbstractVector;
diff --git a/test/examples/helper.jl b/test/examples/helper.jl
index 01122a841..4ff9bd507 100644
--- a/test/examples/helper.jl
+++ b/test/examples/helper.jl
@@ -51,7 +51,7 @@ end
 fitmeasure_names_ml = Dict(
     :AIC => "aic",
     :BIC => "bic",
-    :df => "df",
+    :dof => "df",
     :χ² => "chisq",
     :p_value => "pvalue",
     :nparams => "npar",
@@ -59,7 +59,7 @@ fitmeasure_names_ml = Dict(
 )
 
 fitmeasure_names_ls = Dict(
-    :df => "df",
+    :dof => "df",
     :χ² => "chisq",
     :p_value => "pvalue",
     :nparams => "npar",
diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl
index 6cf10549b..f6a7a230d 100644
--- a/test/examples/multigroup/build_models.jl
+++ b/test/examples/multigroup/build_models.jl
@@ -28,7 +28,7 @@ end
 
 # fit
 @testset "ml_solution_multigroup" begin
-    solution = sem_fit(semoptimizer, model_ml_multigroup)
+    solution = fit(semoptimizer, model_ml_multigroup)
     update_estimate!(partable, solution)
     test_estimates(
         partable,
@@ -36,7 +36,7 @@ end
         atol = 1e-4,
         lav_groups = Dict(:Pasteur => 1, :Grant_White => 2),
     )
-    solution = sem_fit(semoptimizer, model_ml_multigroup2)
+    solution = fit(semoptimizer, model_ml_multigroup2)
     update_estimate!(partable, solution)
     test_estimates(
         partable,
@@ -47,7 +47,7 @@ end
 end
 
 @testset "fitmeasures/se_ml" begin
-    solution_ml = sem_fit(model_ml_multigroup)
+    solution_ml = fit(model_ml_multigroup)
     test_fitmeasures(
         fit_measures(solution_ml),
         solution_lav[:fitmeasures_ml];
@@ -64,7 +64,7 @@ end
         lav_groups = Dict(:Pasteur => 1, :Grant_White => 2),
     )
 
-    solution_ml = sem_fit(model_ml_multigroup2)
+    solution_ml = fit(model_ml_multigroup2)
     test_fitmeasures(
         fit_measures(solution_ml),
         solution_lav[:fitmeasures_ml];
@@ -113,7 +113,7 @@ grad_fd = FiniteDiff.finite_difference_gradient(
 
 # fit
 @testset "ml_solution_multigroup | sorted" begin
-    solution = sem_fit(model_ml_multigroup)
+    solution = fit(model_ml_multigroup)
     update_estimate!(partable_s, solution)
     test_estimates(
         partable_s,
@@ -124,7 +124,7 @@ grad_fd = FiniteDiff.finite_difference_gradient(
 end
 
 @testset "fitmeasures/se_ml | sorted" begin
-    solution_ml = sem_fit(model_ml_multigroup)
+    solution_ml = fit(model_ml_multigroup)
     test_fitmeasures(
         fit_measures(solution_ml),
         solution_lav[:fitmeasures_ml];
@@ -191,7 +191,7 @@ end
 
 # fit
 @testset "solution_user_defined_loss" begin
-    solution = sem_fit(model_ml_multigroup)
+    solution = fit(model_ml_multigroup)
     update_estimate!(partable, solution)
     test_estimates(
         partable,
@@ -226,7 +226,7 @@ model_ls_multigroup = SemEnsemble(model_ls_g1, model_ls_g2; optimizer = semoptim
 end
 
 @testset "ls_solution_multigroup" begin
-    solution = sem_fit(model_ls_multigroup)
+    solution = fit(model_ls_multigroup)
     update_estimate!(partable, solution)
     test_estimates(
         partable,
@@ -237,7 +237,7 @@ end
 end
 
 @testset "fitmeasures/se_ls" begin
-    solution_ls = sem_fit(model_ls_multigroup)
+    solution_ls = fit(model_ls_multigroup)
     test_fitmeasures(
         fit_measures(solution_ls),
         solution_lav[:fitmeasures_ls];
@@ -321,7 +321,7 @@ if !isnothing(specification_miss_g1)
     end
 
     @testset "fiml_solution_multigroup" begin
-        solution = sem_fit(semoptimizer, model_ml_multigroup)
+        solution = fit(semoptimizer, model_ml_multigroup)
         update_estimate!(partable_miss, solution)
         test_estimates(
             partable_miss,
@@ -329,7 +329,7 @@ if !isnothing(specification_miss_g1)
             atol = 1e-4,
             lav_groups = Dict(:Pasteur => 1, :Grant_White => 2),
         )
-        solution = sem_fit(semoptimizer, model_ml_multigroup2)
+        solution = fit(semoptimizer, model_ml_multigroup2)
         update_estimate!(partable_miss, solution)
         test_estimates(
             partable_miss,
@@ -340,7 +340,7 @@ if !isnothing(specification_miss_g1)
     end
 
     @testset "fitmeasures/se_fiml" begin
-        solution = sem_fit(semoptimizer, model_ml_multigroup)
+        solution = fit(semoptimizer, model_ml_multigroup)
         test_fitmeasures(
             fit_measures(solution),
             solution_lav[:fitmeasures_fiml];
@@ -357,7 +357,7 @@ if !isnothing(specification_miss_g1)
             lav_groups = Dict(:Pasteur => 1, :Grant_White => 2),
         )
 
-        solution = sem_fit(semoptimizer, model_ml_multigroup2)
+        solution = fit(semoptimizer, model_ml_multigroup2)
         test_fitmeasures(
             fit_measures(solution),
             solution_lav[:fitmeasures_fiml];
diff --git a/test/examples/political_democracy/by_parts.jl b/test/examples/political_democracy/by_parts.jl
index ddbbfc3fa..3397b5f0a 100644
--- a/test/examples/political_democracy/by_parts.jl
+++ b/test/examples/political_democracy/by_parts.jl
@@ -70,7 +70,7 @@ solution_names = Symbol.("parameter_estimates_" .* ["ml", "ls", "ml", "ml"])
 for (model, name, solution_name) in zip(models, model_names, solution_names)
     try
         @testset "$(name)_solution" begin
-            solution = sem_fit(optimizer_obj, model)
+            solution = fit(optimizer_obj, model)
             update_estimate!(partable, solution)
             test_estimates(partable, solution_lav[solution_name]; atol = 1e-2)
         end
@@ -79,9 +79,9 @@ for (model, name, solution_name) in zip(models, model_names, solution_names)
 end
 
 @testset "ridge_solution" begin
-    solution_ridge = sem_fit(optimizer_obj, model_ridge)
-    solution_ml = sem_fit(optimizer_obj, model_ml)
-    # solution_ridge_id = sem_fit(optimizer_obj, model_ridge_id)
+    solution_ridge = fit(optimizer_obj, model_ridge)
+    solution_ml = fit(optimizer_obj, model_ml)
+    # solution_ridge_id = fit(optimizer_obj, model_ridge_id)
     @test solution_ridge.minimum < solution_ml.minimum + 1
 end
 
@@ -97,8 +97,8 @@ end
 end
 
 @testset "ml_solution_weighted" begin
-    solution_ml = sem_fit(optimizer_obj, model_ml)
-    solution_ml_weighted = sem_fit(optimizer_obj, model_ml_weighted)
+    solution_ml = fit(optimizer_obj, model_ml)
+    solution_ml_weighted = fit(optimizer_obj, model_ml_weighted)
     @test solution(solution_ml) ≈ solution(solution_ml_weighted) rtol = 1e-3
     @test nsamples(model_ml) * StructuralEquationModels.minimum(solution_ml) ≈
           StructuralEquationModels.minimum(solution_ml_weighted) rtol = 1e-6
@@ -109,7 +109,7 @@ end
 ############################################################################################
 
 @testset "fitmeasures/se_ml" begin
-    solution_ml = sem_fit(optimizer_obj, model_ml)
+    solution_ml = fit(optimizer_obj, model_ml)
     test_fitmeasures(fit_measures(solution_ml), solution_lav[:fitmeasures_ml]; atol = 1e-3)
 
     update_se_hessian!(partable, solution_ml)
@@ -123,7 +123,7 @@ end
 end
 
 @testset "fitmeasures/se_ls" begin
-    solution_ls = sem_fit(optimizer_obj, model_ls_sym)
+    solution_ls = fit(optimizer_obj, model_ls_sym)
     fm = fit_measures(solution_ls)
     test_fitmeasures(
         fm,
@@ -176,13 +176,13 @@ if opt_engine == :Optim
     end
 
     @testset "ml_solution_hessian" begin
-        solution = sem_fit(optimizer_obj, model_ml)
+        solution = fit(optimizer_obj, model_ml)
         update_estimate!(partable, solution)
         test_estimates(partable, solution_lav[:parameter_estimates_ml]; atol = 1e-2)
     end
 
     @testset "ls_solution_hessian" begin
-        solution = sem_fit(optimizer_obj, model_ls)
+        solution = fit(optimizer_obj, model_ls)
         update_estimate!(partable, solution)
         test_estimates(
             partable,
@@ -254,7 +254,7 @@ solution_names = Symbol.("parameter_estimates_" .* ["ml", "ls", "ml"] .* "_mean"
 for (model, name, solution_name) in zip(models, model_names, solution_names)
     try
         @testset "$(name)_solution_mean" begin
-            solution = sem_fit(optimizer_obj, model)
+            solution = fit(optimizer_obj, model)
             update_estimate!(partable_mean, solution)
             test_estimates(partable_mean, solution_lav[solution_name]; atol = 1e-2)
         end
@@ -267,7 +267,7 @@ end
 ############################################################################################
 
 @testset "fitmeasures/se_ml_mean" begin
-    solution_ml = sem_fit(optimizer_obj, model_ml)
+    solution_ml = fit(optimizer_obj, model_ml)
     test_fitmeasures(
         fit_measures(solution_ml),
         solution_lav[:fitmeasures_ml_mean];
@@ -285,7 +285,7 @@ end
 end
 
 @testset "fitmeasures/se_ls_mean" begin
-    solution_ls = sem_fit(optimizer_obj, model_ls)
+    solution_ls = fit(optimizer_obj, model_ls)
     fm = fit_measures(solution_ls)
     test_fitmeasures(
         fm,
@@ -336,13 +336,13 @@ end
 ############################################################################################
 
 @testset "fiml_solution" begin
-    solution = sem_fit(optimizer_obj, model_ml)
+    solution = fit(optimizer_obj, model_ml)
     update_estimate!(partable_mean, solution)
     test_estimates(partable_mean, solution_lav[:parameter_estimates_fiml]; atol = 1e-2)
 end
 
 @testset "fiml_solution_symbolic" begin
-    solution = sem_fit(optimizer_obj, model_ml_sym)
+    solution = fit(optimizer_obj, model_ml_sym)
     update_estimate!(partable_mean, solution)
     test_estimates(partable_mean, solution_lav[:parameter_estimates_fiml]; atol = 1e-2)
 end
@@ -352,7 +352,7 @@ end
 ############################################################################################
 
 @testset "fitmeasures/se_fiml" begin
-    solution_ml = sem_fit(optimizer_obj, model_ml)
+    solution_ml = fit(optimizer_obj, model_ml)
     test_fitmeasures(
         fit_measures(solution_ml),
         solution_lav[:fitmeasures_fiml];
diff --git a/test/examples/political_democracy/constraints.jl b/test/examples/political_democracy/constraints.jl
index fb2116023..cc1b0874d 100644
--- a/test/examples/political_democracy/constraints.jl
+++ b/test/examples/political_democracy/constraints.jl
@@ -39,14 +39,14 @@ constrained_optimizer = SemOptimizer(;
 ############################################################################################
 
 @testset "ml_solution_maxeval" begin
-    solution_maxeval = sem_fit(model_ml, engine = :NLopt, options = Dict(:maxeval => 10))
+    solution_maxeval = fit(model_ml, engine = :NLopt, options = Dict(:maxeval => 10))
 
     @test solution_maxeval.optimization_result.problem.numevals == 10
     @test solution_maxeval.optimization_result.result[3] == :MAXEVAL_REACHED
 end
 
 @testset "ml_solution_constrained" begin
-    solution_constrained = sem_fit(constrained_optimizer, model_ml)
+    solution_constrained = fit(constrained_optimizer, model_ml)
 
     @test solution_constrained.solution[31] * solution_constrained.solution[30] >=
           (0.6 - 1e-8)
diff --git a/test/examples/political_democracy/constructor.jl b/test/examples/political_democracy/constructor.jl
index 15999294d..7a8adc72e 100644
--- a/test/examples/political_democracy/constructor.jl
+++ b/test/examples/political_democracy/constructor.jl
@@ -75,7 +75,7 @@ solution_names = Symbol.("parameter_estimates_" .* ["ml", "ml", "ls", "ml", "ml"
 for (model, name, solution_name) in zip(models, model_names, solution_names)
     try
         @testset "$(name)_solution" begin
-            solution = sem_fit(semoptimizer, model)
+            solution = fit(semoptimizer, model)
             update_estimate!(partable, solution)
             test_estimates(partable, solution_lav[solution_name]; atol = 1e-2)
         end
@@ -84,9 +84,9 @@ for (model, name, solution_name) in zip(models, model_names, solution_names)
 end
 
 @testset "ridge_solution" begin
-    solution_ridge = sem_fit(semoptimizer, model_ridge)
-    solution_ml = sem_fit(semoptimizer, model_ml)
-    # solution_ridge_id = sem_fit(semoptimizer, model_ridge_id)
+    solution_ridge = fit(semoptimizer, model_ridge)
+    solution_ml = fit(semoptimizer, model_ml)
+    # solution_ridge_id = fit(semoptimizer, model_ridge_id)
     @test abs(solution_ridge.minimum - solution_ml.minimum) < 1
 end
 
@@ -102,8 +102,8 @@ end
 end
 
 @testset "ml_solution_weighted" begin
-    solution_ml = sem_fit(semoptimizer, model_ml)
-    solution_ml_weighted = sem_fit(semoptimizer, model_ml_weighted)
+    solution_ml = fit(semoptimizer, model_ml)
+    solution_ml_weighted = fit(semoptimizer, model_ml_weighted)
     @test isapprox(solution(solution_ml), solution(solution_ml_weighted), rtol = 1e-3)
     @test isapprox(
         nsamples(model_ml) * StructuralEquationModels.minimum(solution_ml),
@@ -117,7 +117,7 @@ end
 ############################################################################################
 
 @testset "fitmeasures/se_ml" begin
-    solution_ml = sem_fit(semoptimizer, model_ml)
+    solution_ml = fit(semoptimizer, model_ml)
     test_fitmeasures(fit_measures(solution_ml), solution_lav[:fitmeasures_ml]; atol = 1e-3)
 
     update_se_hessian!(partable, solution_ml)
@@ -131,7 +131,7 @@ end
 end
 
 @testset "fitmeasures/se_ls" begin
-    solution_ls = sem_fit(semoptimizer, model_ls_sym)
+    solution_ls = fit(semoptimizer, model_ls_sym)
     fm = fit_measures(solution_ls)
     test_fitmeasures(
         fm,
@@ -182,8 +182,8 @@ end
         obs_colnames = colnames,
     )
     # fit models
-    sol_ml = solution(sem_fit(semoptimizer, model_ml_new))
-    sol_ml_sym = solution(sem_fit(semoptimizer, model_ml_sym_new))
+    sol_ml = solution(fit(semoptimizer, model_ml_new))
+    sol_ml_sym = solution(fit(semoptimizer, model_ml_sym_new))
     # check solution
     @test maximum(abs.(sol_ml - params)) < 0.01
     @test maximum(abs.(sol_ml_sym - params)) < 0.01
@@ -225,13 +225,13 @@ if opt_engine == :Optim
     end
 
     @testset "ml_solution_hessian" begin
-        solution = sem_fit(semoptimizer, model_ml)
+        solution = fit(semoptimizer, model_ml)
         update_estimate!(partable, solution)
         test_estimates(partable, solution_lav[:parameter_estimates_ml]; atol = 1e-2)
     end
 
     @testset "ls_solution_hessian" begin
-        solution = sem_fit(semoptimizer, model_ls)
+        solution = fit(semoptimizer, model_ls)
         update_estimate!(partable, solution)
         test_estimates(
             partable,
@@ -296,7 +296,7 @@ solution_names = Symbol.("parameter_estimates_" .* ["ml", "ml", "ls", "ml"] .* "
 for (model, name, solution_name) in zip(models, model_names, solution_names)
     try
         @testset "$(name)_solution_mean" begin
-            solution = sem_fit(semoptimizer, model)
+            solution = fit(semoptimizer, model)
             update_estimate!(partable_mean, solution)
             test_estimates(partable_mean, solution_lav[solution_name]; atol = 1e-2)
         end
@@ -309,7 +309,7 @@ end
 ############################################################################################
 
 @testset "fitmeasures/se_ml_mean" begin
-    solution_ml = sem_fit(semoptimizer, model_ml)
+    solution_ml = fit(semoptimizer, model_ml)
     test_fitmeasures(
         fit_measures(solution_ml),
         solution_lav[:fitmeasures_ml_mean];
@@ -327,7 +327,7 @@ end
 end
 
 @testset "fitmeasures/se_ls_mean" begin
-    solution_ls = sem_fit(semoptimizer, model_ls)
+    solution_ls = fit(semoptimizer, model_ls)
     fm = fit_measures(solution_ls)
     test_fitmeasures(
         fm,
@@ -381,8 +381,8 @@ end
         meanstructure = true,
     )
     # fit models
-    sol_ml = solution(sem_fit(semoptimizer, model_ml_new))
-    sol_ml_sym = solution(sem_fit(semoptimizer, model_ml_sym_new))
+    sol_ml = solution(fit(semoptimizer, model_ml_new))
+    sol_ml_sym = solution(fit(semoptimizer, model_ml_sym_new))
     # check solution
     @test maximum(abs.(sol_ml - params)) < 0.01
     @test maximum(abs.(sol_ml_sym - params)) < 0.01
@@ -427,13 +427,13 @@ end
 ############################################################################################
 
 @testset "fiml_solution" begin
-    solution = sem_fit(semoptimizer, model_ml)
+    solution = fit(semoptimizer, model_ml)
     update_estimate!(partable_mean, solution)
     test_estimates(partable_mean, solution_lav[:parameter_estimates_fiml]; atol = 1e-2)
 end
 
 @testset "fiml_solution_symbolic" begin
-    solution = sem_fit(semoptimizer, model_ml_sym)
+    solution = fit(semoptimizer, model_ml_sym)
     update_estimate!(partable_mean, solution)
     test_estimates(partable_mean, solution_lav[:parameter_estimates_fiml]; atol = 1e-2)
 end
@@ -443,7 +443,7 @@ end
 ############################################################################################
 
 @testset "fitmeasures/se_fiml" begin
-    solution_ml = sem_fit(semoptimizer, model_ml)
+    solution_ml = fit(semoptimizer, model_ml)
     test_fitmeasures(
         fit_measures(solution_ml),
         solution_lav[:fitmeasures_fiml];
diff --git a/test/examples/proximal/l0.jl b/test/examples/proximal/l0.jl
index da20f3901..374f8e58a 100644
--- a/test/examples/proximal/l0.jl
+++ b/test/examples/proximal/l0.jl
@@ -35,7 +35,7 @@ ram_mat = RAMMatrices(partable)
 
 model = Sem(specification = partable, data = dat, loss = SemML)
 
-fit = sem_fit(model)
+sem_fit = fit(model)
 
 # use l0 from ProximalSEM
 # regularized
@@ -44,11 +44,11 @@ prox_operator =
 
 model_prox = Sem(specification = partable, data = dat, loss = SemML)
 
-fit_prox = sem_fit(model_prox, engine = :Proximal, operator_g = prox_operator)
+fit_prox = fit(model_prox, engine = :Proximal, operator_g = prox_operator)
 
 @testset "l0 | solution_unregularized" begin
     @test fit_prox.optimization_result.result[:iterations] < 1000
-    @test maximum(abs.(solution(fit) - solution(fit_prox))) < 0.002
+    @test maximum(abs.(solution(sem_fit) - solution(fit_prox))) < 0.002
 end
 
 # regularized
@@ -56,12 +56,12 @@ prox_operator = SlicedSeparableSum((NormL0(0.0), NormL0(100.0)), ([1:30], [31]))
 
 model_prox = Sem(specification = partable, data = dat, loss = SemML)
 
-fit_prox = sem_fit(model_prox, engine = :Proximal, operator_g = prox_operator)
+fit_prox = fit(model_prox, engine = :Proximal, operator_g = prox_operator)
 
 @testset "l0 | solution_regularized" begin
     @test fit_prox.optimization_result.result[:iterations] < 1000
     @test solution(fit_prox)[31] == 0.0
     @test abs(
-        StructuralEquationModels.minimum(fit_prox) - StructuralEquationModels.minimum(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 314453df4..beb5cf529 100644
--- a/test/examples/proximal/lasso.jl
+++ b/test/examples/proximal/lasso.jl
@@ -35,18 +35,18 @@ ram_mat = RAMMatrices(partable)
 
 model = Sem(specification = partable, data = dat, loss = SemML)
 
-fit = sem_fit(model)
+sem_fit = fit(model)
 
 # use lasso from ProximalSEM
 λ = zeros(31)
 
 model_prox = Sem(specification = partable, data = dat, loss = SemML)
 
-fit_prox = sem_fit(model_prox, engine = :Proximal, operator_g = NormL1(λ))
+fit_prox = fit(model_prox, engine = :Proximal, operator_g = NormL1(λ))
 
 @testset "lasso | solution_unregularized" begin
     @test fit_prox.optimization_result.result[:iterations] < 1000
-    @test maximum(abs.(solution(fit) - solution(fit_prox))) < 0.002
+    @test maximum(abs.(solution(sem_fit) - solution(fit_prox))) < 0.002
 end
 
 λ = zeros(31);
@@ -54,11 +54,11 @@ end
 
 model_prox = Sem(specification = partable, data = dat, loss = SemML)
 
-fit_prox = sem_fit(model_prox, engine = :Proximal, operator_g = NormL1(λ))
+fit_prox = fit(model_prox, engine = :Proximal, operator_g = NormL1(λ))
 
 @testset "lasso | solution_regularized" begin
     @test fit_prox.optimization_result.result[:iterations] < 1000
-    @test all(solution(fit_prox)[16:20] .< solution(fit)[16:20])
+    @test all(solution(fit_prox)[16:20] .< solution(sem_fit)[16:20])
     @test StructuralEquationModels.minimum(fit_prox) -
-          StructuralEquationModels.minimum(fit) < 0.03
+          StructuralEquationModels.minimum(sem_fit) < 0.03
 end
diff --git a/test/examples/proximal/ridge.jl b/test/examples/proximal/ridge.jl
index 3d116dcd4..fd7ae113d 100644
--- a/test/examples/proximal/ridge.jl
+++ b/test/examples/proximal/ridge.jl
@@ -35,7 +35,7 @@ ram_mat = RAMMatrices(partable)
 
 model = Sem(specification = partable, data = dat, loss = SemML)
 
-fit = sem_fit(model)
+sem_fit = fit(model)
 
 # use ridge from StructuralEquationModels
 model_ridge = Sem(
@@ -46,7 +46,7 @@ model_ridge = Sem(
     which_ridge = 16:20,
 )
 
-solution_ridge = sem_fit(model_ridge)
+solution_ridge = fit(model_ridge)
 
 # use ridge from ProximalSEM; SqrNormL2 uses λ/2 as penalty
 λ = zeros(31);
@@ -54,7 +54,7 @@ solution_ridge = sem_fit(model_ridge)
 
 model_prox = Sem(specification = partable, data = dat, loss = SemML)
 
-solution_prox = @suppress sem_fit(model_prox, engine = :Proximal, operator_g = SqrNormL2(λ))
+solution_prox = @suppress fit(model_prox, engine = :Proximal, operator_g = SqrNormL2(λ))
 
 @testset "ridge_solution" begin
     @test isapprox(solution_prox.solution, solution_ridge.solution; rtol = 1e-3)
diff --git a/test/examples/recover_parameters/recover_parameters_twofact.jl b/test/examples/recover_parameters/recover_parameters_twofact.jl
index fd27dafd2..a3e426cbc 100644
--- a/test/examples/recover_parameters/recover_parameters_twofact.jl
+++ b/test/examples/recover_parameters/recover_parameters_twofact.jl
@@ -73,6 +73,6 @@ optimizer = SemOptimizerOptim(
     Optim.Options(; f_tol = 1e-10, x_tol = 1.5e-8),
 )
 
-solution_ml = sem_fit(optimizer, model_ml)
+solution_ml = fit(optimizer, model_ml)
 
 @test true_val ≈ solution(solution_ml) atol = 0.05
diff --git a/test/unit_tests/StatsAPI.jl b/test/unit_tests/StatsAPI.jl
index 18811dcc2..8648fc363 100644
--- a/test/unit_tests/StatsAPI.jl
+++ b/test/unit_tests/StatsAPI.jl
@@ -9,7 +9,7 @@ model = Sem(
     specification = partable,
     data = data
 )
-model_fit = sem_fit(model)
+model_fit = fit(model)
 
 @testset "params" begin
     out = [NaN]
@@ -20,3 +20,10 @@ end
     @test param_labels(partable) == [:θ_1] == coefnames(partable)
 end
 
+@testset "nobs" begin
+    @test nobs(model) == nsamples(model)
+end
+
+@testset "coeftable" begin
+    @test_throws "StructuralEquationModels does not support" coeftable(model) 
+end
\ No newline at end of file
diff --git a/test/unit_tests/bootstrap.jl b/test/unit_tests/bootstrap.jl
index f30092865..a2d5b6832 100644
--- a/test/unit_tests/bootstrap.jl
+++ b/test/unit_tests/bootstrap.jl
@@ -1,4 +1,4 @@
-solution_ml = sem_fit(model_ml)
+solution_ml = fit(model_ml)
 bs = se_bootstrap(solution_ml; n_boot = 20)
 
 update_se_hessian!(partable, solution_ml)
diff --git a/test/unit_tests/sorting.jl b/test/unit_tests/sorting.jl
index 0908a6497..3c61e13c4 100644
--- a/test/unit_tests/sorting.jl
+++ b/test/unit_tests/sorting.jl
@@ -11,7 +11,7 @@ model_ml_sorted = Sem(specification = partable, data = dat)
 end
 
 @testset "ml_solution_sorted" begin
-    solution_ml_sorted = sem_fit(model_ml_sorted)
+    solution_ml_sorted = fit(model_ml_sorted)
     update_estimate!(partable, solution_ml_sorted)
     @test test_estimates(par_ml, partable, 0.01)
 end