Одна из полезных функций пакета mgcv в R - возможность подгонять отдельные сглаживания для разных уровней фактора с помощью аргумента s(,by). После этого часто возникает следующий вопрос: нужен ли термин «взаимодействие» или достаточно простой аддитивной модели с одним сглаживанием для всех уровней. Однако мне до сих пор неясно, как лучше всего это проверить, особенно при использовании расширенных семейств с автоматической c оценкой параметра дисперсии (например, family = nb() или family = tw()), и иногда получаю странные результаты. Надеюсь, это простое представление, основанное на наборе данных, который я недавно проанализировал, поможет проиллюстрировать, где я запутался.
# Construct a simple dataset
library(mgcv);library(ggplot2);library(dplyr)
data = data.frame(
N =c(10,9,9,14,0,0,0,0,1,6,0,0,0,4,1,0,0,5,2,0,0,0,11,6,6,1,0,2,1,0,0,2,0,0,0,1,0,0,0,9,0,0,0,5),
d =c(0.7,0,0,0,9.1,5,9.4,18.5,36.5,1,3.7,9,18,0.8,1.9,8.9,18.4,0.1,3.5,8.9,18.7,38.5,1.1,0.2,0.1,0,17.2,1.5,3.9,13.7,33.7,0.2,3.3,7.9,19,0.8,2.9,8.3,18,0.9,3.2,5.5,17,0.9),
loc =as.factor(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2))
)
Данные состоят из количества животных (N) как функции расстояния (d) из двух мест ( ло c). График показывает, что численность резко возрастает на меньших расстояниях.
ggplot(data,aes(x=d,y=N,colour=loc))+geom_point()
Данные выглядят немного странно, но симуляции показывают, что, хотя модель Пуассона сильно диспергирована, GAM соответствуют отрицательной биномиальной ошибке распространение имеют приемлемые остаточные свойства. Первоначально я хотел бы проверить гипотезу о том, что численность меняется по-разному с расстоянием от каждого места, поэтому я подбираю модель с сглаженным по факторам членом. Обычно я бы использовал метод «двойного штрафа» с select=T для эффективного удаления неважных терминов:
m = gam(N ~ s(d,by=loc) + loc,family = nb(),data=data,method='REML',select=T)
summary(m)
Family: Negative Binomial(3.988)
Link function: log
Formula:
N ~ s(d, by = loc) + loc
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.9677 0.7263 -1.332 0.183
loc2 -2.9670 1.6949 -1.751 0.080 .
---
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(d):loc1 2.1981 9 24.23 3.19e-06 ***
s(d):loc2 0.9286 9 11.08 0.00055 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.684 Deviance explained = 82.4%
-REML = 60.086 Scale est. = 1 n = 44
Однако, хотя это говорит мне, что N значительно увеличивается с d, он не сообщает мне, улучшает ли соответствие модели включение другого сглаживания для каждого уровня loc.
Тест отношения правдоподобия
Одним из вариантов было бы использовать вероятность Ratio Test или AI C для сравнения вложенных моделей с термином 'by=' и без него:
m1 = update(m,,method='ML',select=F)
m2 = update(m1,~.-s(d,by=loc)+s(d))
anova(m1,m2,test='Chisq')
Analysis of Deviance Table
Model 1: N ~ s(d, by = loc) + loc
Model 2: N ~ loc + s(d)
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 37.866 27.113
2 38.675 27.423 -0.80826 -0.30984 0.492
Это указывает на то, что термин «взаимодействие» не добавляет много, а AI C в целом соглашается
AIC(m1,m2)
df AIC
m1 6.702825 118.0077
m2 5.841031 116.2276
Однако, поскольку я подбираю модели с использованием функции nb() (тета неизвестна), они имеют разные тета-параметры в их распределениях ошибок. Следовательно, можно ли их считать действительно вложенными моделями для целей теста отношения правдоподобия? p-значение, потому что более простая модель также имеет более низкое остаточное отклонение.
m3 = update(m2,,family = negbin(theta = m1$family$getTheta(TRUE)))
anova(m1,m3,test='Chisq')
Analysis of Deviance Table
Model 1: N ~ s(d, by = loc) + loc
Model 2: N ~ loc + s(d)
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 37.866 27.113
2 38.699 27.088 -0.83225 0.02493
Этот пример - всего лишь один вид, взятый из гораздо большего набора данных. В других случаях верно обратное: сравнение m1 и m2 (оба подходят с использованием функции nb()) с использованием теста отношения правдоподобия приводит к неоцениваемому значению p, которое можно исправить, переустановив m2 на параметр тета из m1. Мы можем показать это, подставив счетчики другого вида и снова запустив предыдущий код.
data$y = c(22,2,11,14,0,0,0,0,0,0,0,0,0,0,0,0,0,6,0,0,0,0,3,8,3,4,0,0,0,0,0,0,0,0,0,0,0,0,0,4,0,0,0,7)
Упорядоченные факторы
Еще один подход - преобразовать loc к заказанному коэффициенту и проверьте, отличается ли сглаживание для loc = 2 от контрольного уровня loc = 1. Однако это дает совершенно другую интерпретацию (т.е. сглаживание для местоположения 2 значительно отличается от местоположения 1):
m4 = update(m,data = mutate(data,loc = factor(loc,ordered=T)))
summary(m4)
Family: Negative Binomial(0.439)
Link function: log
Formula:
N ~ s(d, by = loc) + loc
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.527 1.104 -1.383 0.1667
loc.L -3.773 1.561 -2.417 0.0157 *
---
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(d):loc2 0.8777 9 5.367 0.0134 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.106 Deviance explained = 35.5%
-REML = 75.897 Scale est. = 1 n = 44
Хотя я отмечаю, что объясненное отклонение в этом случае намного ниже, чем в оригинальная модель ...
ВОПРОСЫ
Итак, мои вопросы к списку тройные:
- Что люди считают наиболее подходящий способ оценки значимости фактора с помощью переменной сглаживания?
- Можно ли надежно сравнить две модели с использованием
nb() (или других автоматизированных семейств) с разными параметрами дисперсии, используя LRT? - Почему в последнем случае более сложная модель может также иметь более высокое остаточное отклонение и как вы интерпретируете результат LRT? Можем ли мы просто сделать вывод, что более простая модель предпочтительна с p = 1, потому что она также имеет меньшее отклонение?
Любые подсказки, полученные с благодарностью!