Pourquoi je reçois des prédictions différentes pour l'expansion polynomiale manuelle et j'utilise la fonction R `poly`?

10

Pourquoi j'obtiens des prévisions différentes pour l'expansion polynomiale manuelle et l'utilisation de la polyfonction R ?

set.seed(0)
x <- rnorm(10)
y <- runif(10)
plot(x,y,ylim=c(-0.5,1.5))
grid()

# xp is a grid variable for ploting
xp <- seq(-3,3,by=0.01)
x_exp <- data.frame(f1=x,f2=x^2)
fit <- lm(y~.-1,data=x_exp)
xp_exp <- data.frame(f1=xp,f2=xp^2)
yp <- predict(fit,xp_exp)
lines(xp,yp)

# using poly function
fit2 <- lm(y~ poly(x,degree=2) -1)
yp <- predict(fit2,data.frame(x=xp))
lines(xp,yp,col=2)

entrez la description de l'image ici

Ma tentative:

  • Il semble que ce soit un problème d'interception, lorsque j'adapte le modèle avec interception, c'est-à-dire, pas -1de modèle formula, les deux lignes sont identiques. Mais pourquoi sans l'interception les deux lignes sont différentes?

  • Un autre «correctif» utilise l' rawexpansion polynomiale au lieu du polynôme orthogonal. Si nous changeons le code en fit2 = lm(y~ poly(x,degree=2, raw=T) -1), cela fera 2 lignes identiques. Mais pourquoi?

Haitao Du
la source
4
C'est hors sujet de votre question, mais vous êtes souvent très ouvert aux commentaires. Lors de la lecture de votre code, la première chose que je remarque est que vous utilisez =et <-pour l'affectation de manière incohérente. Je ne ferais vraiment pas cela, ce n'est pas vraiment déroutant, mais cela ajoute beaucoup de bruit visuel à votre code sans aucun avantage. Vous devez vous contenter de l'un ou de l'autre à utiliser dans votre code personnel et vous y tenir.
Matthew Drury
merci de m'avoir aidé à coder! question corrigée. @MatthewDrury
Haitao Du
3
Conseil de suivi au hasard pour faire <-moins de tracas à taper: alt+-.
JAD
@JarkoDubbeldam merci pour le conseil de codage. J'adore les raccourcis clavier
Haitao Du

Réponses:

12

Comme vous le constatez correctement, la différence d'origine est due au fait que dans le premier cas vous utilisez les polynômes "bruts" tandis que dans le second cas vous utilisez les polynômes orthogonaux. Par conséquent, si l' lmappel ultérieur était modifié en: fit3<-lm(y~ poly(x,degree=2, raw = TRUE) -1)nous obtiendrions les mêmes résultats entre fitet fit3. La raison pour laquelle nous obtenons les mêmes résultats dans ce cas est "triviale"; nous adaptons exactement le même modèle que nous avons équipé fit<-lm(y~.-1,data=x_exp), pas de surprise là-bas.

On peut facilement vérifier que les matrices de modèles des deux modèles sont identiques all.equal( model.matrix(fit), model.matrix(fit3) , check.attributes= FALSE) # TRUE).


Ce qui est plus intéressant, c'est pourquoi vous obtiendrez les mêmes tracés lorsque vous utiliserez une interception. La première chose à noter est que, lors de l'ajustement d'un modèle avec une interception

  • Dans le cas de fit2nous déplaçons simplement les prédictions du modèle verticalement; la forme réelle de la courbe est la même.

  • D'autre part, inclure une interception dans le cas des fitrésultats non seulement dans une ligne différente en termes de placement vertical, mais avec une forme globale complètement différente.

Nous pouvons facilement voir cela en ajoutant simplement les ajustements suivants sur le tracé existant.

fit_b<-lm(y~. ,data=x_exp)
yp=predict(fit_b,xp_exp)
lines(xp,yp, col='green', lwd = 2)

fit2_b<-lm(y~ poly(x,degree=2, raw = FALSE) )
yp=predict(fit2_b,data.frame(x=xp))
lines(xp,yp,col='blue')

entrez la description de l'image ici

OK ... Pourquoi les ajustements sans interception étaient-ils différents alors que les ajustements avec interception étaient les mêmes? Le hic est de nouveau à la condition d'orthogonalité.

Dans le cas où fit_bla matrice modèle utilisée contient des éléments non orthogonaux, la matrice Gram crossprod( model.matrix(fit_b) )est loin d'être diagonale; dans le cas des fit2_béléments orthogonaux ( crossprod( model.matrix(fit2_b) )est effectivement diagonale).

fitfit_b XTXfitfit2fit2_b

La question secondaire intéressante est pourquoi les fit_bet fit2_bsont les mêmes; après tout, les matrices de modèle de fit_bet fit2_bne sont pas les mêmes en valeur nominale . Ici, nous devons simplement nous souvenir de cela en fin de compte fit_bet fit2_bavoir les mêmes informations. fit2_best juste une combinaison linéaire des fit_bajustements si essentiellement résultants seront les mêmes. Les différences observées dans le coefficient ajusté reflètent la recombinaison linéaire des valeurs de fit_bafin de les rendre orthogonales. (Voir la réponse de G. Grothendieck ici aussi pour un exemple différent.)

usεr11852
la source
+2,5 merci pour la bonne réponse. Pour le graphique final, j'ai appris de @kjetilb halvorsen: Une façon plus abstraite de décrire cela est que le modèle lui-même ne dépend que d'un certain sous-espace linéaire, à savoir l'espace de colonne défini par la matrice de conception. Mais les paramètres dépendent non seulement de ce sous-espace, mais de la base de ce sous-espace, donnée par les variables spécifiques utilisées, c'est-à-dire les colonnes elles-mêmes. Les prédictions du modèle, par exemple, ne dépendront que du sous-espace linéaire et non de la base choisie.
Haitao Du
j'espère que cela ne vous dérange pas, j'ai reformaté un peu ..
Haitao Du
@ hxd1011: Pas de problème du tout, merci d'avoir pris le temps de "peigner" un peu.
usεr11852