David Lillis, Ph. D.

a 3.részben használtuk az lm() parancsot, hogy végre legkisebb négyzetek regressziók. A 4. részben megvizsgáljuk a regressziós modellek fejlettebb aspektusait,és megnézzük, mit kínál R.

az adatok nemlinearitásának ellenőrzésének egyik módja egy polinommodell illesztése, és annak ellenőrzése, hogy a polinommodell jobban illeszkedik-e az adatokhoz, mint egy lineáris modell. Előfordulhat azonban, hogy is szeretném, hogy illeszkedjen a másodfokú vagy magasabb modell, mert oka van, hogy nem hiszem, hogy a kapcsolat a változók eredendően polinom jellegű. Lássuk, hogyan illeszkedik a másodfokú modell R.

Fogunk használni, adatok halmaza számít egy változó, hogy csökken az idő múlásával. Vágja be az alábbi adatokat a R munkaterület.

A <- structure(list(Time = c(0, 1, 2, 4, 6, 8, 9, 10, 11, 12, 13, 14, 15, 16, 18, 19, 20, 21, 22, 24, 25, 26, 27, 28, 29, 30), Counts = c(126.6, 101.8, 71.6, 101.6, 68.1, 62.9, 45.5, 41.9, 46.3, 34.1, 38.2, 41.7, 24.7, 41.5, 36.6, 19.6, 22.8, 29.6, 23.5, 15.3, 13.4, 26.8, 9.8, 18.8, 25.9, 19.3)), .Names = c("Time", "Counts"),row.names = c(1L, 2L, 3L, 5L, 7L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 19L, 20L, 21L, 22L, 23L, 25L, 26L, 27L, 28L, 29L, 30L, 31L),class = "data.frame")

nézzük, csatolja a teljes adatkészlet úgy, hogy tudjuk: minden változók közvetlenül a név.

attach(A)
names(A)

először állítsunk fel egy lineáris modellt, bár valójában először ábrázolnunk kell, majd csak a regressziót kell végrehajtanunk.

linear.model <-lm(Counts ~ Time)

most részletes információkat kapunk a regressziónkról az összefoglaló() paranccsal.

summary(linear.model)Call:lm(formula = Counts ~ Time)Residuals: Min 1Q Median 3Q Max -20.084 -9.875 -1.882 8.494 39.445 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 87.1550 6.0186 14.481 2.33e-13 ***Time -2.8247 0.3318 -8.513 1.03e-08 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 15.16 on 24 degrees of freedomMultiple R-squared: 0.7512, Adjusted R-squared: 0.7408 F-statistic: 72.47 on 1 and 24 DF, p-value: 1.033e-08

a modell a variancia több mint 74%-át magyarázza, és rendkívül jelentős együtthatókkal rendelkezik a metszéspontra és a független változóra, valamint rendkívül jelentős általános modell p-értékkel. Azonban ábrázoljuk a számlálásokat az idő múlásával, és tegyük fel a lineáris modellünket.

plot(Time, Counts, pch=16, ylab = "Counts ", cex.lab = 1.3, col = "red" )abline(lm(Counts ~ Time), col = "blue")

tn_image001

tn_image001

itt a szintaxis cex.lab = 1.3 gyártott tengely címkék egy szép méretű.

a modell jól néz ki, de láthatjuk, hogy a telek görbülete nem magyarázható jól egy lineáris modellel. Most egy olyan modellt illesztünk be, amely időben másodfokú. Létrehozunk egy Time2 nevű változót, amely a változó idő négyzete.

Time2 <- Time^2
quadratic.model <-lm(Counts ~ Time + Time2)

vegye figyelembe a lineáris modell két vagy több prediktorral történő illesztésének szintaxisát. Minden előrejelzőt beleszámítunk, és pluszjelet teszünk közéjük.

summary(quadratic.model)Call:lm(formula = Counts ~ Time + Time2)Residuals: Min 1Q Median 3Q Max -24.2649 -4.9206 -0.9519 5.5860 18.7728 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 110.10749 5.48026 20.092 4.38e-16 ***Time -7.42253 0.80583 -9.211 3.52e-09 ***Time2 0.15061 0.02545 5.917 4.95e-06 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 9.754 on 23 degrees of freedomMultiple R-squared: 0.9014, Adjusted R-squared: 0.8928F-statistic: 105.1 on 2 and 23 DF, p-value: 2.701e-12

másodfokú modellünk lényegében lineáris modell két változóban, amelyek közül az egyik a másik négyzete. Látjuk, hogy bármennyire is jó volt a lineáris modell, a másodfokú modell még jobban teljesít, a variancia további 15% – át magyarázva. Most ábrázoljuk a másodfokú modellt úgy, hogy beállítunk egy 0-30 másodperces időérték-rácsot 0,1 másodperces lépésekben.

timevalues <- seq(0, 30, 0.1)
predictedcounts <- predict(quadratic.model,list(Time=timevalues, Time2=timevalues^2))
plot(Time, Counts, pch=16, xlab = "Time (s)", ylab = "Counts", cex.lab = 1.3, col = "blue")

most a másodfokú modellt a lines() paranccsal felvesszük a cselekménybe.

lines(timevalues, predictedcounts, col = "darkgreen", lwd = 3)

tn_image002

úgy tűnik, hogy a másodfokú modell jobban illeszkedik az adatokhoz, mint a lineáris modell. A következő blogbejegyzésünkben újra megvizsgáljuk az ívelt modellek illesztését.

lásd a teljes R bemutató sorozat és egyéb blogbejegyzések kapcsolatos R programozás.

A szerzőről: David Lillis tanított R sok kutató és statisztikus. Cége, a Sigma Statistics and Research Limited on-line oktatást és szemtől-szembe workshopokat biztosít az R-ről,valamint a kódolási szolgáltatásokat R-ben.

 könyvjelző és megosztás

Vélemény, hozzászólás?

Az e-mail-címet nem tesszük közzé.