Как в R извлечь матрицу шляпы/проекции/влияния или значения из объекта модели `nls`?

Для объектов типа lm или glm или даже для объектов типа lmer можно извлечь значения шляпы из модели с помощью функции R hatvalues(). Однако, по-видимому, это не работает с объектами nls. Я гуглил по-всякому, но не могу найти способ получить эти значения. Разве nls просто не создает матрицу шляпы, или значения шляпы, полученные из нелинейной модели наименьших квадратов, просто как-то ненадежны?

Воспроизводимый пример:

xs = rep(1:10, times = 10)
ys = 3 + 2*exp(-0.5*xs)
for (i in 1:100) {
xs[i] = rnorm(1, xs[i], 2)
}
df1 = data.frame(xs, ys)
nls1 = nls(ys ~ a + b*exp(d*xs), data=df1, start=c(a=3, b=2, d=-0.5))

r nls non-linear-regression
person Bajcz    schedule 26.08.2016    source источник
comment
Привет @lmo, я исправил пример, чтобы он теперь правильно соответствовал nls, что является базой R. Извините за это.   -  person Bajcz    schedule 26.08.2016
comment
Хотя это далеко не окончательный ответ, в этом сообщении в Википедии нелинейный метод наименьших квадратов не упоминается как оценка поддается построению шляпной матрицы. Я подозреваю, что это, вместе с отсутствием метода hatvalues для объектов nls, означает, что этот расчет не является правильным методом / невозможен для оценки nls. crossValidated больше подходит для дальнейшего продвижения по этому пути.   -  person lmo    schedule 26.08.2016
comment
Спасибо за комментарий - я не против перенести этот вопрос в crossValidated, если модераторы сочтут это целесообразным (я не думаю, что смогу сделать это сам?).   -  person Bajcz    schedule 26.08.2016
comment
Самый быстрый способ выполнить миграцию — скопировать текст уценки в текстовый редактор, удалить здесь вопрос, а затем повторно ввести его в резюме. Если вы это сделаете, мне потребуется немного времени, чтобы добавить начальный текст, чтобы мотивировать вопрос с точки зрения того, что вы пытаетесь сделать аналитически.   -  person lmo    schedule 26.08.2016
comment
Для нелинейного метода наименьших квадратов не существует матрицы шляпы в традиционном смысле, потому что подобранные значения нелинейно зависят от значений отклика, поэтому связь между ними не может быть охарактеризована матрицей. Вы, вероятно, можете вычислить локальную матрицу шляпы, вычисленную с использованием производных, которые были бы действительны в окрестности решения, но, по-видимому, R этого не делает.   -  person mrip    schedule 26.08.2016

Ответы (1)


arrow_upward
3
arrow_downward

Есть хорошая статья (о выбросах Обнаружение в нелинейной регрессии), где матрица шляпы аппроксимируется матрицей градиента, вычисленной в оценочной точке.

В твоем случае:

# gradient of the model function at the current parameter values
V <- nls1$m$gradient()

# tangent plane leverage matrix (it plays a similar role as the Hat matrix)
H <- V %*% solve(t(V) %*% V) %*% t(V)

# 'hat' values for nls
nls1.hat_values <- diag(H)

И если вы будете следовать этой статье можно вычислить H немного быстрее:

Q1 <- qr.Q(qr(V)) # V is the same matrix as above
H <- Q1 %*% t(Q1)

Поскольку H может быть довольно большим, и если вам нужны только значения шляпы, вы можете вообще пропустить умножение матриц. Нам нужна только диагональ матрицы H.

###
#' Approximation of hat values for nls.
#'
#' @param model An 'nls' object
#' @param ... Additional parameters (ignored)
#' @return Vector of approximated hat values
###
hatvalues.nls <- function(model, ...) {
  stopifnot(is(model, 'nls'))
  list(...) # ignore additional parameters
  V <- model$m$gradient()
  Q1 <- qr.Q(qr(V))
  rowSums(Q1*Q1)
}
person trope    schedule 20.01.2017