Команда glm в R (или функция GLM из библиотеки statsmodel.api в Python) может соответствовать модели GLM без необходимости знать все сложные уравнения и методы оптимизации, используемые «за кулисами». Но умение что-то кодировать вручную обычно является признаком глубокого понимания.
Итак, в этой статье мы собираемся установить GLM вручную в R и сравнить его со встроенным glm.
Уравнения вероятности
Уравнения правдоподобия относятся к первой производной логарифмического правдоподобия модели GLM. ОМЛ имеют общетеоретический анализ — в нем рассматривается общий случай экспоненциального семейства (точнее, экспоненциальной дисперсионной модели [ЭДМ], являющейся подмножеством экспоненциального семейства), который выглядит следующим образом:
Где тета — «естественный параметр», который является функцией среднего мю; b(theta) — логарифмическая статистическая сумма; среднее значение связано с линейным предиктором эта через некоторую функцию связи; а линейный предиктор — это скалярное произведение x на бета (интересующие нас коэффициенты). Итак, взяв производную логарифмической вероятности относительно бета, мы можем использовать цепное правило, чтобы получить:
Вывод дель-тета / дель-мю использует то, что называется отношением средней дисперсии. В двух словах - с помощью этого анализа общей формы можно показать, что первая производная функции b (тета) (относительно тета) дает среднее значение распределения. А вторая производная дает функцию V(mu), которая фиксирует взаимосвязь между средним значением и дисперсией.
Мы можем записать это в матричной записи:
Где V — диагональная матрица функции среднего отклонения V(mu), а D — диагональная матрица для выбранной функции связи del-mu / del-eta. Обратите внимание, что если мы выбираем функцию связи так, что тета = эта, две частные производные сокращаются, и поэтому матрицы также сокращаются. Это называется «канонической» функцией связи, и она упрощает вывод:
Подсчет очков Фишера
GLM использует метод оптимизации 2-го порядка, называемый подсчетом Фишера. Он очень похож на метод Ньютона-Рафсона, только вместо гессиана (второй производной) логарифмического правдоподобия мы берем его ожидаемое значение, равное обратная) Информационная матрица Фишера.
Можно показать, что информационная матрица Фишера для ОЛМ равна:
Где W — диагональная матрица с этими элементами по диагонали.
Внедрение оценки Фишера
Теперь у нас есть все необходимое для реализации алгоритма оценки Фишера, поэтому давайте сделаем это для регрессии Пуассона. Функция средней дисперсии для распределения Пуассона имеет вид V(mu)=mu. В регрессии Пуассона мы видим, что каноническая функция связи — это логарифмическая функция (ее обратная — экспоненциальная), поэтому давайте используем каноническую функцию связи в качестве нашей функции связи. Это упростит вывод и для матрицы информации Фишера:
Давайте сначала сгенерируем данные:
set.seed(247) n=100 x = runif(n,-3,3) lp = 3+0.4*x mu = exp(lp) y = rpois(n,mu)
Таким образом, истинные коэффициенты бета0 = 3, бета1 = 0,4. Давайте построим это:
library(ggplot2) # for plotting ggplot(data.frame(x=x, y=y)) + geom_point(aes(x=x,y=y), color="#2E9FDF", fill="#2E9FDF", alpha=0.3)
Таким образом, алгоритм следующий:
Или в коде:
library(pracma) # for matrix inversion
set.seed(247)
# let's create the design matrix
X = cbind(rep(1,n), x)
# start with random values for beta's
beta_0 = c(rnorm(1),rnorm(1))
# max number of iteration for the Fisher-scoring algorithm
max_iter = 100
# Do Fisher-scoring
beta_old = beta_0
for (i in 1:max_iter) {
mu = c(exp(X%*%beta_old))
V = diag(mu)
beta_t = beta_old + inv(t(X)%*%V%*%X) %*% t(X) %*% (y-mu)
# if the coefficients change is below some threshold, we stop
nrm = norm(beta_t-beta_old, type="2")
if (nrm < 0.05)
break
beta_old = beta_t
}
print(i)
# [1] 21
Потребовалась 21 итерация, и мы получили довольно близкие результаты:
beta_t # [,1] # 2.9960474 # x 0.3926822
Давайте сравним это с glm:
fit.poisson = glm(y~x, family="poisson")
summary(fit.poisson)
Call:
glm(formula = y ~ x, family = "poisson")
Deviance Residuals:
Min 1Q Median 3Q Max
-2.7848 -0.7181 -0.0102 0.5873 3.2298
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.99601 0.02533 118.26 <2e-16 ***
x 0.39270 0.01443 27.22 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 959.35 on 99 degrees of freedom
Residual deviance: 106.49 on 98 degrees of freedom
AIC: 601.1
Number of Fisher Scoring iterations: 4
Здесь потребовалось всего 4 итерации, и мы получили почти идентичное решение — но не совсем идентичное. Почему это?
IRLS — Итерированные перевзвешенные наименьшие квадраты
R использует вариант оценки Фишера с повторным взвешиванием наименьших квадратов (IRLS). Кроме того, R проверяет не сходимость бета, а сходимость остаточного отклонения. И R также использует другую начальную точку — она инициирует мю, а не бета (бета — это функция мю).
IRLS был известным алгоритмом в 70-х годах, поэтому Нелдер и Веддерберн, которые изобрели структуру GLM, хотели внедрить оценку Фишера, используя уже существующее программное обеспечение. Это потребовало еще немного вывода и небольшой хитрости. Обратите внимание на взаимосвязь между матрицами D и V в градиенте и матрицей W в (ожидаемо-отрицательном) гессиане:
Таким образом, мы можем записать оценку Фишера как:
Где мы использовали хитрость и умножили исходную бету на слагаемые, отменяющие друг друга. Z определяется как «рабочий отклик» (рабочие y), поэтому каждую итерацию оценки Фишера можно рассматривать как выполнение (обобщенного) метода наименьших квадратов для рабочего отклика.
Для модели Пуассона с канонической функцией зацепления имеем
Остаточное отклонение для модели Пуассона:
И начальными значениями для мю являются сами у + 0,1 (во избежание ошибок, поскольку у может быть 0, а мю должны быть строго положительными).
Давайте закодируем это:
set.seed(247)
n=100
x = runif(n,-3,3)
X = cbind(rep(1,n), x)
beta_true = c(3,0.4)
lp = X%*%beta_true
mu = exp(lp)
y = rpois(n,mu)
# start with
mu0 = y+0.1
# max number of iteration for the IRLS algorithm
max_iter = 100
# Do IRLS
mu_t = mu0
for (i in 1:max_iter) {
zt = log(mu_t)+y/mu_t-1 # calculate the working response
V = diag(mu_t)
beta_hat = solve(t(X)%*%V%*%X) %*% t(X)%*%V%*% zt # GLS on the working response
# calculate the residual deviance before and after
deviance_t = 2*sum(y*log(y/mu_t)-y+mu_t)
mu_t = c(exp(X%*%beta_hat)) # update mu
deviance_t_1 = 2*sum(y*log(y/mu_t)-y+mu_t)
# if the Deviance change is below some threshold, we stop
change = abs(deviance_t_1-deviance_t)/(abs(deviance_t_1)+0.1)
if (change < 1e-8)
break
}
print(i)
# [1] 4
beta_hat
# [,1]
# 2.9960065
# x 0.3927027
deviance_t_1
# [1] 106.4857
Теперь у нас получилось точно так же, как glm!
Если вы хотите узнать больше, загляните на мой канал YouTube series и на мой онлайн-курс GLM.
