Блог

У меня есть вопрос относительно относительной важности переменных в GLM, который содержит взаимодействие (непрерывный * фактор).

Я экспериментирую с подходом, основанным на разбиении объясненной вариации, аппроксимируемой через (псевдо) — R-квадрат. Но я не уверен, как это сделать (1) в GLM и (2) с моделью, которая содержит взаимодействие.

Для простоты я подготовил пример модели с Guassian GLM с одним взаимодействием (используя набор данных mtcars, см. Код в конце поста). Но на самом деле я заинтересован в применении метода к обобщенному пуассоновскому GLM, который может содержать несколько взаимодействий. Из тестовой модели возникает несколько вопросов:

1. Как правильно разделить R-квадрат? Я попытался создать раздел, но я не уверен, что это правильный путь.
2. r-квадрат каждого члена не суммируется с r-квадратом полной модели (даже близко). Это также происходит с моделью, которая не содержит взаимодействия. Помимо ошибок в разбиении r-квадрата (я все еще считаю себя новичком в статистике :P); может ли на это также повлиять коллинеарность? Коэффициенты инфляции дисперсии ниже 3 после масштабирования непрерывных предикторов (модель без масштабирования имеет самый высокий VIF = 5,7).

Любая помощь очень истощена!

  library(tidyverse) library(rsq) library(car)  data lt;- mtcars %gt;%  # scale reduces collinearity: without standardizing, the variance inflation factor for the factor is 5.7  mutate(disp = scale(disp)) data$am lt;- factor(data$am)  summary(data)  # test model, continuous response (miles per gallon), type of transmission (automatic/manual) as factor, displacement as continuous model lt;-  glm(mpg ~ am   disp   am:disp,  data = data,  family = gaussian(link = "identity")) drop1(model, test = "F")  # graph the data ggplot(data = data, aes(x = disp, y = mpg, col = am))   geom_jitter()   geom_smooth(method = "glm")  # Attempted partitioning (rsq_full lt;- rsq::rsq(model, adj = TRUE, type = "v"))  (rsq_int lt;- rsq_full - rsq::rsq(update(model, . ~ . - am:disp), adj = TRUE, type = "v"))  (rsq_factor lt;- rsq_full - rsq::rsq(update(model, . ~ . - am - am:disp), adj = TRUE, type = "v"))  (rsq_cont lt;- rsq_full - rsq::rsq(update(model, . ~ . - disp - am:disp), adj = TRUE, type = "v"))  c(rsq_full, rsq_int   rsq_factor   rsq_cont)  car::vif(model)   # A simpler model with no interaction model2 lt;- glm(mpg ~ am   disp, data = data, family = gaussian(link = "identity")) drop1(model2, test = "F")  (rsq_full2 lt;- rsq::rsq(model2, adj = TRUE, type = "v")) (rsq_factor2 lt;- rsq_full2 - rsq::rsq(update(model2, . ~ . - am), adj = TRUE, type = "v")) (rsq_cont2 lt;- rsq_full2 - rsq::rsq(update(model2, . ~ . - disp), adj = TRUE,type = "v"))  c(rsq_full2, rsq_factor2   rsq_cont2)  car::vif(model2)    

Дано:

1. y = A B A * B

Я бы сравнил его значение в квадрате R с его более простыми версиями:

2. y = A B
3. y = A
4. y = B

Если нет взаимодействия, я ожидаю

 r-squared(model1) = r-squared(model2)  

Это должно относиться к любой линейной модели. Это также должно быть полезно для сравнения основного эффекта предикторов, даже если существует взаимодействие. Я знаю, что это спорно, но если вы посмотрите на сценарий, представленный на рисунке ниже, предиктор A является информативным только в том случае, если учитывается предиктор B; и наоборот, предиктор B обладает некоторой прогностической силой даже сам по себе (y для B1 выше, чем y для B2, независимо от уровня A, к которому они принадлежат).

Вот пример с моделируемыми данными (чтобы избежать проблем коллинеарности и ненормальности):

 # simulate data: df lt;- data.frame(Species = as.factor(c(rep("Species A", 200),  rep("Species B", 200)  )),  Treatment = as.factor(rep(c("diet 1", "diet 2","diet 1", "diet 2"), each=100)),  body.weight = c(rnorm(n=100, 30, 5),  rnorm(n=100, 29.9, 5),  rnorm(n=100, 55, 5),  rnorm(n=100, 90, 5)  )  )  

 # Let's fit and compare the alternative models: lm.interactive lt;- lm(body.weight ~ Species * Treatment, data=df) lm.additive lt;- lm(body.weight ~ Species   Treatment, data=df) lm.only.species lt;- lm(body.weight ~ Species, data=df) lm.only.Treatment lt;- lm(body.weight ~ Treatment, data=df) lm.null lt;- lm(body.weight ~ 1, data=df)  # obtain R^2:   summary(lm.only.Treatment)$adj.r.squared # main effect of Treatment summary(lm.only.species)$adj.r.squared # main effect of species ID.  # As the figure suggests, it's larger than the main effect of Treatment  # (species identity affects body weight regardless of treatment) summary(lm.additive)$adj.r.squared # sum of the main effects summary(lm.interactive)$adj.r.squared # main effects   interaction  # fraction of variance explained by the interaction alone: summary(lm.interactive)$adj.r.squared - summary(lm.additive)$adj.r.squared  

Я не уверен, что мы действительно можем говорить о «доле дисперсии, объясняемой только взаимодействием», хотя. Возможно, более уместным будет говорить об увеличении объясненной дисперсии благодаря включению термина взаимодействия.

Я не уверен в том, насколько статистически обоснован подход, который я предлагаю, в его ограничениях или в том, надежно ли он работает для несбалансированных наборов данных. Одна из проблем этого подхода заключается в том, что различия в R-квадрате не могут быть проверены статистически, учитывая, что у нас есть только одно значение R-квадрата для каждой модели. Одним из способов обойти это было бы получить распределение значений в квадрате R для каждой модели с помощью начальной загрузки.