Подгонка многих линейных моделей в R с идентичными матрицами проектирования [дубликаты]

Для приложения нейровизуализации я пытаюсь подогнать многие линейные модели методом наименьших квадратов в R (стандартный вызов lm). Представьте, что у меня есть матрица дизайна X. Эта матрица дизайна будет одинаковой для всех моделей. Подходящие данные (Y) изменятся, и в результате изменятся все подходящие параметры (например, бета, p-значения, остатки и т. д.).

В настоящее время я просто втыкаю его в цикл for, поэтому он выполняет сотни тысяч вызовов lm. Кажется, должен быть лучший способ.

Я считаю, что самая затратная в вычислительном отношении часть — это обращение матрицы. Похоже, это обрабатывается вызовом Fortran в lm.fit.

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

Есть ли способ получить мой торт и съесть его тоже? А именно, получить удобство lm, но использовать эту мощь для вычислительной эффективности множества моделей с идентичными матрицами проектирования?


r parallel-processing regression mlm
person Daniel Kessler    schedule 28.01.2013    source источник
comment
@gung Да, на этот вопрос, я думаю, вы правы. Я пойду вперед и флаг. Спасибо за вашу помощь!   -  person Daniel Kessler    schedule 29.01.2013
comment
Какой бы алгоритм вы ни использовали, после вычисления QR-разложения матрицы модели его не следует пересчитывать для нового y. Единственная часть вычислений, которую следует повторить, — это часть, зависящая от y. Таким образом, ответы, которые избегают этого пересчета, при прочих равных условиях гораздо более эффективны, чем те, которые этого не делают. Для любого предложения, которое не предлагает избегать его, будет эквивалентный ответ, который действительно избегает его, что должно быть быстрее. Если вы хотите выполнять регрессии отдельно, вы можете повторно использовать уже вычисленные ранее части; есть несколько способов сделать это   -  person Glen_b    schedule 29.01.2013

Ответы (3)


arrow_upward
8
arrow_downward

На странице справки для lm:

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

Таким образом, может показаться, что простым подходом будет объединение всех различных векторов y в матрицу и передача ее в качестве ответа в одном вызове lm. Например:

(fit <- lm( cbind(Sepal.Width,Sepal.Length) ~ Petal.Width+Petal.Length+Species, data=iris))
summary(fit)
summary(fit)[2]
coef(summary(fit)[2])
coef(summary(fit))[2]
sapply( summary(fit), function(x) x$r.squared )
person Greg Snow    schedule 28.01.2013
comment
Кажется, это действительно эффективно вычисляет, ура! Однако для действительно больших моделей (ncol(response) ~= 600k) вычисление сводки занимает очень много времени. Мне нужно получить t-значения и p-значения, соответствующие бета-коэффициентам, и единственный известный мне способ сделать это — с помощью сводки. Любые идеи о более эффективном способе их извлечения, или я застрял, переписывая резюме, чтобы быть более эффективным? - person Daniel Kessler; 30.01.2013
comment
Вы можете посмотреть на такие функции, как estVar или vcov, они могут быть быстрее (затем рассчитайте статистику теста на основе этой информации). - person Greg Snow; 30.01.2013

arrow_upward
11
arrow_downward

Да, есть лучший способ. Мы писали примеры функций замены fastLm() на основе использования внешнего кода C/C++ от Armadillo, GSL и Eigen в пакетах RcppArmadillo, RcppGSL и RcppEigen.

На сегодняшний день наибольшее количество времени тратится на настройку матрицы модели и разбор формулы. Вы можете прочитать исходный код lm() или, может быть, наш в fastLm() и посмотреть, как сделать этот синтаксический анализ всего один раз. Сохраняйте правую сторону, а затем перебирайте различные y вектора. Какая функция подгонки вы используете, имеет меньшее значение. Мне нравится fastLm() от RcppArmadillo, но я тоже написал его :)

person Dirk Eddelbuettel    schedule 28.01.2013
comment
Недавно я наткнулся на fastLm, и да, это действительно быстро!! Я пытаюсь заменить вызов lm.wfit() на fastLm(). Поэтому я называю это так: fastLm(X * wts, y * wts, method = 0L). Коэффициенты совпали с коэффициентами из lm.wfit(), но подогнанные значения — нет. У вас есть идеи, что я мог упустить? . Я не могу задать вопрос из-за плохой работы на SO, поэтому задаю его здесь. Сердечно извиняюсь! - person joel.wilson; 14.10.2016
comment
Ошибка мышления с вашей стороны: используйте sqrt(w), а не w. - person Dirk Eddelbuettel; 14.10.2016
comment
Я включил это... на самом деле wts ‹- sqrt(w); Извините, я не упомянул об этом. Однако я получил решение. на самом деле подогнанные значения из fastLm также необходимо разделить на wts fastLm$fitted.values <- fastLm$fitted.values/wts fastLm$residuals <- fastLm$residuals/wts - person joel.wilson; 14.10.2016
comment
Я не уверен, что вы пытаетесь сделать: ни одна из трех реализаций fastLm() не претендует на работу с весами... - person Dirk Eddelbuettel; 14.10.2016

arrow_upward
5
arrow_downward

Я не знаю лучшего способа использования lm; но вы можете рассмотреть возможность использования функции lsfit. Хотя синтаксис lsfit(X,y) проще и с меньшим количеством наворотов, он позволяет y быть не только вектором со значениями переменной ответа, но и матрицей. Затем один вызов lsfit соответствует всем столбцам y путем их регрессии в той же матрице проекта X. Довольно быстро и удобно.

person F. Tusell    schedule 28.01.2013