Warum erhalte ich für poly (raw = T) vs. poly () völlig unterschiedliche Ergebnisse?

10

Ich möchte zwei verschiedene Zeitvariablen modellieren, von denen einige in meinen Daten stark kollinear sind (Alter + Kohorte = Periode). Dabei hatte ich einige Probleme mit lmerund und Interaktionen von poly(), aber es ist wahrscheinlich nicht darauf beschränkt lmer, dass ich mit nlmeIIRC die gleichen Ergebnisse erzielt habe.

Offensichtlich fehlt mein Verständnis dafür, was die poly () - Funktion bewirkt. Ich verstehe, was poly(x,d,raw=T)funktioniert, und ich dachte, ohne raw=Tes entstehen orthogonale Polynome (ich kann nicht sagen, dass ich wirklich verstehe, was das bedeutet), was die Anpassung erleichtert, aber Sie die Koeffizienten nicht direkt interpretieren lässt.
Ich habe gelesen, dass die Vorhersagen gleich sein sollten, da ich die Vorhersagefunktion verwende.

Dies ist jedoch nicht der Fall, selbst wenn die Modelle normal konvergieren. Ich verwende zentrierte Variablen und dachte zuerst, dass das orthogonale Polynom möglicherweise zu einer höheren Korrelation fester Effekte mit dem kollinearen Interaktionsterm führt, aber es scheint vergleichbar. Ich habe hier zwei Modellzusammenfassungen eingefügt .

Diese Darstellungen veranschaulichen hoffentlich das Ausmaß des Unterschieds. Ich habe die Vorhersagefunktion verwendet, die nur im Entwickler verfügbar ist. Version von lme4 ( hier davon gehört ), aber die festen Effekte sind in der CRAN-Version gleich (und sie scheinen auch von selbst aus zu sein, z. B. ~ 5 für die Interaktion, wenn mein DV einen Bereich von 0-4 hat).

Der letzte Anruf war

cohort2_age =lmer(churchattendance ~ 
poly(cohort_c,2,raw=T) * age_c + 
ctd_c + dropoutalive + obs_c + (1+ age_c |PERSNR), data=long.kg)

Die Vorhersage war nur ein fester Effekt auf gefälschte Daten (alle anderen Prädiktoren = 0), wobei ich den in den Originaldaten vorhandenen Bereich als Extrapolation = F markierte.

predict(cohort2_age,REform=NA,newdata=cohort.moderates.age)

Ich kann bei Bedarf mehr Kontext bereitstellen (ich habe es nicht leicht geschafft, ein reproduzierbares Beispiel zu erstellen, kann mich aber natürlich mehr anstrengen), aber ich denke, dies ist eine grundlegendere Bitte: Erklären Sie poly()mir die Funktion, bitte schön.

Rohe Polynome

Rohe Polynome

Orthogonale Polynome (abgeschnitten, bei Imgur nicht abgeschnitten )

Orthogonale Polynome

Ruben
quelle

Antworten:

10

Ich denke, dies ist ein Fehler in der Vorhersagefunktion (und damit meine Schuld), den nlme tatsächlich nicht teilt. ( Bearbeiten : sollte in der neuesten R-Forge-Version von behoben sein lme4.) Ein Beispiel finden Sie unten ...

Ich denke, Ihr Verständnis von orthogonalen Polynomen ist wahrscheinlich in Ordnung. Das Schwierige, was Sie über sie wissen müssen, wenn Sie versuchen, eine Vorhersagemethode für eine Klasse von Modellen zu schreiben, ist, dass die Basis für die orthogonalen Polynome basierend auf einem bestimmten Datensatz definiert wird. Wenn Sie also naiv sind (wie ich! ) Wenn Sie model.matrixversuchen, die Entwurfsmatrix für einen neuen Datensatz zu generieren, erhalten Sie eine neue Basis - was mit den alten Parametern keinen Sinn mehr ergibt. Bis ich dies behoben habe, muss ich möglicherweise eine Falle einfügen, die Personen mitteilt, predictdie nicht mit orthogonalen Polynombasen (oder Spline-Basen mit derselben Eigenschaft) arbeiten.

d <- expand.grid(x=seq(0,1,length=50),f=LETTERS[1:10])
set.seed(1001)
u.int <- rnorm(10,sd=0.5)
u.slope <- rnorm(10,sd=0.2)
u.quad <- rnorm(10,sd=0.1)
d <- transform(d,
               ypred = (1+u.int[f])+
               (2+u.slope[f])*x-
               (1+u.quad[f])*x^2)
d$y <- rnorm(nrow(d),mean=d$ypred,sd=0.2)
ggplot(d,aes(x=x,y=y,colour=f))+geom_line()+
    geom_line(aes(y=ypred),linetype=2)

library(lme4)
fm1 <- lmer(y~poly(x,2,raw=TRUE)+(1|f)+(0+x|f)+(0+I(x^2)|f),
            data=d)


fm2 <- lmer(y~poly(x,2)+(1|f)+(0+x|f)+(0+I(x^2)|f),
            data=d)
newdat <- data.frame(x=unique(d$x))
plot(predict(fm1,newdata=newdat,REform=NA))
lines(predict(fm2,newdata=newdat,REform=NA),col=2)
detach("package:lme4")

library(nlme)
fm3 <- lme(y~poly(x,2,raw=TRUE),
           random=list(~1|f,~0+x|f,~0+I(x^2)|f),
            data=d)
VarCorr(fm3)

fm4 <- lme(y~poly(x,2),
           random=list(~1|f,~0+x|f,~0+I(x^2)|f),
            data=d)

newdat <- data.frame(x=unique(d$x))
lines(predict(fm3,newdata=newdat,level=0),col=4)
lines(predict(fm4,newdata=newdat,level=0),col=5)
Ben Bolker
quelle
Danke, das ist beruhigend. Um es noch einmal zu wiederholen: Ich habe gelesen, dass man die orthogonalen Polynom-Fixeffekte nicht zum Nennwert nehmen kann, aber manchmal scheinen sie wahnsinnig groß zu sein. Wenn ich zum Beispiel eine Wechselwirkung von zwei kubischen Polynomen durchführe, erhalte ich feste Effekte für die Polynome und ihre Wechselwirkungen im Bereich von -22 bis -127400. Das scheint mir einfach weit weg zu sein, besonders wenn man bedenkt, dass alle festen Effekte negativ sind. Würde eine überarbeitete Vorhersagefunktion für diese festen Effekte Sinn machen oder konvergierten die Modelle fälschlicherweise oder stimmt etwas in lmer doch nicht?
Ruben
Wieder vermute ich (aber weiß es offensichtlich nicht genau), dass alles in Ordnung ist. Orth. Polynome eignen sich gut für numerische Stabilitäts- und Hypothesentests, aber (wie Sie herausfinden) können die tatsächlichen Parameterwerte schwieriger zu interpretieren sein. Die aktuelle Version von lme4-devel (ich habe gerade eine Version veröffentlicht, die Tests bestehen sollte. Die Wiederherstellung auf r-forge kann ~ 24 Stunden dauern, es sei denn, Sie können selbst aus SVN erstellen) sollte Ihnen übereinstimmende Vorhersagen zwischen Roh- / Orthopolynomen liefern. Eine Alternative besteht darin, kontinuierliche Prädiktoren à la Schielzeth 2010 Methoden in Ökologie und Evolution zu zentrieren und zu skalieren ...
Ben Bolker
Ja, die beiden Polynome stimmen jetzt vollkommen überein. Vielen Dank! Ich hatte meine Prädiktoren skaliert und zentriert, aber einige Modelle passten nicht zu rohen Polynomen.
Ruben