Трансплантация коэффициентов уравнения линейной регрессии от одной модели другой, в R Опции

[ИНО]

Задача в том, чтобы обучить модель линейной регрессии (см. тему соседнюю) на одних данных, а потом применить к другим, при этом не меняя коэффициентов. Извлечь коэффициенты из объекта класса "lm" легко, но как правильно передать другой модели? С некоторым удивлением не обнаружил в синтаксисе функции lm() возможности задачи пользовательских коэффициентов для каждого члена уравнения. На буржуйских формах что-то пишут о применении служебной функции offset(), но я так и не понял как ее запрограммировать в случае множественной регрессии со взаимодействиями категориальных переменных с числовыми. Поэтому возникла идея тупой "пересадки органов". Логика подсказывает, что помимо $coefficients требуется пересадить еще и $fitted.values (в моем случае все наблюдения во модели-реципиенте являются подмножеством наблюдений на которых построена модель-донор, так что можно просто выкинуть лишние значения). Если бы это было не так, можно было бы подогнать новые, применив к модели-донору функцию predict() с наблюдениями, используемыми в модели-реципиенте, в качестве аргумента newdata. Далее следует поменять остатки ($residuals). Их можно не пересаживать, а вычислить, вычтя пересаженные $fitted.values из родных для реципиента $model$имя_зависимой_переменной. Достаточно ли перечисленных операций для того, чтобы модель-реципиент стала вести себя так, будто была построена на априорно заданных коэффициентах в таких задачах как построение доверительных интервалов, вычисление F-статистики и R2? Или надо пересаживать что-то еще? Просто внутри объекта "lm", помимо вышеназванного, напихано много всякого, что недоступно моему разумению. Быть может, без редактирования чего-нибудь из этого в дополнение к проделанным трансплантациям органов, чье назначение мне ясно, организм донора будет функционировать не совcем правильно, выдавая вместо ожидаемых от него результатов погоду на Луне?

Сообщение отредактировал ИНО - 28.04.2025 - 19:05

[comisora]

Задача в том, чтобы обучить модель линейной регрессии (см. тему соседнюю) на одних данных, а потом применить к другим, при этом не меняя коэффициентов. Извлечь коэффициенты из объекта класса "lm" легко, но как правильно передать другой модели?

Я считаю, что Вам достаточно будед записать код вида

Код

predict(fit, newdata = newdf)

. Коэффициенты модели меняться не будут, они просто будут использоваться на новых данных.

https://rpubs.com/eR_ic/regression - инструменты по автоматизации сходной с Вашей (в моем прочтении) задачи.

[ИНО]

Но постойте, predict() не принимает новые данные по зависимой переменной, а лишь предсказывает их! Это будет предсказание модели-донора. А мне требуется получить предсказание модели-реципиента после пересадки уравнения. Да, регрессионная кривая в обоих случаях будет одинаковой, но вот ошибки и, соответственно, ДИ, F-статистика, R-квадрат и т. д. - совсем разными. Вот мне нужно это все иметь, основанное на остатках модели-реципиента от уравнения модели-донора. Да, это уже будет не подгонка по методу наименьших квадратов, а сравнение наблюдений с неким априори заданным нерушимым законом, но в том и суть задачи.

К сожалению статью по ссылке не смог прочесть - Гуль-переводчик не захотел этот сайт кушать. Завтра попробую Яндекс. Но таки, Вы, вероятно, правы в том, что примочки кросс-валидаторов потенциально могут помочь в моей задаче, хотя она и немного иного рода. Вопрос только в том, реализована ли там полноценная переработка модели, чтобы ее можно было потом использовать как обычную lm, в частности, скармливать anova(), экстрагировать R2, строить ДИ и т. д. Обычно ведь при проверке на новых данных достаточно получить одно число в качестве усредненного показателя ошибки, и на том успокоиться.

Еще вопрос, связанный с препарированием: как узнать внутреннее устройство функции predict() из конкретного пакета?

Сообщение отредактировал ИНО - 1.05.2025 - 12:50

[ИНО]

Эксперименты показали, что пациенты класса "lm" переносят трансплантацию хорошо и ведут себя после нее адекватно. Например, если заменить все остатки нулями, при вызове summary() получается бесконечная F-статистика, R2, равный единицы и нулевые стандартные ошибки коэффициентов, также имеют место выдаваемые функцией predict() доверительные интервалы нулевой ширины. То есть все так, как и ожидалось, без неприятных сюрпризов. Правда, среди органов-кандидатов для пересадки обнаружились еще парочка, а именно $effects и $qr, но я покамест не разобрался с их устройством, назначением, а главное - не знаю функций, которые их используют. Рад буду подсказкам на этот счет. Хотя, похоже, конкретно для моих задач уже пересаженного достаточно.

Увы, далеко не все классы регрессионных моделей в R имеют такое же удобное для трансплантации устройство как "lm". Например, организм объектов класса "lm_robust" из пакета estimatr устроен по совершенно иным принципам, и его так просто не взломаешь, видимо, препарировать придется саму функцию эти объекты порождающую.

[ИНО]

Разобрался с $qr - надо оставлять донорский, там инфа только о предикторах. Все, что касается соответствия независимой переменной модели, содержится в $fitted.values и $residuals. Статистики, появляющиеся при вызове summary(), anova() и др., вычисляются на лету. Это касается сугубо lm, в объектах других классов, даже похожих внешне, в нутре может быть все совсем по-другому. Осталось разобраться с $effects.

[Игорь]

Задача в том, чтобы обучить модель линейной регрессии (см. тему соседнюю) на одних данных

Это стадия обучения. Получится обученная модель.

а потом применить к другим, при этом не меняя коэффициентов.

Это стандартная стадия распознавания.

Задача в любых электронных таблицах решается на "ура".

В чем сложность? Предположу, сложность в выбранной реализации. Когда-то тоже думал, что было бы неплохо иллюстрировать алгоритмы программными текстами на чем-либо типа R в дополнение к теории с формулами. Оказалось всё гораздо печальнее. Книжный рынок заполнили издания по анализу данных, во многих из которых вообще нет формул, только вызовы R или Python. Нескольких лет хватило для выработки иммунитета: знакомство с книгами или статьями с упоминанием R или Python - пустая трата времени.

Сообщение отредактировал Игорь - 4.05.2025 - 13:25

[ИНО]

Да, проблема именно в конкретной реализации. Окзалось, что нельзя вот-так просто вставить в готовую функцию пользовательские коэффициенты, и заставить программу посчитать статистики по остаткам. Просто не предусмотрены такие аргументы. В случае простой линейной регрессии, конечно, несложно свой код написать, или даже в электронной таблице посчитать, но вот если предикторов много, уже сложнее. По крайней мере для меня. Да еще и хочется, чтобы все другие функции для работы с объектами класса "lm", работали коорректно.

ИМХО в случае с lm трансплантация - рабочий лайфхак. Даже если у нас нет донора, можно вырастить органы in vitro! Вот алгоритм:
1. Строим модель обычным образом.
2. Редактируем в полученном объекте $coefficients на свой вкус, не трогая имен перменных.
3. Делаем predict(отредактированная_модель) и заменяем его выдачей значения $fitted.values.
4. Заменяем значения $residuals на разность $model$имя_зависимой_переменной - $fitted.values (замененные ранее).

Вроде бы, это все, что требуется для корректной работы summary() на перепрошитой модели. Осталось разобраться, то за зверь $effects, и кто к нему обращается при вызове.

Но это все касалось функции lm() и ее порождений.

А вот в lm_robust из пакета estimatr все иначе, я даже не смог найти, где там хранится инфа об остатках либо о значениях зависимой переменной. Кажется, вообще не хранится, зато хранится готовый список всех возможных статистик. Такой объект вышеописанным способом не хакнешь.

[comisora]

Окзалось, что нельзя вот-так просто вставить в готовую функцию пользовательские коэффициенты, и заставить программу посчитать статистики по остаткам.

Я предположил, что функция nls для Вас может подойти лучше, так как присутствует ручная установка коэффициентов. Попробовал смастерить (правда для lm) эквивалентные реализации, можете ознакомиться.

CODE

testdata <- data.frame(
group = factor(
c(rep("A", 7), rep("B", 6))
),
time = c(
15, 17, 18, 19, 20, 21, 22,
15, 18, 19, 20, 21, 22
),
response=c(
1.078, 0.949, 0.793, -2.201, -1.181, -1.95, -2.391,
2.628, 2.525, 2.852, 2.811, 3.484, 2.206
)
)

fit <- lm(response~I(time^2):group,
data = testdata)

fit1 <- fit
fit1$coefficients <- c(1, 2, 3)
fit1$fitted.values <- predict(fit1)
fit1$residuals <- fit$fitted.values-fit1$fitted.values

fit2 <- lm(
response ~ -1 +
offset(I(rep(1, 13))) +
offset(
I(
ifelse(group == 'A', I(time^2)*2, 0)
)
) +
offset(
I(
ifelse(group == 'B', I(time^2)*3, 0)
)
)
,
data = testdata
)

fit3 <- nls(
response ~ c +
ifelse(group == 'A', I(time^2)*a, 0) +
ifelse(group == 'B', I(time^2)*b, 0),
start = list(c = 1, a = 2, b = 3),
lower = list(b = 1, a = 2, b = 3),
upper = list(b = 1, a = 2, b = 3),
data = testdata, algorithm = 'port'
)

rsq <- function(actual, predicted) {
ss_total <- sum((actual - mean(actual)) ^ 2)
ss_residual <- sum((actual - predicted) ^ 2)
r_squared <- 1 - (ss_residual / ss_total)

return(r_squared)
}

Также бегло посмотрел stackoverflow по схожим вопросам. все упирается в умение решать задачи по (не)линейному программированию, самому писать функции и т.д., что неудивительно.

[ИНО]

Спасибо, но что-то Ваш самопальный R2 на данном синтетическом примере выдает нечто удивительное, а именно -215031.1. Я даже не представляю, какой статистический смысл может нести это странное число. Хотя в нормальных ситуациях, где подгонка модели осуществлялась по методу наименьших квадратов, вроде все фурычет, как надо, и сама формула соответствует прописям. Выходит что ж R2, в случае с моделями с кожффициентами, взятыми с потолка, теряет всякий смысл? А как же квадрат корреляции? Вот альтернативная формула, которая гарантирует иyтервал [0; 1]:

Код

R2=cor(fit1$model$response, fit1$fitted.values)^2

0.1524337

Но ни Ваш R2, ни мой, не соответствуют

Код

summary(fit1)$r.squared

0.08443

Вот-это уже совсем непонятно, поскольку данное несоответствие наблюдается исключительно в "перепрошитых" моделях, в обычных все три цифры сходятся. Сначала я грешил на так и не понятый мною компонент $effects, который оставлял неизменным, однако даже полное его удаление никак не влияет на работу функции summary(). А все прочее, что относилось к переменной отклика, я отредактировал, как мне кажется, вполне корректно. Проверить соответствие штатного R2 c самопальными для модели fit2 оказалось невозможно - summary() с ней не работает толком, поскольку коэффициентов внутри нет, и предсказания с аргументом newdata она делать не умеет. Т. е. это не полноценная регрессионная модель, а выкидыш какой-то. Так что этот хитровыделанный offset() - не для меня. А для объектов класса "nls" свой штатный R2 не предусмотрен в принципе, как и много чего еще. Так что кроме насилия над объектами "lm" выхода по-прежнему не вижу.

Еще пара слов об $effects: оказывается, без него не работает anova(). Причем суммы квадратов в таблице определяются именно им, хотя на значения F-статистики из той же таблицы влияет и что-то из мною замененного. В итоге получается некая гибридная таблица, корректность которой весьма сомнительна. Так что при трансплантации коэффициентов $effects по-хорошему тоже надо редактировать, но я понятия не имею как, поскольку даже не представляю что оно такое, и откуда берется.

Сообщение отредактировал ИНО - 10.05.2025 - 23:01

[ИНО]

В общем написал я и функцию для F-стактистики, пос кольку на моих "перепрошитых" линейных моделях, summary() привирает в третьем знаке после запятой, точно так же как и в случае с R2. Почему и по каким данным, хранящемся в объекте, она вообще их считает, - тайна велика есть, учитывая что я заменил и $fitted.values и $residuals, а $effects вообще выкинул. И похоже никто не ответит здесь на этот вопрос
Но вот другой: почему даже в девственно чистой модели $fitted.values и выдача predict() без аргумента newdata не совпадают полностью? Расхождение в х. з. каком знаке, но оно достаточно для того, чтобы model$fitted.values==predict(model) время от времени выдавало FALSE. Видимо, имеет место отличие в точности вычислений, но непонятно, где она выше, а где ниже. Не то чтобы, это сколь-нибудь существенно влияло на практические задачи, но перфекционист внутри меня негодует