diff --git a/docs/src/HowTo/custom_solvers.md b/docs/src/HowTo/custom_solvers.md
index ec675db8..0be9a454 100644
--- a/docs/src/HowTo/custom_solvers.md
+++ b/docs/src/HowTo/custom_solvers.md
@@ -1,48 +1,46 @@
-# [Custom Solvers](@id custom_solvers)
+# [Alternative Solvers](@id custom_solvers)
 
 ### Setup some data
 
 ```@Example main
 using Unfold
 using UnfoldMakie, CairoMakie
-using UnfoldSim,
-dat,evts = UnfoldSim.predef_eeg(;noiselevel=10,return_epoched=true)
-
-f  = @formula 0~1+condition+continuous
-designDict = Dict(Any=>(f,range(0,1,length=size(dat,1))))
+using UnfoldSim
+dat, evts = UnfoldSim.predef_eeg(; noiselevel = 10, return_epoched = true)
 
+f = @formula 0 ~ 1 + condition + continuous
+designDict = Dict(Any => (f, range(0, 1, length = size(dat, 1))))
 ```
 
 ### GPU Solvers
-GPU solvers can speed up your modelfit drastically! up to factor of 30 has been observed already
+GPU solvers can significantly speed up your model fitting, with observed improvements of up to a factor of 30!
+
 ```julia
-using Krylov,CUDA # necessary to load the right package extension
-gpu_solver =(x,y)->Unfold.solver_krylov(x,y;GPU=true)
-m = Unfold.fit(UnfoldModel,designDict,evts,dat,solver=gpu_solver)
+using Krylov, CUDA # necessary to load the right package extension
+gpu_solver =(x, y) -> Unfold.solver_krylov(x, y; GPU = true)
+m = Unfold.fit(UnfoldModel, designDict, evts, dat, solver = gpu_solver)
 ```
-We can't run it on the docs though, so try it yourself! If you need something else than CUDA, write an issue, we cant test it with something else right now...
-
+To test it, you will need to run it yourself as we cannot run it on the docs. If you require a different graphicscard vendor than NVIDA/CUDA, please create an issue. Currently, we are unable to test it due to lack of hardware.
 
 ### Robust Solvers
-Robust solvers automatically account for outlying trials. They come at a severe computational cost though!
+Robust solvers automatically account for outlier trials, but they come at a significant computational cost.
 ```@Example main
-using RobustModels # necessary to load the Unfold st
-package extension
-se_solver =(x,y)->Unfold.solver_robust(x,y)
-m = Unfold.fit(UnfoldModel,designDict,evts,dat,solver=se_solver)
-results =coeftable(m)
-plot_erp(results;stderror=true)
+using RobustModels # necessary to load the Unfold package extension
+se_solver = (x, y) -> Unfold.solver_robust(x, y)
+m = Unfold.fit(UnfoldModel, designDict, evts, dat, solver = se_solver)
+results = coeftable(m)
+plot_erp(results; stderror = true)
 ```
 
 ### Back2Back regression
 ```@Example main
-b2b_solver = (x, y) -> Unfold.solver_b2b(x, y;ross_val_reps = 5)
-dat_3d = permutedims(repeat(dat,1,1,20),[3 1 2])
-m = Unfold.fit(UnfoldModel, designDict, evts,dat_3d; solver=b2b_solver)
+b2b_solver = (x, y) -> Unfold.solver_b2b(x, y; ross_val_reps = 5)
+dat_3d = permutedims(repeat(dat, 1, 1, 20), [3 1 2])
+m = Unfold.fit(UnfoldModel, designDict, evts, dat_3d; solver = b2b_solver)
 results = coeftable(m)
 
 plot_erp(results)
 ```
-These are the decoding-results for `conditionA` while taking into account `conditionB` - and vice versa. 
+These are the decoding results for `conditionA` while considering `conditionB`, and vice versa. 
 
 
diff --git a/docs/src/HowTo/multiple_events.md b/docs/src/HowTo/multiple_events.md
index 13258952..cad8e5a8 100644
--- a/docs/src/HowTo/multiple_events.md
+++ b/docs/src/HowTo/multiple_events.md
@@ -1,7 +1,6 @@
 # How to model multiple events
 
-In case of overlapping data, we often need to plot multiple events.
-
+When dealing with overlapping data, it is often necessary to plot multiple events.
 
 ### Load Example Data
 ```@example main
@@ -11,34 +10,42 @@ using DataFrames
 using StatsModels
 using MixedModels
 
-include(joinpath(dirname(pathof(Unfold)), "../test/test_utilities.jl") ) # to load data
+include(joinpath(dirname(pathof(Unfold)), "../test/test_utilities.jl")) # to load data
 dat, evts = loadtestdata("test_case_4b");
 
 evts[1:5,:]
 ```
-As you can see, there are two events here. EventA and EventB. Both are in the column `type` in which the toolbox looks for different events by default.
+The `type` column of table `evts` contains two conditions: EventA and EventB. By default, the toolbox will search for these conditions.
 
 ### Specify formulas and basisfunctions
 
 ```@example main
-bf1 = firbasis(τ=(-0.4,.8),sfreq=50,name="stimulusA")
-bf2 = firbasis(τ=(-0.2,1.2),sfreq=50,name="stimulusB")
+bf1 = firbasis(τ = (-0.4, 0.8), sfreq = 50, name = "stimulusA")
+bf2 = firbasis(τ = (-0.2, 1.2), sfreq = 50, name = "stimulusB")
 ```
-For each event, we have to specify a basisfunction and a formula. We could use the same basis and the same formulas though
+For each event, a basis function and formula must be specified. The same basis and formulas may be used.
 ```@example main
-f  = @formula 0~1
+f  = @formula 0 ~ 1
 ```
 
-Next, we have to specify for each event, what is the formula and what is the basisfunction we want to use
+For each event, we must specify the formula and basis function to be used. 
 ```@example main
-bfDict = ["eventA"=>(f,bf1),
-              "eventB"=>(f,bf2)]
+
+bfDict = [ "eventA" => (f, bf1),
+           "eventB" => (f, bf2) ]
 ```
 
 Finally, fitting & plotting works the same way as always
-```@example main
 
-m = Unfold.fit(UnfoldModel,bfDict,evts,dat;solver=(x,y) -> Unfold.solver_default(x,y;stderror=true),eventcolumn="type")
+```@example main
+m = Unfold.fit(
+    UnfoldModel,
+    bfDict,
+    evts,
+    dat,
+    solver = (x, y) -> Unfold.solver_default(x, y; stderror = true),
+    eventcolumn = "type",
+)
 results = coeftable(m)
-plot_erp(results;stderror=true,mapping=(;col=:group))
+plot_erp(results; stderror = true, mapping = (; col = :group))
 ``` 
diff --git a/docs/src/HowTo/standarderrors.md b/docs/src/HowTo/standarderrors.md
index 3721ea47..411b0ba0 100644
--- a/docs/src/HowTo/standarderrors.md
+++ b/docs/src/HowTo/standarderrors.md
@@ -5,22 +5,21 @@
 ```@Example main
 using Unfold
 using UnfoldMakie, CairoMakie
-using UnfoldSim,
-dat,evts = UnfoldSim.predef_eeg(;noiselevel=10,return_epoched=true)
-
-f  = @formula 0~1+condition+continuous
-designDict = Dict(Any=>(f,range(0,1,length=size(dat,1))))
+using UnfoldSim
+dat, evts = UnfoldSim.predef_eeg(; noiselevel = 10, return_epoched = true)
 
+f = @formula 0 ~ 1 + condition + continuous
+designDict = Dict(Any => (f, range(0, 1, length = size(dat, 1))))
 ```
 
-Similar to @Ref(custom_solver), it is possible to specify a solver than also calculates the standard errors of the (single-subject) estimates.
+It is possible to specify a solver that calculates the standard errors of the estimates for a single subject as it possible for [custom solvers](@ref custom_solvers).
 
 ```@Example main
-se_solver =(x,y)->Unfold.solver_default(x,y,stderror=true)
-m = Unfold.fit(UnfoldModel,designDict,evts,dat,solver=se_solver)
-results =coeftable(m)
-plot_erp(results;stderror=true)
+se_solver = (x, y) -> Unfold.solver_default(x, y, stderror = true)
+m = Unfold.fit(UnfoldModel, designDict, evts, dat, solver = se_solver)
+results = coeftable(m)
+plot_erp(results; stderror = true)
 ```
 !!! warning
-    **In case of overlap-correction:** Use single-subject SE on your own risk. Because EEG data are autocorrelated your standarderrors will be typically too small.
+    **In case of overlap-correction:** Use single-subject standard errors on your own risk. EEG data is autocorrelated, which means that standard errors are typically too small.
 
diff --git a/docs/src/tutorials/lm_mu.md b/docs/src/tutorials/lm_mu.md
index e406d671..95b5eabc 100644
--- a/docs/src/tutorials/lm_mu.md
+++ b/docs/src/tutorials/lm_mu.md
@@ -1,6 +1,6 @@
 # [Mass Univariate Linear Models (no overlap correction)](@id lm_massunivariate)
 
-In this notebook we will fit regression models to (simulated) EEG data. We will see that we need some type of overlap correction, as the events are close in time to each other, so that the respective brain responses overlap.
+In this notebook we will fit regression models to simulated EEG data. We will see that we need some type of overlap correction, as the events are close in time to each other, so that the respective brain responses overlap.
 If you want more detailed introduction to this topic check out [our paper](https://peerj.com/articles/7838/).
 
 
@@ -9,7 +9,7 @@ If you want more detailed introduction to this topic check out [our paper](https
 ```@example Main
 using DataFrames
 using Unfold
-using UnfoldMakie,CairoMakie # for plotting
+using UnfoldMakie, CairoMakie # for plotting
 using UnfoldSim
 
 
@@ -23,91 +23,93 @@ data, evts = UnfoldSim.predef_eeg()
 nothing # hide
 ```
 ## Inspection
-The data has little noise and the underlying signal is a pos-neg-pos spike pattern
+The data has only little noise. The underlying signal pattern is a positive-negative-positive spike.
 ```@example Main
-times = range(1/50,length=200,step=1/50)
-plot(times,data[1:200])
-vlines!(current_axis(),evts[evts.latency.<=200,:latency]./50) # show events, latency in samples!
+times = range(1/50, length=200, step=1/50)
+f = Figure()
+plot(f[1, 1], data[1:200], times)
+vlines!(evts[evts.latency .<= 200, :latency] ./ 50) # show events, latency in samples!
+f
 ```
 
 To inspect the event dataframe we use
 ```@example Main
-show(first(evts,6,),allcols=true)
+show(first(evts, 6), allcols = true)
 ```
-Every row is an event. Note that `:latency` is commonly the timestamp in samples, whereas `:onset` would typically refer to seconds.
+Every row is an experimental event. Note that `:latency` refers to time in samples, whereas `:onset` would typically refer to seconds.
 
 
 ## Traditional Mass Univariate Analysis
-For a mass univariate analysis, will 
-1. epoch the data
-2. specify a formula 
-3. fit a linear model to each time point & channel
-4. visualize the results.
+To perform a mass univariate analysis, you must complete the following steps:
 
+1. Split data into epochs 
+2. Specify a formula 
+3. Fit a linear model to each time point & channel
+4. Visualize the results.
 
-#### 1. Epoch the data
-We have to cut the data into small regular trials.
+
+#### 1. Split data into epochs 
+Initially, you have data with a duration that represents the whole experimental trial. You need to cut the data into small regular epochs related to the some event, e.g. start of fixation.
 
 ```@example Main
 # we have multi channel support
-data_r = reshape(data,(1,:))
+data_r = reshape(data, (1,:))
 # cut the data into epochs
-data_epochs,times = Unfold.epoch(data=data_r,tbl=evts,τ=(-0.4,0.8),sfreq=50);
+data_epochs, times = Unfold.epoch(data = data_r, tbl = evts, τ = (-0.4, 0.8), sfreq = 50);
 size(data_epochs)
 ```
-
-`τ` specifies the epoch size. To convert `τ` to samples, we need to specify the sampling rate.
+- `τ` specifies the epoch size.
+- `sfreq` - sampling rate, converts `τ` to samples.
 
 
 ```@example Main
 typeof(data_epochs)
 ```
 !!! note
-    In julia, `missing` is supported throughout the ecosystem. Thus, we can have partial trials and they will be incorporated / ignored at the respective functions. Helpful functions are the julia-base `disallowmissing` and the internal `Unfold.dropMissingEpochs` functions
-
-
+   Partial trials could be included or ignored by the corresponding functions due to `missing` data type. Check functions such as the Julia-based `disallowmissing` and the internal `Unfold.drop_missing_epochs`.
 
 #### 2. Specify a formula
-We define a formula that we want to apply to each point in time
+Define a formula to be applied to each time point.
+
 ```@example Main
-f  = @formula 0~1+condition+continuous # 0 as a dummy, we will combine wit data later
+f = @formula 0 ~ 1 + condition + continuous # 0 as a dummy, we will combine with data later
 nothing # hide
 ```
 
-#### 3. fit a linear model to each time-point & channel
+#### 3. Fit a linear model to each time point & channel
 
-We fit the `UnfoldLinearModel`.
+Fit the `UnfoldLinearModel`.
 ```@example Main
-m = fit(UnfoldModel,f,evts,data_epochs,times); 
+m = fit(UnfoldModel, f, evts, data_epochs, times); 
 nothing #hide
 ```
 
-Alternatively, we could also call it using:
+Alternative way to call this model is below. This syntax allows you to fit multiple events at once. For example, replacing `Any` with `:fixation =>...` will fit this model specifically to the fixation event type.
 ```@example Main
-m = fit(UnfoldModel,Dict(Any=>(f,times)),evts,data_epochs); 
+m = fit(UnfoldModel, Dict(Any=>(f, times)), evts, data_epochs); 
 nothing #hide
 ```
-A syntax which allows for fitting multiple events at the same time (by replacing Any with e.g. `:fixation =>...` it will fit that model specifically to the fixation event type)
 
-We can inspect the object
+Inspect the fitted model:
 ```@example Main
 m
 ```
-And see that there are several helpful functions to recover `design`, `designmatrix`, raw-`modelfit` and most importantly, `coeftable`. 
+Note these functions to discover the model: `design`, `designmatrix`, `modelfit` and most importantly, `coeftable`. 
 
 !!! info
-        There are of course further methods, e.g. `coef`, `ranef`, `Unfold.formula`, `modelmatrix` which might be helpful at some point, but not important now
+        There are of course further methods, e.g. `coef`, `ranef`, `Unfold.formula`, `modelmatrix` which might be helpful at some point, but not important now.
 
-Using `coeftable`, we can get a *tidy*-dataframe, very useful for your further analysis.
+Using `coeftable`, we can get a *tidy* DataFrames, very useful for your further analysis.
 
 ```@example Main
-first(coeftable(m),6)
+first(coeftable(m), 6)
 ```
 
 #### 4. Visualize the results
-Tidy-Dataframes make them easy to visualize using e.g. AlgebraOfGraphics.jl. We simpliefied this even more, by implementing a wrapper in `UnfoldMakie`
+Tidy DataFrames are easy to visualize using e.g. AlgebraOfGraphics.jl. Function `plot_erp` from `UnfoldMakie`makes it even easier.  
+
 ```@example Main
 results = coeftable(m)
 plot_erp(results)
 ```
-As you can see here, a lot is going on, even in the baseline-period! This is because the signal was simulated with overlapping events. Head over to the next tutorial to find out how to remedy this.
+As you can see, there is a lot going on, even in the baseline period! This is because the signal was simulated with overlapping events. In the next tutorial you will learn how to fix this.
diff --git a/docs/src/tutorials/lmm_mu.md b/docs/src/tutorials/lmm_mu.md
index c7dfcc73..1ecabcc8 100644
--- a/docs/src/tutorials/lmm_mu.md
+++ b/docs/src/tutorials/lmm_mu.md
@@ -5,7 +5,7 @@
 using Unfold
 using UnfoldSim
 using MixedModels # important to load to activate the UnfoldMixedModelsExtension
-using UnfoldMakie,CairoMakie # plotting
+using UnfoldMakie, CairoMakie # plotting
 using DataFrames
 using CategoricalArrays
 nothing;#hide
@@ -13,72 +13,61 @@ nothing;#hide
 !!! important
     You have to run `using MixedModels` before or after loading Unfold to activate the MixedModels abilities!
 
-This notebook is similar to the [Mass Univariate Linear Models (no overlap correction) tutorial](@ref lm_massunivariate) , but fits mass-univariate *mixed* models - that is, one model over all subjects, instead one model per subject. This allows to incorporate e.g. Item-effects.
-
-
+This notebook is similar to the [Mass Univariate Linear Models (no overlap correction) tutorial](@ref lm_massunivariate), but fits mass-univariate *mixed* models - that is, one model over all subjects, instead of one model per subject. This allows to include item effects, for example.
 
 ## Mass Univariate **Mixed** Models
 Again we have 4 steps:
-1. epoch the data
-2. specify a formula 
-3. fit a linear model to each time point & channel
-4. visualize the results.
-
+1. Split data into epochs 
+2. Specify a formula 
+3. Fit a linear model to each time point & channel
+4. Visualize the results.
 
 #### 1. Epoching
 
 ```@example Main
-
-data, evts = UnfoldSim.predef_eeg(10;return_epoched=true) # simulate 10 subjects
-data = reshape(data,1,size(data, 1), :) # we need to concatenate the data to one long EEG dataset
-times = range(0,length=size(data,2),step=1/100)
-transform!(evts,:subject=>categorical=>:subject); # has to be categorical, else MixedModels.jl complains
+data, evts = UnfoldSim.predef_eeg(10; return_epoched = true) # simulate 10 subjects
+data = reshape(data, 1, size(data, 1), :) # concatenate the data into a long EEG dataset
+times = range(0, length = size(data, 2), step = 1 / 100)
+transform!(evts, :subject => categorical => :subject); # :subject must be categorical, otherwise MixedModels.jl complains
 nothing #hide
 ```
 
 The `events` dataFrame has an additional column (besides being much taller): `subject`
 ```@example Main
-first(evts,6)
-```        
-
-#### 2. Specify the formula
-We define the formula. Importantly we need to specify a random effect. We are using `zerocorr` to speed up the calculation and show off that we can use it.
+first(evts, 6)
+```
 
+#### 2. Formula specification
+We define the formula. Importantly, we need to specify a random effect. We use `zerocorr` to speed up the calculation.
 ```@example Main
-f  = @formula 0~1+condition*continuous+zerocorr(1+condition*continuous|subject);
+f = @formula 0 ~ 1 + condition * continuous + zerocorr(1 + condition * continuous | subject);
 nothing #hide
 ```
 
-
-#### 3. Fit the model
-We can now run the LinearMixedModel on each time point
+#### 3. Model fitting
+We can now run the LinearMixedModel at each time point.
 ```@example Main
-
-m = fit(UnfoldModel,f,evts,data,times)
+m = fit(UnfoldModel, f, evts, data, times)
 nothing #hide
 ```
 
+#### 4. Visualization of results
 
-#### 4. Visualize results
+Let's start with the **fixed** effects. 
+We see the condition effects and some residual overlap activity in the fixed effects.
 
-Let's start with the **fixed** Effects
 ```@example Main
 results = coeftable(m)
 
-res_fixef = results[isnothing.(results.group),:]
+res_fixef = results[isnothing.(results.group), :]
 plot_erp(res_fixef)
 ```
 
-
-We see the condition effects and some residual overlap activity in the fixed effects
-
-
-And now the **random** effect results
+And now comes the **random** effect:
 ```@example Main
-res_ranef = results[results.group.==:subject,:]
+res_ranef = results[results.group .== :subject, :]
 plot_erp(res_ranef)
 ```
 
-
 ### Statistics
-Check out the [LMM p-value tutorial](@ref lmm_pvalues)
\ No newline at end of file
+Check out the [LMM p-value tutorial](@ref lmm_pvalues)
diff --git a/docs/src/tutorials/lmm_overlap.md b/docs/src/tutorials/lmm_overlap.md
index bb3e830e..3ab8c154 100644
--- a/docs/src/tutorials/lmm_overlap.md
+++ b/docs/src/tutorials/lmm_overlap.md
@@ -7,7 +7,7 @@ using UnfoldSim
 
 using CategoricalArrays
 using MixedModels
-using UnfoldMakie,CairoMakie
+using UnfoldMakie, CairoMakie
 using DataFrames
 
 nothing;#hide
@@ -17,67 +17,60 @@ nothing;#hide
 This notebook is similar to the Linear Model with Overlap Correction tutorial, but fits **mixed** models with overlap correction
 
 !!! warning 
-    **Limitation**: This is not production ready at all. Still lot's of things to find out and tinker with. Don't use this if you did not look under the hood of the toolbox!
+    **Limitation**: This is not ready for use at all. There are still a lot of things to find out and tinker with. Don't use this if you haven't looked under the hood of the toolbox!
 
-!!! important
-    Even worse - it is right now not working. Some new bug appeared
 
-## Get some  data
+## Get some data
 
 ```@example Main
-dat,evts = UnfoldSim.predef_2x2(;signalsize=20,n_items=16,n_subjects=16)
+dat, evts = UnfoldSim.predef_2x2(; signalsize=20, n_items=16, n_subjects=16)
 
-# we have to fix the latencies as well, they are now relative to 1:size(data,1), but we want one continuous long EEG
-subj_idx = [parse(Int,split(string(s),'S')[2]) for s in evts.subject]
-evts.latency .+= size(dat,1) .* (subj_idx.-1)
+# We also need to fix the latencies, they are now relative to 1:size(data, 1), but we want a continuous long EEG.
+subj_idx = [parse(Int, split(string(s), 'S')[2]) for s in evts.subject]
+evts.latency .+= size(dat, 1) .* (subj_idx .- 1)
 
 dat = dat[:] # we need all data concatenated over subjects
 evts.subject  = categorical(Array(evts.subject))
 nothing #hide
 ```
 
-
 ## Linear **Mixed** Model Continuous Time
 Again we have 4 steps:
-1. specify a temporal basisfunction
-2. specify a formula
-3. fit a linear model for each channel (one for all timepoints!)
-4. visualize the results.
+1. Specify a temporal basisfunction
+2. Specify a formula
+3. Fit a linear model for each channel (one model for all timepoints!)
+4. Visualize the results.
 
-#### 1. specify a temporal basisfunction
-By default, we would want to use a FIR basisfunction. See [Basis Functions](@ref) for more details.
+#### 1. Specify a temporal basisfunction
+By default, we would want to use a FIR basis function. See [Basis Functions](@ref) for more details.
 ```@example Main
-basisfunction = firbasis(τ=(-0.4,.8),sfreq=20,name="stimulus")
+basisfunction = firbasis(τ=(-0.4, .8), sfreq=20, name="stimulus")
 nothing #hide
 ```
 
-
-
-
 #### 2. Specify the formula
-We define the formula. Importantly we need to specify a random effect. 
+Define the formula and specify a random effect. 
 
 !!! note
-    We are using `zerocorr` because we need it here, else the model will try to model all correlations between all timepoints and all factors!
+    We use `zerocorr` to prevent the model from computing all correlations between all timepoints and factors.
 
 ```@example Main
-f  = @formula 0~1+A*B+zerocorr(1+A*B|subject);
+f  = @formula 0 ~ 1 + A  *B + zerocorr(1 + A*B|subject);
 ```
 
 
 #### 3. Fit the model
 ```@example Main
-bfDict = Dict(Any=>(f,basisfunction))
-# for some reason this results in a big error. Skipping this tutorial right now
-m = fit(UnfoldModel,bfDict,evts,dat) 
+bfDict = Dict(Any=>(f, basisfunction))
+# Skipping this tutorial for now due to a significant error.
+m = fit(UnfoldModel, bfDict, evts, dat) 
 
 results = coeftable(m)
-first(results,6)
+first(results, 6)
 ```
 
-
 #### 4. Visualize results
 
 ```@example Main
-plot_erp(results;mapping=(;col=:group))
+plot_erp(results; mapping=(; col = :group))
 ```