Construção do modelo: como criar um modelo gam significativo? (modelo aditivo generalizado)

9

Eu já vi várias questões relacionadas à interpretação e construção de jogos, o que parece ilustrar a dificuldade dos não-estatísticos em lidar com eles. Infelizmente, a partir de nenhum dos tópicos ou tutoriais que li, pude entender claramente como criar um modelo significativo.

Atualmente, estou estudando o efeito da agricultura orgânica no desempenho das colônias de abelhas. Assim, tento relacionar características da paisagem como a porcentagem de agricultura orgânica em um raio de 500m (bio.percent_b500) a um parâmetro de desenvolvimento de colônia, como a reserva de mel. Primeiro, construí um modelo gam básico (modelo0) com apenas a semana do ano como variável explicativa, pois a quantidade de mel nas colméias varia não linearmente ao longo de um ano.

library("gam")
library("mgcv")

model0 <- gam(honey.mean ~ s(week), data= my.data.frame) 
summary(model0)
plot(model0)

insira a descrição da imagem aqui

Depois, tentei incluir um termo suave que continha a porcentagem de agricultura orgânica. No entanto, isso falhou, eu acho, porque mais de 85% das colônias não tinham campos orgânicos em um raio de 500m.

model1 <- gam(honey.mean ~ s(week) + s(bio.percent_b500),data = my.data.frame)
# Error in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots) : 
# A term has fewer unique covariate combinations than specified maximum 
# degrees of freedom

model2 = gam(honey.mean ~ s(week,bio.percent_b500) , data= my.data.frame)

Fiquei surpreso ao ver que o modelo que incluía uma interação entre o percentual de agricultura orgânica e a semana funcionava. No entanto, li em um livro de estatística alemã que os termos de interação não devem ser incluídos nos modelos sem seus efeitos independentes. O autor se referiu a algo chamado "Marginalitätstheorem" (teorema da marginalidade). Como sabia que, a partir do modelo1, que o termo suave para agricultura orgânica causa problemas, incluí apenas um termo suave adicional para a semana do ano. Esse modelo faz sentido intuitivamente para mim, pois a semana do ano sempre tem um efeito; o efeito da agricultura orgânica é sempre sempre dependente da época do ano. Por exemplo, no verão, deve haver maior disponibilidade de flores para plantas daninhas.

model3 = gam(honey.mean ~ s(week) + s(week, bio.percent_b500) , data= my.data.frame)

Como as reservas de mel nas colméias provavelmente dependem de várias características da paisagem, construí modelos, incluindo a porcentagem de colza (osr.percent_b500).

model4 = gam(honey.mean ~ s(week) + s(osr.percent_b500),data = my.data.frame)
vis.gam(model4, type = "response", plot.type = "persp")   
summary(model4)

model5 = gam(honey.mean ~ s(week,osr.percent_b500) + s(week,bio.percent_b500), data = my.data.frame)
summary(model5)

model6 = gam(honey.mean ~ s(week) + s(week,osr.percent_b500) + s(week,bio.percent_b500), data= my.data.frame)
summary(model6)

model7 = gam(honey.mean ~ s(week) + s(week,osr.percent_b500,bio.percent_b500), data= my.data.frame)
summary(model7)

Os modelos 0, 3 e 6 parecem-me mais significativos pelas razões acima mencionadas. Não tenho certeza se devo considerar os modelos construídos de outra maneira e também aceitá-los e compará-los via AIC.

AIC(model0,model2,model3,model4,model5,model6,model7)

A comparação dos valores da AIC identificou o modelo 7 como o melhor, porque possui menos graus de liberdade do modelo que o modelo 3. Isso é novamente surpreendente para mim, pois o modelo 7 inclui uma interação mais complexa que o modelo 3.

Alguém poderia me dar conselhos sobre como construir modelos gam significativos?

1) Os termos de interação podem aparecer em um modelo (gam) sem seus termos independentes?

2) Por que termos mais complexos de interação gam complexa podem levar a uma diminuição nos graus de liberdade do modelo?

3) Quais dos modelos mencionados acima são significativos?

4) Existem alternativas melhores aos modelos de aditivos generalizados para o que estou tentando fazer?

Abaixo você encontra my.data.frame:

structure(list(year = c(2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 
2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 
2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 
2008L, 2008L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 
2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 
2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 
2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 
2009L, 2009L, 2009L, 2009L, 2009L, 2010L, 2010L, 2010L, 2010L, 
2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 
2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 
2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 
2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 
2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 
2010L, 2010L, 2010L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 
2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 
2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 
2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 
2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 
2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 
2011L, 2011L, 2011L, 2011L, 2012L, 2012L, 2012L, 2012L, 2012L, 
2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 
2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 
2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 
2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 
2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 
2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2013L, 
2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
2013L, 2013L), apiary = c(4L, 8L, 8L, 8L, 18L, 18L, 18L, 19L, 
19L, 19L, 23L, 23L, 23L, 23L, 34L, 34L, 34L, 45L, 45L, 45L, 46L, 
46L, 46L, 49L, 49L, 49L, 3L, 3L, 3L, 3L, 9L, 9L, 9L, 9L, 14L, 
14L, 14L, 14L, 17L, 17L, 17L, 17L, 20L, 20L, 20L, 28L, 28L, 28L, 
28L, 31L, 31L, 31L, 31L, 33L, 33L, 33L, 33L, 33L, 35L, 35L, 35L, 
44L, 44L, 44L, 44L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 
12L, 12L, 26L, 26L, 26L, 26L, 26L, 30L, 30L, 30L, 30L, 30L, 32L, 
32L, 32L, 32L, 32L, 37L, 37L, 37L, 37L, 37L, 42L, 42L, 42L, 42L, 
42L, 47L, 47L, 47L, 47L, 47L, 47L, 47L, 48L, 48L, 48L, 48L, 48L, 
50L, 50L, 50L, 50L, 1L, 1L, 1L, 1L, 1L, 6L, 6L, 6L, 6L, 6L, 6L, 
7L, 7L, 7L, 7L, 7L, 7L, 22L, 22L, 22L, 22L, 24L, 24L, 24L, 24L, 
24L, 24L, 27L, 27L, 27L, 27L, 27L, 27L, 36L, 36L, 36L, 36L, 36L, 
40L, 40L, 40L, 40L, 40L, 41L, 41L, 41L, 41L, 41L, 43L, 43L, 43L, 
43L, 43L, 43L, 43L, 2L, 2L, 2L, 2L, 2L, 2L, 5L, 5L, 5L, 5L, 5L, 
10L, 10L, 10L, 10L, 10L, 10L, 13L, 13L, 13L, 13L, 13L, 15L, 15L, 
15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L, 16L, 16L, 21L, 21L, 21L, 
21L, 21L, 21L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 29L, 29L, 29L, 
29L, 29L, 29L, 29L, 39L, 39L, 39L, 39L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 18L, 18L, 18L, 18L, 18L, 18L, 
18L, 19L, 19L, 19L, 19L, 19L, 19L, 23L, 23L, 23L, 23L, 23L, 23L, 
23L, 34L, 34L, 34L, 34L, 34L, 34L, 38L, 38L, 38L, 38L, 38L, 38L, 
38L, 45L, 45L, 45L, 45L, 45L, 45L, 46L, 46L, 46L, 46L, 46L, 46L, 
46L, 49L, 49L, 49L, 49L, 49L, 49L), week = c(26L, 24L, 26L, 28L, 
23L, 28L, 31L, 23L, 24L, 28L, 24L, 26L, 28L, 29L, 23L, 26L, 28L, 
24L, 26L, 29L, 23L, 28L, 29L, 23L, 28L, 31L, 18L, 20L, 22L, 32L, 
18L, 20L, 30L, 32L, 16L, 22L, 26L, 32L, 16L, 18L, 24L, 28L, 16L, 
24L, 32L, 16L, 24L, 28L, 30L, 18L, 20L, 22L, 26L, 16L, 20L, 22L, 
26L, 30L, 16L, 24L, 28L, 18L, 26L, 28L, 32L, 20L, 21L, 33L, 35L, 
39L, 21L, 25L, 27L, 29L, 31L, 35L, 21L, 25L, 27L, 31L, 35L, 21L, 
23L, 29L, 35L, 39L, 17L, 27L, 33L, 35L, 39L, 17L, 20L, 27L, 35L, 
39L, 17L, 21L, 23L, 25L, 35L, 17L, 20L, 21L, 25L, 27L, 31L, 33L, 
17L, 21L, 23L, 29L, 39L, 20L, 31L, 33L, 39L, 19L, 21L, 23L, 29L, 
37L, 19L, 21L, 23L, 29L, 33L, 39L, 17L, 19L, 25L, 29L, 31L, 35L, 
19L, 33L, 37L, 39L, 15L, 19L, 23L, 35L, 37L, 39L, 15L, 17L, 21L, 
29L, 33L, 35L, 17L, 23L, 25L, 29L, 39L, 17L, 19L, 21L, 29L, 35L, 
17L, 19L, 21L, 25L, 39L, 15L, 19L, 27L, 31L, 33L, 37L, 39L, 13L, 
23L, 27L, 33L, 35L, 39L, 23L, 25L, 27L, 31L, 37L, 13L, 15L, 19L, 
23L, 29L, 37L, 29L, 33L, 35L, 37L, 39L, 13L, 21L, 25L, 27L, 29L, 
35L, 23L, 29L, 31L, 35L, 37L, 39L, 15L, 19L, 21L, 27L, 33L, 39L, 
13L, 15L, 23L, 27L, 29L, 35L, 39L, 13L, 15L, 23L, 27L, 29L, 31L, 
35L, 13L, 31L, 35L, 37L, 16L, 20L, 26L, 38L, 40L, 42L, 44L, 16L, 
24L, 32L, 34L, 38L, 40L, 44L, 18L, 20L, 24L, 34L, 38L, 42L, 44L, 
24L, 28L, 32L, 40L, 42L, 44L, 16L, 20L, 26L, 38L, 40L, 42L, 44L, 
18L, 20L, 22L, 32L, 38L, 44L, 16L, 20L, 22L, 28L, 30L, 34L, 38L, 
18L, 20L, 22L, 28L, 32L, 44L, 16L, 22L, 24L, 28L, 32L, 34L, 38L, 
22L, 28L, 32L, 34L, 38L, 40L), bio.percent_b500 = c(0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 16.13, 16.13, 16.13, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15.73, 15.73, 15.73, 
15.73, 15.73, 0, 0, 0, 0, 0, 0.75, 0.75, 0.75, 0.75, 0.75, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 2.14, 2.14, 2.14, 2.14, 2.14, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 13.69, 13.69, 13.69, 13.69, 
13.69, 13.69, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 6.47, 6.47, 6.47, 6.47, 6.47, 6.47, 6.47, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.68, 5.68, 
5.68, 5.68, 5.68, 5.68, 5.68, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 44.93, 44.93, 44.93, 44.93, 44.93, 
44.93, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0), osr.percent_b500 = c(10.12, 1.51, 1.51, 
1.51, 0, 0, 0, 4.85, 4.85, 4.85, 0, 0, 0, 0, 8.94, 8.94, 8.94, 
0, 0, 0, 0, 0, 0, 1.2, 1.2, 1.2, 6.41, 6.41, 6.41, 6.41, 0, 0, 
0, 0, 8.27, 8.27, 8.27, 8.27, 4.67, 4.67, 4.67, 4.67, 7.2, 7.2, 
7.2, 5.84, 5.84, 5.84, 5.84, 20.51, 20.51, 20.51, 20.51, 10.22, 
10.22, 10.22, 10.22, 10.22, 9.85, 9.85, 9.85, 0.02, 0.02, 0.02, 
0.02, 14.33, 14.33, 14.33, 14.33, 14.33, 21.6, 21.6, 21.6, 21.6, 
21.6, 21.6, 0, 0, 0, 0, 0, 6.1, 6.1, 6.1, 6.1, 6.1, 3.18, 3.18, 
3.18, 3.18, 3.18, 5.45, 5.45, 5.45, 5.45, 5.45, 0, 0, 0, 0, 0, 
22.65, 22.65, 22.65, 22.65, 22.65, 22.65, 22.65, 0.52, 0.52, 
0.52, 0.52, 0.52, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
5.59, 5.59, 5.59, 5.59, 5.59, 5.59, 7.41, 7.41, 7.41, 7.41, 4.13, 
4.13, 4.13, 4.13, 4.13, 4.13, 21.77, 21.77, 21.77, 21.77, 21.77, 
21.77, 3.58, 3.58, 3.58, 3.58, 3.58, 7.09, 7.09, 7.09, 7.09, 
7.09, 18.35, 18.35, 18.35, 18.35, 18.35, 0.78, 0.78, 0.78, 0.78, 
0.78, 0.78, 0.78, 0.41, 0.41, 0.41, 0.41, 0.41, 0.41, 12.2, 12.2, 
12.2, 12.2, 12.2, 0.26, 0.26, 0.26, 0.26, 0.26, 0.26, 7.57, 7.57, 
7.57, 7.57, 7.57, 12.8, 12.8, 12.8, 12.8, 12.8, 12.8, 34.1, 34.1, 
34.1, 34.1, 34.1, 34.1, 18.33, 18.33, 18.33, 18.33, 18.33, 18.33, 
12.44, 12.44, 12.44, 12.44, 12.44, 12.44, 12.44, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 1.97, 1.97, 1.97, 1.97, 1.97, 1.97, 1.97, 
18.06, 18.06, 18.06, 18.06, 18.06, 18.06, 18.06, 0, 0, 0, 0, 
0, 0, 0, 16.76, 16.76, 16.76, 16.76, 16.76, 16.76, 0, 0, 0, 0, 
0, 0, 0, 4.99, 4.99, 4.99, 4.99, 4.99, 4.99, 0, 0, 0, 0, 0, 0, 
0, 5.28, 5.28, 5.28, 5.28, 5.28, 5.28, 7.99, 7.99, 7.99, 7.99, 
7.99, 7.99, 7.99, 18.09, 18.09, 18.09, 18.09, 18.09, 18.09), 
    honey.mean = c(2.48, 3.99666666666667, 2.36, 2.94, 3.42, 
    3.71, 4.09, 2.12, 3.92, 4.145, 6.27, 6.92, 9.16, 6.75, 6.8, 
    1.07, 6.06, 1.7, 3.4, 5.805, 4.45, 4.19, 13.61, 3.695, 2.86, 
    8.32, 7.67, 6.81, 3.68, 14.335, 2.78, 3.62, 19.035, 12.77, 
    5.81, 3.05, 10.22, 10.44, 4.43, 8.64, 2.4, 16.41, 2.9, 7.175, 
    15.735, 3.16, 1.49, 5.48, 18.95, 6.885, 4.46, 7.9, 0.68, 
    1.4, 2.5, 8.12, 3.09, 14.72, 5.85, 1.885, 16.44, 8.055, 6.68, 
    8.58, 24.7, 8.135, 8.43, 26.08, 16.83, 9.72, 5.24, 5.65, 
    5.19, 7.35, 17.25, 8.82, 14.95, 12.05, 7.3, 62.4, 16.68, 
    1, 10.65, 10.28, 19.65, 17.26, 6.64, 9.94, 65.15, 12.07, 
    20.62, 7.7, 6.31, 1.68, 20.97, 23.825, 6.5, 6.14, 4.22, 2.47, 
    17.97, 2.61, 3.17, 3.24, 0.57, 0.54, 33.07, 49.8, 9.1, 8.41, 
    7.29, 10.61, 19.67, 3.09, 37.125, 24.99, 18.62, 24.15, 17.96, 
    16.61, 28.86, 7.74, 18.95, 18.45, 15.56, 48.35, 16.045, 8.37, 
    23.47, 5.44, 1.8, 64.27, 17.08, 20.62, 18.465, 18.255, 16.5, 
    23.17, 7.49, 12.55, 7.45, 16.72, 23.29, 7.965, 9.83, 15.39, 
    11.19, 35.85, 16.755, 18.8, 19.51, 10.39, 14.02, 32.82, 12.9466666666667, 
    14.68, 15.79, 12.8, 40.37, 22.27, 14.63, 16.9, 6.65, 2.42, 
    18.24, 9.3, 23.08, 17.94, 57.78, 24.34, 20.06, 18.2, 3.99, 
    6.465, 2.93, 25.98, 19.87, 17.25, 13.21, 9.07, 5.21, 9.48, 
    11.825, 7.58, 3.41, 12.56, 13.58, 22.17, 19.43, 11.7, 36.5, 
    18, 12.675, 5.8, 7.72, 4.41, 1.96, 2.83, 12.04, 17.24, 15.77, 
    17.655, 40.15, 21.87, 17.42, 19.16, 8.91, 5.41, 19.91, 9.65, 
    43.54, 17.72, 2.85, 3.41, 7.4, 7.38, 13.73, 14.16, 20.25, 
    2.77, 5.93, 11.185, 2.36, 12.62, 30.24, 13.97, 9.11, 13.985, 
    12.54, 11.13, 1.54, 8.91, 1.3, 4.03, 9.2, 8.86, 9.12, 1.11, 
    7.83, 17.985, 0.86, 14.5, 4.17, 5.18, 5.76, 6.22, 3.79, 17.18, 
    15.83, 11.195, 9.99, 12.395, 7.42, 26.15, 18.29, 15.955, 
    14.76, 2.18, 4.41, 3.53, 11.77, 10.1, 12.81, 20.25, 4.9, 
    10.43, 0.84, 8.81, 19.59, 24.94, 1.42, 6.57, 11.38, 1.92, 
    6.97, 19.31, 17.885, 8.07, 11.25, 6.05, 5.55, 30.23, 9.82, 
    4.8, 4.94, 3.835, 2.54, 21.73, 20.84, 19.02, 5.62, 0.72, 
    23.335, 10.745, 10.43, 7.34)), .Names = c("year", "apiary", 
"week", "bio.percent_b500", "osr.percent_b500", "honey.mean"), row.names = c(NA, 
296L), class = "data.frame")
r model-selection gam cara de abelha
fonte

Respostas:

10

A razão do primeiro erro é que bio.percent_b500não possui k- 1 valores únicos. Se você definir kalgo mais baixo para esse spline, o modelo será ajustado. Por que a versão em chapa fina bidimensional funciona, eu acho, IIRC, é por causa da maneira como calcula um padrão para kquando não fornecido e deve fazer isso a partir dos dados. Então model1, definindo k = 9funciona:

> model1 <- gam(honey.mean ~ s(week) + s(bio.percent_b500, k = 9), data = df)

Em mgcv, um modelo s(x1, x2) inclui os principais efeitos e a interação. Observe que s()assume uma escala semelhante nas variáveis, então duvido que você realmente queira s()- tente te()um produto tensor suave. Se você deseja uma decomposição nos principais efeitos e interação, tente:

model2 <- gam(honey.mean ~ ti(week) + ti(bio.percent_b500) + 
                ti(week, bio.percent_b500), data = df)

Pelo summary()que sugere, a interação não é necessária:

> summary(model2)

Family: gaussian 
Link function: identity 

Formula:
honey.mean ~ ti(week) + ti(bio.percent_b500) + ti(week, bio.percent_b500)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  12.2910     0.5263   23.35   <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    
ti(week)                  3.753  3.963 22.572  <2e-16 ***
ti(bio.percent_b500)      3.833  3.974  2.250  0.0461 *  
ti(week,bio.percent_b500) 1.246  1.448  1.299  0.2036    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =   0.25   Deviance explained = 27.3%
GCV = 84.697  Scale est. = 81.884    n = 296

E você obtém informações semelhantes no teste da razão de verossimilhança generalizada "

> model1 <- gam(honey.mean ~ ti(week) + ti(bio.percent_b500, k = 9), data = df)
> anova(model1, model2, test = "LRT")
Analysis of Deviance Table

Model 1: honey.mean ~ ti(week) + ti(bio.percent_b500, k = 9)
Model 2: honey.mean ~ ti(week) + ti(bio.percent_b500) + ti(week, bio.percent_b500)
  Resid. Df Resid. Dev       Df Deviance Pr(>Chi)
1    285.65      23363                           
2    286.17      23433 -0.51811  -69.675   0.1831

Aqui, os ti()termos são termos de interação tensor-produto nos quais os efeitos principais / marginais de outros termos no modelo são removidos. Você não estritamente precisa ti()acho que em model1, te()ou mesmo s()deve funcionar, mas os exemplos mgcv todos agora utilizar este formulário, por isso estou indo com isso.

model3não faz muito sentido para mim, especialmente se você usa produtos tensores; o s(week)termo está incluído no s(week, bio.percent_b500), mas, como mencionei, os s()termos bivariados assumem isotropia; portanto, isso pode ter restringido demais o weekcomponente, deixando espaço para s(week)que ele entre e explique alguma coisa. Em geral, você não deveria fazer isso, pois obtém o termo bivariado corretamente. Não é claro se você pode obter um termo bivariado verdadeiro com sua variável de 500m.

Questões:

Q1

Duvido que você queira uma interação sem efeitos principais / marginais. O seu model2inclui os dois.

Q2

O modelo pode ser mais complexo em termos dos conjuntos de funções que você imagina, mas é necessário considerar as bases produzidas para cada termo suave, mais o mgcv usa splines com penalidades que fazem a seleção da suavidade e, portanto, você pode acabar com um modelo que esteja em À primeira vista, é mais complexo, mas os termos suaves foram encolhidos devido às penalidades, de modo que eles usam menos graus de liberdade depois de ajustados.

Q3

Tecnicamente, somente model1é realmente especificado corretamente. Eu não acho que você queira assumir a isotropia model2. Eu acho que você vai correr em montagem problemas com model3, model6e, model7em geral; se você configurar corretamente as bases para os produtos tensores, não precisará de s(week)termos separados nesses modelos.

Gavin Simpson
fonte
Muito obrigado pela resposta incrível. :) Enquanto isso, passo a escala de todas as variáveis ​​e o uso do gamm4 para incluir efeitos aleatórios. É adequado usar s () após o dimensionamento das variáveis? gamm4 não suporta ti (). Você recomendaria que eu usasse o pacote itsadug em vez de gamm4?
bee guy
1
A menos que você realmente precisa das capacidades de manipulação de esparsidade de gamm()ou gamm4(), você pode usar ti()com gam()e usar a base ranhura efeito aleatório bs = "re"para adicionar efeitos aleatórios. itsadug é apenas uma funcionalidade extra criada sobre os modelos gamm (). Se você está usando modelos gaussianos, eu usaria gamm()ou mesmo apenas gam()splines de efeito aleatório.
Gavin Simpson