Brz početak rada sa R: Složeniji regresioni modeli (4. Deo)

U trećem delu koristili smo komandu lm() da ocenimo regresioni model. U 4. Delu ćemo pogledati naprednije aspekte regresionih modela i videti šta R nudi.

Jedan od načina provere nelinearnosti u vašim podacima jeste da prilagodite polinomni model i proverite da li on odgovara podacima bolje od linearnog modela. Međutim, možda ćete takođe hteti da prilagodite kvadratni ili model višeg reda jer imate razloga da verujete da je odnos između promenljivih po prirodi polinomijalni. Da vidimo prvo kako prilagoditi kvadratni model u R.

Koristićemo datoteku sa brojevima (vreme raspada koji se odvija unutar izvora zračenja), merenih Gajgerovim brojačem u jednoj nuklearnoj elektrani. Tačke su registrovane tokom perioda od 30 sekundi za kratkotrajno veštačko radioaktivno jedinjenje. Učitavamo podatke i oduzmemo pozadinski broj od 623.4 brojanja u sekundi kako bismo dobili brojke koje se odnose na radioaktivni izvor. Isecite i zalepite sledeće podatke u radni prostor R.

gajger <- structure(list(Vreme = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30), Brojevi = c(750, 725.2, 695, 708.5, 725, 690.9, 691.5, 678.6, 686.3, 668.9, 665.3, 669.7, 657.5, 661.6, 665.1, 648.1, 664.9, 647.9, 660, 643, 646.2, 653, 646.9, 639.3, 638.7, 636.8, 650.2, 633.2, 642.2, 649.3, 642.7)), .Names = c("Vreme", "Brojevi"), row.names = c(NA, -31L), class = "data.frame")
gajger

Sada uzimamo u obzir pozadinski broj i oduzmimo ga od skupa podataka.

pozadinskibroj <- 623.4
gajger$Brojevi <- gajger$Brojevi - pozadinskibroj

Pridružićemo celu datoteku kako bismo se mogli pozivati na sve promenljive direktno po imenu.

attach(gajger)
names(gajger)

Prvo, specificirajmo linearni model, iako bi zaista trebalo prvo da se grafički prikažu podaci, a tek potom da se specificira odgovarajući regresioni model.

linearni.model <- lm(Brojevi ~ Vreme)
summary(linearni.model)


Model objašnjava preko 76% varijacija i ima veoma statistički značajne koeficijente. Hajde da prikažemo brojeve u odnosu na vreme na grafikonu i ucrtajmo i liniju koja predstavlja ocenjeni regresioni model.

plot(Vreme, Brojevi, pch = 16, ylab = "Brojevi u sekundi", cex.lab = 1.3)
abline(lm(Brojevi ~ Vreme))

Ovde je sintaksa cex.lab = 1.3 proizvela oznake osi odgovarajuće veličine.

Model izgleda dobro, ali možemo videti zakrivljenost na dijagramu rasturanja koja nije dobro objašnjena linearnim modelom. Sada ćemo prilagoditi kvadratni model. Prethodno kreiramo promenljivu koja se zove Vreme2 koja je kvadrat promenljive Vreme.

Vreme2 <- Vreme^2
kvadratni.model <- lm(Brojevi ~ Vreme + Vreme2)

Uočite sintaksu koja je korišćena u postavljanje linearnog modela sa dve ili više promenljive. Uključujemo svaku promenljivu i stavimo znak “+” između njih.
summary(kvadratni.model)

Naš kvadratni model je u osnovi linearni model sa dve promenljive, od kojih je jedan kvadrat drugog. Vidimo da je, koliko god je linearni model dobar, kvadratni model još bolje odgovara podacima, objašnjavajući dodatnih 14% varijacija. Sada ćemo grafički prikazati kvadratni model tako što ćemo na horizontalnoj osi uneti vrednosti vremena koje se kreću od 0 do 30 sekundi u koracima od 0.1s.

vrednostivremena <- seq(0, 30, 0.1)
prognoziranibrojevi <- predict(kvadratni.model, list(Vreme= vrednostivremena, Vreme2= vrednostivremena^2))
plot(Vreme, Brojevi, pch=16, xlab = "Vreme (s)", ylab = "Brojevi po sekundi", cex.lab = 1.3)

Dodajmo sada kvadratni model na grafikon pomoću naredbe lines().

lines(vrednostivremena, prognoziranibrojevi, col = "red", lwd = 3)

Izgleda da kvadratni model odgovara podacima bolje od linearnog modela. Međutim, možemo ići još i dalje. Pošto fizika atomskog raspada uključuje eksponencijalne odnose, treba proveriti da li eksponencijalna funkcija odgovara podacima čak i bolje od kvadratnih. To ćemo uraditi sledeći put.

Sve najbolje za sada!
David

Dodatak: Korišćeni R kodovi

# Učitavanje gajger datoteke u R radni prostor
gajger <- structure(list(Vreme = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30), Brojevi = c(750, 725.2, 695, 708.5, 725, 690.9, 691.5, 678.6, 686.3, 668.9, 665.3, 669.7, 657.5, 661.6, 665.1, 648.1, 664.9, 647.9, 660, 643, 646.2, 653, 646.9, 639.3, 638.7, 636.8, 650.2, 633.2, 642.2, 649.3, 642.7)), .Names = c("Vreme", "Brojevi"), row.names = c(NA, -31L), class = "data.frame")
gajger

# Dobijanje brojeva koji se odnose na izvor radioaktivnosti
# tako što će se oduzeti pozadinski broj
pozadinskibroj <- 623.4 
gajger$Brojevi <- gajger$Brojevi - pozadinskibroj

# Pridruživanje datoteke i pozivanje promenljivih po imenu
attach(gajger)
names(gajger) 

# Ocenjivanje linearnog regresionog modela
linearni.model <- lm(Brojevi ~ Vreme)
summary(linearni.model)

# Dijagram rasturanja i linearni regresioni model 
plot(Vreme, Brojevi, pch=16, ylab = "Brojevi po sekundi", cex.lab = 1.3)
abline(lm(Brojevi ~ Vreme))

# Ocenjivanje kvadratnog regresionog modela
Vreme2 <- Vreme^2
kvadratni.model <- lm(Brojevi ~ Vreme + Vreme2)
summary(kvadratni.model)

# Dijagram rasturanja i kvadratni regresioni model 
vrednostivremena <- seq(0, 30, 0.1)
prognoziranibrojevi <- predict(kvadratni.model, list(Vreme= vrednostivremena, Vreme2= vrednostivremena^2))
plot(Vreme, Brojevi, pch=16, xlab = "Vreme (s)", ylab = "Brojevi po sekundi", cex.lab = 1.3)
lines(vrednostivremena, prognoziranibrojevi, col = "red", lwd = 3)