diff --git a/src/StructuralEquationModels.jl b/src/StructuralEquationModels.jl
index e9d529ded..547549ccb 100644
--- a/src/StructuralEquationModels.jl
+++ b/src/StructuralEquationModels.jl
@@ -4,6 +4,7 @@ using LinearAlgebra,
     Optim,
     NLSolversBase,
     Statistics,
+    StatsBase,
     SparseArrays,
     Symbolics,
     NLopt,
@@ -18,6 +19,8 @@ using LinearAlgebra,
 import DataFrames: DataFrame
 export StenoGraphs, @StenoGraph, meld
 
+const SEM = StructuralEquationModels
+
 # type hierarchy
 include("types.jl")
 include("objective_gradient_hessian.jl")
diff --git a/src/additional_functions/helper.jl b/src/additional_functions/helper.jl
index 62100907d..05253cd8c 100644
--- a/src/additional_functions/helper.jl
+++ b/src/additional_functions/helper.jl
@@ -10,7 +10,7 @@ function neumann_series(mat::SparseMatrixCSC)
     return inverse
 end
 
-#= 
+#=
 function make_onelement_array(A)
     isa(A, Array) ? nothing : (A = [A])
     return A
@@ -108,11 +108,8 @@ function sparse_outer_mul!(C, A, B::Vector, ind) #computes A*S*B -> C, where ind
 end
 
 function cov_and_mean(rows; corrected = false)
-    data = transpose(reduce(hcat, rows))
-    size(rows, 1) > 1 ? obs_cov = Statistics.cov(data; corrected = corrected) :
-    obs_cov = reshape([0.0], 1, 1)
-    obs_mean = vec(Statistics.mean(data, dims = 1))
-    return obs_cov, obs_mean
+    obs_mean, obs_cov = StatsBase.mean_and_cov(reduce(hcat, rows), 2, corrected = corrected)
+    return obs_cov, vec(obs_mean)
 end
 
 function duplication_matrix(nobs)
diff --git a/src/additional_functions/parameters.jl b/src/additional_functions/parameters.jl
index 4305d56e7..1b7d4643f 100644
--- a/src/additional_functions/parameters.jl
+++ b/src/additional_functions/parameters.jl
@@ -1,5 +1,14 @@
-function fill_A_S_M(A, S, M, A_indices, S_indices, M_indices, parameters)
-    for (iA, iS, par) in zip(A_indices, S_indices, parameters)
+# fill A, S, and M matrices with the parameter values according to the parameters map
+function fill_A_S_M!(
+    A::AbstractMatrix,
+    S::AbstractMatrix,
+    M::Union{AbstractVector, Nothing},
+    A_indices::AbstractArrayParamsMap,
+    S_indices::AbstractArrayParamsMap,
+    M_indices::Union{AbstractArrayParamsMap, Nothing},
+    parameters::AbstractVector,
+)
+    @inbounds for (iA, iS, par) in zip(A_indices, S_indices, parameters)
         for index_A in iA
             A[index_A] = par
         end
@@ -10,7 +19,7 @@ function fill_A_S_M(A, S, M, A_indices, S_indices, M_indices, parameters)
     end
 
     if !isnothing(M)
-        for (iM, par) in zip(M_indices, parameters)
+        @inbounds for (iM, par) in zip(M_indices, parameters)
             for index_M in iM
                 M[index_M] = par
             end
@@ -18,14 +27,20 @@ function fill_A_S_M(A, S, M, A_indices, S_indices, M_indices, parameters)
     end
 end
 
-function get_parameter_indices(parameters, M; linear = true, kwargs...)
-    M_indices = [findall(x -> (x == par), M) for par in parameters]
-
-    if linear
-        M_indices = cartesian2linear.(M_indices, [M])
+# build the map from the index of the parameter to the linear indices
+# of this parameter occurences in M
+# returns ArrayParamsMap object
+function array_parameters_map(parameters::AbstractVector, M::AbstractArray)
+    params_index = Dict(param => i for (i, param) in enumerate(parameters))
+    T = Base.eltype(eachindex(M))
+    res = [Vector{T}() for _ in eachindex(parameters)]
+    for (i, val) in enumerate(M)
+        par_ind = get(params_index, val, nothing)
+        if !isnothing(par_ind)
+            push!(res[par_ind], i)
+        end
     end
-
-    return M_indices
+    return res
 end
 
 function eachindex_lower(M; linear_indices = false, kwargs...)
@@ -49,9 +64,6 @@ function linear2cartesian(ind_lin, dims)
     return ind_cart
 end
 
-cartesian2linear(ind_cart, A::AbstractArray) = cartesian2linear(ind_cart, size(A))
-linear2cartesian(ind_linear, A::AbstractArray) = linear2cartesian(ind_linear, size(A))
-
 function set_constants!(M, M_pre)
     for index in eachindex(M)
         δ = tryparse(Float64, string(M[index]))
@@ -85,12 +97,18 @@ function get_matrix_derivative(M_indices, parameters, n_long)
     return ∇M
 end
 
-function fill_matrix(M, M_indices, parameters)
+# fill M with parameters
+function fill_matrix!(
+    M::AbstractMatrix,
+    M_indices::AbstractArrayParamsMap,
+    parameters::AbstractVector,
+)
     for (iM, par) in zip(M_indices, parameters)
         for index_M in iM
             M[index_M] = par
         end
     end
+    return M
 end
 
 function get_partition(A_indices, S_indices)
diff --git a/src/additional_functions/start_val/start_fabin3.jl b/src/additional_functions/start_val/start_fabin3.jl
index eaf899127..bce6471c7 100644
--- a/src/additional_functions/start_val/start_fabin3.jl
+++ b/src/additional_functions/start_val/start_fabin3.jl
@@ -1,6 +1,6 @@
 """
     start_fabin3(model)
-    
+
 Return a vector of FABIN 3 starting values (see Hägglund 1982).
 Not available for ensemble models.
 """
@@ -58,8 +58,8 @@ function start_fabin3(ram_matrices::RAMMatrices, Σ, μ)
         if !isnothing(M)
             in_M = length.(M_ind) .!= 0
             in_any = in_A .| in_S .| in_M
-        else 
-            in_any = in_A .| in_S 
+        else
+            in_any = in_A .| in_S
         end
 
         if !all(in_any)
diff --git a/src/diff/Empty.jl b/src/diff/Empty.jl
index 03923a851..57fa9ee98 100644
--- a/src/diff/Empty.jl
+++ b/src/diff/Empty.jl
@@ -34,5 +34,5 @@ update_observed(optimizer::SemOptimizerEmpty, observed::SemObserved; kwargs...)
 ############################################################################################
 
 function Base.show(io::IO, struct_inst::SemOptimizerEmpty)
-    StructuralEquationModels.print_type_name(io, struct_inst)
+    print_type_name(io, struct_inst)
 end
diff --git a/src/frontend/fit/fitmeasures/n_obs.jl b/src/frontend/fit/fitmeasures/n_obs.jl
index 3a403ff3d..cd4bdca30 100644
--- a/src/frontend/fit/fitmeasures/n_obs.jl
+++ b/src/frontend/fit/fitmeasures/n_obs.jl
@@ -13,4 +13,4 @@ n_obs(sem_fit::SemFit) = n_obs(sem_fit.model)
 
 n_obs(model::AbstractSemSingle) = n_obs(model.observed)
 
-n_obs(model::SemEnsemble) = sum(n_obs.(model.sems))
+n_obs(model::SemEnsemble) = sum(n_obs, model.sems)
diff --git a/src/frontend/fit/standard_errors/hessian.jl b/src/frontend/fit/standard_errors/hessian.jl
index bb84efe71..396d3b98c 100644
--- a/src/frontend/fit/standard_errors/hessian.jl
+++ b/src/frontend/fit/standard_errors/hessian.jl
@@ -4,9 +4,9 @@
 Return hessian based standard errors.
 
 # Arguments
-- `hessian`: how to compute the hessian. Options are 
+- `hessian`: how to compute the hessian. Options are
     - `:analytic`: (only if an analytic hessian for the model can be computed)
-    - `:finitediff`: for finite difference approximation 
+    - `:finitediff`: for finite difference approximation
 """
 function se_hessian(sem_fit::SemFit; hessian = :finitediff)
     c = H_scaling(sem_fit.model)
@@ -17,7 +17,7 @@ function se_hessian(sem_fit::SemFit; hessian = :finitediff)
         hessian!(H, sem_fit.model, sem_fit.solution)
     elseif hessian == :finitediff
         H = FiniteDiff.finite_difference_hessian(
-            x -> objective!(sem_fit.model, x),
+            Base.Fix1(objective!, sem_fit.model),
             sem_fit.solution,
         )
     elseif hessian == :optimizer
@@ -33,7 +33,7 @@ function se_hessian(sem_fit::SemFit; hessian = :finitediff)
             ),
         )
     else
-        throw(ArgumentError("I dont know how to compute `$hessian` standard-errors"))
+        throw(ArgumentError("I don't know how to compute `$hessian` standard-errors"))
     end
 
     invH = c * inv(H)
diff --git a/src/frontend/specification/EnsembleParameterTable.jl b/src/frontend/specification/EnsembleParameterTable.jl
index fcbf188fa..79283953f 100644
--- a/src/frontend/specification/EnsembleParameterTable.jl
+++ b/src/frontend/specification/EnsembleParameterTable.jl
@@ -27,7 +27,7 @@ function Dict(partable::EnsembleParameterTable)
 end
 
 #= function DataFrame(
-        partable::ParameterTable; 
+        partable::ParameterTable;
         columns = nothing)
     if isnothing(columns) columns = keys(partable.columns) end
     out = DataFrame([key => partable.columns[key] for key in columns])
diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl
index fd3ed71da..1910d666e 100644
--- a/src/frontend/specification/ParameterTable.jl
+++ b/src/frontend/specification/ParameterTable.jl
@@ -224,9 +224,9 @@ update_partable!(partable::AbstractParameterTable, sem_fit::SemFit, vec, column)
 # update estimates -------------------------------------------------------------------------
 """
     update_estimate!(
-        partable::AbstractParameterTable, 
+        partable::AbstractParameterTable,
         sem_fit::SemFit)
-    
+
 Write parameter estimates from `sem_fit` to the `:estimate` column of `partable`
 """
 update_estimate!(partable::AbstractParameterTable, sem_fit::SemFit) =
@@ -236,7 +236,7 @@ update_estimate!(partable::AbstractParameterTable, sem_fit::SemFit) =
 """
     update_start!(partable::AbstractParameterTable, sem_fit::SemFit)
     update_start!(partable::AbstractParameterTable, model::AbstractSem, start_val; kwargs...)
-    
+
 Write starting values from `sem_fit` or `start_val` to the `:estimate` column of `partable`.
 
 # Arguments
@@ -262,10 +262,10 @@ end
 # update partable standard errors ----------------------------------------------------------
 """
     update_se_hessian!(
-        partable::AbstractParameterTable, 
-        sem_fit::SemFit; 
+        partable::AbstractParameterTable,
+        sem_fit::SemFit;
         hessian = :finitediff)
-    
+
 Write hessian standard errors computed for `sem_fit` to the `:se` column of `partable`
 
 # Arguments
diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl
index 8124595ed..e0fcc575c 100644
--- a/src/frontend/specification/RAMMatrices.jl
+++ b/src/frontend/specification/RAMMatrices.jl
@@ -2,11 +2,15 @@
 ### Type
 ############################################################################################
 
+# map from parameter index to linear indices of matrix/vector positions where it occurs
+AbstractArrayParamsMap = AbstractVector{<:AbstractVector{<:Integer}}
+ArrayParamsMap = Vector{Vector{Int}}
+
 struct RAMMatrices
-    A_ind::Any
-    S_ind::Any
-    F_ind::Any
-    M_ind::Any
+    A_ind::ArrayParamsMap
+    S_ind::ArrayParamsMap
+    F_ind::Vector{Int}
+    M_ind::Union{ArrayParamsMap, Nothing}
     parameters::Any
     colnames::Any
     constants::Any
@@ -18,9 +22,9 @@ end
 ############################################################################################
 
 function RAMMatrices(; A, S, F, M = nothing, parameters, colnames)
-    A_indices = get_parameter_indices(parameters, A)
-    S_indices = get_parameter_indices(parameters, S)
-    isnothing(M) ? M_indices = nothing : M_indices = get_parameter_indices(parameters, M)
+    A_indices = array_parameters_map(parameters, A)
+    S_indices = array_parameters_map(parameters, S)
+    M_indices = !isnothing(M) ? array_parameters_map(parameters, M) : nothing
     F_indices = findall([any(isone.(col)) for col in eachcol(F)])
     constants = get_RAMConstants(A, S, M)
     return RAMMatrices(
@@ -50,7 +54,7 @@ end
 import Base.==
 
 function ==(c1::RAMConstant, c2::RAMConstant)
-    res = ((c1.matrix == c2.matrix) & (c1.index == c2.index) & (c1.value == c2.value))
+    res = ((c1.matrix == c2.matrix) && (c1.index == c2.index) && (c1.value == c2.value))
     return res
 end
 
@@ -425,13 +429,13 @@ end
 
 function ==(mat1::RAMMatrices, mat2::RAMMatrices)
     res = (
-        (mat1.A_ind == mat2.A_ind) &
-        (mat1.S_ind == mat2.S_ind) &
-        (mat1.F_ind == mat2.F_ind) &
-        (mat1.M_ind == mat2.M_ind) &
-        (mat1.parameters == mat2.parameters) &
-        (mat1.colnames == mat2.colnames) &
-        (mat1.size_F == mat2.size_F) &
+        (mat1.A_ind == mat2.A_ind) &&
+        (mat1.S_ind == mat2.S_ind) &&
+        (mat1.F_ind == mat2.F_ind) &&
+        (mat1.M_ind == mat2.M_ind) &&
+        (mat1.parameters == mat2.parameters) &&
+        (mat1.colnames == mat2.colnames) &&
+        (mat1.size_F == mat2.size_F) &&
         (mat1.constants == mat2.constants)
     )
     return res
diff --git a/src/frontend/specification/Sem.jl b/src/frontend/specification/Sem.jl
index a627b5f09..208ef3000 100644
--- a/src/frontend/specification/Sem.jl
+++ b/src/frontend/specification/Sem.jl
@@ -9,11 +9,11 @@ function Sem(;
     optimizer::D = SemOptimizerOptim,
     kwargs...,
 ) where {O, I, L, D}
-    kwargs = Dict{Symbol, Any}(kwargs...)
+    kwdict = Dict{Symbol, Any}(kwargs...)
 
-    set_field_type_kwargs!(kwargs, observed, imply, loss, optimizer, O, I, D)
+    set_field_type_kwargs!(kwdict, observed, imply, loss, optimizer, O, I, D)
 
-    observed, imply, loss, optimizer = get_fields!(kwargs, observed, imply, loss, optimizer)
+    observed, imply, loss, optimizer = get_fields!(kwdict, observed, imply, loss, optimizer)
 
     sem = Sem(observed, imply, loss, optimizer)
 
@@ -27,11 +27,11 @@ function SemFiniteDiff(;
     optimizer::D = SemOptimizerOptim,
     kwargs...,
 ) where {O, I, L, D}
-    kwargs = Dict{Symbol, Any}(kwargs...)
+    kwdict = Dict{Symbol, Any}(kwargs...)
 
-    set_field_type_kwargs!(kwargs, observed, imply, loss, optimizer, O, I, D)
+    set_field_type_kwargs!(kwdict, observed, imply, loss, optimizer, O, I, D)
 
-    observed, imply, loss, optimizer = get_fields!(kwargs, observed, imply, loss, optimizer)
+    observed, imply, loss, optimizer = get_fields!(kwdict, observed, imply, loss, optimizer)
 
     sem = SemFiniteDiff(observed, imply, loss, optimizer)
 
diff --git a/src/imply/RAM/generic.jl b/src/imply/RAM/generic.jl
index 7bba224eb..b9d6fad69 100644
--- a/src/imply/RAM/generic.jl
+++ b/src/imply/RAM/generic.jl
@@ -222,7 +222,7 @@ gradient!(imply::RAM, par, model::AbstractSemSingle) =
 
 # objective and gradient
 function objective!(imply::RAM, parameters, model, has_meanstructure::Val{T}) where {T}
-    fill_A_S_M(
+    fill_A_S_M!(
         imply.A,
         imply.S,
         imply.M,
@@ -250,7 +250,7 @@ function gradient!(
     model::AbstractSemSingle,
     has_meanstructure::Val{T},
 ) where {T}
-    fill_A_S_M(
+    fill_A_S_M!(
         imply.A,
         imply.S,
         imply.M,
@@ -346,7 +346,7 @@ function check_acyclic(A_pre, n_par, A_indices)
     A_rand = copy(A_pre)
     randpar = rand(n_par)
 
-    fill_matrix(A_rand, A_indices, randpar)
+    fill_matrix!(A_rand, A_indices, randpar)
 
     # check if the model is acyclic
     acyclic = isone(det(I - A_rand))
diff --git a/src/imply/RAM/symbolic.jl b/src/imply/RAM/symbolic.jl
index 79faf509f..6ce7b6ae5 100644
--- a/src/imply/RAM/symbolic.jl
+++ b/src/imply/RAM/symbolic.jl
@@ -112,7 +112,7 @@ function RAMSymbolic(;
     F[CartesianIndex.(1:n_var, ram_matrices.F_ind)] .= 1.0
 
     set_RAMConstants!(A, S, M, ram_matrices.constants)
-    fill_A_S_M(A, S, M, ram_matrices.A_ind, ram_matrices.S_ind, ram_matrices.M_ind, par)
+    fill_A_S_M!(A, S, M, ram_matrices.A_ind, ram_matrices.S_ind, ram_matrices.M_ind, par)
 
     A, S, F = sparse(A), sparse(S), sparse(F)
 
diff --git a/src/imply/empty.jl b/src/imply/empty.jl
index a29fb3cac..ba8580d16 100644
--- a/src/imply/empty.jl
+++ b/src/imply/empty.jl
@@ -36,11 +36,10 @@ end
 
 function ImplyEmpty(; specification, kwargs...)
     ram_matrices = RAMMatrices(specification)
-    identifier = StructuralEquationModels.identifier(ram_matrices)
 
     n_par = length(ram_matrices.parameters)
 
-    return ImplyEmpty(identifier, n_par)
+    return ImplyEmpty(identifier(ram_matrices), n_par)
 end
 
 ############################################################################################
diff --git a/src/loss/WLS/WLS.jl b/src/loss/WLS/WLS.jl
index 3260c8579..18649ab68 100644
--- a/src/loss/WLS/WLS.jl
+++ b/src/loss/WLS/WLS.jl
@@ -10,20 +10,20 @@ Weighted least squares estimation.
 
     SemWLS(;
         observed,
-        meanstructure = false, 
-        wls_weight_matrix = nothing, 
-        wls_weight_matrix_mean = nothing, 
-        approximate_hessian = false, 
+        meanstructure = false,
+        wls_weight_matrix = nothing,
+        wls_weight_matrix_mean = nothing,
+        approximate_hessian = false,
         kwargs...)
 
 # Arguments
 - `observed`: the `SemObserved` part of the model
 - `meanstructure::Bool`: does the model have a meanstructure?
 - `approximate_hessian::Bool`: should the hessian be swapped for an approximation
-- `wls_weight_matrix`: the weight matrix for weighted least squares. 
-    Defaults to GLS estimation (``0.5*(D^T*kron(S,S)*D)`` where D is the duplication matrix 
-    and S is the inverse ob the observed covariance matrix)
-- `wls_weight_matrix_mean`: the weight matrix for the mean part of weighted least squares. 
+- `wls_weight_matrix`: the weight matrix for weighted least squares.
+    Defaults to GLS estimation (``0.5*(D^T*kron(S,S)*D)`` where D is the duplication matrix
+    and S is the inverse of the observed covariance matrix)
+- `wls_weight_matrix_mean`: the weight matrix for the mean part of weighted least squares.
     Defaults to GLS estimation (the inverse of the observed covariance matrix)
 
 # Examples
diff --git a/src/objective_gradient_hessian.jl b/src/objective_gradient_hessian.jl
index 11bc5669f..53d68ec2c 100644
--- a/src/objective_gradient_hessian.jl
+++ b/src/objective_gradient_hessian.jl
@@ -87,9 +87,10 @@ end
 
 function objective!(loss::SemLoss, par, model)
     return mapreduce(
-        fun_weight -> fun_weight[2] * objective!(fun_weight[1], par, model),
+        (fun, weight) -> weight * objective!(fun, par, model),
         +,
-        zip(loss.functions, loss.weights),
+        loss.functions,
+        loss.weights,
     )
 end
 
@@ -108,19 +109,19 @@ end
 
 function objective_gradient!(gradient, loss::SemLoss, par, model)
     return mapreduce(
-        fun_weight ->
-            objective_gradient_wrap_(gradient, fun_weight[1], par, model, fun_weight[2]),
+        (fun, weight) -> objective_gradient_wrap_(gradient, fun, par, model, weight),
         +,
-        zip(loss.functions, loss.weights),
+        loss.functions,
+        loss.weights,
     )
 end
 
 function objective_hessian!(hessian, loss::SemLoss, par, model)
     return mapreduce(
-        fun_weight ->
-            objective_hessian_wrap_(hessian, fun_weight[1], par, model, fun_weight[2]),
+        (fun, weight) -> objective_hessian_wrap_(hessian, fun, par, model, weight),
         +,
-        zip(loss.functions, loss.weights),
+        loss.functions,
+        loss.weights,
     )
 end
 
@@ -134,16 +135,11 @@ end
 
 function objective_gradient_hessian!(gradient, hessian, loss::SemLoss, par, model)
     return mapreduce(
-        fun_weight -> objective_gradient_hessian_wrap_(
-            gradient,
-            hessian,
-            fun_weight[1],
-            par,
-            model,
-            fun_weight[2],
-        ),
+        (fun, weight) ->
+            objective_gradient_hessian_wrap_(gradient, hessian, fun, par, model, weight),
         +,
-        zip(loss.functions, loss.weights),
+        loss.functions,
+        loss.weights,
     )
 end
 
@@ -174,9 +170,10 @@ end
 
 function objective!(ensemble::SemEnsemble, par)
     return mapreduce(
-        model_weight -> model_weight[2] * objective!(model_weight[1], par),
+        (model, weight) -> weight * objective!(model, par),
         +,
-        zip(ensemble.sems, ensemble.weights),
+        ensemble.sems,
+        ensemble.weights,
     )
 end
 
@@ -201,20 +198,20 @@ end
 function objective_gradient!(gradient, ensemble::SemEnsemble, par)
     fill!(gradient, zero(eltype(gradient)))
     return mapreduce(
-        model_weight ->
-            objective_gradient_wrap_(gradient, model_weight[1], par, model_weight[2]),
+        (model, weight) -> objective_gradient_wrap_(gradient, model, par, weight),
         +,
-        zip(ensemble.sems, ensemble.weights),
+        ensemble.sems,
+        ensemble.weights,
     )
 end
 
 function objective_hessian!(hessian, ensemble::SemEnsemble, par)
     fill!(hessian, zero(eltype(hessian)))
     return mapreduce(
-        model_weight ->
-            objective_hessian_wrap_(hessian, model_weight[1], par, model_weight[2]),
+        (model, weight) -> objective_hessian_wrap_(hessian, model, par, weight),
         +,
-        zip(ensemble.sems, ensemble.weights),
+        ensemble.sems,
+        ensemble.weights,
     )
 end
 
@@ -236,16 +233,11 @@ function objective_gradient_hessian!(gradient, hessian, ensemble::SemEnsemble, p
     fill!(gradient, zero(eltype(gradient)))
     fill!(hessian, zero(eltype(hessian)))
     return mapreduce(
-        model_weight -> objective_gradient_hessian_wrap_(
-            gradient,
-            hessian,
-            model_weight[1],
-            par,
-            model,
-            model_weight[2],
-        ),
+        (model, weight) ->
+            objective_gradient_hessian_wrap_(gradient, hessian, model, par, model, weight),
         +,
-        zip(ensemble.sems, ensemble.weights),
+        ensemble.sems,
+        ensemble.weights,
     )
 end
 
diff --git a/src/observed/covariance.jl b/src/observed/covariance.jl
index aba1ef629..5785ba404 100644
--- a/src/observed/covariance.jl
+++ b/src/observed/covariance.jl
@@ -45,7 +45,7 @@ struct SemObservedCovariance{B, C, D, O} <: SemObserved
     n_man::D
     n_obs::O
 end
-using StructuralEquationModels
+
 function SemObservedCovariance(;
     specification,
     obs_cov,
@@ -76,7 +76,8 @@ function SemObservedCovariance(;
 
     if !isnothing(spec_colnames)
         obs_cov = reorder_obs_cov(obs_cov, spec_colnames, obs_colnames)
-        obs_mean = reorder_obs_mean(obs_mean, spec_colnames, obs_colnames)
+        isnothing(obs_mean) ||
+            (obs_mean = reorder_obs_mean(obs_mean, spec_colnames, obs_colnames))
     end
 
     n_man = Float64(size(obs_cov, 1))
@@ -107,25 +108,19 @@ function reorder_obs_cov(obs_cov, spec_colnames, obs_colnames)
     if spec_colnames == obs_colnames
         return obs_cov
     else
-        new_position = [findall(x .== obs_colnames)[1] for x in spec_colnames]
-        indices = reshape(
-            [CartesianIndex(i, j) for j in new_position for i in new_position],
-            size(obs_cov, 1),
-            size(obs_cov, 1),
-        )
-        obs_cov = obs_cov[indices]
+        new_position = [findfirst(==(x), obs_colnames) for x in spec_colnames]
+        obs_cov = obs_cov[new_position, new_position]
         return obs_cov
     end
 end
 
 # reorder means ----------------------------------------------------------------------------
-reorder_obs_mean(obs_mean::Nothing, spec_colnames, obs_colnames) = nothing
 
 function reorder_obs_mean(obs_mean, spec_colnames, obs_colnames)
     if spec_colnames == obs_colnames
         return obs_mean
     else
-        new_position = [findall(x .== obs_colnames)[1] for x in spec_colnames]
+        new_position = [findfirst(==(x), obs_colnames) for x in spec_colnames]
         obs_mean = obs_mean[new_position]
         return obs_mean
     end
diff --git a/src/observed/data.jl b/src/observed/data.jl
index 0dc4b07fb..77bbc8e8b 100644
--- a/src/observed/data.jl
+++ b/src/observed/data.jl
@@ -103,28 +103,14 @@ function SemObservedData(;
         data = Matrix(data)
     end
 
-    n_obs, n_man = Float64.(size(data))
-
-    if compute_covariance
-        obs_cov = Statistics.cov(data)
-    else
-        obs_cov = nothing
-    end
-
-    # if a meanstructure is needed, compute observed means
-    if meanstructure
-        obs_mean = vcat(Statistics.mean(data, dims = 1)...)
-    else
-        obs_mean = nothing
-    end
-
-    if rowwise
-        data_rowwise = [data[i, :] for i in 1:convert(Int64, n_obs)]
-    else
-        data_rowwise = nothing
-    end
-
-    return SemObservedData(data, obs_cov, obs_mean, n_man, n_obs, data_rowwise)
+    return SemObservedData(
+        data,
+        compute_covariance ? Statistics.cov(data) : nothing,
+        meanstructure ? vec(Statistics.mean(data, dims = 1)) : nothing,
+        Float64.(size(data, 2)),
+        Float64.(size(data, 1)),
+        rowwise ? [data[i, :] for i in axes(data, 1)] : nothing,
+    )
 end
 
 ############################################################################################
@@ -152,8 +138,8 @@ function reorder_data(data::AbstractArray, spec_colnames, obs_colnames)
     if spec_colnames == obs_colnames
         return data
     else
-        new_position = [findall(x .== obs_colnames)[1] for x in spec_colnames]
-        data = data[:, new_position]
-        return data
+        obs_positions = Dict(col => i for (i, col) in enumerate(obs_colnames))
+        new_positions = [obs_positions[col] for col in spec_colnames]
+        return data[:, new_positions]
     end
 end
diff --git a/src/observed/missing.jl b/src/observed/missing.jl
index cf1620ad9..6cfd09391 100644
--- a/src/observed/missing.jl
+++ b/src/observed/missing.jl
@@ -33,7 +33,7 @@ For observed data with missing values.
 - `get_data(::SemObservedMissing)` -> observed data
 - `data_rowwise(::SemObservedMissing)` -> observed data as vector per observation, with missing values deleted
 
-- `patterns(::SemObservedMissing)` -> indices of non-missing variables per missing patterns 
+- `patterns(::SemObservedMissing)` -> indices of non-missing variables per missing patterns
 - `patterns_not(::SemObservedMissing)` -> indices of missing variables per missing pattern
 - `rows(::SemObservedMissing)` -> row indices of observed data points that belong to each pattern
 - `pattern_n_obs(::SemObservedMissing)` -> number of data points per pattern
diff --git a/src/types.jl b/src/types.jl
index 641b34c60..803bc733a 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -45,7 +45,7 @@ function SemLoss(functions...; loss_weights = nothing, kwargs...)
     return SemLoss(functions, loss_weights)
 end
 
-# weights for loss functions or models. If the weight is nothing, multiplication returs second argument
+# weights for loss functions or models. If the weight is nothing, multiplication returns the second argument
 struct SemWeight{T}
     w::T
 end
@@ -147,7 +147,7 @@ Constructor for ensemble models.
 - `weights::Vector`:  Weights for each model. Defaults to the number of observed data points.
 
 All additional kwargs are passed down to the constructor for the optimizer field.
-    
+
 Returns a SemEnsemble with fields
 - `n::Int`: Number of models.
 - `sems::Tuple`: `AbstractSem`s.
@@ -170,7 +170,7 @@ function SemEnsemble(models...; optimizer = SemOptimizerOptim, weights = nothing
     # default weights
 
     if isnothing(weights)
-        nobs_total = sum(n_obs.(models))
+        nobs_total = sum(n_obs, models)
         weights = [n_obs(model) / nobs_total for model in models]
     end
 
diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl
index 01200f406..2f0fc749f 100644
--- a/test/examples/multigroup/build_models.jl
+++ b/test/examples/multigroup/build_models.jl
@@ -86,7 +86,7 @@ grad_fd = FiniteDiff.finite_difference_gradient(
     solution = sem_fit(model_ml_multigroup)
     update_estimate!(partable_s, solution)
     @test compare_estimates(
-        partable,
+        partable_s,
         solution_lav[:parameter_estimates_ml];
         atol = 1e-4,
         lav_groups = Dict(:Pasteur => 1, :Grant_White => 2),
@@ -134,16 +134,15 @@ struct UserSemML <: SemLossFunction end
 ### functors
 ############################################################################################
 
-import LinearAlgebra: isposdef, logdet, tr, inv
-import StructuralEquationModels: Σ, obs_cov, objective!
+using LinearAlgebra: isposdef, logdet, tr, inv
 
-function objective!(semml::UserSemML, parameters, model::AbstractSem)
-    let Σ = Σ(imply(model)), Σₒ = obs_cov(observed(model))
-        if !isposdef(Σ)
-            return Inf
-        else
-            return logdet(Σ) + tr(inv(Σ) * Σₒ)
-        end
+function SEM.objective!(semml::UserSemML, parameters, model::AbstractSem)
+    Σ = imply(model).Σ
+    Σₒ = SEM.obs_cov(observed(model))
+    if !isposdef(Σ)
+        return Inf
+    else
+        return logdet(Σ) + tr(inv(Σ) * Σₒ)
     end
 end
 
diff --git a/test/examples/multigroup/multigroup.jl b/test/examples/multigroup/multigroup.jl
index d5309d40c..759c24eda 100644
--- a/test/examples/multigroup/multigroup.jl
+++ b/test/examples/multigroup/multigroup.jl
@@ -1,12 +1,9 @@
 using StructuralEquationModels, Test, FiniteDiff
-import LinearAlgebra: diagind, LowerTriangular
-# import StructuralEquationModels as SEM
-include(
-    joinpath(
-        chop(dirname(pathof(StructuralEquationModels)), tail = 3),
-        "test/examples/helper.jl",
-    ),
-)
+using LinearAlgebra: diagind, LowerTriangular
+
+const SEM = StructuralEquationModels
+
+include(joinpath(chop(dirname(pathof(SEM)), tail = 3), "test/examples/helper.jl"))
 
 dat = example_data("holzinger_swineford")
 dat_missing = example_data("holzinger_swineford_missing")
diff --git a/test/examples/political_democracy/constructor.jl b/test/examples/political_democracy/constructor.jl
index 417e199cd..e13bc3816 100644
--- a/test/examples/political_democracy/constructor.jl
+++ b/test/examples/political_democracy/constructor.jl
@@ -1,4 +1,4 @@
-import Statistics: cov, mean
+using Statistics: cov, mean
 
 ############################################################################################
 ### models w.o. meanstructure
diff --git a/test/examples/political_democracy/political_democracy.jl b/test/examples/political_democracy/political_democracy.jl
index 093b07588..389800745 100644
--- a/test/examples/political_democracy/political_democracy.jl
+++ b/test/examples/political_democracy/political_democracy.jl
@@ -1,5 +1,5 @@
 using StructuralEquationModels, Test, FiniteDiff
-# import StructuralEquationModels as SEM
+
 include(
     joinpath(
         chop(dirname(pathof(StructuralEquationModels)), tail = 3),
diff --git a/test/examples/recover_parameters/recover_parameters_twofact.jl b/test/examples/recover_parameters/recover_parameters_twofact.jl
index f7e555813..dc47a9055 100644
--- a/test/examples/recover_parameters/recover_parameters_twofact.jl
+++ b/test/examples/recover_parameters/recover_parameters_twofact.jl
@@ -1,5 +1,5 @@
 using StructuralEquationModels, Distributions, Random, Optim, LineSearches
-import StructuralEquationModels as SEM
+
 include(
     joinpath(
         chop(dirname(pathof(StructuralEquationModels)), tail = 3),
@@ -63,10 +63,10 @@ Random.seed!(1234)
 x = transpose(rand(true_dist, 100000))
 semobserved = SemObservedData(data = x, specification = nothing)
 
-loss_ml = SemLoss(SEM.SemML(; observed = semobserved, n_par = length(start)))
+loss_ml = SemLoss(SemML(; observed = semobserved, n_par = length(start)))
 
 optimizer = SemOptimizerOptim(
-    BFGS(; linesearch = BackTracking(order = 3), alphaguess = InitialHagerZhang()),# m = 100), 
+    BFGS(; linesearch = BackTracking(order = 3), alphaguess = InitialHagerZhang()),# m = 100),
     Optim.Options(; f_tol = 1e-10, x_tol = 1.5e-8),
 )
 
diff --git a/test/runtests.jl b/test/runtests.jl
index a21ff9b87..c3b15475f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,7 +1,9 @@
 using Test, SafeTestsets, JuliaFormatter, StructuralEquationModels
 
 @testset "JuliaFormatter.jl" begin
-    @test format(StructuralEquationModels; verbose = false, overwrite = false)
+    if !format(StructuralEquationModels; verbose = false, overwrite = false)
+        @warn "Julia code formatting style inconsistencies detected."
+    end
 end
 @time @safetestset "Unit Tests" begin
     include("unit_tests/unit_tests.jl")
diff --git a/test/unit_tests/data_input_formats.jl b/test/unit_tests/data_input_formats.jl
index 540e43c25..bef7ba8df 100644
--- a/test/unit_tests/data_input_formats.jl
+++ b/test/unit_tests/data_input_formats.jl
@@ -1,6 +1,5 @@
 using StructuralEquationModels, Test, Statistics
-import StructuralEquationModels: obs_cov, obs_mean, get_data
-
+using StructuralEquationModels: obs_cov, obs_mean, get_data
 ### model specification --------------------------------------------------------------------
 
 spec = ParameterTable(nothing)