Mitt detaljerade svar är nedan, men det allmänna (dvs. verkliga) svaret på denna typ av fråga är: 1) experimentera, skruva runt, titta på data, du kan inte bryta datorn oavsett vad du gör, så . . . experimentera; eller 2) RTFM.

Här är någon R -kod som replikerar problemet som identifierats i denna fråga, mer eller mindre:

  # Detta program skrivet som svar på en korsvaliderad fråga # http://stats.stackexchange.com/questions/95939/# # Det är en undersökning av varför resultatet från lm (y_x + I (x ^ 2 )) # ser så annorlunda ut än resultatet från lm (y ~ poly (x, 2)) bibliotek (ggplot2) epsilon <- 0,25 * rnorm (100) x <- seq (från = 1, till = 5, längd. ut = 100) y <- 4 - 0,6 * x + 0,1 * x ^ 2 + epsilon # Minimum är vid x = 3, det förväntade y-värdet är4 - 0,6 * 3 + 0,1 * 3 ^ 2ggplot (data = NULL, aes ( x, y)) + geom_point () + geom_smooth (metod = "lm", formel = y ~ poly (x, 2)) sammanfattning (lm (y ~ x + I (x ^ 2))) # Ser rättighetsmätning (lm (y ~ poly (x, 2))) # Ser ut som sopor # Vad hände? # Hur ser x och x ^ 2 ut: huvud (cbind (x, x ^ 2)) # Hur ser poly (x, 2) ser ut som: head (poly (x, 2))  

Den första lm returnerar det förväntade svaret:

  Ring : lm (formel = y ~ x + I (x ^ 2)) Rester: Min 1Q Median 3Q Max -0,53815 -0,13465 -0,01262 0,15369 0,61645 Koefficienter: Uppskattning Std. Fel t värde 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. koder: 0 '***' 0,001 '**' 0,01 '*' 0,05 '.' 0,1 '' 1 Restfel kvarvarande: 0,2241 på 97 frihetsgrader Flera R-kvadrat: 0,1985, Justerat R-kvadrat: 0,182 F-statistik: 12.01 på 2 och 97 DF, p-värde: 2.181e-05  

Den andra lm returnerar något udda:

  Samtal: lm (formel = y ~ poly (x, 2)) Rester: Min 1Q Median 3Q Max
-0.53815 -0.13465 -0.01262 0.15369 0.61645 Koefficienter: Uppskattning Std. Fel t värde 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. koder: 0 '***' 0,001 '**' 0,01 '*' 0,05 '.' 0,1 '' 1 Restfel kvarvarande: 0,2241 på 97 frihetsgrader Flera R-kvadrat: 0,1985, Justerat R-kvadrat: 0,182 F-statistik: 12.01 på 2 och 97 DF, p-värde: 2.181e-05  

Eftersom lm är densamma i de två samtalen måste det vara argumenten för lm som är olika. Så, låt oss titta på argumenten. Uppenbarligen är y samma. Det är de andra delarna. Låt oss titta på de första observationerna på högervariablerna i det första anropet av lm . Returen för -huvud (cbind (x, x ^ 2)) ser ut som:

  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  

Detta är som förväntat. Första kolumnen är x och den andra kolumnen är x ^ 2 . Vad sägs om det andra samtalet av lm , det med poly? Returen av -huvud (poly (x, 2)) ser ut som:

  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, det är verkligen annorlunda. Den första kolumnen är inte x , och den andra kolumnen är inte x ^ 2 . Så vad poly (x, 2) gör, returnerar det inte x och x ^ 2 . Om vi ​​vill veta vad poly gör kan vi börja med att läsa dess hjälpfil. Så vi säger hjälp (poly) . Beskrivningen säger:

Returnerar eller utvärderar ortogonala polynom av grad 1 till grad över den angivna uppsättningen punkter x. Dessa är alla ortogonala mot den konstanta polynomgraden 0. Alternativt kan du utvärdera råa polynomier.

Nu vet du antingen vad "ortogonala polynomier" är eller inte. Om du inte gör det, använd sedan Wikipedia eller Bing (naturligtvis inte Google för att Google är ont --- inte lika illa som Apple, naturligtvis, men fortfarande dåligt). Eller så kan du bestämma att du inte bryr dig om vad ortogonala polynom är. Du kanske märker frasen "råa polynom" och du kanske märker lite längre ner i hjälpfilen att poly har ett alternativ råkod som standard är lika med FALSE . Dessa två överväganden kan inspirera dig att testa head (poly (x, 2, raw = TRUE)) som returnerar:

  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  

Upphetsad av denna upptäckt (det ser rätt ut, nu, ja?), kanske du försöker sammanfattning (lm (y ~ poly (x, 2, rå = SANT)) Detta returnerar:

  Samtal: lm (formel = y ~ poly (x, 2, rå = SANT)) Rester: Min 1Q Median 3Q Max -0,53815 -0,13465 -0,01262 0,15369 0,61645 Koefficienter: Uppskattning Std. Fel t värde 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, rå = SANT) 2 0,09029 0,01843 4,900 3,84e-06 *** --- Signif. koder: 0 '***' 0,001 '**' 0,01 '*' 0,05 '.' 0,1 '' 1 Restfel kvarvarande: 0,2241 på 97 frihetsgrader Flera R-kvadrat: 0,1985, Justerat R-kvadrat: 0,182 F-statistik: 12.01 på 2 och 97 DF, p-värde: 2.181e-05  

Det finns minst två nivåer i ovanstående svar. Först svarade jag på din fråga. För det andra, och mycket viktigare, illustrerade jag hur du ska gå att svara på frågor som detta själv. Varje enskild person som "vet hur man programmerar" har gått igenom en sekvens som den ovanför sextio miljoner gånger. Även människor som är så deprimerande dåliga att programmera som jag går igenom denna sekvens hela tiden. Det är normalt att koden inte fungerar. Det är normalt att missförstå vad funktioner gör. Sättet att hantera det är att skruva runt, experimentera, titta på data och RTFM. Gå ut ur "tanklöst att följa ett recept" -läge och in i "detektiv" -läge.

Om du bara vill ha en knuff i rätt riktning utan riktigt så mycket omdöme: poly () skapar ortogonala (inte korrelerade) polynom, i motsats till I () , som helt ignorerar korrelationen mellan de resulterande polynomema. Korrelation mellan prediktorvariabler kan vara ett problem i linjära modeller (se här för mer information om varför korrelation kan vara problematisk), så det är nog bättre (i allmänhet) att använda poly () kod> istället för I () .

Varför ser resultaten så annorlunda ut nu? Tja, både poly () och I () tar x och konverterar det till en ny x. När det gäller I () är den nya x bara x ^ 1 eller x ^ 2. När det gäller poly () är de nya x: erna mycket mer komplicerade. Om du vill veta var de kommer ifrån (och du antagligen inte gör det) kan du komma igång här eller den ovan nämnda Wikipedia-sidan eller en lärobok.

Poängen är att när du beräknar (förutsäger) y baserat på en viss uppsättning x-värden, måste du använda de konverterade x-värdena producerade av antingen poly () eller I () (beroende på vilken som fanns i din linjära modell). Så:

  -bibliotek (ggplot2)

# ställ in fröet så att resultaten kan reproduceras.
set.seed (3)

#### simulera lite data ####
# epsilon = slumpmässigt feluttryck
epsilon <- 0,25 * rnorm (100)
# x-värden är bara en sekvens från 1 till 5
x <- seq (från = 1, till = 5, längd ut = 100)
# y är en polynomfunktion av x (plus något fel)
y <- 4 - 0,6 * x + 0,1 * x ^ 2 + epsilon

# Minimum är vid x = 3, det förväntade y-värdet finns
4 - 0,6 * 3 + 0,1 * 3 ^ 2

# visualisera data (med en polynom-anpassad linje)
ggplot (data = NULL, aes (x, y)) + geom_point () +
   geom_smooth (metod = "lm", formel = y ~ poly (x, 2))
 

  #### Modellera data ####
# först försöker vi modellera data med bara jag ()
modI <- lm (y ~ x + I (x ^ 2))
# modellöversikten ser rätt ut
sammanfattning (modI)
# nästa försöker vi med poly ()
modp <- lm (y ~ poly (x, 2))
# modellöversikten ser konstig ut
sammanfattning (modp)

#### gör förutsägelser vid x = 3 baserat på varje modell ####
# förutsäga y med modI
# allt vi behöver göra är att ta modellkoefficienterna och ansluta dem till formeln: fånga + beta1 * x ^ 1 + beta2 * x ^ 2
coef (modI) [1] + coef (modI) [2] * 3 ^ 1 + coef (modI) [3] * 3 ^ 2
 

(Intercept)
3.122988

  # förutsäga y med modp
# detta tar ett extra steg.
# först, beräkna de nya x-värdena med hjälp av predict.poly ()
x_poly <- stats ::: förutsäga.poly (objekt = poly (x, 2), newdata = 3)
# använd sedan samma formel som ovan, men den här gången istället för det råa x-värdet (3), använd det polynomiska x-värdet (x_poly)
coef (modp) [1] + coef (modp) [2] * x_poly [1] + coef (modp) [3] * x_poly [2]
 

(Intercept)
3.122988

I det här fallet returnerar båda modellerna samma svar, vilket antyder att korrelation mellan prediktorvariabler inte påverkar dina resultat. Om korrelation var ett problem skulle de två metoderna förutsäga olika värden.

Det finns ett intressant synsätt på tolkning av polynomregression av Stimson et al.(1978).Det handlar om omskrivning

$ Y = \ beta_ {0} + \ beta_ {1} X + \ beta_ {2} X ^ {2} + u $

som

$ Y = m + \ beta_ {2} \ left (f - X \ right) ^ {2} + u $

där $ m = \ beta_ {0} - \ vänster.\ beta_ {1} ^ {2} \ right / 4 \ beta_ {2} $ är lägsta eller högsta (beroende på tecknet på $ \ beta_ {2} $) och $ f = \ left.- \ beta_ {1} \ right / 2 \ beta_ {2} $ är brännvärdet.Det förvandlar i grunden den tredimensionella kombinationen av backar till en parabel i två dimensioner.Deras artikel ger ett exempel från statsvetenskap.

'poly' utför Graham-Schmidt-orto-normalisering på polynomerna 1, x, x ^ 2, ..., x ^ deg Till exempel gör denna funktion samma sak som "poly" utan att självklart returnera "coef" -attribut.

  MyPoly <-
funktion (x, deg)
{
    n <- längd (x)
    ans <- NULL
    för (k i 1: deg)
    {
        v <- x ^ k
        cmps <- rep (0, n)
        om (k>0) för (j i 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]
}
 

Jag landade på den här tråden för att jag var intresserad av den funktionella formen.Så hur uttrycker vi resultatet av "poly" som ett uttryck?Invertera bara Graham-Schmidt-förfarandet.Du får en röra!