Обратите внимание, что для этой модели может иметь смысл использовать ответ Tweedie через семейство tw с gam() и bam(), которое нельзя использовать с gamm(). Фактически Саймон Вуд и Маттео Фазиоло сопоставляют эти данные с масштабом местоположения Tweedie GAM (в котором они моделируют средние, дисперсионные и энергетические параметры распределения Твиди, каждый из которых имеет отдельный линейный предиктор [модель]).
По предложению @ BenBolker: я бы даже не стал проверять случайный эффект в этой модели, и часто мне все равно, значительный он или нет. Это зависит от вопроса или гипотезы, над которой я работаю. Часто это требуется в модели из-за некоторой кластеризации данных, которые я хочу включить в модель, независимо от их значимости.
Однако я не уверен, что теория (Обобщенного) теста отношения правдоподобия (GLRT) не применима к использованию квази-правдоподобия в данном случае . Саймон Вуд представляет выводы в Приложении A 2-го издания своего учебника по GAMS, которые показывают, что ранее полученные результаты для оценки максимального правдоподобия (которые включают в себя результаты для GLRT) имеют место, если мы заменим логарифмическую правдоподобие логарифмической квази-правдоподобием. Это, по крайней мере, кажется, что Саймон спорит, предполагает, что интерпретация теста, который я упоминаю ниже и который реализован в summary.gam() для случайных эффектов, настолько надежна, как если бы она основывалась на правильной вероятности.
Если мне действительно не нужно, я подгоню эту модель к gam() или bam(), а затем gamm4() (последний из пакета gamm4 ), до gamm(), особенно для негауссовы модели, так как функция gamm() должна соответствовать этой модели как модели смешанных эффектов с использованием штрафованного квази-правдоподобия, что не обязательно является лучшим способом оценки этих моделей.
library(mgcv)
library(gamair)
devtools::install_github('gavinsimpson/gratia')
library(gratia)
data(sole)
sole$off <- log(sole$a.1-sole$a.0)
sole$a<-(sole$a.1+sole$a.0)/2
solr <- sole
solr$t <- solr$t-mean(sole$t)
solr$t <- solr$t/var(sole$t)^0.5
solr$la <- solr$la-mean(sole$la)
solr$lo <- solr$lo-mean(sole$lo)
solr$station <- factor(with(solr,paste(-la,-lo,-t,sep="")))
som <- gam(eggs ~ te(lo, la, t, bs = c('tp','tp'), k = c(25, 5), d = c(2,1)) +
s(t, k = 5, by = a) + s(station, bs = 're') + offset(off),
family = quasipoisson, data = solr, method = 'REML')
Затем summary(som) дает тест, основанный на тесте отношения правдоподобия, как предложено @BenBolker, но эталонное распределение корректируется для тестирования на границе пространства параметров.
> summary(som)
Family: quasipoisson
Link function: log
Formula:
eggs ~ te(lo, la, t, bs = c("tp", "tp"), k = c(25, 5), d = c(2,
1)) + s(t, k = 5, by = a) + s(station, bs = "re") + offset(off)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.4016 0.3061 -11.11 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
te(lo,la,t) 56.025 65.456 2.547 4.62e-10 ***
s(t):a 4.535 4.886 54.790 < 2e-16 ***
s(station) 128.563 388.000 1.175 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.833 Deviance explained = 88%
-REML = -7.9014 Scale est. = 0.58148 n = 1575
У меня возникли проблемы при сближении модели без случайного эффекта с использованием gamm(), поэтому я не смог проверить член случайного эффекта и даже столкнулся с ошибкой при попытке использовать мультимодальную форму anova().
Если вы хотите получить случайные эффекты, используя модель gam(), вы можете использовать мой пакет gratia (возможно, через несколько дней на CRAN, но который можно установить из github, как показано выше) и затем:
> evaluate_smooth(som, 's(station)')
# A tibble: 394 x 5
smooth by_variable station est se
<chr> <fct> <chr> <dbl> <dbl>
1 s(statio… NA -0.0004304761904734280.419685714285714-… -0.0396 2.55
2 s(statio… NA -0.0004304761904734280.6586857142857140… 1.48 1.20
3 s(statio… NA -0.0004304761904734281.15968571428571-1… -0.00606 2.63
4 s(statio… NA -0.0004304761904734281.176685714285710.… -0.0767 2.48
5 s(statio… NA -0.002430476190475870.9096857142857141.… -0.00654 2.63
6 s(statio… NA -0.01243047619047390.4106857142857140.0… -0.802 1.61
7 s(statio… NA -0.0154304761904740.631685714285714-0.4… -0.138 2.35
8 s(statio… NA -0.02043047619047660.375685714285714-0.… -0.426 1.94
9 s(statio… NA -0.02543047619047911.14668571428571-0.4… -0.0333 2.57
10 s(statio… NA -0.02743047619047450.875685714285714-0.… -0.0673 2.49
# … with 384 more rows
и вам нужен столбец est.
Смещение - это термин в модели, который имеет фиксированный эффект, равный 1. В этом случае он используется для стандартизации ответа на счетчик, так что каждый из них сравнивается как для; в этом случае он используется для интеграции возрастов яиц, обнаруженных в этом образце. Читать стр. 143 из 2-го издания книги Симона о GAM, чтобы узнать больше о том, что делается для этой модели и что означает смещение.
В общем, скажем, вы пробовали реку с двумя сетями; одна сеть имеет вдвое большую площадь, чем другая. У вас больше шансов собрать больше вещей в более крупной сети, и, следовательно, подсчет из более крупной сети будет выше из-за больших усилий по сэмплированию & mdash; Вы охватили большую часть реки более крупной сетью (при условии, что вы взяли пробы за то же время). Чтобы убедиться, что вы учитываете эту разницу в усилиях, вы можете включить смещение в модель. Смещение будет (для модели Пуассона с лог-ссылкой) offset(log(net_area)). Мы должны включить смещение в масштаб ссылки, следовательно, log(). Теперь то, что мы моделируем, это количество на единицу площади нетто.