x = c(8.3, 9.3, 12.1, 12.3, 17.0, 17.3, 24.3, 24.5, 33.6) #protok saobracaja [1000 automobila/dan]
y = c(227, 312, 362, 521, 539, 728, 945, 1000, 1263) #prisustvo olova u kori drveta
n = length(x)

plot(x,y,xlim = c(-5, 40), ylim = c(-100, 1300)) #xlim i ylim odredjuju raspone osa

# Pravljenje modela linearne regresije je jako jednostavno
model = lm(y ~ x)
# U pozadini se optimizuju koeficijenti modeli, minimozovanjem srednjekvadratne greske.

summary(model) #ispis koeficijenata, R^2, kao i t-statistike za hipotezu beta_i = 0
# Interpretacija: Slobodan clan je -73.36, sto znaci da u startu imamo toliku
# koncentraciju olova u kori, a koeficijent uz x je 41.32, pa se za
# svakih 1000 automobila koncentracija toliko poveca.

# Napomena: R^2 je mera kvaliteta regresionog modela koja pripada intervalu
# [0,1]. Vece vrednosti su bolje - detaljnije na ESU/LSM.


abline(model, col = "red") #dodavanje modela na grafik, abline crta prave

beta_0 = model$coefficients[1]
beta_1 = model$coefficients[2]
plot(x,y)
abline(beta_0, beta_1, col = "red") #standardna sintaksa za abline

fitted_values = model$fitted.values #pristup predvidjenim vrednostima modela
plot(x,y)
points(x, fitted_values, col = "red") #sa points mozemo da ih docrtamo na grafik
points(x, fitted_values, col = "red", type = "l") #type = "l", povezuje tacke linijama
#dobar nacin da nacrtamo model nezavisno od tipa, dok abline moze samo prave da crta

# Napomena: Nazive parametara modela, coefficients i fitted.values mozemo da
# nadjemo preko ugradjenog helpa, ?lm


# Predvidjanje nad novim vrednostima:
# Koristimo funckiju predict() za skoro sve modele u R-u.
# Kljucno je da se podaci spakuju u data.frame i da se kolone isto zovu kao pri
# pozivu modela, posto R prediktore pronalazi po nazivu.

x_new = c(7, 15, 40)
predict(model, data.frame(x = x_new))

# Primer: Dodajemo kolonu koja nicemu ne sluzi.
z_new = c(100,3123,615)
predict(model, data.frame(x = x_new, z = z_new)) 
# z_new nicemu ne sluzi, ali ne smeta, model prepoznaje prediktore po nazivu


plot(x,y,xlim = c(0, 40), ylim = c(0, 1300))
# Najbolji nacin da se docrta model:
# Pravimo 'resetku' vrednosti od -10 do 40 i vrsimo predvidjanje modela za sve
# vrednosti:
xs = seq(0, 40, 0.01)
y_pred = predict(model, data.frame(x = xs))
lines(xs, y_pred, col = 'red')


# Pri kreaciji modela, izraz y~x se u R-u naziva formulom.
# Sa leve strane formule je ciljan promenljiva, a sa desne stavljamo prediktore
# Ako zelimo da napravimo model sa vise prediktora, odvajamo ih oznakom +.
# Npr. ako imamo prediktore x1 i x2, i zelimo oba da ukljucimo u model, pisali
# bismo y ~ x1+x2 pri pozivu modela.
# Oznakom - mozemo eksplicitno da izbacujemo prediktore:

# Nije idealno sto model predvidja negativne vrednosti za koncentraciju olova.
# Ovo mozemo da izbegnemo tako sto napravimo model bez slobodnog clana, tj.
# model sa regresionom funkcijom f(x) = beta_1*x.
# Model bez slobodnog clana pravimo tako sto ga 'oduzmemo' iz spiska prediktora
# sa desne strane formule dodavanjem '-1'.
model1 = lm(y ~ x - 1)
summary(model1)
abline(model1, col = 'blue') # Prolazi kroz koordinatni pocetak



# Sta kada zavisnost izmedju ciljne promenljive i prediktora nije linearna?
q = seq(-1, 1, 0.02)
w = q^2 + rnorm(length(q), 0, 0.02)
plot(q,w)

model = lm(w ~ q)
summary(model) #jako los model, visoka p-vrednost t-testa, R^2 prakticno 0
abline(model, col = "red")

# Nove prediktore mozemo i da pravimo izmenom postojecim. Kada zelimo da vrsimo
# matematicke operacije nad prediktorima unutar izraza za formulu modela, te
# operacije moraju da budu sadrzane unutar takozvanog operatora izolacije - I().
# Pravimo model koji uzima u obzir i kvadrat prediktora q:
model2 = lm(w ~ q + I(q^2)) #dodajemo kvadratni clan u model
summary(model2) #sad mnogo bolje
abline(model2, col = "blue") # ne moze ovako, izbacuje i warning
fitted_values = model2$fitted.values
plot(q,w)
points(q, fitted_values, type = "l", col = "blue")



# Zadatak:
# Baza E1.19 iz paketa SenSrivastava sadrzi podatke o ceni, broju stranica, i
# ukoricenosti raznih knjiga. Cilj je da napravimo model koji sto bolje opisuje,
# odnosno predvidja, cenu knjige na osnovu ostalih parametara.
library(SenSrivastava)
data(E1.19)

# Razliciti nacini da se dobiju razni podaci o bazi podataka:
?E1.19 # Dokumentacija
head(E1.19) # Ispis prvih nekoliko redova
str(E1.19) # Struktura baze - Dobar nacin da se vide tipovi obelezja 
summary(E1.19) # Ispis sumarnih statistika za svako obelezje
View(E1.19) # Prikaz cele baze u novom prozoru

# Pravimo model cene u odnosu na broj stranica (P) i ukoricenost (B).
model_books = lm(Price ~ P + B, data = E1.19)
summary(model_books)
# Vidimo prediktor koji se naziva Bp - on ustvari oznacava indikator B==b, tj.
# da li je knjiga neukoricena.
# Interpretacija modela: Slobodan clan je jednak 20.6 (dolara) - to je osnovna
# cena knjige. Za svaku novu stranicu cena raste za 0.027. Ako knjiga nije
# ukoricena, cena se smanji za 13.04.

# Napomena: B je faktor, ali ga lm u pozadini tretira isto kao numericku
# promenljivu sa vrednostima 0 i 1, za svrhe racuna.
# Pretvoricemo B u numericki tip da se uverimo:
as.numeric(E1.19$B) - 1 # -1 jer inace pretvara 'c' u 1 i 'p' u 2.
# Drugi nacin koji je citljiviji:
ifelse(E1.19$B == 'p', 1 , 0)

E1.19$B = ifelse(E1.19$B == 'p', 1 , 0)
model_books_num = lm(Price ~ B + P, data = E1.19)
summary(model_books_num) # isti model

# Zelimo da rekodiramo nule i jedinice - Prirodnije je da ukoricenost bude
# obelezeno sa 1, jer je dodatna usluga koja dize cenu knjige.
E1.19$B = 1 - E1.19$B

model_books_num_naopako = lm(Price ~ B + P, data = E1.19)
summary(model_books_num_naopako)
# Isti model sustinski, ali drugacije zapisan. Sada je osnovna cena 7.56, a
# ukoricenje dodaje 13.04 na cenu.


# Kako da vizuelizujemo model sa vise promenljivih? Obicno ne mozemo, ali
# ovaj mozemo.
# Posto postoji samo jedna numericka promenljiva - Nju cemo da nacrtamo na x
# osi, a regresionu pravu cemo da nacrtamo posebno za svaku kategoriju.
plot(E1.19$P, E1.19$Price)
#lines(E1.19$P, model$fitted.values, col = 'red') lose - razmisliti zasto
xs = seq(100, 500, 1)
df_paperback = data.frame(P = xs, B = rep(0, length(xs)))
y_pred_paperback = predict(model_books_num_naopako, df_paperback)
lines(xs, y_pred_paperback, col = 'red')
df_cloth = data.frame(P = xs, B = rep(1, length(xs)))
y_pred_cloth = predict(model_books_num_naopako, df_cloth)
lines(xs, y_pred_cloth, col = 'blue')
# Plava prava predstavlja predvidjanje za ukoricene, a crvena za neukoricene
# knjige.

# Napravljeni model predvidja isti rast cene po stranici za ukoricene i
# neukoricene knjige. Mozda je realnije da cena za ukoricene raste brze, ili
# uopste pod drugacijim koeficijentom. Ovo mozemo da postignemo tako sto
# dodajemo novi prediktor koji predstavlja proizvod indikatora ukoricenosti
# i broja stranica. On ce aktivan samo kada je knjiga ukoricena, pa ce za
# ukoricene knjige dodatno da se naplacuje po stranici.
model_int = lm(Price ~ P + B + I(B*P), data = E1.19)
summary(model_int) # dodatni koeficijent

plot(E1.19$P, E1.19$Price)
#lines(E1.19$P, model$fitted.values, col = 'red')
xs = seq(100, 500, 1)
df_paperback = data.frame(P = xs, B = rep(0, length(xs)))
y_pred_paperback = predict(model_int, df_paperback)
lines(xs, y_pred_paperback, col = 'red')
df_cloth = data.frame(P = xs, B = rep(1, length(xs)))
y_pred_cloth = predict(model_int, df_cloth)
lines(xs, y_pred_cloth, col = 'blue')
# Plava prava ima drugaciji nagib.

# Laksi nacin da se postigne ista stvar je da se u formuli modela izvrsi
# 'mnozenje' predikora umesto sabiranja:
model_int = lm(Price ~ P*B, data = E1.19)
summary(model_int) # Isti model