mgcv: Fix Glättungsparameter in GAMM und Gültigkeit der Modellverschachtelung

Question

Ich versuche derzeit, die Interaktion zwischen einer Kovariaten x (Alter) und einem Geschlecht 0-1 Variable mit gamm aus dem mgcv Paket zu modellieren. Nach Angabe eines Hauptmodells (wir nennen es M0) mit einem glatten Begriff für jedes Geschlecht, möchte ich die einfachere Hypothese testen, dass der Unterschied zwischen den Geschlechtern linear ist (und nicht beliebig glatt). Ich habe die folgenden zwei Fragen:

1. Wenn nisten versuchen, die Modelle richtig würde Ich mag von M0 den Glättungsparameter für das 0-Geschlecht glatt bekommen, und es in der einfacheren Modellspezifikation zu verwenden. Aber ich bekomme die folgende Fehlermeldung:

Fehler in gamm.setup (gp, PTERMs = PTERMs, data = mf, Knoten Knoten =, parametric.only = FALSCH:

GAMM kann nicht Griff verbunden Glättungsparameter (wahrscheinlich aus der Benutzung `id‘ oder adaptiven glättend)

1. Zweite Frage ist es, die dümmer ein. Sind die Modelle sogar verschachtelte, wenn ich von einer glatten gehen für jedes Geschlecht ein 0-Geschlecht glatt und ein linearer Unterschied zu 1-Geschlecht?

Unten folgt ein Beispiel. Ich habe einige zufällige Daten simuliert, daher zeigen die Daten nicht das Verhalten meiner tatsächlichen Daten, aber die Probleme bleiben gleich.

library(mgcv) 

### Simulate random data 
x <- rnorm(100, mean = 10, sd = 1.5) 
y <- rnorm(100, mean = 1, sd = 0.025) 

id <- sample(1:10, size = 100, replace = T) 
id <- as.factor(id) 

gender <- sample(c(0,1), size = 100, replace = T) 


### Specify main model, M0 
ctrl <- list(niterEM=0,optimMethod="L-BFGS-B", msMaxIter = 100) 
M0 <- gamm(y~s(x, by = as.factor(gender)) + gender, 
      random=list(id=~1+x), control=ctrl) 

### Specify model with linear difference between gender0 and gender1 
M1 <- gamm(y~s(x) + gender:x + gender, 
      random=list(id=~1+x), control=ctrl) 

### Testing 
anova(M0$lme, M1$lme) 

### Problems: 
sp0 <- M0$gam$sp[1] 
M1 <- gamm(y~s(x, sp = sp0) + gender:x + gender, 
      random=list(id=~1+x), control=ctrl) 

Fehler in gamm.setup (gp, PTERMs = PTERMs, data = mf, Knoten = Knoten, parametric.only = FALSCH:

GAMM keine Parameter verarbeiten kann (wahrscheinlich verknüpft Glättung von Verwendung von `id‘ oder adaptive glättend)

Irgendwelche Gedanken? Vielen Dank im Voraus.

Answer 1

in Bezug auf die gamm Fehler

Das ist eine sehr interessante Sache ... Nun, ich sollte zuerst die Logik erklären.

Grundsätzlich ist es illegal, einen Glättungsparameter in gamm zu beheben, weil gamm wird die Wackel Komponenten der glatten als zufällige Effekte behandeln, die Varianz von denen von lme geschätzt werden (wie Sie Gaussian Daten). Wenn Sie versuchen, den Glättungsparameter zu korrigieren, heißt das im Wesentlichen, dass Sie die Varianz des Zufallseffekts korrigieren möchten. Nun, lme erlaubt es Ihnen nicht, dies zu tun (und ich bezweifle, dass ein solcher Versuch auch in der gemischten Modellierung legal ist).)

Daher gamm mögliche Einschränkungen Glättungsparameter deaktivieren würde, einschließlich:

1. untere Grenze des Glättungsparameters, über min.sp;
2. verknüpfte glatte teilen die gleiche Glättungsparameter, über id in s();
3. direkt angegebenen Glättungsparameter über sp, über s().

Die ersten beiden sind perfekt überprüft. gamm hat kein min.sp Argument wie gam; Selbst wenn Sie es über ... angeben, gibt es keine Chance, dass es verwendet wird (wie später ist es NULL, die gam.setup während gamm.setup übergeben wird, so dass Ihre angegebene min.sp ignoriert wird). Die Angabe id wird auch durch die Fehlermeldung, die Sie gesehen haben, erkannt, aber Sie haben id natürlich nicht angegeben, so dass der obige Fehler nicht das korrekte Problem meldet, daher ein Fehler.

Die dritte wurde tatsächlich nicht direkt über gamm überprüft. Idealerweise sollte, sobald die Formel gamm/gam interpretiert wurde (durch interpret.gam), sp auf -1 zurückgesetzt werden, wenn dies nicht sofort der Fall ist, und eine Warnmeldung darüber ausgegeben werden sollte. Dieser Teil fehlt jedoch. Daher ist das Beste, was Sie jetzt tun können, einfach sp nicht anzugeben.

---

Haben Sie eine gültige Modellverschachtelung?

Kommen wir nun zu Ihrer Sorge über die Verschachtelung. Die Verschachtelung bezieht sich auf die Basiseinstellung und nicht auf die Wahl des Glättungsparameters. Solange Sie den gleichen Basissatz haben (gleicher Basistyp, gleiche Anzahl und/oder gleiche Stelle von "Knoten"), ist die Modellmatrix identisch. Zum Beispiel haben Sie in Ihrem Modell M0 und M1 die gleiche Konfiguration des s() mit mgcv Standard bs = 'tp', k = 10. So ist die Design-Matrix für s() die gleiche in Ihren beiden Modellen. Die by = factor(gender) repliziert nur diese s() auf alle Ebenen von gender in Ihrem Hauptmodell M0. Vielleicht ist es nicht leicht zu sehen, aber tatsächlich ist Ihr M1 in M0 verschachtelt.

Betrachten wir ein kleines Beispiel, um dies zu überprüfen. Zur Vereinfachung verwende ich nicht s(x) von mgcv, aber verwenden Sie poly(x, degree = 2) (vorstellen, es ist s(x)). Lassen Sie uns einige Spielzeug-Daten erzeugen zuerst:

x <- 1:10 
f <- gl(2, 5, labels = c("M", "F")) 

Da f ist kein Faktor bestellt, s(x, by = factor(f)) erzeugt Designmatrix durch s(x) für alle Ebenen der f replizierende:

## original design matrix for `s(x)` 
X0 <- poly(x, 2) 
## design matrix for `f`, without contrasting 
Xf <- model.matrix(~ f + 0) 
## design matrix for level `M` 
X1 <- X0 * Xf[, 1] 
## design matrix for level `F` 
X2 <- X0 * Xf[, 2] 
## design matrix for `s(x, by = f)` "please, imagine it as `poly`" 
X <- cbind(X1, X2) 

#    1   2   1   2 
# [1,] -0.49543369 0.52223297 0.00000000 0.00000000 
# [2,] -0.38533732 0.17407766 0.00000000 0.00000000 
# [3,] -0.27524094 -0.08703883 0.00000000 0.00000000 
# [4,] -0.16514456 -0.26111648 0.00000000 0.00000000 
# [5,] -0.05504819 -0.34815531 0.00000000 0.00000000 
# [6,] 0.00000000 0.00000000 0.05504819 -0.34815531 
# [7,] 0.00000000 0.00000000 0.16514456 -0.26111648 
# [8,] 0.00000000 0.00000000 0.27524094 -0.08703883 
# [9,] 0.00000000 0.00000000 0.38533732 0.17407766 
#[10,] 0.00000000 0.00000000 0.49543369 0.52223297 

Ihr zweites Modell M1, hat nur glatte Bezeichnung s(x), deren Design-Matrix X0 ist.

Hier ist, wie wir sehen können, dass Ihr M1 in M0 verschachtelt ist:

1. Als X1 + X2 = X0, Design-Matrix von s(x) und s(x, by = f) die gleiche Spaltenspanne haben, so ist s(x) in s(x, by = f) verschachtelt;
2. seit x ist verschachtelt in s(x), x:f ist in s(x, by = f) verschachtelt.

---

Was würde ich

tun Obwohl Sie Modelle sind bereits schön verschachtelt, das Hauptmodell M0 keine vergleichbare Interpretation mit Ihrem M1 hat. Ihr Hauptmodell M0 wird mit einem unabhängigen Glatt für jedes Level enden, während Ihr M1 auf den Unterschied zwischen zwei Gruppen konzentrieren.

Es wäre gut, wenn wir M0 steuern könnten, um eine Form von "Referenz glatt + Unterschied glatt" zuzulassen. Dann, wenn die Differenz glatt eine Linie ausmacht, ohne tatsächlich zu passen M1 wir bereits wissen, gibt es keinen Beweis für nicht-lineare Unterschied zwischen den Gruppen.

In mgcv, Unterschied glatt wird konstruiert, wenn Ihr Faktor bestellt wird. Also schlage ich vor, Sie Ihr Haupt-Modell passen von:

gender1 <- ordered(gender) ## create an ordered factor 
s(x) + s(x, by = gender1) + gender 

Wenn Schätzergebnis der Differenz glatt s(x, by = gender1) als Linie zeigt, wissen Sie, Sie stattdessen ein einfacheres Modell passen

s(x) + gender:x + gender 

auch ohne F-test zu verwenden.

Hinweis, es ist sehr wichtig, einen geordneten Faktor by zu haben, um "Unterschied" glatt zu konstruieren. Wenn Sie das tun

s(x) + s(x, by = gender) + gender ## note, it is "gender" in "by" 

s(x) und s(x, by = gender) sind vollständig linear abhängig sind. Die resultierende Modellmatrix ist Rang-defizient.

---

Einige Follow-up

Ich habe vergessen, in meinem Beispiel enthalten, die ich anfangs das gleiche Modell wie s(x, by = as.factor(gender)) und s(x) + s(x, by = gender) von AIC parametrisiert verglichen (man erinnere sich gender 0-1 numerische Variable). Diese Modelle sind statistisch äquivalent, aber die Glättungsparameter werden in den Fällen offensichtlich unterschiedlich geschätzt und die AIC unterscheiden sich daher ein wenig.

Oh, ja. Ihr gender ist binär, also ist eine numerische by auch eine gute Idee für das Konstruieren eines glatten Unterschieds. Aber mach das vorsichtig. Numerisch by ergibt keine zentrierte Glättung. Daher wird s(x) + s(x, by = gender) implizit eine Intercept-Spalte haben, die mit dem Modellabschnitt verwechselt wird.Sie sollten mit s(x) + s(x, by = gender) - 1 gehen.