Variance analysis of mixed effects based on the example of data on the specific range of cars

Introduction

Imagine that you are a test engineer, and you need to find out whether different car models really differ in specific power reserve, or whether the observed differences are just an accident related, for example, to the specifics of a particular factory. The answer to this question is provided by analysis of variance, a statistical method that decomposes the overall variability of data into parts caused by different sources. The key idea is as follows: if the differences between the groups significantly exceed the level of random "noise", the influence of the factor is recognized as significant.

This example uses a mixed linear model. Each car model acts as a fixed effect, and the manufacturer and the interaction of the factory with the model act as random effects.

The result of training the model is the variance components, which quantify the contribution of each source to the overall spread.:

(Manufacturer) — how different are the cars from different factories on average;
(Manufacturer × Model) — the degree of inconsistency of the effect of the model from plant to plant (interaction of factors);
(Remainder) is an unexplained variation related to the individual characteristics of cars and measurement errors.

To visually represent the data, we use boxplots (boxplot, or "box with a mustache"). This graph compactly displays the median, quartiles, range, and potential outliers, allowing you to instantly estimate the center of the distribution, its spread, and symmetry in each group.

Formal hypothesis testing relies on the likelihood ratio test: it compares the complete model with a reduced (no fixed effect model) and determines whether the exclusion of the car model leads to a statistically significant deterioration in the description of the data.

By combining the capabilities of variance analysis of mixed effects with expressive statistical graphics, we reconstruct the objective picture hidden behind the columns of numbers in the test reports.

Importing libraries

We will add the libraries necessary for working with variance analysis.

neededLibs = ["AnovaGLM", "CSV", "CategoricalArrays", "MixedModels", "Statistics", "Distributions", "GLM", "StatsPlots"]
for lib in neededLibs
    try
        eval(Meta.parse("using $lib"))
    catch ex
        Pkg.add(lib)
        eval(Meta.parse("using $lib"))
    end
end

Initial data

Let's determine the initial data on the specific power reserve of each car model (пробег). We transform the data into a table in which the data is aligned.

# Initial data
исходные = [
    33.3  34.5  37.4
    33.4  34.8  36.8
    32.9  33.8  37.6
    32.6  33.4  36.6
    32.5  33.7  37.0
    33.0  33.9  36.7
]

# Parameters of the experiment plan
заводов = 3
строк = 6
моделей = 2
повторов = div(строк, моделей)  
пробег = Float64[]
завод = Int[]
модель = Int[]

for з in 1:заводов
    for стр in 1:строк
        push!(пробег, исходные[стр, з])
        push!(завод, з)
        push!(модель, стр <= 3 ? 1 : 2)
    end
end

# Data table
таблица = DataFrame(
    Пробег = пробег,
    Завод = string.(завод),
    Модель = string.(модель)
)
таблица.Завод = categorical(таблица.Завод)
таблица.Модель = categorical(таблица.Модель)

println("Specific power reserve of cars")
println(first(таблица, 18))

Удельный запас хода автомобилей
18×3 DataFrame
 Row │ Пробег   Завод  Модель
     │ Float64  Cat…   Cat…
─────┼────────────────────────
   1 │    33.3  1      1
   2 │    33.4  1      1
   3 │    32.9  1      1
   4 │    32.6  1      2
   5 │    32.5  1      2
   6 │    33.0  1      2
   7 │    34.5  2      1
   8 │    34.8  2      1
   9 │    33.8  2      1
  10 │    33.4  2      2
  11 │    33.7  2      2
  12 │    33.9  2      2
  13 │    37.4  3      1
  14 │    36.8  3      1
  15 │    37.6  3      1
  16 │    36.6  3      2
  17 │    37.0  3      2
  18 │    36.7  3      2

Building models

Let's build a mixed linear model in which mileage is explained by the fixed effect of the car model and the random effects of the factory and the interaction of the factory with the model. Let's select the model parameters using the maximum likelihood method and analyze the obtained estimates of the variance components and fixed effects.

формула1 = @formula(Пробег ~ 1 + Модель + (1 | Завод) + (1 | Завод & Модель))
подбор1 = fit(MixedModel, формула1, таблица)
println(подбор1)

Minimizing 2    Time: 0:00:00 ( 0.13  s/it)
   objective: 42.098336394833446

┌ Warning: ProgressMeter by default refresh meters with additional information in IJulia via `IJulia.clear_output`, which clears all outputs in the cell. 
│  - To prevent this behaviour, do `ProgressMeter.ijulia_behavior(:append)`. 
│  - To disable this warning message, do `ProgressMeter.ijulia_behavior(:clear)`.
└ @ ProgressMeter /usr/local/ijulia-core/packages/ProgressMeter/N660J/src/ProgressMeter.jl:607

Minimizing 34    Time: 0:00:00 (19.87 ms/it)

Linear mixed model fit by maximum likelihood
 Пробег ~ 1 + Модель + (1 | Завод) + (1 | Завод & Модель)
   logLik   -2 logLik     AIC       AICc        BIC    
   -12.1071    24.2142    34.2142    39.2142    38.6661

Variance components:
                  Column   Variance Std.Dev. 
Завод & Модель (Intercept)  0.000000 0.000000
Завод          (Intercept)  2.948321 1.717068
Residual                    0.093778 0.306232
 Number of obs: 18; levels of grouping factors: 6, 3

  Fixed-effects parameters:
───────────────────────────────────────────────────
                 Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────
(Intercept)  34.9444      0.996591  35.06    <1e-99
Модель: 2    -0.566667    0.144359  -3.93    <1e-04
───────────────────────────────────────────────────

Let's build a reduced mixed model without a fixed effect of the car model, leaving only the random effects of the factory and interaction. Let's perform a likelihood ratio test between the full and reduced models in order to statistically verify the significance of the influence of the car model on mileage.

формула0 = @formula(Пробег ~ 1 + (1 | Завод) + (1 | Завод & Модель))
подбор0 = fit(MixedModel, формула0, таблица)
тест_правд = MixedModels.likelihoodratiotest(подбор0, подбор1)
println("Likelihood ratio test results:")
println(тест_правд)

Результаты теста отношения правдоподобия:
Model Formulae
1: Пробег ~ 1 + (1 | Завод) + (1 | Завод & Модель)
2: Пробег ~ 1 + Модель + (1 | Завод) + (1 | Завод & Модель)
─────────────────────────────────────────────────
     model-dof  -2 logLik      χ²  χ²-dof  P(>χ²)
─────────────────────────────────────────────────
[1]          4    31.5367                        
[2]          5    24.2142  7.3224       1  0.0068
─────────────────────────────────────────────────

Estimates of variance components

Let's extract the variance components from the trained model using the function VarCorr Let's store the numerical variance values for each random effect in a convenient data structure.

компоненты = VarCorr(подбор1)
println(компоненты)
# A function for extracting variances
function извлечь_дисперсии(vc)
    строка_vc = string(vc)
    дисперсии = Dict{String, Float64}()
    for line in split(строка_vc, "\n")
        очищенная = strip(line)
        if isempty(очищенная)
            continue
        end
        if occursin("Variance components:", очищенная) || 
           occursin("Column", очищенная) || 
           (occursin("Variance", очищенная) && occursin("Std.Dev.", очищенная)) ||
           startswith(очищенная, "──")
            continue
        end
        
        части = split(очищенная, r"\s+")
        if length(части) >= 3
            предпоследнее = части[end-1]
            if occursin(r"^\d+\.?\d*$", предпоследнее) || occursin(r"^\d*\.?\d+$", предпоследнее)
                зн = tryparse(Float64, предпоследнее)
                if зн !== nothing
                    имя = join(части[1:end-2], " ")
                    дисперсии[имя] = зн
                end
            end
        end
    end
    return дисперсии
end

дисп = извлечь_дисперсии(компоненты)
σ²_завод = 0.0
σ²_взаимодействие = 0.0
σ²_остаток = 0.0

for (k, v) in дисп
    if occursin("Производитель & Модель", k)
        σ²_взаимодействие = v
    elseif occursin("Производитель", k) && !occursin("Модель", k)
        σ²_завод = v
    elseif occursin("Residual", k)
        σ²_остаток = v
    end
end

println("Estimates of random effect variances:")
println("  σ2(Manufacturer)                = ", round(σ²_завод, digits=4))
println("  σ2(Manufacturer:The model)         = ", round(σ²_взаимодействие, digits=4))
println("  σ2(Error)               = ", round(σ²_остаток, digits=4))

Variance components:
                  Column   Variance Std.Dev. 
Завод & Модель (Intercept)  0.000000 0.000000
Завод          (Intercept)  2.948321 1.717068
Residual                    0.093778 0.306232

Оценки дисперсий случайных эффектов:
  σ²(Производитель)                = 0.0
  σ²(Производитель:Модель)         = 0.0
  σ²(Ошибка)               = 0.0938

Visualization of results

Let's set up the visual design of the graphs.

theme(:default)
default(fontfamily="sans-serif", titlefontsize=11, guidefontsize=9, tickfontsize=7)
цвета_заводов = [:steelblue, :darkorange, :forestgreen]
цвета_комп = [:steelblue, :darkorange, :forestgreen]

3-element Vector{Symbol}:
 :steelblue
 :darkorange
 :forestgreen

Using boxplot, we will plot the distribution of the specific power reserve of cars for each manufacturer.

gr()
гр1 = @df таблица boxplot(:Завод, :Пробег,
    legend=false,
    title="Distribution of specific power reserve",
    xlabel="Manufacturer",
    ylabel="Specific power reserve",
    color=цвета_заводов,
    linewidth=2,
    size=(500, 400))
display(гр1)

Let's plot the distribution of the specific power reserve of cars for each model.

гр2 = @df таблица boxplot(:Модель, :Пробег,
    legend=false,
    title="Distribution of specific power reserve",
    xlabel="Car model",
    ylabel="Specific power reserve",
    color=[:salmon, :lightblue],
    linewidth=2,
    size=(500, 400))
display(гр2)

Let's build a graph of Manufacturer-Model interaction.

средние_зм = combine(groupby(таблица, [:Завод, :Модель]), :Пробег => mean => :Средний)
средние_зм.Завод = string.(средние_зм.Завод)
средние_зм.Модель = string.(средние_зм.Модель)
гр3 = plot(legend=:bottomright, size=(500, 400),
    title="Manufacturer-Model Interaction Schedule",
    xlabel="Manufacturer",
    ylabel="Specific power reserve")
цвета_моделей = [:blue, :red]
for (i, мод) in enumerate(["1", "2"])
    данные = средние_зм[средние_зм.Модель .== мод, :]
    plot!(гр3, данные.Завод, данные.Средний,
        marker=:circle, markersize=8, linewidth=2,
        color=цвета_моделей[i],
        label="Модель $мод")
end
display(гр3)

Let's group the data by manufacturer and calculate the specific power reserve, standard deviation, and number of observations for each one. Calculate the standard error of the mean and 95% confidence intervals. Let's build a bar chart with average values and error ranges in order to visually compare manufacturers by specific power reserve and evaluate the accuracy of the estimates obtained.

стат_завод = combine(groupby(таблица, :Завод), 
    :Пробег => mean => :Среднее,
    :Пробег => std => :СтОткл,
    :Пробег => length => :N)
стат_завод.Завод = string.(стат_завод.Завод)
стат_завод.СтОшибка = стат_завод.СтОткл ./ sqrt.(стат_завод.N)
стат_завод.ДИ = 1.96 .* стат_завод.СтОшибка

гр4 = bar(стат_завод.Завод, стат_завод.Среднее,
    yerror=стат_завод.ДИ,
    legend=false,
    title="Distributions with 95% confidence intervals",
    xlabel="Manufacturer",
    ylabel="Specific power reserve",
    color=цвета_заводов,
    linewidth=2,
    size=(500, 400))
display(гр4)
println()

Let's construct a scattering diagram of all observations

таблица.Группа = string.(таблица.Завод) .* "-" .* string.(таблица.Модель)
таблица.Группа = categorical(таблица.Группа)
гр5 = plot(legend=:topleft, size=(600, 400),
    title="All observations of specific power reserve",
    xlabel="Manufacturer-Model Combination",
    ylabel="Specific power reserve")
for (i, з) in enumerate(["1", "2", "3"])
    данные = таблица[таблица.Завод .== з, :]
    scatter!(гр5, данные.Группа, данные.Пробег,
        label="Производитель $з",
        color=цвета_заводов[i],
        markersize=8,
        markerstrokewidth=1.5)
end
display(гр5)

Let's build a heat map of the distribution of specific power reserve by models and manufacturers.

матрица_средних = zeros(заводов, моделей)
for з in 1:заводов
    for м in 1:моделей
        маска = (таблица.Завод .== string(з)) .& (таблица.Модель .== string(м))
        матрица_средних[з, м] = mean(таблица.Пробег[маска])
    end
end
гр6 = heatmap(["Модель 1", "Модель 2"], 
    ["Производитель 1", "Производитель 2", "Производитель 3"],
    матрица_средних, 
    title="Specific power reserve",
    xlabel="Model",
    ylabel="Manufacturer",
    color=:viridis,
    size=(400, 300))
display(гр6)

Conclusion

The analysis made it possible to quantify the variability of the specific power reserve into components related to the car model, the manufacturer and their interaction. The likelihood ratio test provided a formal answer to the question of the significance of the fixed effect of the model, and graphical support — from box diagrams to a heat map — made these conclusions clear and interpretable. Now a test engineer can rely not on intuition, but on statistically sound estimates: whether it is worth upgrading the design of the model or directing efforts to unify production processes between plants.

The considered method is universal and goes far beyond the automotive industry. In pharmaceuticals, it is used to assess the reproducibility of a drug's effect between different production sites: in agriculture, to separate the effect of plant variety, soil type and their interaction on yield; in psychology and education, when it is necessary to take into account both the individual differences of the subjects and the effect of experimental exposure. Wherever the data has a hierarchical structure (samples within batches, patients within clinics, students within classes), mixed models become an indispensable tool for correctly assessing the effects of interest and avoiding false conclusions due to ignoring the group structure of variability.

Thus, mastering the technique demonstrated in this example provides a universal opportunity for analyzing complexly organized data in a wide variety of subject areas.