# Linearna regresija


# (1) 
# Prosta linearna regresija i ocene parametara.


# (2)
# # protok saobracaja u hiljadama automobila u toku dana
x=c(8.3, 9.3, 12.1, 12.3, 17.0, 17.3, 24.3, 24.5, 33.6) 
# prisustva olova u kori drveca:
y=c(227, 312, 362, 521, 539, 728, 945, 1000, 1263)
# Oceniti koeficijente modela y=a*x+b

# primenom linearne regresije sa srednje kvadratnom greskom:
model=lm(y~x)
summary(model)
plot(model)

a=41.318
b=-73.356

plot(x, y)
# regresiona prava
abline(b,a, col="red")
new=28
# predvidjanje na osnovu ovog modela za vrednosti x=5, 13.3, 28, 30 i 100.
predict(model, newdata=data.frame(x=c(5, 13.3, 28, 30, 100)))

# (3)
# broj sati ucenja
x=c(4,9, 10, 14, 4, 7, 12, 22, 1, 3, 8, 11, 5, 6, 10, 11, 16, 13, 13, 10)
# rezultati testa.
y=c(390, 580, 650, 730, 410, 530, 600, 790, 350, 400, 590, 640,450, 520, 690, 690, 770,700, 730, 640)


model=lm(y~x)
summary(model)
plot(x, y)
# regresiona prava
a=25.326
b=353.165
abline(b,a, col="red")
# Vidimo da postoji skoro da linearna zavisnost medju brojem sati ucenja i rezultatima.

reziduali=model$residuals # zelimo da ovo pripada nekoj normalnoj raspodeli.

# (4)
q=seq(-1, 1, 0.05)
w=q^2+rnorm(length(q),0, 0.05)
w
[1]  1.078499039  0.884586722  0.854204947  0.756695964  0.693174977  0.586613777
[7]  0.426887412  0.493764092  0.304910978  0.332567489  0.350705839  0.185258685
[13]  0.178543423  0.116121110  0.201922674  0.095142984 -0.007614178  0.078469689
[19]  0.020015986 -0.017067335  0.128702875  0.010107586  0.014829872  0.056055063
[25] -0.080430093  0.076724353  0.075218465  0.044692828  0.172234174  0.271440680
[31]  0.258209703  0.293491675  0.343640829  0.434966424  0.477094938  0.625986227
[37]  0.604295805  0.713034165  0.770949824  0.856426849  1.063928001
# naci regresionu pravu za w u odnosu na q.

# 
model=lm(w~q)
summary(model)
# Vidimo da nemamo "zvezdice" za q - nije znacajna, i p-vrednost je blizu 1
# a u ovom slucaju sto je manja to je bolja. Dok vrednosti za R-squared treba 
# da budu sto blizi jedinici.
plot(q, w)
a=   model$coefficients[2]
b= model$coefficients[1]
abline(b,a, col="red")
# Ovakve stvari se desavaju iz prostog razloga sto veza nije linearna.

# Probamo kvadratni model:
model2=lm(w~I(q)+I(q^2))
summary(model2)
# Ovde je od znacaja samo parametar uz q^2, sto i ima smisla skroz.
plot(w~q)
points(q, fitted(model2), col="red",pch=20)


# (4)
# # Nadalje bice nam potreban paket MASS, inace je od velikog znacaja, mnogo razlicitih testova, nacina ocenjivanja i slicnog se nalazi bas tu.
library(MASS)

# I pokusacemo da instaliramo paket ISLR (koji detaljno prati knjigu Hastie-ja i Tibshirani-ja).
install.packages("ISLR")
library(ISLR)

# Obicna linearna regresija na primeru baze- Boston.
Boston

# funkcija "names" vraca nazive promenljivih u datoj bazi.
names(Boston)

# ?Boston=help(Boston), tj opisuje sta predstavlja ta baza. Vezana je za ucestalost krsenja zakona u Bostonu. 
?Boston

# plotujemo sledece i dobijamo grafik raspsenosti promenljive medv u odnosu na lstat. 
plot(medv~lstat,Boston)

# Sledeci poziv se vec odnosi na linearni model, funkcija lm = linear model vraca najbolji linearni model koji odgovara datim podacima.
fit1=lm(medv~lstat,data=Boston)
fit1
# Od znacaja su nam ocene koeficijenata.
# Da bismo detaljnije videli kako i sta se desava, pozovemo summary().
summary(fit1)
# I vidimo da su obe promenljive znacajne.

# Sledecim pozivom mozemo da dodamo regresionu pravu na grafik, a ne samo kao sto to radili na prvom casu.
abline(fit1,col="red")
names(fit1)

# Ovde mozemo da vidimo intervali poverenja za ocene nasih koeficijenata modela.
confint(fit1)

# funkcija pairs() vraca sve moguce grafike rasprsenosti nase baze. Tj. pomocu nje mozemo ponekad i da ustanovimo gde odmah imamo zgodnu linearnu zavisnost
# i samim tim da primenimo linearnu regresiju. Jedini je problem, sto je zgodno predstaviti samo po 2 promenljive na grafiku,
# vec zavinost od 3 promenljive
# nije lako i lepo uocljiva na grafiku, a vece dimenzije cak i ne mozemo da nacrtamo. No, svejedno mozemo pokretati linearne regresije.

### Nadalje, linearna regresija se moze koristiti i na vise promenljivih.
fit2=lm(medv~lstat+age,data=Boston)
summary(fit2)

# Pri ovom zapisu znaci da smo uzeli u obzir sve moguce promenljive koje ima baza.
fit3=lm(medv~.,Boston)
summary(fit3)

# funkcija update(model, ...) menja model na zadati nacin. Prosto da ne bismo uvek pravili od nule.
fit4=update(fit3,~.-age-indus)
summary(fit4)

# I ne mora samo plus da bude izmedju prediktora, moze se javiti, recimo i proizvod: 
### NAPOMENA: ako stavimo x*y, tada se u model ukljucuje sledeci prediktori: x, y i x*y(koji se u pozvanom summary oznacava kao x:y).
# Ako bismo zeleli da pozovemo samo za prediktor x*y, tada bismo morali da pokrenemo:
# lm(z~ I(x*y)) ili u nasem slucaju fit=lm(medv~I(lstat*age))
fit5=lm(medv~lstat*age,Boston)
summary(fit5)

# Kad stavimo I(x^2) - time trazimo da x^2 posmatra totalno nezavisno od x. 
fit6=lm(medv~lstat +I(lstat^2),Boston); summary(fit6)
attach(Boston)
plot(medv~lstat)
points(lstat,fitted(fit6),col="red",pch=20)
# vidimo da je kvadratna regresija ovde bila preciznija od obicne. 

# Hocemo da probamo sa vecim stepenom mogli smo u fit6 da stavimo poly(lstat, 2) dobicemo isto.
fit7=lm(medv~poly(lstat,4))
points(lstat,fitted(fit7),col="blue",pch=20)
# Vidimo da je neznacajno bolja, ali zato mozda overfituje (preprilagodjava) model, sto svakako bolje da izbegnemo.
summary(fit7)