Derivatives_HelperFunctions.jl

# # Helper functions for the sensitivity analysis
# ## Load additional packages
using FiniteDifferences
using DataFrames
using CategoricalArrays
using CSV
using MixedModels
using StatsPlots
using Printf
using Combinatorics

# ## Analyze fixed models with Shapley regression
function shapley_regression(df, lhs, p)
    Δr² = similar(p, Float64)
    Δr² .= 0

    for ipick ∈ eachindex(p)
        p_pick = p[ipick]
        p_rest = filter(x -> x != p_pick, p)
        r_terms = term(1) + (term(1) | term(:seq_type)) + (term(1) | term(:seq_name)) + (term(1) | term(:ROI))

        for nterms = 0:length(p_rest)
            weight = 1 / length(combinations(p_rest, nterms)) / (length(p_rest) + 1)
            for combination ∈ combinations(p_rest, nterms)
                rhs = combination == Symbol[] ? r_terms : sum(t -> term(t), combination) + r_terms
                r²_rest = _r2(df, term(lhs), rhs)

                rhs += term(p_pick)
                r²_pick = _r2(df, term(lhs), rhs)

                Δr²[ipick] += weight * (r²_pick - r²_rest)
            end
        end
    end
    return Δr²
end

function _r2(df, lhs, rhs)
    model = fit(MixedModel, lhs ~ rhs, df)
    fixedeffects = Symbol.(rhs[collect(typeof.(rhs) .== Term)])

    σ²_fixed = isempty(fixedeffects) ? 0 : var(sum(param -> model.βs[param] .* df[!, param], fixedeffects); corrected=false)
    σ²_ROI = model.σs[:ROI][1]^2
    σ²_seqname = model.σs[:seq_name][1]^2
    σ²_seqtype = model.σs[:seq_type][1]^2
    σ²_ϵ = model.σ^2
    σ²_all = σ²_fixed + σ²_ROI + σ²_seqtype + σ²_seqname + σ²_ϵ
    r² = σ²_fixed / σ²_all
    return r²
end
nothing #hide #md

# ## Helper function:
function get_main_path()
    path_parts = splitpath(pwd())
    i_docs = findfirst(==("docs"), path_parts)
    main_path = i_docs === nothing ? pwd() : joinpath(path_parts[1:i_docs-1]...)
    return main_path
end
nothing #hide #md