Руководство по GLM.jl
Установка
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 возвращает обновленное отклонение.