Comment interpréter les coefficients d'un ajustement de modèle polynomial?

36

J'essaie de créer un ajustement polynomial du second ordre à certaines données que j'ai. Disons que je trace cette correspondance avec ggplot():

ggplot(data, aes(foo, bar)) + geom_point() + 
       geom_smooth(method="lm", formula=y~poly(x, 2))

Je reçois:

tracé d'ajustement parabolique avec bande de confiance sur le diagramme de dispersion

Ainsi, un ajustement de deuxième ordre fonctionne assez bien. Je le calcule avec R:

summary(lm(data$bar ~ poly(data$foo, 2)))

Et je reçois:

lm(formula = data$bar ~ poly(data$foo, 2))
# ...
# Coefficients:
#                     Estimate Std. Error t value Pr(>|t|)    
# (Intercept)         3.268162   0.008282 394.623   <2e-16 ***
# poly(data$foo, 2)1 -0.122391   0.096225  -1.272    0.206
# poly(data$foo, 2)2  1.575391   0.096225  16.372   <2e-16 ***
# ....

Maintenant, je suppose que la formule pour mon ajustement est la suivante:

bar=3.268-0,122foo+1,575foo2

Mais cela ne me donne que de mauvaises valeurs. Par exemple, avec étant 3, je pense que deviendrait quelque chose autour de 3.15. Cependant, en insérant dans la formule ci-dessus, je reçois: barfoobar

bar=3.268-0,1223+1,57532=17.077

Ce qui donne? Est-ce que j'interprète mal les coefficients du modèle?

utilisateur13907
la source
2
Plusieurs réponses permettent de répondre à cette question en recherchant sur notre site des polynômes orthogonaux
whuber
6
@ whuber Si j'avais su que le problème était avec les "polynômes orthogonaux", j'aurais probablement trouvé une réponse. Mais si vous ne savez pas quoi chercher, c'est un peu difficile.
user13907
2
Vous pouvez également trouver des réponses en effectuant une recherche sur poly , qui apparaît en évidence dans votre code. Je mets ces informations dans les commentaires pour deux raisons: (1) les liens peuvent aider les lecteurs à venir ainsi que vous-même et (2) ils peuvent vous aider à vous montrer comment exploiter notre système de recherche (quelque peu idiosyncratique).
whuber
7
Vous avez posté une question relative à votre utilisation de polysans taper d'abord ?polydans R? Cela dit " Calculez les polynômes orthogonaux " en haut en grosses lettres amicales.
Glen_b -Reinstate Monica
4
@Glen_b Ouais, eh bien, je l'ai fait taper ?polypour comprendre la syntaxe. Certes, je ne connais que très peu les concepts qui la sous-tendent Je ne savais pas qu'il y avait autre chose (ou une différence si importante entre les polynômes "normaux" et les polynômes orthogonaux), et les exemples que j'ai vus en ligne étaient tous utilisés poly()pour le montage, en particulier avec ggplot- alors pourquoi ne pas simplement utiliser cela et être confus si le résultat était "faux"? Remarquez que je ne suis pas doué en mathématiques - je ne fais que mettre en pratique ce que j'ai vu faire et que j'essaie de comprendre.
user13907

Réponses:

55

Ma réponse détaillée est ci-dessous, mais la réponse générale (c'est-à-dire réelle) à ce type de question est la suivante: 1) expérimentez, vissez, regardez les données, vous ne pouvez pas casser l'ordinateur, peu importe ce que vous faites, donc. . . expérience; ou 2) RTFM .

Voici un Rcode qui reproduit plus ou moins le problème identifié dans cette question:

# This program written in response to a Cross Validated question
# http://stats.stackexchange.com/questions/95939/
# 
# It is an exploration of why the result from lm(y_x+I(x^2))
# looks so different from the result from lm(y~poly(x,2))

library(ggplot2)


epsilon <- 0.25*rnorm(100)
x       <- seq(from=1, to=5, length.out=100)
y       <- 4 - 0.6*x + 0.1*x^2 + epsilon

# Minimum is at x=3, the expected y value there is
4 - 0.6*3 + 0.1*3^2

ggplot(data=NULL,aes(x, y)) + geom_point() + 
       geom_smooth(method = "lm", formula = y ~ poly(x, 2))

summary(lm(y~x+I(x^2)))       # Looks right
summary(lm(y ~ poly(x, 2)))   # Looks like garbage

# What happened?
# What do x and x^2 look like:
head(cbind(x,x^2))

#What does poly(x,2) look like:
head(poly(x,2))

Le premier lmretourne la réponse attendue:

Call:
lm(formula = y ~ x + I(x^2))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.53815 -0.13465 -0.01262  0.15369  0.61645 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.92734    0.15376  25.542  < 2e-16 ***
x           -0.53929    0.11221  -4.806 5.62e-06 ***
I(x^2)       0.09029    0.01843   4.900 3.84e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

Residual standard error: 0.2241 on 97 degrees of freedom
Multiple R-squared:  0.1985,    Adjusted R-squared:  0.182 
F-statistic: 12.01 on 2 and 97 DF,  p-value: 2.181e-05

La seconde lmretourne quelque chose d'étrange:

Call:
lm(formula = y ~ poly(x, 2))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.53815 -0.13465 -0.01262  0.15369  0.61645 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.24489    0.02241 144.765  < 2e-16 ***
poly(x, 2)1  0.02853    0.22415   0.127    0.899    
poly(x, 2)2  1.09835    0.22415   4.900 3.84e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

Residual standard error: 0.2241 on 97 degrees of freedom
Multiple R-squared:  0.1985,    Adjusted R-squared:  0.182 
F-statistic: 12.01 on 2 and 97 DF,  p-value: 2.181e-05

Puisque lmest le même dans les deux appels, il faut que les arguments lmsoient différents. Alors, regardons les arguments. Evidemment, yc'est pareil. Ce sont les autres parties. Regardons les premières observations sur les variables de droite lors du premier appel de lm. Le retour de head(cbind(x,x^2))ressemble à:

            x         
[1,] 1.000000 1.000000
[2,] 1.040404 1.082441
[3,] 1.080808 1.168146
[4,] 1.121212 1.257117
[5,] 1.161616 1.349352
[6,] 1.202020 1.444853

C'est comme prévu. La première colonne est xet la deuxième colonne est x^2. Que diriez-vous du deuxième appel de lm, celui avec poly? Le retour de head(poly(x,2))ressemble à:

              1         2
[1,] -0.1714816 0.2169976
[2,] -0.1680173 0.2038462
[3,] -0.1645531 0.1909632
[4,] -0.1610888 0.1783486
[5,] -0.1576245 0.1660025
[6,] -0.1541602 0.1539247

OK, c'est vraiment différent. La première colonne ne l'est pas xet la deuxième colonne ne l'est pas x^2. Donc, quoi poly(x,2)qu’il en soit , il ne revient pas xet x^2. Si nous voulons savoir ce qui se polypasse, nous pourrions commencer par lire son fichier d’aide. Alors on dit help(poly). La description dit:

Retourne ou évalue les polynômes orthogonaux de degré 1 à degré sur l'ensemble de points spécifié x. Ceux-ci sont tous orthogonaux au polynôme constant du degré 0. Vous pouvez également évaluer les polynômes bruts.

Maintenant, vous savez ce que sont les "polynômes orthogonaux" ou vous ne le savez pas. Si vous ne le faites pas, utilisez alors Wikipedia ou Bing (pas Google, bien sûr, car Google est diabolique - pas aussi grave qu'Apple, bien sûr, mais toujours aussi mauvais). Ou bien, vous pouvez décider que vous ne vous souciez pas de ce que sont les polynômes orthogonaux. Vous remarquerez peut-être l'expression "polynômes bruts" et vous remarquerez peut-être un peu plus loin dans le fichier d'aide qui polya une option rawqui est, par défaut, égale à FALSE. Ces deux considérations peuvent vous inciter à choisir les head(poly(x, 2, raw=TRUE))retours:

            1        2
[1,] 1.000000 1.000000
[2,] 1.040404 1.082441
[3,] 1.080808 1.168146
[4,] 1.121212 1.257117
[5,] 1.161616 1.349352
[6,] 1.202020 1.444853

Excité par cette découverte (ça a l'air bien, maintenant, oui?), Vous pouvez continuer d'essayer summary(lm(y ~ poly(x, 2, raw=TRUE))) Cela revient:

Call:
lm(formula = y ~ poly(x, 2, raw = TRUE))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.53815 -0.13465 -0.01262  0.15369  0.61645 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)              3.92734    0.15376  25.542  < 2e-16 ***
poly(x, 2, raw = TRUE)1 -0.53929    0.11221  -4.806 5.62e-06 ***
poly(x, 2, raw = TRUE)2  0.09029    0.01843   4.900 3.84e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

Residual standard error: 0.2241 on 97 degrees of freedom
Multiple R-squared:  0.1985,    Adjusted R-squared:  0.182 
F-statistic: 12.01 on 2 and 97 DF,  p-value: 2.181e-05

La réponse ci-dessus comporte au moins deux niveaux. Tout d'abord, j'ai répondu à votre question. Deuxièmement, et bien plus important encore, j’ai illustré comment vous êtes censé répondre vous-même à de telles questions. Chaque personne qui "sait programmer" a suivi une séquence semblable à celle décrite plus de soixante millions de fois. Même des gens aussi décevants en matière de programmation que moi parcourent cette séquence tout le temps. Il est normal que le code ne fonctionne pas. Il est normal de mal comprendre ce que font les fonctions. La façon de le gérer consiste à faire le tour, à expérimenter, à regarder les données et à RTFM. Sortez du mode "suivre une recette sans réfléchir" et passez en mode "détective".

Facture
la source
7
Je pense que cela mérite un +6. Je vais essayer de me rappeler dans quelques jours quand cela devient possible. FTR, je pense que cela n’a pas besoin d’être aussi sarcastique, mais cela montre bien ce que sont les polynômes orthogonaux / comment ils fonctionnent et montre le processus que vous utilisez pour résoudre de tels problèmes.
gung - Réintégrer Monica
13
Super réponse, merci. Bien que je sois un peu offensé par un "RTFM" (mais c'est peut-être juste moi): Le problème est que dans tout ce que j'ai lu, du moins en ce qui concerne la régression linéaire en R, les gens font parfois cela, d'autres le font. Franchement, je ne comprends pas l'entrée de Wikipedia sur les polynômes orthogonaux. Je ne vois pas pourquoi on utiliserait ceci pour la régression si les coefficients que vous obtenez sont "faux". Je ne suis pas mathématicien - j'essaie de suivre les recettes parce que je ne suis pas un cuisinier averti, mais je dois néanmoins manger quelque chose.
user13907
12
@ user13907, ce n'est pas que vous. C’est une bonne réponse qui mérite d’être votée à la hausse, mais il serait bon d’avoir un ton plus doux.
Waldir Leoncio
8
Vous n'avez pas vraiment besoin de comprendre quels polynômes orthogonaux sont ici - vous avez juste besoin de comprendre qu'ils ne sont pas ce que vous voulez. Pourquoi quelqu'un pourrait-il vouloir des polynômes orthogonaux? Soumettez cov (poly (x, 2)) pour constater que la covariance entre les deux termes du polynôme est égale à zéro (erreur maximale). C'est la propriété clé des polynômes orthogonaux - leurs termes ont une covariance nulle les uns avec les autres. Il est parfois pratique que vos variables RHS aient une corrélation nulle entre elles. Leurs coefficients ne sont pas faux, en réalité, ils doivent simplement être interprétés différemment.
Bill
2
Oh, d'accord, cette explication en anglais clair a maintenant du sens. Merci.
user13907
5

Il existe une approche intéressante pour l’interprétation de la régression polynomiale par Stimson et al. (1978) . Cela implique la réécriture

Y=β0+β1X+β2X2+vous

comme

Y=m+β2(F-X)2+vous

m=β0-β12/4β2β2F=-β1/2β2

Durden
la source
2
+1 Pour des analyses connexes, veuillez consulter stats.stackexchange.com/questions/28730 et stats.stackexchange.com/questions/157629 .
whuber
4

Si vous voulez juste un coup de pouce dans la bonne direction sans autant de jugement: poly()créez des polynômes orthogonaux (non corrélés), par opposition à I(), ce qui ignore complètement la corrélation entre les polynômes résultants. La corrélation entre les variables prédictives peut être un problème dans les modèles linéaires (voir ici pour plus d'informations sur les raisons pour lesquelles la corrélation peut être problématique), il est donc probablement préférable (en général) d'utiliser poly()plutôt que d' utiliser I(). Maintenant, pourquoi les résultats sont-ils si différents? Eh bien, à la fois poly()et I()prendre x et le convertir en un nouveau x (dans le cas I(), la nouvelle x est juste x ^ 1 ou x ^ 2, dans le cas poly(), sont beaucoup plus compliquées sont les nouveaux x (si vous voulez savoir d'où ils viennent (et vous n'avez probablement pas), vous pouvez commencerici ou la page Wikipedia susmentionnée ou un manuel). Le fait est que, lorsque vous calculez (prédisez) y en fonction d'un ensemble particulier de valeurs x, vous devez utiliser les valeurs x converties générées par l'un poly()ou l' autre I()(en fonction de celle qui figurait dans votre modèle linéaire). Alors:

library(ggplot2)    

set.seed(3)
epsilon <- 0.25*rnorm(100)
x       <- seq(from=1, to=5, length.out=100)
y       <- 4 - 0.6*x + 0.1*x^2 + epsilon

# Minimum is at x=3, the expected y value there is
4 - 0.6*3 + 0.1*3^2

ggplot(data=NULL,aes(x, y)) + geom_point() + 
   geom_smooth(method = "lm", formula = y ~ poly(x, 2))

modI <- lm(y~x+I(x^2)) 
summary(modI) # Looks right
modp <- lm(y ~ poly(x, 2))
summary(modp)  # Looks like garbage

# predict y using modI
coef(modI)[1] + coef(modI)[2] * 3^1 + coef(modI)[3] * 3^2

# predict y using modp
# calculate the new x values using predict.poly()
x_poly <- stats:::predict.poly(object = poly(x,2), newdata = 3)
coef(modp)[1] + coef(modp)[2] * x_poly[1] + coef(modp)[3] * x_poly[2]

Dans ce cas, les deux modèles renvoient la même réponse, ce qui suggère que la corrélation entre les variables de prédicteur n'influence pas vos résultats. Si la corrélation posait problème, les deux méthodes prédiraient des valeurs différentes.

filups21
la source
1

'poly' exécute l'ortho-normalisation de Graham-Schmidt sur les polynômes 1, x, x ^ 2, ..., x ^ deg Par exemple, cette fonction fait la même chose que 'poly' sans renvoyer bien sûr d'attributs 'coef'.

MyPoly <- 
function(x, deg)
{
    n <- length(x)
    ans <- NULL
    for(k in 1:deg)
    {
        v <- x^k
        cmps <- rep(0, n)
        if(k>0) for(j in 0:(k-1)) cmps <- cmps + c(v%*%ans[,j+1])*ans[,j+1]
        p <- v - cmps
        p <- p/sum(p^2)^0.5
        ans <- cbind(ans, p)
    }
    ans[,-1]
}

J'ai atterri sur ce fil parce que la forme fonctionnelle m'intéressait. Alors, comment pouvons-nous exprimer le résultat de 'poly' en tant qu'expression? Il suffit d'inverser la procédure Graham-Schmidt. Vous allez vous retrouver avec un désordre!

izmirlig
la source