параметр весов в GAM

Question

В моем наборе данных много избыточных наблюдений (но каждое наблюдение должно учитываться). Поэтому я рассматриваю использование опции «весов» в GAM, потому что это значительно сокращает время вычисления.

gam функция (в mgcv пакете) объясняет, что они являются «эквивалент» (от ?gam аргументам weights):

«Обратите внимание, что вес 2, например, эквивалентно тому, что сделал именно такое же наблюдение дважды ".

Но это не кажется правильным.

yy = c(5,2,8,9) 
xx = 1:4 
wgts = c(3,2,4,1) 
yy2 = rep(yy, wgts) 
xx2 = rep(xx, wgts) 
mod1 = gam(yy2 ~ xx2) 
mod2 = gam(yy ~ xx, weights = wgts) 
mod3 = gam(yy ~ xx, weights = wgts/mean(wgts)) 

predict(mod1,data.frame(xx2=1:4)) 
predict(mod2,data.frame(xx=1:4)) 
predict(mod3,data.frame(xx=1:4)) 

Оценки идентичны во всех трех моделях. Стандартная ошибка одинакова для моделей 2 и 3, но различна в модели 1. GCV отличается во всех трех моделях.

Я понимаю, что GCV могут быть разными. Но как мы можем сказать, что модели идентичны, если стандартные ошибки различны? Это ошибка, или есть ли хорошее объяснение этому?

Answer 1

Проблемы, которые вы видели, не касаются GAM. Вы использовали gam, чтобы соответствовать параметрической модели, и в этом случае gam ведет себя почти так же, как lm. Чтобы ответить на ваши вопросы, достаточно сосредоточиться на случае линейной регрессии. Что произойдет с линейной моделью, произойдет и с GLM и GAM. Вот как мы можем воспроизвести проблему с lm:

yy <- c(5,2,8,9) 
xx <- 1:4 
wgts <- c(3,2,4,1) 
yy2 <- rep(yy,wgts) 
xx2 <- rep(xx,wgts) 
fit1 <- lm(yy2 ~ xx2) 
fit2 <- lm(yy ~ xx, weights = wgts) 
fit3 <- lm(yy ~ xx, weights = wgts/mean(wgts)) 
summary1 <- summary(fit1) 
summary2 <- summary(fit2) 
summary3 <- summary(fit3) 
pred1 <- predict(fit1, list(xx2 = xx), interval = "confidence", se.fit = TRUE) 
pred2 <- predict(fit2, list(xx = xx), interval = "confidence", se.fit = TRUE) 
pred3 <- predict(fit3, list(xx = xx), interval = "confidence", se.fit = TRUE) 

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

1. Для взвешенной регрессии fit2 и fit3, почему почти все то же самое, за исключением остаточной стандартной ошибки?
2. Почему взвешенная регрессия (fit2 или fit3) не эквивалентна обычной регрессии со связями?

---

Ваш первый вопрос о масштабной инвариантности весовых наименьших квадратов с весами. Вот краткое резюме я сделал:

Если мы отмасштабировать W произвольным положительным значением, только остаточная ошибка стандартной и немасштабированная ковариация изменится. Такое изменение не подразумевает другую, неэквивалентную модель. На самом деле все, что связано с предсказанием, не затрагивается. При взвешенной регрессии не просто смотрите на sigma2; это лишь незначительная дисперсия. Что действительно представляет интерес, так это валовая разница после умножения веса. Если вы разделите свои веса на 2, вы найдете sigma2 удвоений, но вы все равно получите тот же результат при их умножении.

summary2$coef 
summary3$coef 

#   Estimate Std. Error t value Pr(>|t|) 
#(Intercept) 2.128713 3.128697 0.6803832 0.5664609 
#xx   1.683168 1.246503 1.3503125 0.3094222 

pred2 
pred3 

#$fit 
#  fit  lwr  upr 
#1 3.811881 -5.0008685 12.62463 
#2 5.495050 -0.1299942 11.12009 
#3 7.178218 0.6095820 13.74685 
#4 8.861386 -1.7302209 19.45299 
# 
#$se.fit 
#  1  2  3  4 
#2.048213 1.307343 1.526648 2.461646 
# 
#$df 
#[1] 2 
# 
#$residual.scale ## for `pred2` 
#[1] 3.961448 
# 
#$residual.scale ## for `pred3` 
#[1] 2.50544 

---

Ваш второй вопрос о значение веса. Массы используются для моделирования гетеросекастической реакции для преодоления эффекта плеча при обычной наименьшей квадратичной регрессии. Массы пропорциональны взаимной дисперсии: вы даете больший вес данным с меньшими ожидаемыми ошибками. Весы могут быть нецелыми, поэтому он не имеет естественного объяснения с точки зрения повторяющихся данных. Таким образом, то, что написано в пакете mgcv, не совсем корректно.

Реальная разница между fit1 и fit2? это степень свободы. Проверьте приведенную выше таблицу на (n - p). n - это количество данных, которое у вас есть, а p - количество коэффициентов не NA, поэтому n - p - это остаточная степень свободы. Для обеих моделей есть p = 2 (перехват и уклон), но для fit1 у нас есть n = 10, а для fit2 у нас есть n = 4. Это оказывает существенное влияние на вывод, поскольку теперь стандартные ошибки для коэффициентов и прогнозов (следовательно, доверительные интервалы) будут отличаться. Эти две модели далеко не эквивалентны.

summary1$coef 
#   Estimate Std. Error t value Pr(>|t|) 
#(Intercept) 2.128713 1.5643486 1.360766 0.21068210 
#xx2   1.683168 0.6232514 2.700625 0.02704784 

summary2$coef 

#   Estimate Std. Error t value Pr(>|t|) 
#(Intercept) 2.128713 3.128697 0.6803832 0.5664609 
#xx   1.683168 1.246503 1.3503125 0.3094222 

pred1 

#$fit 
#  fit  lwr  upr 
#1 3.811881 1.450287 6.173475 
#2 5.495050 3.987680 7.002419 
#3 7.178218 5.417990 8.938446 
#4 8.861386 6.023103 11.699669 
# 
#$se.fit 
#  1   2   3   4 
#1.0241066 0.6536716 0.7633240 1.2308229 
# 
#$df # note, this is `10 - 2 = 8` 
#[1] 8 
# 
#$residual.scale 
#[1] 1.980724 

pred2 

#$fit 
#  fit  lwr  upr 
#1 3.811881 -5.0008685 12.62463 
#2 5.495050 -0.1299942 11.12009 
#3 7.178218 0.6095820 13.74685 
#4 8.861386 -1.7302209 19.45299 
# 
#$se.fit 
#  1  2  3  4 
#2.048213 1.307343 1.526648 2.461646 
# 
#$df # note, this is `4 - 2 = 2` 
#[1] 2 
# 
#$residual.scale ## for `pred2` 
#[1] 3.961448