Значение фиксированных эффектов в отрицательной биномиальной смешанной модели с использованием MGCV гам - PullRequest
Новая
       105

Значение фиксированных эффектов в отрицательной биномиальной смешанной модели с использованием MGCV гам

0 голосов
/ 19 октября 2018

Я использую gam из пакета mgcv для анализа набора данных с 24 записями: 

ran  f1     f2   y
1   3000    5   545
1   3000    10  1045
1   10000   5   536
1   10000   10  770
2   3000    5   842
2   3000    10  2042
2   10000   5   615
2   10000   10  1361
3   3000    5   328
3   3000    10  1028
3   10000   5   262
3   10000   10  722
4   3000    5   349
4   3000    10  665
4   10000   5   255
4   10000   10  470
5   3000    5   680
5   3000    10  1510
5   10000   5   499
5   10000   10  1422
6   3000    5   628
6   3000    10  2062
6   10000   5   499
6   10000   10  2158

Данные имеют два фиксированных эффекта (f1 и f2) и один случайныйэффект (ran).Зависимые данные y.Поскольку зависимые данные y представляют счет и имеют избыточную дисперсию, я использую отрицательную биномиальную модель.

Модель gam и ее вывод summary выглядят следующим образом: 

library(mgcv)
summary(gam(y ~ f1 * f2 + s(ran, bs = "re"), data = df2, family = nb, method = "REML"))

Family: Negative Binomial(27.376) 
Link function: log 

Formula:
y ~ f1 * f2 + s(ran, bs = "re")

Parametric coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  5.500e+00  3.137e-01  17.533  < 2e-16 ***
f1          -3.421e-05  3.619e-05  -0.945    0.345    
f2           1.760e-01  3.355e-02   5.247 1.55e-07 ***
f1:f2        2.665e-07  4.554e-06   0.059    0.953    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
         edf Ref.df Chi.sq p-value    
s(ran) 4.726      5  85.66  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.866   Deviance explained = 93.6%
-REML = 185.96  Scale est. = 1         n = 24

Тест Вальда из summary дает очень высокое значение для f2 (P = 1.55e-07),Однако, когда я проверяю значимость f2, сравнивая две разные модели с использованием anova, я получаю совершенно разные результаты: 

anova(gam(y ~ f1 * f2 + s(ran, bs = "re"), data = df2, family = nb, method = "ML"),
gam(y ~ f1 + s(ran, bs = "re"), data = df2, family = nb, method = "ML"),
test="Chisq")

Analysis of Deviance Table

Model 1: y ~ f1 * f2 + s(ran, bs = "re")
Model 2: y ~ f1 + s(ran, bs = "re")
  Resid. Df Resid. Dev      Df Deviance Pr(>Chi)
1    14.843     18.340                          
2    16.652     21.529 -1.8091   -3.188   0.1752

f2 больше не имеет значения.Модели были изменены с REML на ML, как рекомендовано для оценки фиксированных эффектов.

Если взаимодействие сохраняется, f2 все еще остается незначительным при использовании anova: 

anova(gam(y ~ f1 + f2 + f1:f2 + s(ran, bs = "re"), data = df2, family = nb, method = "ML"),
gam(y ~ f1 + f1:f2 + s(ran, bs = "re"), data = df2, family = nb, method = "ML"),
test="Chisq")
Analysis of Deviance Table

Model 1: y ~ f1 + f2 + f1:f2 + s(ran, bs = "re")
Model 2: y ~ f1 + f1:f2 + s(ran, bs = "re")
  Resid. Df Resid. Dev       Df Deviance Pr(>Chi)
1    14.843     18.340                           
2    15.645     19.194 -0.80159 -0.85391   0.2855

Я был бы очень признателен за совет, какой из этих подходов является более подходящим.Большое спасибо!

1 Ответ

0 голосов
/ 23 октября 2018

В разделе ПРЕДУПРЕЖДЕНИЕ ?anova.gam говорится:

Если модели a и b отличаются только в терминах без компонентов без штрафа (таких как случайные эффекты), тогда значения p отanova(a,b) ненадежны и, как правило, слишком низки.

Я бы предположил, что значение p ненадежно, но в этом случае у вас есть ситуация, когда наблюдается случай, противоположный ожидаемому -Значения р намного больше.

Однако, я думаю, вы не сравниваете правильные модели.Если вы не знаете, что делаете, при сравнении моделей с взаимодействиями следует соблюдать принцип маржинальности.

Поэтому я бы сравнил модель с основными эффектами f1 и f2 с моделью, включающей этиОсновной эффект и их взаимодействие.

  • модель 1: y ~ f1 * f2 + s(ran, bs = "re")
  • модель 2: y ~ f1 + f2 + s(ran, bs = "re")

Разве чтоЕсли вы не рассказываете о том, как настроены ваши модели, вы не должны включать термин более высокого порядка, не включая термины более низкого порядка, которые встречаются в терминах более высокого порядка.Например, у вас есть f1 + f1:f2 и f2 найдено во втором члене заказа, но не найдено в качестве первого порядка в модели.

...