x <- 2*rnorm(50) y <- x^2 - 2*x + rnorm(50) plot(x,y) reg1 <- lm(y ~ x) abline(reg1$coef) plot(x,reg1$resid) plot(reg1$fitted,reg1$resid) # a hibák parabolát követnek! # regresszióba kéne négyzetes tag! reg2 <- lm(y ~ x + I(x^2)) # I(): identitás, azért kell, mert lm() # hasában mást jelentenek a műveleti jelek! reg2 plot(x,y) points(x,reg2$fitted,col="red") lines(x,reg2$fitted,col="red") # ha x nem növekvő sorrendben van, keszekusza... # növekvő rendben kéne x, de ezzel együtt át kell # rendezni a fitted vector-t is! sort(x) order(x) x[order(x)] plot(x,y) points(x[order(x)],y[order(x)],col=2,pch=2) lines(x[order(x)],y[order(x)],col=2) plot(x,y) lines(x[order(x)],reg2$fitted[order(x)],col=2) # ha új alappontokon akarjuk ezt a becsült függvényt: t <- seq(-6,6,0.2) # predict függvény: illesztett regressziós függvényt # új alappontokban képes kiszámolni # háklis: új pontok data.frame-ben kellenek neki, # és az oszlopok neve a becslő változók: most x f <- predict(reg2,data.frame(x=t)) plot(x,y) lines(t,f) plot(x,reg2$resid) plot(reg2$fitted,reg2$resid) summary(reg2) # t-teszt szerint: konstans tag nagy eséllyel 0 # új modell, konstans nélkül reg3 <- lm(y ~ x + I(x^2) + 0) # y = bx + cx^2 + hiba plot(x,y) lines(x[order(x)],reg3$fitted[order(x)]) lines(x[order(x)],reg2$fitted[order(x)],col="blue") # grfikusan: nagyon közel vna eymáshoz a két modell/ # két függvény # statisztikai összehasonlítás: anova # analysis of variance anova(reg2,reg3) # RSS: residual sum of squares: hibák négyzetösszege # F = F-statisztika, Pr(..)=p-value: F-teszt: # H0: gyakorlatilag ekvivalens a két modell # H1: lényeges a különbség # nagy p-value: H0 -- ilyenkor az egyszerűbb modellt # használjuk/preferáljuk; kis p-value: H1: megéri # a komplexebb modellt használni anova(reg1,reg2) # anova() csak akkor van értelme, ha az egyik modell # részmodellje a másiknak: néhány változó # elhagyásával kaptuk # reg1, reg3 nem hasonlíthatók össze: mindkettőben # van olyan változó, ami a másikban nincs. # több változós előrejelzés: a <- rnorm(30) b <- rnorm(30) c <- 2*a - b + rnorm(30) reg4 <- lm(c ~ a + b) reg4 # itt is lehet szorzat, magasabb rendű tag(ok), ... reg5 <- lm(c ~ a + b + I(a*b) + I(a^2))