# import Pkg
# Pkg.add(["AnovaGLM", "CSV", "DataFrames", "CategoricalArrays", "MixedModels", "Statistics", "Distributions", "GLM", "StatsPlots"])
using AnovaGLM, CSV, DataFrames, CategoricalArrays, MixedModels, Statistics, Distributions, GLM, StatsPlots

# 初始数据
来源=[
    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[]

力在1:工厂
    对于第1页：行
        推！(里程,初始[页,s])
        推！(工厂,z)
        推！（模型，页<=3？ 1 : 2)
    end
end

# 数据表
表=DataFrame（
    哩程=哩程,
    工厂=字符串。(工厂),
    模型=字符串。(型号)
)
表。工厂=分类（表。工厂)
表。模型=分类（表。的模型）

println("汽车的特定动力储备")
println(第一(表,18))

公式1=@公式（里程~1+模型+（1/工厂）+（1/工厂和模型））
fit1=fit（MixedModel，公式1，表）
println(选择1)

公式0=@公式（里程~1+（1/工厂）+（1/工厂和模型））
fit0=fit（MixedModel，formula0，table）
test_pravd=MixedModels。likelihoodratiotest（选择0，选择1）
println("似然比测试结果:")
println(test_pravd)

组件=VarCorr(选择1)
println(组件)
# 提取方差的函数
函数提取色散(vc)
    string_vc=字符串（vc）
    方差=Dict{String,Float64}()
    for line in split(строка_vc, "\n")
        清洁=条带(线)
        如果isempty（清除）
            continue
        end
        if occursin("Variance components:", очищенная) || 
           occursin("Column", очищенная) || 
           (occursin("Variance", очищенная) && occursin("Std.Dev.", очищенная)) ||
           startswith(очищенная, "──")
            continue
        end
        
        части = split(очищенная, r"\s+")
        如果长度（部分）>=3
            倒数第二=零件[end-1]
            if occursin(r"^\d+\.?\d*$", предпоследнее) || occursin(r"^\d*\.?\d+$", предпоследнее)
                zn=tryparse（Float64，倒数第二）
                如果zn！==没什么
                    имя = join(части[1:end-2], " ")
                    方差[名称]=zn
                end
            end
        end
    end
    方差的返回
end

disp=提取分散体（组分）
σ2_zavod=0.0
σ2_interaction=0.0
σ2_state=0.0

for(k,v)int
    if occursin("制造商及型号", k)
        σ2_interaction=v
    elseif occursin("制造商", k) && !occursin("模型", k)
        σ2_zavod=v
    elseif occursin("Residual", k)
        σ2_state=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)
tsveta_zavodov=[:steelblue,:darkorange,:forestgreen]
color_comp=[:steelblue,:darkorange,:forestgreen]

gr()
gr1=@df箱线图表（：工厂，：里程,
    legend=false,
    title="比功率储备的分配",
    xlabel="制造商",
    ylabel="特定动力储备",
    颜色=color_zavodov,
    linewidth=2,
    size=(500, 400))
显示(gr1)

gr2=@df箱线图表（：型号，：里程,
    legend=false,
    title="比功率储备的分配",
    xlabel="汽车模型",
    ylabel="特定动力储备",
    color=[:salmon, :lightblue],
    linewidth=2,
    size=(500, 400))
显示(gr2)

mean_zm=combine(groupby(table,[:Factory,:Model]),:Mileage=>mean=>:Average)
中间主义。工厂=字符串。（平均年龄。工厂)
中间主义。模型=字符串。（平均年龄。的模型）
gr3=plot(图例=:bottomright,size=(500,400),
    title="制造商-模型交互时间表",
    xlabel="制造商",
    ylabel="特定动力储备")
colormodels=[:蓝色,:红色]
for (i, мод) in enumerate(["1", "2"])
    数据=平均值[平均值。模特。==国防部,:]
    阴谋！（gr3，数据。工厂，数据。平均,
        marker=:circle, markersize=8, linewidth=2,
        颜色=颜色模型[i],
        label="$Mod模型")
end
显示(gr3)

stat_zavod=combine(groupby(table,:Factory), 
    :里程=>均值=>:平均值,
    ：里程=>std=>：StOtkl,
    ：里程=>长度=>：N）
stat_zavod。工厂=字符串。（stat_zavod。工厂)
stat_zavod。错误=stat_zavod。号/sqrt.（stat_zavod。N）
stat_zavod。DI=1.96。*stat_zavod。一百个错误

gr4=bar(stat_zavod.工厂，stat_zavod。平均,
    yerror=stat_zavod。迪伊,
    legend=false,
    title="置信区间为95%的分布",
    xlabel="制造商",
    ylabel="特定动力储备",
    颜色=color_zavodov,
    linewidth=2,
    size=(500, 400))
显示器(gr4)
println()

таблица.Группа = string.(таблица.Завод) .* "-" .* string.(таблица.Модель)
表。组=分类（表。组）
gr5=plot(图例=:topleft,尺寸=(600,400),
    title="所有关于特定功率储备的观察",
    xlabel="制造商-型号组合",
    ylabel="特定动力储备")
for (i, з) in enumerate(["1", "2", "3"])
    数据=表[表。工厂。==z，：]
    散开！（gr5，数据。组，数据。哩程数,
        label="制造商№z",
        颜色=color_zavodov[i],
        markersize=8,
        markerstrokewidth=1.5)
end
显示(gr5)

平均值矩阵=零（工厂，模型）
力在1:工厂
    对于m in1：模型
        掩模=（表。工厂。==字符串（h））。&(表.模特。==字符串(m))
        平均值矩阵[z,m]=平均值(表).里程[面具])
    end
end
гр6 = heatmap(["模型1", "模型2"], 
    ["制造商1", "制造商2", "制造商3"],
    平均矩阵, 
    title="特定动力储备",
    xlabel="模型",
    ylabel="制造商",
    color=:viridis,
    size=(400, 300))
显示(gr6)