Мне интересно, как начальные значения по умолчанию указаны в glm
.
Этот пост предполагает, что значения по умолчанию установлены в виде нулей. Это один говорит , что существует алгоритм позади него, однако соответствующая связь нарушена.
Я попытался согласовать простую модель логистической регрессии с алгоритмом трассировки:
set.seed(123)
x <- rnorm(100)
p <- 1/(1 + exp(-x))
y <- rbinom(100, size = 1, prob = p)
# to see parameter estimates in each step
trace(glm.fit, quote(print(coefold)), at = list(c(22, 4, 8, 4, 19, 3)))
Во-первых, без указания начальных значений:
glm(y ~ x, family = "binomial")
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
NULL
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.386379 1.106234
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3991135 1.1653971
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3995188 1.1669508
На первом этапе начальные значения NULL
.
Во-вторых, я устанавливаю начальные значения равными нулю:
glm(y ~ x, family = "binomial", start = c(0, 0))
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0 0
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3177530 0.9097521
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3909975 1.1397163
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3994147 1.1666173
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3995191 1.1669518
И мы видим, что итерации между первым и вторым подходом различаются.
Чтобы увидеть начальные значения, указанные glm
мной, я попытался согласовать модель только с одной итерацией:
glm(y ~ x, family = "binomial", control = list(maxit = 1))
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
NULL
Call: glm(formula = y ~ x, family = "binomial", control = list(maxit = 1))
Coefficients:
(Intercept) x
0.3864 1.1062
Degrees of Freedom: 99 Total (i.e. Null); 98 Residual
Null Deviance: 134.6
Residual Deviance: 115 AIC: 119
Оценки параметров (что неудивительно) соответствуют оценкам первого подхода во второй итерации, т. Е. [1] 0.386379 1.106234
Установка этих значений в качестве начальных значений приводит к той же последовательности итераций, что и в первом подходе:
glm(y ~ x, family = "binomial", start = c(0.386379, 1.106234))
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.386379 1.106234
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3991135 1.1653971
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3995188 1.1669508
Итак, вопрос в том, как эти значения рассчитываются?
источник
start
значения, они используются при расчете того, что передается вC_Cdqrls
рутину. В противном случае переданные значения вычисляются (включая вызовeval(binomial()$initialize)
), ноglm.fit
никогда явно не вычисляют значения дляstart
. Потратьте час или два и изучитеglm.fit
код.glm.fit
код, но до сих пор не знаю, как рассчитываются начальные значения.Ответы:
TL; DR
start=c(b0,b1)
инициализирует eta вb0+x*b1
(mu в 1 / (1 + exp (-eta)))start=c(0,0)
инициализирует eta к 0 (mu к 0.5) независимо от значения y или x.start=NULL
инициализирует eta = 1.098612 (mu = 0.75), если y = 1, независимо от значения x.start=NULL
инициализирует eta = -1.098612 (mu = 0.25), если y = 0, независимо от значения x.После того, как ETA (и , следовательно му и вар (MU)) были вычислены,
w
иz
вычисляются и отправляются в QR решатель, в духеqr.solve(cbind(1,x) * w, z*w)
.полной формы
Строительство от комментариев Роланда: я сделал
glm.fit.truncated()
, когда я взялglm.fit
доC_Cdqrls
звонка, а затем прокомментировал это.glm.fit.truncated
выводитz
иw
значение (а также значения величин , используемых для расчетаz
иw
) , которые затем быть переданыC_Cdqrls
вызовом:Подробнее можно прочитать
C_Cdqrls
здесь . К счастью, функцияqr.solve
в base R напрямую подключается к версиям LINPACK, вызываемым вglm.fit()
.Таким образом, мы запускаем
glm.fit.truncated
различные спецификации начальных значений, а затем выполняем вызовqr.solve
со значениями w и z и видим, как вычисляются «начальные значения» (или первые отображаемые значения итерации). Как указал Роланд, указаниеstart=NULL
илиstart=c(0,0)
в glm () влияет на вычисления для w и z, а не дляstart
.Для начала = NULL:
z
это вектор, где элементы имеют значение 2.431946 или -2.431946, иw
это вектор, в котором все элементы имеют значение 0,4330127:Для начала = c (0,0):
z
это вектор, где элементы имеют значение 2 или -2, иw
это вектор, где все элементы равны 0,5:Так что это все хорошо, но как рассчитать
w
иz
? В нижней частиglm.fit.truncated()
мы видимПосмотрите на следующие сравнения между выведенными значениями величин, используемых для расчета,
z
иw
:Обратите внимание, что
start.is.00
вектор будет иметьmu
только значения 0,5, потому что значение eta равно 0, а mu (eta) = 1 / (1 + exp (-0)) = 0,5.start.is.null
устанавливает те, у которых y = 1, равным mu = 0,75 (что соответствует eta = 1,098612), а те, у которых y = 0, равны mu = 0,25 (что соответствует eta = -1,098612), и, таким образом,var_mu
= 0,75 * 0,25 = 0,1875.Тем не менее, интересно отметить, что я изменил начальное число и перезапустил все, и mu = 0,75 для y = 1 и mu = 0,25 для y = 0 (и, таким образом, остальные величины остались прежними). То есть start = NULL приводит к тому же
w
иz
независимо от того, чтоy
иx
есть, потому что они инициализируют eta = 1.098612 (mu = 0.75), если y = 1, и eta = -1.098612 (mu = 0.25), если y = 0.Таким образом, представляется, что начальное значение для коэффициента перехвата и для X-коэффициента не установлено для start = NULL, а вместо этого начальные значения задаются для eta в зависимости от значения y и не зависят от значения x. Оттуда
w
иz
рассчитываются, затем отправляются вместе сx
qr.solver.Код для запуска перед фрагментами выше:
источник