Modellbau: Wie baue ich ein aussagekräftiges Spielmodell? (verallgemeinertes additives Modell)

9

Ich habe gesehen, dass es verschiedene Fragen bezüglich der Interpretation und Konstruktion von Gams gibt, was die Schwierigkeit für Nicht-Statistiker zu veranschaulichen scheint, mit diesen umzugehen. Leider konnte ich aus keinem der von mir gelesenen Threads oder Tutorials ein klares Verständnis dafür gewinnen, wie man ein aussagekräftiges Modell erstellt.

Derzeit untersuche ich die Auswirkungen des ökologischen Landbaus auf die Leistung von Honigbienenkolonien. Ich versuche dabei, Landschaftsmerkmale wie den Prozentsatz des ökologischen Landbaus in einem Radius von 500 m (bio.percent_b500) mit einem Kolonieentwicklungsparameter wie dem Honigreservat in Beziehung zu setzen. Ich habe zuerst ein grundlegendes Spielmodell (Modell 0) mit nur der Woche des Jahres als erklärende Variable erstellt, da die Menge an Honig in den Bienenstöcken im Laufe eines Jahres nicht linear variiert.

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

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

Geben Sie hier die Bildbeschreibung ein

Dann habe ich versucht, einen glatten Begriff aufzunehmen, der den Prozentsatz des ökologischen Landbaus enthält. Dies schlug jedoch fehl, da mehr als 85% der Kolonien in einem Radius von 500 m keine organischen Felder hatten.

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)

Ich war dann überrascht zu sehen, dass das Modell, das eine Wechselwirkung zwischen dem Prozentsatz des ökologischen Landbaus und der Woche beinhaltete, funktionierte. Ich habe jedoch in einem deutschen Statistikbuch gelesen, dass Interaktionsterme nicht in Modelle ohne deren unabhängige Auswirkungen aufgenommen werden sollten. Der Autor bezog sich auf etwas, das als "Marginalitätsstheorem" bezeichnet wird. Da ich wusste, dass aus Modell 1, dass der glatte Begriff für den ökologischen Landbau Probleme verursacht, habe ich nur einen zusätzlichen glatten Begriff für die Woche des Jahres aufgenommen. Dieses Modell macht für mich intuitiv Sinn, da die Woche des Jahres immer einen Effekt hat; Die Wirkung des ökologischen Landbaus hängt jedoch immer von der Jahreszeit ab. Zum Beispiel sollte es im Sommer eine höhere Verfügbarkeit von Unkrautblumen geben.

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

Da die Honigreserven in den Bienenstöcken wahrscheinlich von verschiedenen Landschaftsmerkmalen abhängen, habe ich Modelle gebaut, einschließlich des Prozentsatzes an Raps (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)

Die Modelle 0, 3 und 6 erscheinen mir aus den oben genannten Gründen am aussagekräftigsten. Ich bin mir nicht sicher, ob ich Modelle, die auf andere Weise konstruiert wurden, in Betracht ziehen und sie auch über AIC akzeptieren und vergleichen soll.

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

Der Vergleich der AIC-Werte ergab, dass Modell 7 das beste ist, da es weniger Modellfreiheitsgrade als Modell 3 aufweist. Dies ist für mich erneut überraschend, da Modell 7 eine komplexere Interaktion als Modell 3 enthält.

Kann mir jemand einen Rat geben, wie man sinnvolle Spielmodelle konstruiert?

1) Können Interaktionsbegriffe in einem (Spiel-) Modell ohne ihre unabhängigen Begriffe erscheinen?

2) Warum können komplexere Begriffe der Spielinteraktion zu einer Verringerung der Modellfreiheitsgrade führen?

3) Welche der oben genannten Modelle sind sinnvoll?

4) Gibt es bessere Alternativen zu gerneralisierten additiven Modellen für das, was ich versuche zu tun?

Unten finden Sie 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 Bienenmensch
quelle

Antworten:

10

Der Grund für den ersten Fehler ist, dass bio.percent_b500es keine keindeutigen Werte gibt. Wenn Sie kfür diesen Spline einen niedrigeren Wert festlegen , passt das Modell. Warum die 2-D-Dünnplattenversion funktioniert, denke ich, IIRC, liegt an der Art und Weise, wie sie einen Standard berechnet, kwenn sie nicht geliefert wird, und dies muss aus den Daten erfolgen. Die model1Einstellung k = 9funktioniert also wie folgt :

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

In mgcv, ein Modell mit s(x1, x2) umfasst die wesentlichen Effekte und die Wechselwirkungen. Beachten s()Sie, dass in den Variablen eine ähnliche Skalierung angenommen wird. Ich bezweifle, dass Sie dies wirklich wollen. s()Versuchen Sie te()stattdessen, ein Tensorprodukt glatt zu machen. Wenn Sie eine Zerlegung in Haupteffekte und Interaktion wünschen, versuchen Sie:

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

Aus dem summary()geht hervor, dass die Interaktion nicht erforderlich ist:

> 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

Und ähnliche Informationen erhalten Sie aus dem allgemeinen Likelihood-Ratio-Test "

> 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

Hier sind ti()Terme Tensor-Produkt-Interaktionsterme, bei denen die Haupt- / Randeffekte für andere Terme im Modell entfernt werden. Sie müssen nicht unbedingt brauchen , ti()ich glaube , in model1, te()oder sogar s()sollte funktionieren, aber die Beispiele in mgcv alle jetzt dieses Formular verwenden, so dass ich damit gehe.

model3macht für mich nicht viel Sinn, besonders wenn Sie Tensorprodukte verwenden; Der s(week)Begriff ist in der enthalten s(week, bio.percent_b500), aber wie ich bereits erwähnt habe, gehen bivariate s()Begriffe von Isotropie aus, so dass dies die weekKomponente möglicherweise übermäßig eingeschränkt hat und somit Raum für das s(week)Eintreten und Erklären von etwas lässt. Im Allgemeinen sollten Sie dies nicht tun, da Sie den bivariaten Begriff richtig verstehen. Ob Sie mit Ihrer 500-m-Variablen einen echten bivariaten Term richtig machen können, ist mir nicht klar.

Fragen:

Q1

Ich bezweifle, dass Sie eine Interaktion ohne Haupt- / Randeffekte wollen. Ihr model2beinhaltet beides.

Q2

Das Modell mag in Bezug auf die von Ihnen geplanten Funktionssätze komplexer sein, aber Sie müssen die für jeden glatten Term erzeugten Basen berücksichtigen. Außerdem verwendet mgcv Splines mit Strafen, die die Auswahl der Glätte ermöglichen, und daher könnten Sie ein Modell erhalten, das sich in befindet Der erste Blick ist komplexer, aber die glatten Begriffe wurden aufgrund der Strafen so geschrumpft, dass sie nach der Anpassung weniger Freiheitsgrade verbrauchen.

Q3

Technisch ist nur model1wirklich richtig spezifiziert. Ich glaube nicht, dass Sie Isotropie für annehmen wollen model2. Ich glaube , Sie mit Problemen in pass laufen werden model3, model6und model7im Allgemeinen; Wenn Sie die Grundlagen für die Tensorprodukte richtig eingerichtet haben, sollten Sie die separaten s(week)Begriffe in diesen Modellen nicht benötigen .

Gavin Simpson
quelle
Vielen Dank für die tolle Antwort. :) In der Zwischenzeit habe ich alle Variablen skaliert und gamm4 verwendet, um zufällige Effekte einzuschließen. Ist s () nach der Skalierung der Variablen geeignet? gamm4 unterstützt ti () nicht. Würden Sie mir empfehlen, das itsadug-Paket anstelle von gamm4 zu verwenden?
Bienenmensch
1
Es sei denn , Sie wirklich die sparsity-Handling Fähigkeiten benötigen gamm()oder gamm4(), können Sie ti()mit gam()und die Zufallseffekt Spline - Basis verwenden , bs = "re"um zufällige Effekte hinzufügen. itsadug ist nur eine zusätzliche Funktionalität, die auf gamm () -Modellen basiert. Wenn Sie Gaußsche Modelle anpassen, würde ich mich an gamm()oder sogar nur gam()an Splines mit zufälligen Effekten halten.
Gavin Simpson