基于汽车特定范围数据示例的混合效应方差分析
导言
想象一下，你是一名测试工程师，你需要找出不同的车型是否真的在特定的动力储备上有所不同，或者观察到的差异是否只是一个事故，例如与特定工厂 这个问题的答案由方差分析提供，方差分析是一种统计方法，将数据的整体变异性分解为由不同来源引起的部分。 关键思想如下：如果组之间的差异显着超过随机"噪声"的水平，则该因子的影响被识别为显着。
此示例使用混合线性模型。 每个汽车模型充当固定效应，制造商和工厂与模型的相互作用充当随机效应。
训练模型的结果是方差分量，它量化了每个源对整体传播的贡献。:
-
（制造商）-来自不同工厂的汽车平均有多大不同;
-
（制造商×模型）-模型从植物到植物的效果不一致程度（因素的相互作用）;
-
（余数）是与汽车的个体特性和测量误差有关的无法解释的变化。
为了直观地表示数据，我们使用箱线图（boxplot，或"带胡子的盒子"）。 此图紧凑地显示中位数，四分位数，范围和潜在异常值，使您能够立即估计每个组中的分布中心，其分布和对称性。
.png)
正式假设检验依赖于似然比检验：它将完整模型与减少（无固定效应模型）进行比较，并确定汽车模型的排除是否导致数据描述的统计显着恶化。
通过将混合效应方差分析的能力与表现性统计图形相结合，我们重建了隐藏在测试报告中数字列后面的客观图像。
导入库
我们将添加使用方差分析所需的库。
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
初始数据
让我们确定每个车型的特定动力储备的初始数据（пробег). 我们将数据转换为数据对齐的表。
# 初始数据
исходные = [
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
]
# 实验计划的参数
заводов = 3
строк = 6
моделей = 2
повторов = div(строк, моделей)
пробег = Float64[]
завод = Int[]
модель = Int[]
for з in 1:заводов
for стр in 1:строк
push!(пробег, исходные[стр, з])
push!(завод, з)
push!(модель, стр <= 3 ? 1 : 2)
end
end
# 数据表
таблица = DataFrame(
Пробег = пробег,
Завод = string.(завод),
Модель = string.(модель)
)
таблица.Завод = categorical(таблица.Завод)
таблица.Модель = categorical(таблица.Модель)
println("汽车的特定动力储备")
println(first(таблица, 18))
建筑模型
让我们建立一个混合线性模型，其中里程由汽车模型的固定效应和工厂的随机效应以及工厂与模型的相互作用来解释。 让我们使用最大似然法选择模型参数，并分析获得的方差分量和固定效应的估计值。
формула1 = @formula(Пробег ~ 1 + Модель + (1 | Завод) + (1 | Завод & Модель))
подбор1 = fit(MixedModel, формула1, таблица)
println(подбор1)
让我们构建一个没有汽车模型固定效果的缩减混合模型，只留下工厂和交互的随机效应。 让我们在完整模型和缩减模型之间进行似然比测试，以便统计验证汽车模型对里程的影响的显着性。
формула0 = @formula(Пробег ~ 1 + (1 | Завод) + (1 | Завод & Модель))
подбор0 = fit(MixedModel, формула0, таблица)
тест_правд = MixedModels.likelihoodratiotest(подбор0, подбор1)
println("似然比测试结果:")
println(тест_правд)
方差分量的估计
让我们使用函数从训练好的模型中提取方差分量 VarCorr 让我们将每个随机效应的数值方差值存储在一个方便的数据结构中。
компоненты = VarCorr(подбор1)
println(компоненты)
# 提取方差的函数
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("随机效应差异的估计:")
println(" σ2(制造商)= ", round(σ²_завод, digits=4))
println(" σ2(制造商:型号)= ", round(σ²_взаимодействие, digits=4))
println(" σ2(误差)= ", round(σ²_остаток, digits=4))
结果可视化
让我们设置图形的视觉设计。
theme(:default)
default(fontfamily="sans-serif", titlefontsize=11, guidefontsize=9, tickfontsize=7)
цвета_заводов = [:steelblue, :darkorange, :forestgreen]
цвета_комп = [:steelblue, :darkorange, :forestgreen]
使用boxplot，我们将绘制每个制造商的汽车特定动力储备的分布。
gr()
гр1 = @df таблица boxplot(:Завод, :Пробег,
legend=false,
title="比功率储备的分配",
xlabel="制造商",
ylabel="特定动力储备",
color=цвета_заводов,
linewidth=2,
size=(500, 400))
display(гр1)
.png)
让我们绘制每个模型的汽车特定动力储备的分布。
гр2 = @df таблица boxplot(:Модель, :Пробег,
legend=false,
title="比功率储备的分配",
xlabel="汽车模型",
ylabel="特定动力储备",
color=[:salmon, :lightblue],
linewidth=2,
size=(500, 400))
display(гр2)
.png)
让我们构建一个制造商-模型交互图。
средние_зм = combine(groupby(таблица, [:Завод, :Модель]), :Пробег => mean => :Средний)
средние_зм.Завод = string.(средние_зм.Завод)
средние_зм.Модель = string.(средние_зм.Модель)
гр3 = plot(legend=:bottomright, size=(500, 400),
title="制造商-模型交互时间表",
xlabel="制造商",
ylabel="特定动力储备")
цвета_моделей = [:blue, :red]
for (i, мод) in enumerate(["1", "2"])
данные = средние_зм[средние_зм.Модель .== мод, :]
plot!(гр3, данные.Завод, данные.Средний,
marker=:circle, markersize=8, linewidth=2,
color=цвета_моделей[i],
label="Модель $мод")
end
display(гр3)
.png)
让我们按制造商对数据进行分组，并计算每个数据的比功率储备、标准偏差和观测值数。 计算均值和95%置信区间的标准误差。 让我们建立一个具有平均值和误差范围的条形图，以便按特定功率储备直观地比较制造商，并评估获得的估计的准确性。
стат_завод = combine(groupby(таблица, :Завод),
:Пробег => mean => :Среднее,
:Пробег => std => :СтОткл,
:Пробег => length => :N)
стат_завод.Завод = string.(стат_завод.Завод)
стат_завод.СтОшибка = стат_завод.СтОткл ./ sqrt.(стат_завод.N)
стат_завод.ДИ = 1.96 .* стат_завод.СтОшибка
гр4 = bar(стат_завод.Завод, стат_завод.Среднее,
yerror=стат_завод.ДИ,
legend=false,
title="置信区间为95%的分布",
xlabel="制造商",
ylabel="特定动力储备",
color=цвета_заводов,
linewidth=2,
size=(500, 400))
display(гр4)
println()
.png)
让我们构建所有观测值的散射图
таблица.Группа = string.(таблица.Завод) .* "-" .* string.(таблица.Модель)
таблица.Группа = categorical(таблица.Группа)
гр5 = plot(legend=:topleft, size=(600, 400),
title="所有关于特定动力储备的观察",
xlabel="制造商-型号组合",
ylabel="特定动力储备")
for (i, з) in enumerate(["1", "2", "3"])
данные = таблица[таблица.Завод .== з, :]
scatter!(гр5, данные.Группа, данные.Пробег,
label="Производитель $з",
color=цвета_заводов[i],
markersize=8,
markerstrokewidth=1.5)
end
display(гр5)
.png)
让我们建立一个模型和制造商比功率储备分布的热图.
матрица_средних = zeros(заводов, моделей)
for з in 1:заводов
for м in 1:моделей
маска = (таблица.Завод .== string(з)) .& (таблица.Модель .== string(м))
матрица_средних[з, м] = mean(таблица.Пробег[маска])
end
end
гр6 = heatmap(["Модель 1", "Модель 2"],
["Производитель 1", "Производитель 2", "Производитель 3"],
матрица_средних,
title="特定动力储备",
xlabel="模型",
ylabel="制造商",
color=:viridis,
size=(400, 300))
display(гр6)
.png)
结论
该分析使得将特定动力储备的可变性量化为与汽车模型，制造商及其相互作用相关的组件成为可能。 似然比检验为模型的固定效应的显着性问题提供了正式的答案，图形支持—从箱形图到热图-使这些结论清晰可解释。 现在，测试工程师可以不依赖于直觉，而是依赖于统计上可靠的估计：是否值得升级模型的设计或指导统一工厂之间的生产流程。
所考虑的方法是通用的，并且远远超出了汽车行业。 在药品中，它被用来评估药物在不同生产地点之间作用的可重复性：在农业中，分离植物品种、土壤类型及其对产量的相互作用的影响；在心理学和教育中，当需要考虑受试者的个体差异和实验暴露的影响时。 只要数据具有层次结构（批次内的样本、诊所内的病人、班级内的学生），混合模型就成为正确评估感兴趣的影响和避免因忽略变异性的群体结构而得出错误结论的不可或缺的工具。
因此，掌握本示例中演示的技术为分析各种主题领域中复杂组织的数据提供了一个普遍的机会。