Документация Engee

Руководство по GLM.jl

Линейные и обобщенные линейные модели в Julia

Установка

Pkg.add("GLM")

устанавливает этот пакет и его зависимости, включая пакет Distributions.

Пакет RDatasets полезен для подгонки моделей на основе стандартных наборов данных R с целью сравнения результатов с результатами в R.

Подгонка моделей GLM

Для подгонки обобщенной линейной модели (GLM) можно использовать два метода: glm(formula, data, family, link) и glm(X, y, family, link). Аргументы должны быть следующими.

  • formula: объект Formula из StatsModels.jl, ссылающийся на столбцы в data; например, если столбцы имеют имена :Y, :X1 и :X2, правильной формулой будет @formula(Y ~ X1 + X2)

  • data: таблица из определения Tables.jl, например фрейм данных; строки со значениями missing игнорируются

  • X: матрица, содержащая значения независимых переменных в столбцах

  • y: вектор, содержащий значения зависимых переменных в столбцах (включая, если применимо, свободный коэффициент)

  • family: один из вариантов Bernoulli(), Binomial(), Gamma(), Geometric(), Normal(), Poisson() или NegativeBinomial(θ)

  • link: вариант из приведенного ниже списка, например LogitLink() — допустимая связь для семейства Binomial()

Типичные распределения для использования с glm и их канонические связывающие функции:

       Bernoulli (LogitLink)
        Binomial (LogitLink)
           Gamma (InverseLink)
       Geometric (LogLink)
 InverseGaussian (InverseSquareLink)
NegativeBinomial (NegativeBinomialLink, часто используется с LogLink)
          Normal (IdentityLink)
         Poisson (LogLink)

В настоящее время доступны следующие типы связей:

CauchitLink
CloglogLink
IdentityLink
InverseLink
InverseSquareLink
LogitLink
LogLink
NegativeBinomialLink
PowerLink
ProbitLink
SqrtLink

Обратите внимание, что канонической связью для отрицательной биномиальной регрессии является NegativeBinomialLink, но на практике обычно используется LogLink. Распределение NegativeBinomial относится к экспоненциальному семейству только в том случае, если θ (параметр формы) является фиксированным, поэтому при использовании glm с семейством NegativeBinomial необходимо задать θ. Если необходимо также оценить θ, то вместо этого следует использовать negbin(formula, data, link).

По умолчанию в любую модель GLM включен свободный коэффициент.

Категориальные переменные

Категориальные переменные по умолчанию задаются в коде как фиктивные, если они не являются числовыми или представляют собой векторы CategoricalVector в таблице Tables.jl (DataFrame, таблице JuliaDB, именованном кортеже векторов и т. д.). Вы также можете явно передать аргумент contrasts, если требуется другая система кодирования контрастов или вы не используете фреймы DataFrame.

Переменная реакции (зависимая) не может быть категориальной.

Использование вектора CategoricalVector, созданного с помощью categorical или categorical!:

julia> using CategoricalArrays, DataFrames, GLM, StableRNGs

julia> rng = StableRNG(1); # Обеспечиваем воспроизводимость примера

julia> data = DataFrame(y = rand(rng, 100), x = categorical(repeat([1, 2, 3, 4], 25)));


julia> lm(@formula(y ~ x), data)
StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}}

y ~ 1 + x

Coefficients:
───────────────────────────────────────────────────────────────────────────
                  Coef.  Std. Error      t  Pr(>|t|)   Lower 95%  Upper 95%
───────────────────────────────────────────────────────────────────────────
(Intercept)   0.490985    0.0564176   8.70    <1e-13   0.378997    0.602973
x: 2          0.0527655   0.0797865   0.66    0.5100  -0.105609    0.21114
x: 3          0.0955446   0.0797865   1.20    0.2341  -0.0628303   0.25392
x: 4         -0.032673    0.0797865  -0.41    0.6831  -0.191048    0.125702
───────────────────────────────────────────────────────────────────────────

С использованием contrasts.

julia> using StableRNGs

julia> data = DataFrame(y = rand(StableRNG(1), 100), x = repeat([1, 2, 3, 4], 25));

julia> lm(@formula(y ~ x), data, contrasts = Dict(:x => DummyCoding()))
StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}}

y ~ 1 + x

Coefficients:
───────────────────────────────────────────────────────────────────────────
                  Coef.  Std. Error      t  Pr(>|t|)   Lower 95%  Upper 95%
───────────────────────────────────────────────────────────────────────────
(Intercept)   0.490985    0.0564176   8.70    <1e-13   0.378997    0.602973
x: 2          0.0527655   0.0797865   0.66    0.5100  -0.105609    0.21114
x: 3          0.0955446   0.0797865   1.20    0.2341  -0.0628303   0.25392
x: 4         -0.032673    0.0797865  -0.41    0.6831  -0.191048    0.125702
───────────────────────────────────────────────────────────────────────────

Сравнение моделей с помощью F-теста

Сравнение двух или более линейных моделей можно проводить с помощью функции ftest, которая выполняет F-тест для каждой пары последовательных моделей и сообщает статистику соответствия:

julia> using DataFrames, GLM, StableRNGs

julia> data = DataFrame(y = (1:50).^2 .+ randn(StableRNG(1), 50), x = 1:50);

julia> ols_lin = lm(@formula(y ~ x), data);

julia> ols_sq = lm(@formula(y ~ x + x^2), data);

julia> ftest(ols_lin.model, ols_sq.model)
F-test: 2 models fitted on 50 observations
─────────────────────────────────────────────────────────────────────────────────
     DOF  ΔDOF           SSR           ΔSSR      R²     ΔR²            F*   p(>F)
─────────────────────────────────────────────────────────────────────────────────
[1]    3        1731979.2266                 0.9399
[2]    4     1       40.7581  -1731938.4685  1.0000  0.0601  1997177.0357  <1e-99
─────────────────────────────────────────────────────────────────────────────────

Методы, применимые к подогнанным моделям

Многие методы в этом пакете имеют имена, похожие на используемые в R.

  • adjr2: скорректированный коэффициент детерминации для линейной модели (псевдоним для adjr²)

  • aic: информационный критерий Акаике

  • aicc: скорректированный информационный критерий Акаике для выборок небольшого размера

  • bic: байесовский информационный критерий

  • coef: оценки коэффициентов модели

  • confint: доверительные интервалы для коэффициентов

  • cooksdistance: расстояние Кука для каждого наблюдения

  • deviance: мера соответствия модели, взвешенная остаточная сумма квадратов для линейной модели

  • dispersion: параметр дисперсии (или масштаба) для распределения модели

  • dof: количество степеней свободы, используемых в модели

  • dof_residual: степени свободы для остатков, если имеет смысл

  • fitted: подобранные значения модели

  • glm: подгоняет обобщенную линейную модель (псевдоним для fit(GeneralizedLinearModel, ...))

  • lm: подгоняет линейную модель (псевдоним для fit(LinearModel, ...))

  • loglikelihood: логарифмическое правдоподобие модели

  • modelmatrix: матрица эксперимента

  • nobs: количество строк или сумма весов, если заданы априорные веса

  • nulldeviance: отклонение модели с удаленными предикторами

  • nullloglikelihood: логарифмическое правдоподобие модели с удаленными предикторами

  • predict: прогнозируемые значения зависимой переменной из подогнанной модели

  • r2: коэффициент детерминации линейной модели (псевдоним для r²)

  • residuals: вектор остатков из подогнанной модели

  • response: реакция модели (зависимая переменная)

  • stderror: среднеквадратичные погрешности коэффициентов

  • vcov: матрица дисперсий и ковариаций для оценочных коэффициентов

Обратите внимание, что канонической связью для отрицательной биномиальной регрессии является NegativeBinomialLink, но на практике обычно используется LogLink.

julia> using GLM, DataFrames, StatsBase

julia> data = DataFrame(X=[1,2,3], y=[2,4,7]);

julia> mdl = lm(@formula(y ~ X), data);

julia> round.(coef(mdl); digits=8)
2-element Vector{Float64}:
 -0.66666667
  2.5

julia> round(r2(mdl); digits=8)
0.98684211

julia> round(aic(mdl); digits=8)
5.84251593

Метод predict возвращает прогнозируемые значения переменной реакции на основе ковариантных значений во входном объекте newX. Если аргумент newX не задан, возвращаются подобранные значения реакции из модели.

julia> test_data = DataFrame(X=[4]);

julia> round.(predict(mdl, test_data); digits=8)
1-element Vector{Float64}:
 9.33333333

Метод cooksdistance возвращает расстояние Кука для каждого наблюдения, используемого для подгонки линейной модели, с учетом оценки влияния каждой точки данных. Обратите внимание, что в настоящее время он реализован только для линейных моделей без весов.

julia> round.(cooksdistance(mdl); digits=8)
3-element Vector{Float64}:
 2.5
 0.25
 2.5

Разделение объекта реакции и объекта предиктора

Согласно общему подходу в этом коде функции, связанные с реакцией, отделяются от функций, связанных с линейным предиктором. Это обеспечивает большую универсальность благодаря возможности комбинировать различные подтипы абстрактных типов LinPred и ModResp.

Тип LinPred включает в себя вектор параметров и матрицу модели. Вектор параметров представляет собой плотный числовой вектор, а матрица модели может быть плотной или разреженной. Тип LinPred должен предусматривать какое-либо разложение взвешенной матрицы модели, обеспечивающее решение системы X’W * X * delta=X’wres, где W — диагональная матрица весов X, предоставленных в виде вектора квадратных корней из диагональных элементов, а wres — взвешенный вектор остатков.

В настоящее время есть два типа плотных предикторов: DensePredQR и DensePredChol, с их обычными особенностями. Вариант Холецкого обрабатывается быстрее, но немного менее точно, чем вариант QR. В коде заложена основа для распределенного типа предиктора, но окончательно он еще не воплощен. Так как в Julia по умолчанию используется библиотека OpenBLAS, которая уже является многопоточной на многоядерных компьютерах, применение распределенных типов предикторов может не дать особых преимуществ.

Тип ModResp должен предоставлять методы для универсальных функций wtres и sqrtxwts. Их значения являются аргументами для методов updatebeta типа LinPred. Значение Float64, возвращаемое updatedelta, — это значение критерия сходимости.

Аналогичным образом, типы LinPred должны предоставлять метод для универсальной функции linpred. Как правило, linpred принимает экземпляр типа LinPred и коэффициент шага. Методы, принимающие только экземпляр LinPred, используют коэффициент шага по умолчанию, равный 1. Значение linpred является аргументом метода updatemu для типов ModResp. Метод updatemu возвращает обновленное отклонение.