Извлечение конкретных коэффициентов взаимодействия из той же линейной модели в R

Меня интересуют конкретные эффекты (или параметры) линейной модели. Например, взяв набор данных радужной оболочки, мне было бы интересно:

  1. Как Sepal.Width модулирует эффект versicolor и virginica (по сравнению с эталонным уровнем setosa).
  2. Линейная зависимость между Sepal.Width и Sepal.Length на каждом из трех уровней species.

Единственный способ, который я нашел, это использовать две модели, одну с взаимодействием, а другую с вложенным взаимодействием:

fit <- lm(Sepal.Length ~ Species * Sepal.Width, data=iris)
summary(fit)  # The two last lines

fit <- lm(Sepal.Length ~ Species / Sepal.Width, data=iris)
summary(fit)  # The three last lines

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

Есть ли способ получить из одной модели эквивалент двух последних строк первой модели и трех последних строк второй модели? Спасибо


r interaction lm glm
person Dominique Makowski    schedule 24.06.2018    source источник
comment
Получили ли вы желаемый ответ, выполнив fit <- lm(Sepal.Length ~ Species * Sepal.Width + Species / Sepal.Width, data=irs) ?   -  person Corey Levinson    schedule 24.06.2018
comment
К сожалению, нет, это возвращает тот же результат, что и fit 1 (только с видами * sepal.width)   -  person Dominique Makowski    schedule 24.06.2018

Ответы (1)


arrow_upward
2
arrow_downward

Если вас интересуют оценки коэффициентов из сводки, то из одной модели путем сложения можно получить эквивалент двух последних строк первой модели и трех последних строк второй модели. Позволь мне объяснить.

На самом деле оба резюме дают вам одну и ту же информацию, только по-разному. Для первого резюме, вот мой вывод:

                              Estimate Std. Error t value Pr(>|t|)    
(Intercept)                     2.6390     0.5715   4.618 8.53e-06 ***
Speciesversicolor               0.9007     0.7988   1.128    0.261    
Speciesvirginica                1.2678     0.8162   1.553    0.123    
Sepal.Width                     0.6905     0.1657   4.166 5.31e-05 ***
Speciesversicolor:Sepal.Width   0.1746     0.2599   0.672    0.503    
Speciesvirginica:Sepal.Width    0.2110     0.2558   0.825    0.411

Это результат подгонки, когда setosa удаляется из модели. Это происходит автоматически в R, потому что, когда вы передаете фактор линейной модели, вам нужно только n-1 части для его описания, поскольку часть базового уровня (в данном случае setosa) полностью коллинеарна с другой частью. два фактора.

Во втором резюме подходит, это то, что я вижу.

                              Estimate Std. Error t value Pr(>|t|)    
(Intercept)                     2.6390     0.5715   4.618 8.53e-06 ***
Speciesversicolor               0.9007     0.7988   1.128    0.261    
Speciesvirginica                1.2678     0.8162   1.553    0.123    
Speciessetosa:Sepal.Width       0.6905     0.1657   4.166 5.31e-05 ***
Speciesversicolor:Sepal.Width   0.8651     0.2002   4.321 2.88e-05 ***
Speciesvirginica:Sepal.Width    0.9015     0.1948   4.628 8.16e-06 ***

В этой модели, несмотря на то, что она полностью коллинеарна, базовый уровень setosa сохраняется в модели. Обычно, если в модели оставить полностью коллинеарную переменную, R вернет ошибку, потому что будет жаловаться, что не может инвертировать матрицу, чтобы найти наилучшее соответствие для желаемой формулы (это потому, что матрица не полноценный).

Но с моделью R все в порядке, потому что она успешно ее вернула. Это означает, что матрица была обратимой, потому что, несмотря на наличие всех трех уровней фактора в терминах взаимодействия, она все же была матрицей полного ранга. Используя это знание, я подумал, что информация о том, что «отсутствовало» во второй модели, была просто запечена в первой модели.

И это именно то, что я видел. Посмотрите на следующие уравнения с коэффициентом из модели 1 в левой части уравнения и коэффициентом из модели 2 в правой части уравнения.

Sepal.Width = Speciessetosa:Sepal.Width
Sepal.Width + Speciesversicolo:Sepal.Width = Speciesversicolor:Sepal.Width
Sepal.Width + Speciesvirginica:Sepal.Width = Speciesvirginica:Sepal.Width

Поэтому, чтобы получить оценки коэффициентов из последних двух строк первой модели и последних трех строк второй модели, просто возьмите сводку первой модели, а затем выполните три уравнения I написал выше, чтобы получить «лишнюю» информацию, хранящуюся в последних трех строках.

person Corey Levinson    schedule 24.06.2018
comment
Спасибо за это ясное объяснение! Как вы думаете, можно ли сделать что-то подобное, чтобы получить значимость (значение p) этих реконструированных коэффициентов? Кроме того, я считаю, что это решение будет работать для частотной модели, но не для байесовской модели, где вместо точечной оценки вы получаете распределение. Простое сложение двух апостериорных распределений кажется сложным, но, возможно, я ошибаюсь. В любом случае спасибо, это уже помогает. - person Dominique Makowski; 24.06.2018
comment
Да, я пытался выяснить, как получить стандартное отклонение, t-значение и p-значение из других коэффициентов, но у меня возникли проблемы, поскольку не сразу очевидно, как их можно получить. На этот раз это не было простым умножением или сложением. На самом деле, я просто подумал. Стандартную ошибку и p-значения можно получить, если повторно инвертировать матрицу. Этот пост показывает, что если бы вы изменили матрицу дизайна, добавив, это могло бы выявить - person Corey Levinson; 24.06.2018
comment
Но пройти через все эти проблемы кажется слишком сложным, и не уверен, что это вообще будет работать правильно. В любом случае, по крайней мере, у вас есть коэффициент решения... - person Corey Levinson; 24.06.2018