#####
# INTRO
#####
# vektor je min element u r
x <- 5.5
x
x[1]
class(x)
# numericki vektor
# integer vektor
y <- 5L
y
class(y)
z <- 'HELLO'
class(z)
z == "HELLO"
# logicke vrednosti
l <- TRUE
class(l)
# creating vectors with more elements
vektor1 <- c(1.5, 3, 4.5, 6, 7)
vektor2 <- c(1.5, 18, -41, -181)
class(vektor1)
class(vektor2)
print(vektor1)
summary(vektor1)
summary(vektor2)
# vektorska aritmetika
# +, -, *, /, ^
vektor1 + 1
class(vektor1 + 1)
summary(vektor1 + 1)
# get the max/min value
max(vektor1)
min(vektor1)
# get the range of values
help("range")
range(vektor1)
range(vektor1)[2]
# calculate the sum of all elements
sum(vektor1)
# calculate the product of all elements
prod(vektor1)
###
# Generating regular sequences
###
a = 1:10
a <- 1:10 #ispravnija notacija u R-u
a
b <- rev(a)
b
# generate a sequence from 3.2 to 4.7, with a step 0.2
sekvenca2 <- seq(from = 3.2, to = 4.7, by = 0.2)
sekvenca3 <- seq(1,100,2.5)
table(sekvenca3)
length(sekvenca3)
###
# Factor values
###
# create vector and convert to factor
vektor3 <- c('cold', 'mild', 'mild', 'hot', "cold")
help('c')
help('table')
class(vektor3)
print(vektor3)
faktor1 <- as.factor(vektor3)
print(faktor1)
# print jedinstvene vrednosti tj nivoe faktora1
levels(faktor1)
print(levels(faktor1))
# specify the order of factor levels
help('factor')
faktor2 <- factor(vektor3, levels = c("cold", "mild", "hot"),
labels = c("bolekurcina", "ladokurcina", "kurcina"))
print(faktor2)
print(faktor1)
faktor3 <- factor(vektor3,levels = c("cold", "mild", "hot"))
levels(faktor2)
levels(faktor3)
# napravljen je novi faktor koristeci predefinisanu funkciju factor kojoj za
# za ulazni parametar damo jedan vektor koji u sebi ima vrednosti cold, mild, hot
# i onda oznacimo redosled tih vrednosti kao sto smo i uradili
summary(faktor3)
summary(faktor2)
mode(faktor2)
mode(faktor3)
mode(vektor1)
mode(vektor2)
help('mode')
median(vektor2)
mean(vektor2)
median(vektor1)
mean(vektor1)
mode(vektor1)
###
# Missing values
###
# check which values are NAs
vektor4 <- c(1, NA, 3, 4, NA)
is.na(vektor4)
sum(is.na(vektor4))
# NaN missing values
# nastaje kada se primene pogresna izracunavanja npr 0 deljeno sa 0
is.nan(0/0)
is.na(0/0)
is.nan(vektor4)
###
# Data frame
###
# create a data frame
# najcesce cemo to uraditi tako sto cemo ucitavati podatke iz neke tabele
datafrejm1 <- data.frame(
day = c("Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"),
t = c(18, 20, 16, 20, 17, 21, 22))
View(datafrejm1)
# create data frame from existing vectors
datafrejm2 <- data.frame(colors, emission)
# create dataframe called 'co2.emissions' containing
# two columns: 'years' 2000, 2002, 2004, 2006, 2008, 2010 and
# 'emission' 2.7, 2.9, 4, 4.9, 5.3, 6.2
co2.emissions.novi <- data.frame(
year = c(2000, 2002, 2004, 2006, 2008, 2010),
emission = c(2.7, 2.9, 4, 4.9, 5.3, 6.2)
)
# reading data frame from csv file
getwd()
beatles0 <- read.csv("data/beatles_v1.csv")
beatles0
View(beatles0)
str(beatles0)
summary(beatles0)
head(beatles)
tail(beatles)
head(beatles0, 2)
tail(beatles0, 2)
# print number of rows and cols
nrow(beatles1)
ncol(beatles1)
# retrieve and change column names
colnames(beatles1)
beatles2 <- beatles1
colnames(beatles2) <- c("SongTitle", "ReleaseYear", "SongDuration")
colnames(beatles2)
str(beatles2)
colnames(beatles2)[2] <- "YearPublished"
# remove column duration
beatles2$SongDuration <- NULL
str(beatles2)
###
# Subsetting
###
# get element from row 3, col 1
beatles[3, 1]
# get 3rd row
beatles[3, ]
# get rows from 3 to 6
beatles[3:6, ]
# get rows at positions 3 and 6
beatles[c(3, 6), ]
# get 2nd column
beatles[, 2]
# get column Year
beatles[ , "Year"]
years1 <- beatles1$Year
years1 <- beatles[['Year']]
years1
# retrieve all songs released before year 1965
songsBefore1965_1 <- beatles1[beatles$Year < 1965, ]
songsBefore1965_1
# alternative
songsBefore1965_2 <- subset(beatles, Year < 1965)
# retrieve all songs before 1965 with duration lower than 150 sec
songsBefore1965_3_less150 <- subset(songsBefore1965_2, Duration < 150)
songsBefore1965_3_less150
# alternatives
songsBefore1965_4_less150 <- beatles[beatles$Year < 1965 & beatles$Duration < 150,]
# ! voditi racuna dal stavljas jedan ampersand ili dva
# kada se stave && on radi samo selektivno poredjenje
# isto je i sa uspravnom crticom ili |
# ako koristimo po dva znaka on bi kod prvog elementa koji bi se slagao
# rekao da je sve super da se slaze i izbacio bi ne ono sto zelimo
songsBefore1965_4_less150
# druga alternativa
songsBefore1965_5_less150 <- subset(beatles,Year < 1965 & Duration < 150)
songsBefore1965_5_less150
###
# Task 2
###
# print number of songs from 1963 that last more than 2 minutes (120 sec)
songsFrom1963_less120_1 <- subset(beatles, Year == 1963 & Duration > 120)
songsFrom1963_less120_1
# alternativa
songsFrom1963_less120_2 <- beatles[beatles$Year == 1963 & beatles$Duration > 120,]
songsFrom1963_less120_2
###
# Plotting - skiciranje grafikona
###
# create plot for given vectors
x <- seq(1, 15, length.out = 10)
y <- seq(10,22, length.out = 10) + 4
plot(x,y)
help('plot')
plot(x, y, main = 'Example plot', sub = "Djodjo i Roza")
# ima tih nekih plotova za brz prikaz necega
# npr boxplot za sagledavanje raspodele neke varijable
boxplot(x)
# include ggplot2 library
# cesce cemo koristiti plotove iz biblioteke ggplot2 jer su dobri
library('ggplot2')
# mozemo instaliramo i dole na panelu Packages pa Install pa mozemo pronaci tu
# a ucitavanje koje mora svaki put da se radi moze samo da se stiklira na istom panelu
# render plot for givent data frame (columns Year and Duration)
ggplot(data = beatles,
mapping = aes(x = Year, y = Duration))
# render a scatter plot for Year and Duration
# na ovaj ggplot mozemo dodavati razne stvari na panel
# te dodatne modifikacije mozemo dodati plusicem +
ggplot(data = beatles,
mapping = aes(x = Year, y = Duration)) +
geom_point()
# sad mozemo da se igramo estetikom
ggplot(data = beatles,
mapping = aes(x = Year, y = Duration)) +
geom_col()
# render a scatter plot with custom title and axes labels
# ajmo bolan malo estetike
p1 <- ggplot(data = beatles,
mapping = aes(x = Year, y = Duration)) +
geom_point()
p1
p1 + labs(x = "Publishing Year", y = "Song Duration",
title = "Song duration per year") +
theme_bw()
# render a bar chart displaying number of songs in each year
p2 <- ggplot(data = beatles,
aes(x = Year)) +
geom_bar()
p2 <- p2 + labs(x="Publishing Year",
y = "Number of songs",
title = "Published songs per year")
# improve previous chart so that only years with songs are shown...
# include title and axes labels
# to cemo odradimo tako sto godinu promenimo u faktorsku varijablu a ne vektorsku, tj
# prikazemo je kao faktorsku varijablu...
p3 <- ggplot(data = beatles,
aes(x = as.factor(Year))) +
geom_bar()
p3
p3 <- p3 +
labs(x = "YEAR",
y = "NUMBER",
title = "Published songs per year:")
p3
# uvedemo i neke bojice malo u igru
p4 <- ggplot(data = beatles,
aes(x = as.factor(Year),
fill = as.factor(Year))) +
geom_bar(show.legend = F) +
labs(x = "Year",
y = "Number",
title = "Number of published songs per year:") +
scale_fill_brewer(palette = "Set3")
p4
# redner a line chart (Year, Duration) for the first five songs
# with specific line and points properties
p5 <- ggplot(data = beatles[1:5,],
aes(x = Year, y = Duration)) +
geom_line(linetype = 2, linewidth = 1.5, color = 'red') +
geom_point(size = 3, fill = "blue", color = "black", shape = 3) +
labs(x = "Year",
y = "Number",
title = "First five songs:")
p5
###
# Task 3
###
# Create a line chart from a dataset *co2.emissions* (created in Task 1) with the x-axis representing the years, and the y-axis representing the values of CO2 emissions.
# Chart title should be "China CO2 Emissions" and y-axis should have a label "CO2 Emissions".
p6 <- ggplot(
data = co2.emissions.novi,
mapping = aes(x = year, y = emission, title = "China CO2 Emissions")
) +
geom_line(size = 1, color = "#125215", ) +
labs(y = "CO2 Emissions")
P6
#####
# Kraj prvih vezbi
#####
#####
# INTRO_v2
#####
###
# Data Sampling
###
# create vector with values from 1 to 10
x <- 1:10
x
# create a sample of size 5 from the vector
s <- sample(x, 5)
s
# create sample of size 20 from vector where duplicates allowed
s1 <- sample(x, 20, T)
s1
# !!! delicemo dataSet na trening i test (svaki put kad radimo neki zad iz masinskog ucenja)
# !!! obucavamo taj model na skupu za trening a potom testiramo na skupu za test
# set seed and create two sample of size 20 from the vector, where duplicates are allowed
set.seed(10)
sample(x, 20, replace = T)
###
# Matrices
###
# create a 2 x 4 matrix with values from 8 to 1, filled by rows
m <- matrix(data = 8:1,
nrow = 2,
byrow = T)
m
# get first row
m[1, ]
# get element from row 1, col 2
m[1, 2]
# get num of rows
nrow(m)
# get num of cols
ncol(m)
# create two matrices of same dimension
m2 <- matrix(data = 8:1, nrow = 2)
m2
m
# add matrix2 to matrix1
m+m2
# transpose matrix
m3 <- t(m+m2)
m3
# matricno mnozenje
m %*% m2
m
m2
t(m) %*% m2
###
# Lists
###
# create new list with attributes: passport, age, diplomatic
passenger <- list(passport = "P123456",
age = "35",
diplomatic = F)
passenger
# get 2nd element
passenger[2]
# get the value of 2nd element
passenger[[2]]
# get the value of age element
passenger$age #pristup elementu liste
passenger[['age']] #pristup vrednosti elementa liste
# get list length
length(passenger)
# add new list after 2nd element
passenger1 <- append(passenger, list(country = "Serbia"), after = 2)
passenger1$city <- "Belgrade"
# delete 3rd element
passenger1[3] <- NULL
passenger1
# concatenate 2 lists
passenger2 <- list(passport = "P123455",
age = "15",
diplomatic = F)
passengers <- c(passenger1, passenger2)
passengers
# check if travelers is a list
is.list(passanger)
# get names of all list elements
names(passengers)
# get elements with 'age' in their name
?grepl
grepl("age", x = names(passengers))
passengers[grepl("age", x = names(passengers))]
grepl(pattern = "passport", x = names(passangers))
names(passengers)
passengers
###
# Foreach loop
###
# print all odd numbers from 1 to 10 using for each loop
for (i in 1:10) {
if (i %% 2 != 0)
print(i)
}
###
# While loop
###
# print all odd number from 1 to 10 using while loop
i <- 1
while (i <= 10) {
print(i)
i <- i + 2
}
###
# Task 1
###
# create 2x3 matrix with following elements: 3, 9, -1, 4, 2, 6 (by row)
# print only positive values from the first row
matrica1 <- matrix(
data = c(3, 9, -1, 4, 2, 6),
nrow = 2,
ncol = 3,
byrow = T
)
matrica1
for(i in matrica1[1, ]) {
if (i > 0)
print(i)
}
###
# if-else
###
# use if-else function to create new attribute called request with the
# value 'assistance required' if a traveler is younger than 10 years,
# and value 'no special requests' otherwise
passenger1$request <- ifelse(test = passenger1$age < 10,
yes = "assistance required",
no = "no special requests")
passenger1
###
# user defined functions and apply
###
# create function that adds two numbers.
# default value for the second argument is 1
add <- function(a, b = 1) {
#return(a + b)
a + b
}
add(5)
add(5, 6)
# create func returning an absolute value of x
# return result using return() func
abs_val1 <- function(x){
if(is.integer(x) | is.numeric(x)){
abs_x <- ifelse(x >= 0,yes = x,no = -x)
return(abs_x)
}
print("ERROR! Input value is NaN")
}
abs_val1(1)
abs_val1(-47)
abs_val1("jEBOTE")
# samom sebi zadaca
# napravi funkciju koja pretvara dinare (unos kao ulazni parametar defaul vrednost je 117,024)
# ali i proverava da li je ulazni parametar ispravno unet
din2eur <- function(dinari = 117.024, kurs = 117.024){
if (is.numeric(dinari)) {
eur <- dinari / kurs
ifelse(dinari == 117.024, yes = return(dinari), no = return(eur))
}
}
din2eur(60000)
din2eur()
###
# Applying a function over rows and cols in data frame
###
# load data "data/beatles_v2.csv"
beatles_v2 <- read.csv("data/beatles_v2.csv")
str(beatles_v2)
# get number of characters in ths song title "Yellow Submarine"
nchar("Yellow Submarine")
# get number of characters of the first 10 songs
prvih10 <- beatles_v2[1:10,]
prvih10
title_lens3 <- sapply(prvih10$Title, nchar)
head(title_lens3)
# calculate MEAN value of duration and Top.50.Billboard values of all songs from 1963
apply(beatles[beatles$Year == 1963, c("Duration", "Top.50.Billboard")], MARGIN = 2, FUN = mean)
# beatles[beatles$Year == 1963, c("Duration", "Top.50.Billboard")] - OVO SMO DEFINISALI DATA FRAME
# MARGIN = 2 - ako radimo nad KOLONAMA, a 1 stavljamo ako radimo nad REDOVIMA
# FUN = mean - funkcija koju primenjujemo nad unetim datafrejmom!
# ajmo uklonimo NA
apply(beatles[beatles$Year == 1963, c("Duration", "Top.50.Billboard")],
MARGIN = 2,
FUN = function(a) mean(a, na.rm = TRUE))
# evo i alternativa za istu stvar u nastavku:
apply(beatles[beatles$Year == 1963, c("Duration", "Top.50.Billboard")],
MARGIN = 2,
FUN = mean, na.rm = TRUE)
###
# Working with tables
###
# create contingency table of column Year values
year_t1 <- table(beatles$Year)
year_t1
# get 4th element from table
year_t1[4]
# store 4th element from table in a variable
year_1960a <- year_t1[4]
# convert variable to numeric
as.integer(year_1960a)
# sort table in desc order
sort(year_t1, decreasing = T)
# get proportions table for values of the Year column
prop.table(year_t1)
sort(prop.table(year_t1), decreasing = T)
# get proportions table for values of Year column but limiting number
sort(round(prop.table(year_t1), 2), decreasing = T)
# create a contigency table Top.50.Billboard vs. Year
table(beatles$Year, beatles$Top.50.Billboard)
table(beatles$Top.50.Billboard, beatles$Year)
###
# Manipulating data frames
###
###
# Adding new rows and columng
###
# create new column On_album and set FALSE for all songs
beatles$On_album <- FALSE
# create new data frame with 2 cols (with sample data)
set.seed(1)
additional_data <- data.frame(
Rating = sample(
x = 1:5,
size = nrow(beatles_v2),
replace = T
),
Platinum = sample(
x = c(T, F),
size = nrow(beatles_v2),
replace = T
)
)
str(additional_data)
# combine two data frames ***
beatles_v3 <- cbind(beatles_v2, additional_data)
str(beatles_v3)
# get first song
firstSong1 <- beatles_v3[1, ]
# add song to the end of data frame
beatles_v4 <- rbind(beatles_v3, firstSong1)
tail(beatles_v4)
# add song after 3rd song in data frame
beatles_v5 <- rbind(beatles_v4[1:2, ], firstSong1, beatles_v4[3:nrow(beatles_v4), ])
View(head(beatles_v5))
###
# Removing columng and rows
###
# remove attribute On.album
beatles_v5$On_album <- NULL
# remove columns Platinum (at index 10) and Score (at index 11)
beatles_v6 <- beatles_v5[, -c(10,11)]
View(beatles_v6)
# create subset of data frame without songs in wors 2, 4 and 6
# dodela novih indexa :)
rownames(beatles_v6) <- 1:nrow(beatles_v6)
beatlesSubset3 <- beatles_v6[-c(2, 4, 6), ]
View(beatlesSubset3)
# create subset of data frame without songs in rows 1 to 8
beatlesSubset4 <- beatles_v6[-(1:8),]
View(head(beatlesSubset4))
###
# Updating columng and row names
###
# get column names
colnames(beatles_v6)
names(beatles_v6)
# change name of column that starts with 'Genre' to 'Song.genre'
?startsWith
colnames(beatles_v6)[startsWith(colnames(beatles_v6), "Genre")] <- "Song.genre"
# change name of column at index 6 to 'Genre'
colnames(beatles_v6)[6] <- "Genre"
# change row names to string containing word 'song' and song order number
rownames(beatles_v6) <- paste("song", 1:nrow(beatles_v6), sep = "_")
rownames(beatles_v6)
# change row names to string containing order number
rownames(beatles_v6) <- 1:nrow(beatles_v6)
# ove rownames se smatraju labelama a ne stvarnim brojevima koje mozemo koristiti u racunicama
###
# Retrieving and changing values
###
# get songs in rows from 1 to 5, but only attributtes Title and Album.debut
songsSubset4 <- beatles_v6[1:5, c(1,3)]
songsSubset5 <- beatles_v6[1:5, c("Title","Album.debut")]
songsSubset4 == songsSubset5
# get songs from year 1964 not having McCartney as a lead vocal
head(grepl(pattern = "McCartney", x = beatles_v6$Lead.vocal), 10)
songsSubset6 <- beatles_v6[beatles_v6$Year == 1964 &
grepl(pattern = "McCartney", x = beatles_v6$Lead.vocal) == FALSE, ]
View(songsSubset6)
# get songs from year 1958, but only attributes Title and Album.debut
songsSubset7 <- beatles_v6[beatles_v6$Year == 1958, c("Title", "Album.debut")]
View(songsSubset7)
# alternativno moze i ovako:
songsSubset8 <- subset(beatles_v6, Year == 1958, c("Title", "Album.debut"))
songsSubset7 == songsSubset8
# create vector of logical values denoting whether the attribute Album.debut has a value or not
head(beatles_v6$Album.debut)
tail(beatles_v6$Album.debut, 20)
noAlbumDebut1 <- beatles_v6$Album.debut == ""
# compute how many songs lack tha data about the debut album
sum(noAlbumDebut1)
# for songs without debut album data, set the value of the Album.debut attribute empty
beatles_v7 <- beatles_v6
beatles_v7$Album.debut[noAlbumDebut1] <- "empty"
View(beatles_v7)
# set the value beck to empty string
beatles_v7$Album.debut[noAlbumDebut1] <- ""
###
# Saving dataset
###
# save dataset to CSV file, but without row names (row numbers) column
write.csv(beatles_v7, "data/beatles_v7_pro.csv", row.names = F)
# save R object for next session into file "data/beatles_v7.RData"
saveRDS(beatles_v7, "data/beatles_v7.RData")
# restore R object from file "data/beatles_v7.RData" in next session
beatles_v8 <- readRDS("data/beatles_v7.RData")
str(beatles_v8)
###
# Task 2
###
# create new column in 'beatles' data frame called 'Billboard.hit'
# having TRUE for all songs that were in top 50 Billboard
# (songs that have the Top.50.Billboard defined),
# and FALSE for all other songs (not having this value set).
beatles_v8$Billboard.hit <- FALSE
notbillboardHit <- ifelse(is.na(beatles_v8$Top.50.Billboard),
yes = FALSE,
no = TRUE)
beatles_v8$Billboard.hit <- notbillboardHit
#####
# Kraj drugih vezbi
#####
#####
# 03 LINEARNA REGRESIJA
#####
# ucitavanje biblioteka
library(MASS)
library(ggplot2)
install.packages('corrplot')
library(corrplot)
install.packages('MASS')
library(MASS)
# examine docs for Boston dataset
?Boston
# examine structure of dataset, kojeg su tipa atribut
str(Boston)
# za regresiju nam treba num podaci ili int
# kako su svi podaci num onda smo lagani
# svi ovi atributi dolaze u obzir da budu ukljuceni u nase predvidjanje
# druga vazna stvar je da procenimo dal neka od ovih ulaznih varijabli moze
# znacajnije da utice na predvidjanje izlazne varijable
# jedan od nacina je da pokusamo da razumemo znacenje svake varijable
# pa na osnovu toga kazemo ovo ima smisla ukljuciti ili ne npr rm broj soba, logicno da ako ima vise soba bice cena visa
# npr zagadjenost vazduha, ali to je subjektivno, mozda nekom drugom to nije vazno
# da bismo napravili objektivniju procenu moramo racunati korelacija izmedju nezavisnih
# varijabli (atributa itd - ovo su sinonimi) varijabla = atribut itd...
# znaci treba da vidimo postoji li zavisnost odnosno korelisanost izmedju ulaznih varijabli i izlazne varijable
# i koji su to koje se izdvajaju po vecem stepenu zavisnosti tj korelisanosti
# zato pravimo korelacionu matricu
# compute the correlation matrix
# moramo videti koju cemo metodu izabrati za izracunavanje korelacije npr PIRSON, KENDL, SPERMAN
# po defaultu se koristi PIRSONova korelacija
# PIRSONovu mozemo koristiti onda kada znamo da su nam varijable normalno rasporedjene
# sada cemo raditi samo PIRSONovu metodu da racunamo korelaciju
# na drugim casovima cemo uciti i druge metode
bostonCorr <- cor(Boston)
View(bostonCorr)
# simetricna matrica je matrica ciji su redovi i kolone isti
# korelaciona matrica nam pomaze da uocimo koji su to atributi visoko korelisani tj visoko zavisni za izlaznu varijablu
# njih cemo selektovati za nas model
# druga VAZNA stvar jeste da uocimo korelisanost izmedju samih ulaznih varijabli
# postoji nesto sto se zove MULTIKOLINEARNOST a to je kada su ove ulazne varijable visoko medjusobno korelisane
# tada te ulazne varijable (nezavisne varijable) cemo izbegavati da ubacujemo zajedno u model!!! VAZNO
# one option for plotting correlations:
# using colors to represent the extent of correlation
corrplot(bostonCorr)
# s obzirom da je simetricna matrica mozemo pokazati samo jednu polovinu
corrplot(bostonCorr, type = "upper", diag = F)
# kako sad da vidimo visoko korelisane varijable za izlaznu
# posmatramo ovu zadnju kolonu medv jer se ona odnosi na varijablu cije vrednosti predvidjamo
# u toj koloni gledamo gde su nam veliki tamni kruzici
# naravno evo rm, sto je veci broj soba to je i cena veca...
# i ovo lstat lower status of population je visoko korelisano sa y
# sad treba da vidimo kolike su otprilike same vrednosti korelacije, pa da se odlucimo sta cemo izabrati
# another option, with both colors and exact correlation score
corrplot.mixed(bostonCorr)
# mi cemo ovog puta uzeti u razmatranje sledece ulazne varijable:
# 'rm', 'lstat' i 'ptratio' koji ima vrednost -0.51 sto i nije bas visoko zavisni ali ok
# pored ovih varijabli koje su visoko korelisane sa y, treba da povedemo racuna
# i o medjusobnoj korelaciji ovih odabranih ulaznih varijabli!!!
# ako jesu visoko korelisane medjusobno, moracemo izbeci te slucajeve
# npr @lstat i @ptratio nisu visoko korelisani svega 0.37 je vrednost izmedju njih
# ali @lstat i @rm su prilicno korelisani medjusobno tu je vr -0.61 TO JE PROBLEMATICNO
# ipak cemo ih ukljuciti obe jer -0.61 nije bas bas visoko...
# npr ne bismo smeli NIKAKO da zanemarimo varijablu @tax i njenu korelisanost sa @rad
# tu je vrednost 0.91 i njih dve NIKAKO ne bi smeli da koristimo u modelu zajedno! VAZNO!!!
# izabrali bi samo onu varijablu koja je vise korelisana sa y
# create scatterplots for the dependent variable @medv and
# candidate indepenet variables @lstat, @rm and @ptratio
# osnovna pretpostavka lin regresije je linearnost
# linearna relacija izmedju nezavisnih i zavisnih promenjljivih
# mi nismo ispitali dal je relacija y i x-eva linearna ili nije. VAZNO!!!
ggplot(Boston, aes(x = medv, y = lstat)) + geom_point()
# proveravanje da li su @medv i @lstat u medjusobno linearnoj ili nelinearnoj relaciji?
# prema ovom grafu mozemo uociti da i nije bas linearana veza ove dve varijable
# u pocetku tamo negde do vrednosti 30 na x osi bi mogla da se povuce jedna linija
# koja bi otprilike mogla kazati da su vrednosti jedne i druge varijable u nekom lin odnosu
# medjutim tamo posle 30 ka 50 se vidi da bi se ta linija znatno promenila
# otuda mozemo reci da nisu kolinearne varijable
# hajd da proverimo druge neke varijable
ggplot(Boston, aes(x = medv, y = rm)) + geom_point()
# i ovde moze prava linija da se nacrta negde kroz tacke
# ali imamo i ovih koje iskacu
# nije idealno ali nije ni tako lose
###
# PRIPREMA PODATAKA
###
###
# Split the data into training and test sets
###
# install and load caret package
install.packages('caret')
library(caret)
# assure the replicability of the results by setting the seed
# kada kreiramo trening i test set moramo da obezbedimo da se nasumice izaberu podaci
# voditi racuna o tome, ne mozemo samo prvih 400 izabrati da budu za treniranje a
# preostalih 400 da budu za testiranje
# odnosno ne mozemo biti pristrasni 'bias' u podeli podataka na @trening i @test setove
# mozda ima neka pravilnost pa su prvi uneti podaci po necemu slicni, a ovi poslednji imaju neku drugu slicnost
# zato ja VAZNO!!! da obezbedimo NASUMICNOST
# druga BITNA stvar jeste da u trening i test setu imamo priblizno jednaku raspoedelu izlazne promenljive
summary(Boston$medv)
# mi treba da obezbedimo da kad podelimo dataset
# da ce ta raspodela da se zadrzi priblizno i u trening i u test setu podataka
# to je npr kao kad bi vezbali neke zadatke na casovima
# onda izadjemo na ispit i dobijemo zadatke koje nemaju puno veze sa onim sto smo vezbali
# zato je vazno da raspodela izlazne promenljive bude SLICNA i u treningu i u testu
# rezime 1. stvar je obezbedjivanje nasumicnih podataka za trenint i test
# rezime 2. stvar je da se zadrzi 'stratifikacija' tj. raspodela izlazne prom u treningu i testu
# da bismo radili sa istim nasumicnim podacima, radicemo sa isto setovanim seedovima
set.seed(3)
# obicno uzmemo 80% podataka za treniranje a ostalo za testiranje
# set.seed je zapravo fja koja fiksira nacin na koji ce se odabrati nasumicni podaci
# ako stavimo npr set.seed(3), obezbedicemo da bilo ko, ko radi sa istim podacima
# moze na svom modelu da koristi isti sed.seed(3)
# i tada cemo moci nas model i njegov model da valjano uporedjujemo
trainIndices <- createDataPartition(Boston$medv, p = 0.8, list = F)
# list = FALSE, da nam rezultati budu kao dataFrame a ne kao lista
# set.seed i createDataPartition bi trebalo zajedno da selektujemo i izvrsimo
# kako bi inicirani seed imao efekta
trainIndices
# select observ. at the positions defined by trainIndices vector
bostonTrain <- Boston[trainIndices,]
# select observ. at the positions that are NOT in the trainIndices vector
bostonTest <- Boston[-trainIndices,]
summary(bostonTrain$medv)
summary(bostonTest$medv)
# prema rezimeu mozemo videti da smo uspeli podeliti podatke caret paketom, a da
# nismo narusili stratifikaciju odnosno raspoedelu izlazne varijable niti u trening nit u test setu
###
# KRAJ PRVOG KLIPA
###
###
# Create a linear regression model
###
# build an lm model with train dataset using formula: medv ~ lstat + rm + indus + ptratio
# pitanje studenta na vezbama
# na koji nacin smo zakljucili da koristimo model lin regresije za ove podatke?
# odgovor: prvo proveravamo dal su zavisne i nezavisne varijable u linearnoj relaciji? to je pretpostavka lin regresije
# zbog toga smo radili scatterplot i pokusavali da utvrdimo dal postoji linearna povezanost...
# otprilike moze da se povuce linija po onim vrednostima sa scatterplota, ne idealno ali moze
# uglavnom mi sada ucimo da koristimo model lin regresije pa smo i zbog toga odabrali ovaj model
lm1 <- lm(medv ~ lstat + rm + indus + ptratio, data = bostonTrain)
# print the model summary
summary(lm1)
class(lm1)
?lm
# prvo pogledamo F-statistiku koja nam predstavlja osnovu za sagledavanje dal ima smisla da
# dalje razmatramo nas model ili ne
# ova stat vrsi proveru stat testa cija je nulta hipoteza tj pretpostavka da
# nijedna od varijabli nije znacajna za predvidjanje izlazne promenljive, odnosno
# da su svi ovi koeficijenti koje vidimo da su svi nula
# ukoliko verovatnoca f-stat pokaze da je jako mala, onda cemo odbaciti nasu pretpostavku
# u nasem primeru ta verovatnoca je p-value: < 2.2e-16
# ovo nam govori da je verovatnoca veoma mala, te zakljucujemo da su nam podaci ipak okej za nase modeliranje
# odbacujemo hipotezu da je nas model beskoristan i ima dalje smisla razmatrati ovaj nas model...
# R-squared to je kod nas prevodjeno sa Koeficijent Determinacije
# to je indikator koji nam pokazuje koji procenat varijabiliteta izlazne promenljive
# ovaj model opisuje...
# ovde je to skoro 69%, dakle toliko procenata varijabiliteta izlazne prom nas model opisuje...
# kojim god modelom da radimo, treba pokriti sto veci procenat varijabiliteta izlazne prom
# da bismo tacnije utvrdili i znali na koji se nacin menja izlazna prom...
# ovih 69% nije lose
# ova Adjusted R-squared se racuna zbog toga sto ce R-squared rasti sto imamo vise varijabli u modelu
# bez obzira dal su te varijable korisne ili ne korisne
# i zbog toga se radi njegovo priloagodjavanje u modelu
# zato pise ovo ADJUSTED...
# te vrednosti ova originalan i ova prilagodjena su ovde gotovo jednake
# jer imamo samo cetiri varijable u nasem modelu, pa zbog toga nije znacajno uticalo...
# Estimate(procena) pokazuje kako je varijabla povezana sa izlaznom varijablom
# pozitivno ili negativno i u kojoj meri
# Error - neka greska procene koeficijenta
# t statistika - vazna je kao indikacija da se ovde radi jedan statistcki test
# koji za nultu hipotezu (prvobitnu pretpostavku) ima da ova varijabla lstat nije povezana
# sa izlaznom varijablom, tj da je ovaj koeficijent @lstat/Estimate jednak nula, a nama pokazuje
# da je -0.55340, sto znaci da je veoma razlicit od nule!
# sto ima vise zvezdica to mi mozemo prtpostaviti da je neka varijabla znacajnaija i korelisanija sa y
# vidimo da @indus nije dovoljno znacajan za nas model
# kao sto smo i mogli zakljuciti na pocetku u korelacionoj matrici
# znaci on nije prediktor izlazne varijable
# kada uocimo ovo mi obicno korigujemo model i izbacimo iz razmatranja ovu varijablu
lm1 <- lm(medv ~ lstat + rm + ptratio, data = bostonTrain)
summary(lm1)
# na primer kada gledamo Estimate
# ta kolona kazuje koliko ce se y smanjiti ili uvecati kada unesemo neku od ulaznih varijabli
# na primer za uvecanje lstat za 1%, cena kvadrata ce se smanjiti za -0.54546 * 1.000 dolara =>
# => smanjice se za 545,46 dolara!!! i to pod uslovom da ostale varijable ostaju nepromenjene
# tzv old thing equal
# npr imamo @rm, ako je broj soba veci za 1 => tada ce cena biti povecana za skoro 5.106 dolara...
# sumarne statistike vezane za reziduale
# sta su reziduali?
# reziduali su razlike izmedju predvidjenih i stvarnih vrednosti
# nama je cilj da ti reziduali budu sto manji, odnosno da teze nuli
# bilo bi idealno da prosecna vrednost bude nula, ili mediana
# ovde kod nas bas i nije ali boze moj nije ni daleko od nule!!!
# bilo bi dobro da su reziduali normalno raspodeljeni zbog samih karakteristika lin regresije
# ove vrednosti za kvartile nam ne govore bas o normalnoj raspodeli ali pogledacemo posle detaljnije
# intercept je neka predvidjena vrednost izlazne varijable ukolko bi sve ulazne var bile nula
# prakticno ne mozemo to bas zamisliti jer ako bi kuca imala nula soba, kakvi bi mogla biti vrednost izlazne prom??
# ali teroijski ajde.. ovde je to neka vrednost 14.2374515
# intervali poverenja za koeficijente???
# koji je to neki interval za kretanje vrednosti tih koef za koje mozemo smatrati
# sa 99% verovatnoce da ce se stvarne vrednosti koef nalaziti u tom intervalu
# print coefficients
coef(lm1)
# compute 95% confidnece interval for the coefficients
# intervali poverenja za koeficijente???
# koji je to neki interval za kretanje vrednosti tih koef za koje mozemo smatrati
# sa 95% ili npr 99% verovatnoce da ce se stvarne vrednosti koef nalaziti u tom intervalu
confint(lm1, level = 0.95)
###
# Diagnostic Plots
###
# jos jedna stvar vezana za one pretpostavke, ono o pretpos linearnosti, itd.. preko scatterplotova
# postoje jos neke pretpostavke da se provere dal su zadovoljene
# da bi ovaj nas model lin regresijie mogli smatrati pouzdanim
# postoje razliciti nacini da se ovo proveri
# zgodno je preko 4 dijagnosticka plota za ovu svrhu
# da bi ih sagledali istovremeno, podelicemo plot na 4 dela, kvadranta da bi ih lepo sagledali
# split the plotting area into 4 cells
par(mfrow = c(2,2))
# print the diagnostic plots
plot(lm1)
# reset the plotting area
# vratimo zbog kasnijeg prikazivanja plotova
par(mfrow = c(1,1))
# komentari vezani za data 4 dijagnosticka plota
# prvi plot Residuals vs Fitted
# odnosi se na pretpostavku linearnosti, evo jedan drugaciji nacin da se to proveri
# ne kao na pocetku sto smo gledali u scatterplotu
# ovde su prikazane fittovane vrednosti, to su u stvari 'predvidjene' vrednosti
# a na y osi su vrednosti reziduala, rekli smo to su razlike stvarnih i predvidjenih vrednosti
# idealna situacija bi bila da su ti reziduali svi nula
# ima ova jedna mala tanka horizontalna linija koja prolazi kroz nulu
# znaci idealna sit bi bila da ove tackice sve leze na toj horizontalnoj nula liniji
# ono sto bi hteli realno da ocekujemo je da su ove tacke podjednako razbacane iznad i ispod linije
# i da se relativno nalaze blizu nje
# to bi znacilo da pravimo greske u slicnom obimu za sve predvidjene vrednosti
# ne izdvajaju se vrednosti gde negde pravimo vecu a negde manju gresku
# crvena linija pokazuje trend u podacima
# vidimo da je malo zakrivljena, i odudara od nulte linije
# nije u potpunosti zadovoljen zahtev linearnosti
# tj ne mozemo potpuno smatrati da postoji linearna relacija izmedju zavisne i nezavisne prom
# prakticno jos jedna potvrda sto smo vec videli na scatter plotu
# drugi plot ili ovaj QQ plot, normal QQ plot
# njime se proverava da li neka varijabla ima normalnu raspodelu ili nema
# ono sto proveravamo jeste raspored reziduala
# jer je jedna od pretpostavki linearne regresije jeste da su reziduali normalno raspodeljeni
# kako proveravamo?
# ukoliko su normalno raspodeljeni, ovi reziduali bi svi lezali na ovoj isprekidanoj liniji
# u nekoj meri je zadovoljeno ali ne skroz jer ih dosta odstupa
# zato ne mozemo smatrati da je raspodela reziduala normalna
# reziduali su razlika stvarne i predvidjene vrednosti
# ono sto mi predvidjamo je cena kuce
# znaci cena kuca je stvarna vrednost minus cena predvidjena i to je graf
# treci plot
# odnosi se na pretpostavku koja se zove 'homoshedasticnost'
# rec je o varijabilitetu reziduala
# tj da reziduali nekako ujednaceno variraju
# to se strucno zove homoshedasticnost tj ujednaceno variranje
# hocemo da vidimo tackice koje predstavljaju varijabilitet reziduala
# hocemo da vidimo da budu nekako razasute po ovom plotu ovde, onako bez plana i redosleda
# jer to nam onda govori da za razlicite predvidjene vrednosti onako nasumice variraju ti reziduali
# to i jeste ta osnovna pretpostavka
# ova pretpostavka da taj varijabilitet bude ujednacen a ne mora biti mali
# hteli bi da vidimo neku horizontalnu liniju bilo gde povucenu
# oko koje ce i s gornje i s donje strane tackice biti razbacane
# mi vidimo da linije nije bas horizontalna ali ne odstupa mnogo
# nije idealno ali okejjjj
# uslovno receno moglo bi da prodje
# naspricane su tackice svuda i otprilike pola je iznad linije pola je ispod...
# cetvrti plot
# se odnosi na leverage points
# to su obzervacije koje imaju neuobicajeno velike ili male vrednosti za ulazne promenljive
# znaci ne za izlaznu prom nego za vrednosti na osnovu kojih predvidjamo
# to su one obzervacije koje po jednoj varijabli ili po vise varijabli dosta odudaraju od ostalih
# i onda one kao takve mogu dosta da uticu na samu procenu modela lin regresije
# pretpostavka je da takvih obzervacija nema
# mi ovde to proveravamo
# ovde je implementiran pristup tzv Kukova Distanca
# da li ovde ima tackica koje odskacu preko te isprekidane linije
# mi imamo ovde obz pod brojem 369 npr
# ali nisu toliko razlicite da nam prave problem
# jer nisu preko te isprekidane linije
# ima i 365 ali ni ono ne prelazi liniju cookovu distancu
# kukala nama majka
# ako predje distancu onda nam cookala majka
# ovaj 4. plot cu nazvati cookala nam majka zbog cookove distance... asocijacija
# nas neki zakljucak
# ovaj nas model ne ispunjava bas idealno uslove
# ali da je neprihvatljiv tako da idemo dalje sa njim
###
# Make predictions and evaluate the model
###
# calculate the predictions with the fitted model over the test
lm1Pred <- predict(lm1, newdata = bostonTest)
# print out a few predictions
head(lm1Pred)
# calculate the predictions with fitted model over the test data,
# including the confidence interval
lm1Pred <- predict(lm1, newdata = bostonTest, interval = 'conf')
head(lm1Pred)
# sad imamo i gornju i donju granicu za predvidjenu vrednost
# dobijeno je uzimajuci u obzir da koeficijenti modela mogu da variraju od neke gornje do donje granice
# na osnovu toga su odredjene ove predikcije ovde
# combine the test set with the predictions
bostonTest$medvPred <- lm1Pred[,1]
# plot actual (medv) vs predicted values
ggplot(bostonTest) +
geom_density(aes(x = medv, color = 'actual')) +
geom_density(aes(x = medvPred, color = 'predicted')) +
theme_minimal()
# ovaj graf pokazuje funkciju gustine raspodele za predv i stvarnu vrednost
# najveci broj kuca ima cenu oko 22 23 hiljade dolara
# vidi se po crvenoj krivi
# ali u uspon delu, dosta dobro sledimo stvarne vrednosti
# ali u padu, ova nasa plava ide preko crvene krive
# to znaci da predvidjamo da ima vise kuca od 25 do 30ak hiljada dolara nego sto je realno
# tako da nam ovakav grafikon omogucuje da sagledamo gde nas model pravi gresku
# i onda da razmislimo dal ima jos neka varijabla koja bi mogla da nam pomogne,
# da ispravimo nas model, da bude precizniji...
# videli da postoji preklapanja ali i da pravimo gresku...
###
# KRAJ DRUGOG KLIPA LIN REGRESIJE
###
# EVOLUCIONE METRIKE
# metrike za procenu koliko je nas regresioni model dobar
# to su standardne metrike za procenu, mi ih radimo 2-3 ali ih ima mnogo vise...
# jedna od tih metrika je RSS
# calculate RSS (Residual Sum of Squares)
# suma kvadrata reziduala
# reziduali su razlike izmedju stvarne i predv vrednosti
# cilj nam je da napravimo model gde ce ta razlika biti sto manja
# tj da reziduali imaju sto manju vrednost
# ova mera RSS nam ukazuje koliko mi i dalje pored kreiranog reg modela mi neuspevamo da predvidimo
# izlaznu varijablu, tj u kojoj meri ne uspevamo da predvidimo taj varijabilitet izlazne varijable
RSS_lm1 <- sum((bostonTest$medv - bostonTest$medvPred)^2)
# dakle kvadrirali smo reziduale i sumirali, suma kvadrata reziduala
# i dobili smo obim varijabiliteta izlazne promenljive koji nismo uspeli da objasnimo
# nasim regresionim modelom
# calculate TSS (Total Sum of Squares)
# ova metrika ukazuje koliki je generalno varijabilitet izlazne prom
# nju dobijamo tako sto umesto predvidj vrednosti koju smo dobili kroz nas regres model,
# u ovom slucaju ce to biti razlika stvarne vrednosti iz test seta i nekog oblika predvidj te vrednosti
# koji bi mogli da napravimo ukoliko ne znamo sad neku posebnu regresionu metodu ali imamo podatke iz trening seta
# sad neko nam da trening set i ajde ti predvidi kavka ce biti vrednost cene kuce, a mi
# ne znamo da koristimo nijednu od regresionih metoda
# pa onda mozemo da se snadjemo da uzmemo srednju vrednost mean npr
# to je neka ocekivana vrednost varijable, pogotovo kad ta varijabla sledi normalnu raspodelu
# to je ono sto bi trebalo nasa izlazna varijabla da sledi
# tako da cemo mi za ovaj TSS, koristiti kao predvidjenu vrednost
# ustvari srednju vrednost na trening setu!!!
TSS <- sum((bostonTest$medv - mean(bostonTrain$medv))^2)
# calculate R-squared on the test data
# R2 = (TSS-RSS)/TSS = 1 - RSS/TSS
# R2 smo vec videli u summary naseg modela, ali to je R2 trening seta a ne za test set
# Sad ga racunamo na test setu
# R2 nam je bas u fokusu interesovanja to je koeficijent determinacije
# R2 je procenat varijabiliteta izlazne prom koji smo reg modelom uspeli da objasnimo
# dobijmo ga: ukupni varijabilitet TSS - varijabilitet koji nismo uspeli da objasnimo RSS, normalizovano
# tj deljeno sa ukupnim varijabilitetom i to kad se sredi dobijemo 1 - RSS/TSS
R2_lm1 <- 1 - RSS_lm1/TSS
R2_lm1
# ovaj rezultat je manje u odnosu na 0.688 sa trening seta
# ali uvek je slabije u testu nego u treningu, jer su kod treninga uvek najbolji moguci rezultati
# dobro je sto nije znacajno nizi rezultat vec je tu u nekoliko procenata razlike
# calculate RMSE
# RMSE = sqrt(RSS/n)
# ukazuje kolko pravimo prosecne greske
# iskazan je u onoj jedinici u kojoj je i izlazna promenljiva
RMSE_lm1 <- sqrt(RSS_lm1/nrow(bostonTest))
RMSE_lm1
# ispada oko 5250 dolara, sto nije mala greska
# compare medv mean to the RMSE
prosecna <- mean(bostonTest$medv)
prosecna
RMSE_lm1 / prosecna
# greska koju pravimo je oko 23% prosecne vrednosti kuca koju predvidjamo, sto je ogromno
# i ukazuje da bi trebalo neki malo drugacijim nacinom da predvidjamo
###
# kraj ovih gore metrika
###
###
# Create a more complex model and examine multicollinearity
###
# build a new model using all of the variables
lm2 <- lm(medv ~ ., data = bostonTrain)
# print summary for model
summary(lm2)
# pnovo prva stvar koju gledamo je F stat da je ona verovatnoca jako mala
# kako bismo odbacili hipotezu da nijedna od ulaz var nije znacajan prediktor izlazne var
# i usvojiti alternativnu da bar neka od sluc prom jeste prediktor izlazne var
# zvezdice ukazuju na znacaj
# multi kolinearnost moramo odbaciti
# npr @rad i @tax imaju visoku kolinearnost prema korelacionoj matrici
# moze se desi da postoji visoka korelisanost izmedju dve ili vise varijabli
# to bas ne mozemo da otkrijemo samo preko matrice
# da bi procenili to, koristicemo statistiku vif various inflation factor
install.packages('car')
library(car)
vif(lm2)
?vif
# ova fja za svaku varijablu izracuna neku vrednost, to je ovo vif
# pitanje kako koristimo ovo
# postoji neki prag, granica, ukoliko vif vrednost neke varijable prelazi tu granicu
# onda smatramo da ta var uslovljava multikolinearnost i da je treba iskljuciti iz modela
# sta je ta granica ne postoji potpuno saglasnost statisticara, a to je
# da kvadratni koren ove vif vrednosti ne bi trebalo da bude veca od 2
sort(sqrt(x = vif(lm2)))
# sad cemo da iskljucimo jednu od ove dve problematicne varijable pa da vidimo dal je bolje
# build a new model using all vars except those with highest VIF values, problematic ones
# dobra je praksa da iskljucujem minimal broj variajbli, jednu po jednu pa da vidimo dal
# cemo na taj nacin resiti problem
lm3 <- lm(medv ~ . - tax, data = bostonTrain)
sort(sqrt(x = vif(lm3)))
# posto je nox visok preko 2, i njega cemo iskljucimo pa sleduje
lm4 <- lm(medv ~ . - tax - nox, data = bostonTrain)
sort(sqrt(x = vif(lm4)))
# print the model summary
summary(lm4)
# i dalje je R2 visok, skoro 73% varijabiliteta izlazne promenljive objasnjavamo
# primetimo ovaj @rad
# ona se pokazala posebno znacajna
# a sad se pokazuje kao visoko znacajna
# jer je taj rad i onaj tax su bili medjusobno visoko korelisani
# i predimenzionisali su sopstveni uticaj na izlaznu var
# zato je bitno da multikorelacionalnost otklonimo
# calcucalte predictions over test data
lm4Pred <- predict(lm4, newdata = bostonTest)
# combine test set with the predictions
bostonTest$pred_lm4 <- lm4Pred
colnames(bostonTest)[colnames(bostonTest) == 'medv_pred_lm1'] <- 'medvPred_lm4'
colnames(bostonTest)[colnames(bostonTest) == 'medvPred'] <- 'medvPred_lm1'
# plot actual (medv) vs. predicted values
library(ggplot2)
ggplot(bostonTest) +
geom_density(aes(x = medv, color = 'actual')) +
geom_density(aes(x = medvPred_lm1, color = 'lm1-predicted')) +
geom_density(aes(x = medvPred_lm4, color = 'lm4-predicted')) +
theme_minimal()
# bolji je lm4 od lm1
# i dalje od 25k do 30k pravimo gresku u proceni ali je bolje
# calculacte RSS
RSS_lm4 <- sum((bostonTest$medv - bostonTest$medvPred_lm4)^2)
# suma kvadrata reziduala
# calcucalte R2 on test data
R2_lm4 <- 1 - RSS_lm4/TSS
R2_lm4
# bolje nego na prvom modelu
# calculate RMSE
RMSE_lm4 <- sqrt(RSS_lm4/nrow(bostonTest))
RMSE_lm4
# malo smo smanjili gresku, bilo je preko 5k dolara sad je ispod 5k dolara
# i dalje je to velika greska ali boze moj
#####
# KRAJ LINEARNE REGRESIJE
#####
#####
# CLASSIFICATION TREES
#####
# ostajemo u sferi nadgledanog masinskog ucenja
# danasnja tema su klasifikaciona stabla ili stabla odlucivanja
# sa regresije prebacujemo se na zadatak klasifikacije
# kod regresije smo predvidjali neku numericku vrednost
# a sada predvidjamo neku kategoricku vrednost, tj pripadnost
# nekoj od unapred definisanih kategorija
###
# Classification trees
###
# load ISLR package
install.packages('ISLR')
library(ISLR)
# get the Carseats dataset docs
?Carseats
# Sales ce nama predstavljati na neki nacin izlaznu varijablu
# necemo se baviti regresijom vec klasifikacijom pa cemo na osnovu ove varijable napraviti novu varijablu
# koja ce biti binarana tj kategoricka npr visoka prodaja, niska prodaja i to ce nam biti izlaznu var
# predvidjacemo je na osnovu ostalih varijabli
# examine dataset structure
str(Carseats)
# examine Sales var distribution
# mi cemo reci da je visoka prodaja ona koja je veca od 3. kvartila
# one koje su ispod 3. kvartile su niska prodaja
# get 3rd quartile of the Sales var
summary(Carseats$Sales)
# 3. kvartil je 9.32 tj ovo je 75 percentil to je to
# a sta ako je menadzer rekao da je visoka prodaja iznad 65og percentila
# onda nam summary fja ne pomaze previse
# morali bi drugim nacinom da odredimo:
sales_75perc <- quantile(Carseats$Sales, probs = 0.75)
sales_65perc <- quantile(Carseats$Sales, probs = 0.65)
# napomena: percentile, quantile to se sve odnosi na isti pojam
# create new var HighSales based on value of the Sales var
Carseats$HighSales <- ifelse(test = Carseats$Sales > sales_75perc,
yes = "Yes",
no = "No")
# ona je karakter tipa i to nam nije bas odgovarajuce
# treba da je prebacimo u factor
Carseats$HighSales <- as.factor(Carseats$HighSales)
# get distribution of HighSales var
# da vidimo kakva je distribucija
# mozemo ocekivati da je 75% u klasi no a 25% u klasi yes jer smo tako postaivli uslov
# ali generalno je dobro da vidimo kako su distribusane klase
table(Carseats$HighSales)
prop.table(table(Carseats$HighSales))
# remove Sales var
# mi smo ovu var iskoristili za kreiranje izlazne var
# sto znaci da ako bi je sada iskoristili za predvidjanje izlazne var
# tada bismo sto posto predvideli izlaznu var
# jer smo HighSales prosto izveli iz Sales i stopostotno su asocirane jedna na drugu
# zato nikako ne bi trebalo da Sales var koristimo za predvidjanje
# moze da ostane u datasetu
# ali ne sme da bude u modelu!!! VAZNO
# zato je bolje je odmah izbaciti da je ne bismo slucajno uvrstili u model
Carseats$Sales <- NULL
###
# Create train and test datasets
###
# isti nacin kao i prosli put sto smo radili
# mi mozemo odmah da nastavimo sa podelom na test i train zato sto
# je stablo odlucivanja takav tip algoritma da je jako tolerantnatn u smislu da moze da radi
# sa najrazlicitijim vrstama varova i sa faktorima i sa num i sa faktoriskim binarnim var (koja imaju 2 nivoa)
# videcemo kasnije na casovima da cemo na drugim algoritmima morati da radimo transformacije varijabli
# neke cemo var morati iskljucivati jer neki algoritmi rade samo sa num, ili samo sa faktorima itd...
# tako da ovo sa stablima je sto se toga tice laganija, fleksibilnija varijanta
# load caret package
library(caret)
library(ggplot2)
library(lattice)
# create train and test datasets
# postavljamo seed naravno
set.seed(1)
train_indices <- createDataPartition(Carseats$HighSales, p = 0.8, list = F)
# prvo smo nasli nasumicne indexe po nasumicnom inicijatoru 1
# pa cemo te indexe koristiti za kreiranje trening i test seta:
trainData <- Carseats[train_indices,]
testData <- Carseats[-train_indices,]
# print distributions of the outcome var on train and test datasets
prop.table(table(trainData$HighSales))
prop.table(table(testData$HighSales))
# bogme distribuirano je slicno oko 75% No je i u testu i u treningu
# dok je Yes skoro 25% takodje u oba seta
###
# Create a prediction model using Classification Trees
###
# load rpart library
install.packages('rpart')
library(rpart)
# rpart used random sampling, so we have to set the seed value before calling the
# build the model
?rpart
# ukoliko hocemo sve var da ubacimo onda stavimo samo tackicu, da ne navodimo sve
# ekspliciton moze se navede method = class onda ova fja zna da treba da napravi
# klasifikacioni model, mada i ne mora jer ona sama vidi kog je tipa nasa izlazna prom
# u nasem slucaju ona je factor, pa ce on po defaultu staviti method = class
# mogu se praviti i numericka stabla tj regresiona stabla
# delimicno se bazira ova rpart fja na nekim nasumicnim random procesima
# onda treba opet seed da postavimo
set.seed(1)
tree1 <- rpart(HighSales ~ ., data = trainData)
# print the model
print(tree1)
# kaze n je 321
# da smo koristili 321 observ za kreiranje stabla
# to je bas onoliko koliko ima u treningu setu
# mi kazemo stable mada tree je u stvari drvo
# uvek stablo krece od nekog korena
# mada je ovde nama koren na vrhu i ide naopacke
# evo ga na prvom mestu root
# u njega stizu sve observacije
# ovo u zagradi pokazuje koji je odnos klasa medju ovih 321 observacija
# evo 75% pripara klasi No
# i oko 25% pripada klasi Yes
# to smo videli pomocu onog prop.table fje
# posto su dominante ove obser sa No
# onda ce na tom nivou ako se radi predikcija
# predikcija ce biti No
# i u tom slucaju greska koja se pravi ce biti 80
# na 80 observacija ce biti pogresno predvidjena vrednost
# i evo na prvom tom stavu vidimo da pise 80
# to bi bila neka predikcija ukoliko ne bi znali vrednost nijednih ulaznih argumenata
# medjutim mi imamo neke od ovih ulaznim argumenata ovi koji opisuju te predovanice i njihovu prodaju
# i algoritam u stablu odlucivanja je procenio da je ShelveLoc najznacajnija varijabla za klasifikaciju
# i nju je prvu izabrao da na osnovu nje radi granjane
# e sad ukoliko je ShelvLoc Bad ili Medium
# takvih observacija ima 256, medju njima 84% su one kod koje nije ostvarena visoka prodaja
# i samim tim posto je to dominantna klasa, mi predvidjamo da ukolko je ShelvLoc Bad ili Medium
# nece biti ostvarena visoka prodaja i u tom slucaju pravimo gresku na 42 observacije
# dole vidimo na 3) situaciju kada je ShelveLoc=Good to je druga grana
# ovakvih obs ima 65 i ima malo vise ovih koji su ostvarili visoku prodaju oko 58.4%
# i zato predvidjamo da ce sve one kod kojih je ShelvLoc=good ostvariti visoku prodaju
# i u tom slucaju pravimo gresku na 27 observacija...
# mi cemo ovo lepse preko plota sagledati
# load rpart.plot library
# install.packages("rpart.plot")
install.packages("rpart.plot")
library(rpart.plot)
?rpart.plot
# plot the tree
rpart.plot(tree1, extra = 104)
# sta je extra 104, ova stotka koja figurise u brojci obezbedjuje da u okviru
# svakog cvora vidimo koji je to % observac koji
# su stigli do tog cvora... npro na prvom cvoru stiglo je 100% observacija
# pa onda posle prvog granjanja, na no je stiglo 80% observacija a 20% na yes itd...
# a ovo 4 nam daje koliko je procentualno zastupljena svaka klasa
# npr u prvom cvoru 75% obzerv pripada klasi no a 25% klasi yes itd...
# zavisno od brojke 3,5,6,7 druge ce informacije biti prikazane...
rpart.plot(tree1, extra = 103)
# na cvoru pise Yes ili No, i to predstavlja dominantnu klasu
# vidimo da na provm cvoru ima 75% observ u klasi No, a svega 25% u klasi Yes
# zato je on na prvom mestu upisao No, jer je dominantna, tj vise zastupljena
# zelenom bojom je ofarbao one cvorove u kojima se primecuje da su Yes klase zastupljenije,
# tj veci broj observacija ima vrednost Yes!
###
# KRAJ PRVOG KLIPA
###
# make predictions with tree1 over test dataset
# na isti nacin pravimo predikcije
tree1Pred <- predict(tree1, newdata = testData)
head(tree1Pred)
# dobijamo verovatnoce pripadanja jedno klasi i verovatnoce pripadanja drugoj klasi
# modifikujemo gornju fju da nam odmah kaze koja je to klasa kojoj pridruzuje odredjenu observaciju
# po defaultu uzima da je neki prag verovatnoce 0.5 i ako je p pozitivne klase veca od 0.5
# onda ce svrstati, tj pridruziti labelu pozitivnoj klasi
tree1Pred <- predict(tree1, newdata = testData, type = 'class')
head(tree1Pred)
# create the confusion matrix
# treba odradimo evaluaciju modela kojeg smo kreirali
# to radimo nekim evolaluacionim metrikama, kao one RSS, R2, RMSE sto smo radili kod lin regresije
# samo ovde imamo neke druge metrike koji ce nam pomoci prilikom evaluacije
# pomocu njih smo gldali koliko je neki model dobar
# zato sada radimo ovu matricu zbunjenosti koja na neki nacin uksrsta
# stvarne i predvidjene vrednosti i na neki nacin ova zbunjenost dolazi od toga koliko je nas model
# zbunjen ovim predikcijama i samim observacijama i tim koliko dobro njih radi
cm1 <- table(actual = testData$HighSales,
predicted = tree1Pred)
cm1
# po redovima:
# 53 + 7 su observ koje realno pripadaju klasi No
# 8 + 11 su observ koje realno pripadaju klasi Yes
# po kolonama:
# 53 + 8 su observ za koje smo predvideli da pripadaju klasi No
# 7 + 11 su observ za koje smo predvideli da pripadaju klasi Yes
# sta je pozitivna a sta negativna klasa?
# za poz klasu proglasavamo onu koja nam je posebno znacajna da je dobro predvidimo
# npr kod prod sedista, mozemo da smatramo da je menadzmenu kompanije veoma vazno da
# predvide one koje nece da ostvare visoku prodaju, jer npr ako pretpostavljaju to onda mogu
# da rade marketing akcije u toj prodavnici da animiraju kupce i da onda tu poboljsaju prodaju
# eto to je jedan primer recimo
# ili na primer kod medicinskih istrazivanja tj medicinskih primena masinskog ucenja se cesto za
# pozitivnu klasu smatra utvrdjivanje prisutnosti bolesti
# naravno ako je neko bolestan nije pozitivna stvar ali u modelu nam je vazno da to znamo
# stoga ce to obicno biti pozitivna klasa
# npr ako je neko inficiran kovidom, jako je vazno da prepoznamo da je neko bolestan
# zato je to pozitivna klasa!!! ona koja je vazna da se ispituje!!!
# znaci ovde pozitivno i negativno nema veze sa etikom!
# ovde ce nama pozitivna klasa biti klasa No, odnosno prodavnica koja
# nije ostvarila visoku prodaju, tj ne ocekujemo da ce ostvariti visoku prodaju, da mi to
# stavimo za pozitivnu klasu, jer takva informacija je posebno znacajna nekom menadzmentu
# kako bi ranije odreagovala itd...
# na ispitu ce eksplicitno biti receno sta da uzmemo za pozitivnu klasu!!!
# ovih 53 iz nase conf matrice, 53 prodavnice stvarno nisu ostvarile visoku prodaju i
# mi smo predvideli da nece ostvariti visoku prodaju, zato su ovo TRUE POSITIVE observacije
# ova 11 su TRUE NEGATIVE, stvarno su ostvarile visoku prodaju i mi smo to predvideli
# a ovo na gl dijagonali su nasa tacna predvidjanja
# na suprotnoj dijagonli su greske koje smo pravili prilikom predvidjanja
# ovih 8 smo predvideli da te prodavnice pripadaju pozit klasi a one u stvari ne pripadaju pozitivnoj klasi
# pa su one FALSE POSITIVE
# a ovih 7 mi smo predvideli da pripadaju negativnoj klasi a one realno pripadaju pozitivnoj klasi jer
# nisu ostvarile tu visoku prodaju i zato su FALSE NEGATIVE
# ukratko da ponovim na ovom primeru:
# 53 su TRUE POSITIVE
# 11 su TRUE NEGATIVE
# 8 su FALSE POSITIVE
# 7 su FALSE NEGATIVE
# BURAYERUUUUUUU
# OVe 4 stavke, ove cetiri klase obzervacija predstavljaju OSNOVU
# za definisanje evalucionih metrika, i imamo neke 4 metrike
# 1. ACCURACY tj TACNOST to je accuracu = (TP+TN)/M gde je M ukupan broj observacija
# accuracy ona predstavlja odnos tacno predvidjenih i ukupnih broja obzervacija M
# 2. PRECISION tj PRECIZNOST to je precision = TP/(TP+FP)
# precision je odnos tacno predvidjenih kroz ukupan broj svih onih koje smo predvideli da su pozitivne
# medju svim prodavnicama koje smo predvideli da nisu ostvarile visoku prodaju, koji je udeo onih koje stvarno nisu
# ostvarile visoku prodaju
# 3. RECALL tj ODZIV to je recall = TP/(TP + FN)
# recall je odnos tacno predvidjenih pozitivnih u odnosu na sve koje su pozitivne
# kad god pokusamo da povecamo preciznost obicno ce odziv da se umanjuje
# i obratno, nekako su ove dve metrike u antagonizmu
# mera koja ih objedinjuje se zove F1
# ona se definise preko preciznosti i odziva
# 4. F1 = 2*Preciznost*Recall/(preciznost + recall)
# function for computing evaluation measures
# obratiti paznju da ono sto je poslednje u funkciji to ce funkcija vratiti kao return vrednost!!! VAZNO VAZNO
computeEvalMeasures <- function(cm){
#a <- (cm[1,1] + cm[2,2])/nrow(cm)
a <- sum(diag(cm))/sum(cm)
p <- cm[1,1]/sum(cm[,1])
r <- cm[1,1]/sum(cm[1,])
f1 <- 2 * p * r / (p + r)
c(accuracy = a, precision = p, recall = r, F1 = f1)
}
# compute the evaluation metrics
eval1 <- computeEvalMeasures(cm1)
eval1
# nisu nam losi rezultati nasih evaluacionih metrika
# mada ipak pravimo modifikaciju da probamo dobiti jos bolji model
# generalno kad god pravimo model, ne zadrzavamo se samo na prvoj verziji modela
# nego gledamo sta bismo mogli da modifikujemo kako bismo dobili jos bolji model!!!
# mi cemo sada pokusati da malo modifikujemo parametre naseg algoritma stabla odlucivanja
# tj rpart funkcije, kako bi videli mozemo li bolje stablo napraviti
# mi smo ga pustili da se napravi stablo na neki default nacin
# kad smo pozivali rpart nista mu dodatno nismo zadavali
# get the docs for rpart.control function
?rpart.control
# npr minsplit sto je manji minsplit to mi dozvoljavamo stablu da vise raste da bude slozenije
# minbucket ako dozvolimo da ima manje obser u tim terminalnim cvorovima mi dozvoljavamo da se stablo
# nekako vise razgrana itd...
# cp najcesce koriscen parametar za kontrolisanje rasta stabla
# neki parametar kompleksnosti stabla
# koji je to neki minimalni doprinos ukupnoj funkciji cilja (koja kreira model) koji grananje na nekom
# cvoru mora da da, da bi se ono realizovalo!!
# i naravno sto je taj doprinos manji to mi prakticno dozvoljavamo na vise mesta da se stablo grana
# tj da stablo bude sve kompleksnije i kompleksnije
# u suprotnom ako povecamo vrednost min doprinosa cvora za grananje on ce se teze granati i
# stablo ce biti kompaktije
# mojim recima, KOLIKO MINIMALNO RADA TREBA DA ANGAZUJE JEDAN CVOR DA BI SE DALJE DELIO!!! AKO SE MIN
# RAD POVECA TO CE SE STABLO TEZE GRANATI I BICE KOMPAKTINIJE I OBRATNO
# build second model with minsplit = 10 and cp = 0.001
# minsplit je default 20, a cp za red velicine manju vrednost za 0.001
# i time ocekujemo da dobijemo vece i kompleksnije stablo
# postoji potencijalni problem sa kompleksnim stablima i modelima jer oni
# imaju tendenciju da pokazu jako dobre podatke na trening setu
# ali kada im damo nove podatke, slabo se snalaze, i slabo ih uklapaju u svoja definisana pravila
# na adekvatan nacin, nemaju dovoljnu sposobnost generalizacije i
# tada kazemo da je model PRETRENIRAN odnosno OVERFITOVAN
set.seed(1)
tree2 <- rpart(HighSales ~ ., data = trainData, control = rpart.control(minsplit = 10, cp = 0.001))
# print model
print(tree2)
# ovaj bogami deluje dosta vece i kompleksnije
# plot the model
rpart.plot(tree2, extra = 104)
# make predictions with tree2 over test dataset
tree2Pred <- predict(tree2, newdata = testData, type = 'class')
# create confusion matrix for tree2 predictions
cm2 <- table(actual = testData$HighSales,
predicted = tree2Pred)
cm2
# compute eval metrics
eval2 <- computeEvalMeasures(cm2)
eval2
# compare evaluation metrics for tree1 and tree2
# bolje su ove nove metrike
data.frame(rbind(eval1, eval2), row.names = paste0("Tree_", 1:2))
# zasto smo stavili cp 0.001 i minsplit 10 pa da znamo da ce biti bolji rezultat?
# u principu i nismo znali... probali smo hehe
# mada ima neki pristup sistematican kako mozemo da odredimo koje su vrednosti parametara koji
# ce nam doneti bolje rezultate
# to se zove TUNING PARAMETARA, podesavanje parametara i radi se kroz postupak cross-validacije
# cross-validation JAKO JE BITNAAAA
# stalno se koristi
# kros validacija nam omogucuje da taj neki dataset koji nam je na raspolaganju
# iskoristimo max efektno za utvrdjivanje performansi modela
# i za isprobavanje razlicitih vrednosti parametara, kao sto
# mi sada hocemo da uradimo za onaj cp parametar!!!
# sta podrazumeva
# podrazumeva da uzimamo onaj skup podataka koji imamo za trening i delimo ga na neki k delova
# to k na ovoj slici koju profesorka prikazuje je ovde 5, 5 delova
# moze da bude k = 10 moze sta god ali obicno se uzima k = 5 iil k = 10
# sta radimo
# uzmemo k-1, tj 4 dela za treniranje modela
# i uzmemo 1 deo za izracunavanje eval metrika, tj za evaluaciju modela
# najcesce se uzima ona accuracy tacnost za klasifikaciju
# i to se ponovi k puta, tj 5 puta...
# ali u svakom ponavljanju drugi se segment uzima za evaluaciju
# pri svakoj iteraciji zapisujemo performanse koje smo dobili
# kao da smo dobili tacnost 1, pa tacnost 2 pa itd do 5
# i kada ih dobijemo sve mi ih sve zajedno uprosecimo
# zasto?
# pa da bismo dobili stabilnije procene performasni naseg modela...
# zasto je to bitno?
# pa posto se ovi podaci nasumicno biraju na test i trening u ovih 5 iteracija
# moze da se desi da u jednoj iteraciji imamo jako dobre rezultate jer su podaci dobro podeljeni
# ali mogu i lose podeljeni podaci da uticu na lose podatke
# pa da bismo to neutralisali, mi dobijene rezultate uprosecimo...
# prakticno uprosecavanje nam daje stabilniju vrednost
# zbog toga se cesce radi 10 iteracaij tzv 10-fold cross-validation, pa sve to uprosecimo!
# podesavacemo cp parametar pa radimo cross-validaciju pa ponovo
# sad ce pomoci caret paket
# load e1071 library
install.packages('e1071')
# ovo je pomocni paket za kros validaciju
library(e1071)
library(caret)
# define cross validation (cv) parameters;
# we'll do 10-fold cross-validation
tr_ctrl <- trainControl(method = "cv", number = 10)
# then, define range of cp values to examine in cross-validation
# fokusiramo se na cp param jer je kljucan za kontrolisanje rasta stabla
ls()
abe <- 0
rm(abe)
# koje su to vrednosti parametra kojeg hocemo da razmotrimo
# ova tacka pre cp to je prosto zahtev ove funkcije da na taj nacin
# identifikujemo parametre modela
# i kazemo da ovde hocu u nekoj sekvenci vrednosti od 0.001 do 0.02 sa korakom od 0.0005
# to kaze uzimacu cp od 0.001, poslednja vrednost koju cu razmatrati je 0.02 a
# izmedju cu uzimati vrednosti sa korakom od 0.0005, zasto?
# zato sto hocu da uzmem mnogo razlicitih vrednosti za taj cp da sto preciznije odredimo
# koja je to vrednost za cp koja ce mi dati najbolje rezultate
# i to zapisujemo u neku prom cp_grid
cp_grid <- expand.grid( .cp = seq(0.001, 0.02, 0.0005))
# i sada treba jos da pozovemo trening proces sa ovim novim parametrima
# since cross-validation is a probabilistic process, we need to set the seed
# so that the results can be replicated
# tress cross validated tree_cv
# sta su ulazni argumenti, to su sve varijable iz trening seta sem izlazne var HighSales, zasto? Ne znam?
# y sta je to sto predvidjamo? aham ovde ubacujemo 11tu kolonu... sad je jasno heheh
# method tu stavljamo koju metodu hocemo da koristimo za kreiranje klasifikacionog modela
# i to je rpart
# kako hocu da kontrolisem ovaj proces treniranja, to je ovo sto smo prethodno kreirali
# ova kros validacija i 10 iteracija, pa kazem trControl je nas tr_ctrl
# i kazem jos kako hocu da radim TUNING parametera to je ono cp_grid
# i eventualno hocu da zadrzim da mi minsplit bude 10, a ne 20 kako je po defaultu
# inace ovaj minsplit parametar ove train fje, ona kad prepozna da to nije njen parametar
# ona ce da shvati da je to parametar ovoj funkciji rpart koju smo zadali kao method!
# o tome se sve brine ova train funkcija iz careta!
set.seed(1)
tree_cv <- train(x = trainData[,-11],
y = trainData$HighSales,
method = 'rpart',
trControl = tr_ctrl,
tuneGrid = cp_grid,
minsplit = 10
)
# relativno se brzo izvrsilo iako je radio brojna izracunavanja
tree_cv
# on je ovde uzeo cp vrednost, sve one raznorazne opcije koje smo mu zadali da razmatra
# i za svaku od njih je odradi cross-validaciju i izracunao 2 metrike
# metriku Accuracy i Kappa (ovom se necemo baviti)
# mi se fokus na Accuracy, i generalno on ovde bira cp vrednost na osnovu Accuracy metrike
# i kaze nam
# Accuracy was used to select the optimal model using the largest value.
# The final value used for the model was cp = 0.0135.
# to moze i kroz plot da se vidi
# run the cross validation
# plot the cross-validation results
plot(tree_cv)
# kao najbolju vrednost za cp je dao ovu poslednju zasto?
# jer sto je veci cp, to je stablo jednostavnije, manje slozeno
# cim je ono jednostavnije, ono je manje sklono podudaranju sa podacima za trening i
# bolje generalizcuje
# uvek cemo se truditi da imamo sto jednostavnije modele tj jednostavnija stabla
# create tree3 using the new cp value
best_cp <- tree_cv$bestTune$cp
set.seed(7)
tree3 <- rpart(HighSales ~ ., data = trainData, method = 'class',
control = rpart.control(minsplit = 10, cp = best_cp))
# plot new tree
plot(tree3)
rpart.plot(tree3, extra = 104)
# make the predictions with tree3 over the test dataset
tree3Pred <- predict(tree3, newdata = testData, type = 'class')
# create the confusion matrix for tree3 predictions
cm3 <- table(actual = testData$HighSales,
predicted = tree3Pred)
cm3
# compute the evaluation metrics
eval3 <- computeEvalMeasures(cm3)
eval3
# compare the evaluation metrics for tree1, tree2 and tree3
data.frame(rbind(eval1, eval2, eval3), row.names = paste0("Tree_", 1:3))
# pitanje studenta: kako znamo da treba da nastavimo dalje jer je 88% prec sasvim ok
# kaze profesorka da nekad je i 65% sasvim okej uz najbolji moguci model
# zavisi od problema i kakvi su podaci,
# prosto moramo uvek da poredimo sa nekim modelom, znaci po defaulto treba da uradimo nekoliko modela
# pa da ih medjusobno poredimo
# i da svaki model bude po necemu razlicit
# po nekim karakteristikama, dal ce to biti varijable koje su ukljucene u model
# ili ce to biti neki drugi algoritam koji cemo da koristimo
# npr ne stablo odlucivanja nego recimo k najblizih suseda itd...
# ili siti algo ali sa drugim varijablama, parametrima.. itd..
# zato se proba nekoliko razlicitih modela da bismo bili sigurni
# da ono sto smo dobili bude najbolje moguce
# zavisi i od dataseta samog
# bar tri modela isprobamo pa teramo kera
# ili ako resavamo neki siroko poznat zadatak i postoje objavljeni rezultati
# onda vidimo kako se nasi rezultati matchuju u odnosu na te vec objavlejne
# i na osnovu toga mozemo znati dal je nasih 88 posto za tacnost je stvarno dobra vrednost
# ili treba probati dalje
#####
# END CLASSIFICATION TREES
#####
#####
# KNN Neighbours
#####
# knn k najblizih suseda
# ovo isto klasifikacija
# samo ne primenom stabala odlucivanja
# nego k nalbizih suseda knn
# load ISLR
library(ISLR)
# koristicemo Carseats
# pa da uporedimo sa onim rezultatima iz stabla
# pa cemo algoritam naivnog bajesa koristiti kasnije
# pa cemo videti koji algoritam daje najbolje rezultate za ovaj dataset
# obicno se primene nekoliko razlicita i biramo model koji daje najbolje performanse
# ili menjamo parametre...
# print dataset structure
print(Carseats)
?Carseats
# calc 3rd quartile
sales_3Q <- quantile(Carseats$Sales, probs = 0.75)
# create a new var HighSales based on 3rd quartile
Carseats$HighSales <- ifelse(test = Carseats$Sales > sales_3Q,
yes = "Yes",
no = "No")
Carseats$HighSales <- as.factor(Carseats$HighSales)
# jer mora bude faktor ne moze char
# provera distribucije
table(Carseats$HighSales)
# procentualna zastupljenost:
prop.table(table(Carseats$HighSales))
# remove the Sales variable
# posto su 100% asocirane, pa ne mozemo Sales koristiti za predvidjanje, pa je izbacujemo
Carseats$Sales <- NULL
###
# Standardize numerical attributes
###
# za razliku od stabla odluc gde mozemo razne vrste podataka da koristimo
# knn trazi iskljucimo numericke podatke
# knn je krajnje jednostavan
# on odmah razmatra novu obzervaciju u kontekstu trening seta
# racuna se udaljenost nove obser od svih ostalih obser u trening setu
# na osnovu odstojanja pronalazimo one obser koje su najblize
# i pronalazimo k najblizzih zato se tako i zove ovaj algoritam
# to kolko je k to je jedan parametar modela kao onaj cp parametr kod stabla sto smo racunali
# tako cemo i ovde morati to da sracunamo
# uvek se za k uzima neparan broj da bi jedna klasa bila dominantnija u odnosu na drugu
# zbog racunanja odstojanja izmedju obser iz trening seta i obser za koje predvidjamo, moramo imati num
# vrednosti, dakle zbog euklidskog rastojanja koja trazi numericke vrednosti
# standardizacija
# ako ima smisla varijable transformisemo
# ako nema smisla onda ih odbacujemo
# print summary of dataset
summary(Carseats)
# ima 7 num var
# jos jedan uslov za ove var
# to je da se njihove vrednosti nalaze u priblizno jednakim opsezima
# to je uslov za dalji rad
# moramo ih transformisati ako ne zadovoljavaju
# sta je opseg vrednosti, to je ovo od min do max vrednosti
# potrebno da ove var budu priblizno jednakih opsega
# bitno zbog racunanja odstojanja
# jer ako su neke var sa velikim opsezima vrednosti, a druge jako male,
# onda ce ove sa velikim opsezima postati dominantne za racunanje odstojanja
# i onda ce biti bayes prema tim pristrasnost prema tim varijablama
# zato svodimo varijable na priblizno jednake opsege
# compPrice i income imaju priblizno jednake opsege
# ali npr adv i population imaju znatnu razliku ovi od 0 do 29 a ovi od 10 do 509
# to moramo menjati
# moramo SKALIRATI
# za SKALIRANJE imamo 2 pristupa:
# 1. NORMALIZACIJA
# 2. STANDARDIZACIJA
# STANDARDIZAIJA od svake vrednosti oduzimamo prosecnu vrednost i delimo sa merom varijabiliteta
# dal ce ta prosecna vr biti mean ili median i dal ce ono sa cim delimo biti standardna devijac ili
# nesto drugo, to zavisi od toga dal je varijabla normalno raspodeljena ili ne
# NORMALIZACIJA ona podraz da svodimo sve vrednosti varijable na jedan isti opseg vrednosti koji je
# tipicno 0 1, ne mora biti, moze biti -1 1 itd...
# kada koristimo koji pristup? jel postoji neki uslov?
# Prisustvo outliyara
# ako ima outliera u okviru vrednosti varijabla, ne preporucuje se NORMALIZACIJA jer
# ako medju tim outlierima imamo neke velike vrednosti, te velike treba da svedemo na max mogucu
# koja je vrednos 1, ove znacajno manje vrednosti treba da se svode na neke jako male vrednosti,
# al se te male vrednosti toliko smanje da se njihov uticaj skoro potpuno izgubi
# tako da je neka praksa jeste da ukoliko imamo OUTLIERE u varijabli, onda te varijable skaliramo
# STANDARDIZACIJOM, a ako nemamo OUTLIERE onda je svejedno koji cemo pristup al tipicno NORMALIZACIJA
# plot boxplot for CompPrice var
boxplot(Carseats$CompPrice)
# debelea linija je mediana
# gornja linija boxa je 3 kvartil
# donja linija boxa je 1 kvartil
# tacke su outlieri
# gornja crta je max
# donja crta je min
boxplot(Carseats$Income)
boxplot.stats(Carseats$CompPrice)
# 85 je min vr
# 115 je 1q
# 125 je mediana
# 135 je 3q
# 162 je max
# 400 je broj obser
# conf je interval poverenja za vrednost mediane
# out su outlieri
length(boxplot.stats(Carseats$CompPrice)$out)
# pring number of outliers in the CompPrice var
length(boxplot.stats(Carseats$CompPrice)$out)
# select num vars
num_vars <- c(1:5, 7,8)
# apply fja koja return num of outliers to each num column
# margina je 2 znaci primeni ovu funkciju po kolonama
apply(Carseats[,num_vars], 2,
function(x) length(boxplot.stats(x)$out))
# za kompprice i price da primenimo standardizaciju a za ostale normalizaciju
# ali cemo za sve uraditi standardizaciju da ne bi sada komplikovali
# e sad da povedemo racuna dal su var normalno raspodeljene ili nisu i da u skladu sa time
# uradimo na jedan ili drugi nacin standardizaciju
# za utvrdjivanje dal je normalna rasp koristimo shapiro.test
?shapiro.test
# apply shapiro-wilk test to each numerical column (var)
# test koji testira nultu hipotezu da neki uzorak dolazi iz populacije sa normal raspodelom
# i ukoliko p vrednost bude jako mala na testu, onda znaci da je jako mala verovatnosca
# da ta hipoteza moze da zadrzi i smatracemo da var nema normalnu raspod
# a ako je p veca od 0.05 smatracemo da nultu hipotezu ne mozem oda odbacimo i
# smatracemo da var ima normalnu raspodelu
apply(Carseats[,num_vars], 2, shapiro.test)
# za compprice ne mozemo odbaciti nultu hipotezu pa znaci ima norm raspodelu
# z aincome jako mala p vrednost, znaci mala verovatnoca da nulta hipo vazi te da nema norm raspo
# isto vazi i za adverising
# isto i za popul
# price ima normalnu
# Age nema
# Educ nema
# zakljucak
# price i compprice imaju norm raspodelu a ostale nemaju
# zavisno od toga kakva je rasp koristimo medianu ili mean odnosno
# stad devij ili interkvartilni opseg za skaliranje
# korsitimo fju scale iz baznog r paketa
# get the doc for scale fun
?scale
# select non-normally distributed num columns
non_norm_vars <- c(2,3,4,7,8)
# apply the scalling fun to each column
transformed_data <-apply(Carseats[,non_norm_vars], 2,
function(a) scale(x = a, center = median(a), scale = IQR(a)))
# ovo je skaliranje na ne normalne vars
# zato smo koristili medianu i inter quartilni opseg to je razlika izmedju 3rd i 1st kvartila i neka
# mera varijabiliteta u ovom slucaju ne normalne raspodele
# since apply() f returns a list, convert it to a data frame
transformed_data <- as.data.frame(transformed_data)
# stadardize Price var ( and convert to vector)
transformed_data$Price <- as.vector(scale(Carseats$Price,
center = mean(Carseats$Price),
scale = sd(Carseats$Price)))
# stadardize CompPrice var ( and convert to vector)
transformed_data$CompPrice <- as.vector(scale(Carseats$CompPrice,
center = mean(Carseats$CompPrice),
scale = sd(Carseats$CompPrice)))
# ovo center = TRUE znaci da centriranje radi preko mean vrednosti po defaultu
# znaci nismo morali da zadajemo posebno ovo center
# scale = TRUE, znaci da skaliranje radi koristeci stand devijaciju po defaultu ukoliko je
# centriranje bilo sa mean vrednoscu, pa ce gornja fja izgledati vako:
transformed_data$CompPrice <- as.vector(scale(Carseats$CompPrice))
summary(transformed_data)
# ne nalaze se svi u istom opsegu jer nismo radili NORMALIZACIJU, ali se nalaze u
# priblizno jednakom opsegu jer smo radili STANDARDIZACIJU
###
# KRAJ KLIPA
###
###
# Transform factor (binary and categorical) vars
###
# treba biti obazriv jer nekad ne mozemo da transformisemo a nekad mozemo
# kod binarnih faktroskih varijabli moze lagano jer su jednostavne
# transform Urban var to int
transformed_data$Urban <- as.integer(Carseats$Urban) -1
table(transformed_data$Urban)
# transform US var to int
transformed_data$US <- as.integer(Carseats$US) -1
table(transformed_data$US)
# examine levels of ShelveLoc var
levels(Carseats$ShelveLoc)
# postoji gradacija od loseg ka dobrom
# ako imaju upitnici npr dal se slazemo, dal se neslazemo, potpuno se slazemo itd to je fakt var
# i takva bi mogla da se transform u num
# ali ako imamo var koja predstvalja boju
# tu nikakvog smisla ne bi bilo da pridruzujemo numericku vrednost
# da oduzimamo zutu od plave i takve neke stvari
# vrsta materijala sedista isto ne bismo mogli da transformisemo u brojke
# ako na to naidjemo na ispotu na kolokvijumu, nju prosto izostavimo...
# to nije bas sasvim tacno da se ne moze transformisati ali dobro
# postoji nesto sto se zove dummy varijable
# npr boje
# kvalitativne varijable da transformisemo u dummy varijable
# npr crvena zelena plava... drugi put...nikad, ikad kad..
# update order of lvls for the ShelveLoc var to Bad, Med and Good
transformed_data$ShelveLoc <- factor(Carseats$ShelveLoc,
levels = c("Bad", "Medium", "Good"))
levels(transformed_data$ShelveLoc)
# convert tu num var shelvloc
transformed_data$ShelveLoc <- as.integer(transformed_data$ShelveLoc)
table(transformed_data$ShelveLoc)
# dodati outcome HighSales to transformed dataset
transformed_data$HighSales <- Carseats$HighSales
summary(transformed_data)
# svi su num sem ove izlazne koja je factorska
# sve var su u pribliznom opsegu vredsnoti
# to nam daje osnovu za koriscenje KNN algoritma
###
# Create train and test data sets
###
#load caret
library(caret)
#set.seed
set.seed(14)
# create train and test sets
train_indices <- createDataPartition(y = transformed_data$HighSales,
p = 0.8,
list = F)
train_data <- transformed_data[train_indices,]
test_data <- transformed_data[-train_indices,]
###
# Model building
###
# load class package for KNN
library(class)
?knn
# create knn model with k=5
# ipak cemo mi k nasumice izbrati
set.seed(14)
rand_k <- sample(seq(3, 21, 2), 1)
knn1_pred <- knn(train = train_data[,-11],
test = test_data[,-11],
cl = train_data$HighSales,
k = rand_k)
# stavili smo pred zato sto ce nam knn fja vratiti predikcije za test set
# print several predictions
head(knn1_pred)
# create confus matrix
cm1 <- table(actual = test_data$HighSales,
predicted = knn1_pred)
cm1
# ne pravimo false negative gresku :)
# precision ce nam biti max moguca
# fun for computing eval metrics
compute_eval_metrics <- function(cm){
TP <- cm[1,1]
TN <- cm[2,2]
FP <- cm[2,1]
FN <- cm[1,2]
a <- (TP+TN)/sum(cm)
p <- TP/(TP+FP)
r <- TP/(TP+FN)
f1 <- 2*p*r/(p+r)
c(accuracy = a, precision = p, recall = r, F1 = f1)
}
# compute eval metrics
eval1 <- compute_eval_metrics(cm1)
eval1
# find optimal val for k through 10-fold cross-validation
# load 1071 lib
library(e1071)
# define cross valid (cv) parameters, we'll perform 10-fold crossval
train_ctrl <- trainControl(method = "cv", number = 10)
# define range for k values to examine in cross val
k_grid <- expand.grid(.k = seq(3,25,2))
# since cross val is probailastic postavi seed
set.seed(14)
knn_cv <- train(x = train_data[,-11],
y = train_data$HighSales,
method = 'knn',
trControl = train_ctrl,
tuneGrid = k_grid)
# plotuj cross validation results
plot(knn_cv)
# na x broj suseda, na y tacnost koja koriscena kao mera performansi u kros validaciji
# peak za ovaj accuray dostize za vrednost 9
# build a new model with best value for k
best_k <- knn_cv$bestTune$k
knn2_pred <- knn(train = train_data[,-11],
test = test_data[,-11],
cl = train_data$HighSales,
k = best_k)
# create confusion matrix
cm2 <- table(actual = test_data$HighSales,
predicted = knn2_pred)
cm2
# compute eval metrics
eval2 <- compute_eval_metrics(cm2)
eval2
data.frame(rbind(eval1, eval2),
row.names = paste0("K=", c(rand_k, best_k)))
# imamo bolji model!!! opala
# cross validacija se uvek radi na train delu podataka
# i on te podatke podeli na deo za obucavanje i deo za validaciju - izracunavanje eval metrika
# i onda kroz vise iteracija svaki put uzima drugi deo na kome ce raditi tu evaluaciju
#####
# Kraj knn najblizih suseda
#####
#####
# NAIVE BAYES
#####
# klasifikacija primenom naivnog bajesa
# uprkos jednostavnosti pokazuje dosta dobre rezultate
# zato ga dosta koristimo
###
# Naive Bayes Classifier
###
# koristicemo Carseats
# treba da ga pripremimo kao i na prethodnim casovima
# ponavljanje...
# ucitavanje r scripte
# load utility functions
source('util.R')
# get adapted Carseates data set
carseats <- get_adapted_carseats_dataset()
###
# Discretise numerical variables
###
# naivni bajes radi sa faktorima, a sa num moze ako imaju norm rasp
# ako nemaju norm rasp onda ih moramo diskretizovati tj pretvoriti u faktorsku var
# ajmo malo transformisati dataset
# vidimo situaciju sa num var prvo, dal su norm rasp ili ne
# faktorske mogu odmah udju u model lagano
# selektovanje num varova
num_vars <- c(1:5,7,8)
# provera dal su norm raspod
# shapiro testom naranvo
# margina je 2 znaci primenjujemo fju na kolone
apply(carseats[,num_vars], 2, shapiro.test)
# nulta hipoteza je da uzorak dolaze iz norm rasp
# ako je P jako mala onda nije norm rasp
# education nema norm rasp, nema ni age, price ima jer je veca od 0.05, popul nije, income nije itd
# fokus je da transform ove koji nemaju norm rasp
# njih treba da diskretizujemo
# da vektor prebvacimo u faktor
# diskretizacija podrazumeva da num var pretvorimo u fact var na nacin da izdelimo
# ukupan opseg vrednosti var na nekoliko intervala i onda
# svaki od tih intervala ce predstavljati jedan level novog faktora
# moze na puno nacina da se uradi ovo
# nadgledana ili nenadgledana diskretizacija
# nadgledana je kad uzimamo u obzir izlaznu var i na osnovu toga metode rade diskretizaciju
# to je slozenije necemo raditi
# mi radimo samo nenadgledano diskretizaciju
# razlika je sto ovde se ne uzima izlazna varijabla kada se radi ta podela na intervale!!!
# radimo ovo sa paketom bnlearn
install.packages('bnlearn')
library(bnlearn)
?discretize
# trazi metodu, jedna je intervalna, a druga je kvantilna diskretizacija
# u oba slucaja uzimamo ukupan opseg od min do max i delimo na k intervala
# kod intervalne delimo tako da ti k intervali budu jednake duzine
# kod kvantilne delimo tako da u okviru tih k intervala imamo jednak broj opservacija, tj jednaku ucestalsot
# metoda kvantilne se pokazala sa boljim rezultatima pa cemo nju koristiti
# sta je to k to je ovo breaks
# obicno se uzima 5 ili 10 pa vidimo sta dobijemo, koje rezultate, kako izgledaju diskretizovane varijable
# mi cemo uzmemo 5 pa da vidimo dal je adekvatno pa cemo menjati ako treba
# slicno onom cp parametru ili onom k parametru iz knn-a
to_discretize <- c(2,3,4,7,8)
# diskretizujemo:
# quantile method ne moramo da navodimo ekspl jer je to default
transformed_data <- discretize(carseats[,to_discretize],
breaks = 5)
summary(transformed_data)
# advertizing je neujednacen iako ga ja nisam uradio jer baca neku greskku
# starije verzije bnlearna bacaju gresku za var koji su neujednaceni imaju preveliki i preniske val
# bas evo pominje profesorka, znaci ja imam neku stariju verziju bnlearna
update.packages('bnlearn')
# nece radi update
# ali ajd sad cemo Advertising da isprovervaamo by plotting histogram
library(ggplot2)
ggplot(carseats,
aes(x = Advertising)) +
geom_histogram() +
theme_minimal()
# probacemo 3 intervala za Advertising var a za ostale da ostavimo 5
transformed_data <- discretize(carseats[,to_discretize],
breaks = c(5,2,5,5,5))
summary(transformed_data)
# meni nece sa tri intervala pa sam stavio 2 jer jedino to hoce da radi
# create vector of var names to be added to data frame with discretised variables
to_add_col <- setdiff(colnames(carseats),colnames(transformed_data))
to_add_col
# erge discret data frame with other cols from original one
transformed_data <- cbind(transformed_data,
carseats[,to_add_col])
# update order variabli isto kao iriginalno
transformed_data <- transformed_data[,colnames(carseats)]
# JOS JEDNOM REXZIME
# za var koje su numericke i nemaju norm raspodelu radimo diskretizaciju
# dva oblika diskretizacija to su intervalna i kvantilna
# mi smo koristili kvantilnu koja obicno daje bolje rezultate, ali nije nuzno
# podela trening i test kao i inace
library(caret)
set.seed(8421)
train_indices <- createDataPartition(transformed_data$HighSales,
p = 0.8,
list = F)
train_data <- transformed_data[train_indices,]
test_data <- transformed_data[-train_indices,]
# kreiranje modela
###
library(e1071)
?naiveBayes
nb1 <- naiveBayes(HighSales ~ ., data = train_data)
print(nb1)
# a-priori verovatnoce
# to su verovatnoce za da je ostvareno niska i visoka prodaja
# to je ono P(c) u formula bajes
# a kako da izracunamo P(x|c) uslovne verovatnoce
# to radimo preko naivnog bajesa
# uvodi se naivna pretpostavka za izracunavanje uslovne verovatnoce
# zato se zove naivni bajes
# a to je pretpostavka da su var koje opisuju ove nase prodavnice medjusobno nezavisne
# da je cena nezavisna od urbane regije, od advertisinga i tako dalje
# to je malo nerealna pretpostavka ali je neophodna da bi ova verovatnoca
# koja je slozena za izracunavanja mogla da se predstavi kao proizvod ovih pojedinacnih
# uslovnih verovatnoca
# po tabelicama sto imamo to su pojedinacne uslovne verovatnoce!!!
# sve je faktorski sem kod price i compprice gde su num
# ali sa normalnom raspodelom
# ali su date kao mean pod [1] i standardna devijacija pod [2] vredsnoti
# zasto, zato sto su nam te vrednosti potrebne za izracunavanje fje gustine raspodele koja
# prakticno omogucuje izracunavanje verovatnoce
# zato tu ima vrednosti preko 1
# P(x) njega ne razmatramo jer je on konstantna prilikom otpimizacije fje cilja
# make predictions
nb1_pred <- predict(nb1, newdata = test_data, type = 'class')
head(nb1_pred)
nb1_pred <- predict(nb1, newdata = test_data, type = 'raw')
# print
head(nb1_pred)
# conf matrix
cm1 <- table(actual = test_data$HighSales,
predicted = nb1_pred)
cm1
# eval metrike
eval1 <- compute_eval_metrics(cm1)
eval1
# dal mozemo bolji model
# ranije smo koristili cross validaciju
# kod naivnog bajesa nemamo argumente koje bi mogli da podesavamo TUNING
# ali mozemo napravimo drugi model sa drugim varijablama
# build model with vars thet prooved relevant in decision tree
# namelu ShelveLoc, Price, Advertising, Income, Age, US
nb2 <- naiveBayes(HighSales ~ ShelveLoc + Price + Advertising + Income + Age + US,
data = train_data)
# make the predictions with nb2 model over the test dataset
nb2.pred <- predict(nb2, newdata = test_data, type = 'class')
# create the confusion matrix for nb2 predictions
nb2.cm <- table(true = test_data$HighSales, predicted = nb2.pred)
nb2.cm
# compute the evaluation metrics for the nb2 model
nb2.eval <- compute_eval_metrics(nb2.cm)
nb2.eval
# compare the evaluation metrics for nb1 and nb2
data.frame(rbind(nb1.eval, nb2.eval), row.names = c("NB_1", "NB_2"))
###
# ROC curves
###
# compute probabilities for each class value for the observations in the test set
nb2.pred.prob <- predict(nb2, newdata = test_data, type = "raw") # note that the type parameter is now set to 'raw'
head(nb2.pred.prob)
#install.packages('pROC')
# load pROC package
library(pROC)
# create a ROC curve
nb2.roc <- roc(response = as.numeric(test_data$HighSales),
predictor = nb2.pred.prob[,1],
levels = c(2, 1))
# print the Area Under the Curve (AUC) value
nb2.roc$auc
# plot the ROC curve, using the "youden" method
plot.roc(nb2.roc,
print.thres = TRUE,
print.thres.best.method = "youden")
# get the coordinates for all local maximas
nb2.coords <- coords(nb2.roc,
ret = c("accuracy", "spec", "sens", "thr"),
x = "local maximas", transpose = FALSE)
nb2.coords
# choose a threshold that assures a high value for sensitivity
prob.threshold <- nb2.coords[1,4]
# create predictions based on the new threshold
nb2.pred2 <- ifelse(test = nb2.pred.prob[,1] >= prob.threshold, # if probability of the positive class (No) is greater than the chosen probability threshold ...
yes = "No", #... assign the positive class (No)
no = "Yes") #... otherwise, assign the negative class (Yes)
nb2.pred2 <- as.factor(nb2.pred2)
# create the confusion matrix for the new predictions
nb2.cm2 <- table(actual = test_data$HighSales, predicted = nb2.pred2)
nb2.cm2
# compute the evaluation metrics
nb2.eval2 <- compute_eval_metrics(nb2.cm2)
nb2.eval2
# compare the evaluation metrics for all three models
data.frame(rbind(nb1.eval, nb2.eval, nb2.eval2),
row.names = c(paste("NB_", 1:3, sep = "")))
#####
# Kraj naive Bayes
#####
#####
# 8. i 9. cas
#####
###
# Data preparation and feature engineering
###
###
# Titanic data set
###
# load Titanic train ("data/train.csv") and test sets ("data/test.csv")
titanic.train <- read.csv("data/train.csv", stringsAsFactors = F)
titanic.test <- read.csv("data/test.csv", stringsAsFactors = F)
# print the structure of the train set
str(titanic.train)
# print the structure of the test set
str(titanic.test)
###
# Detecting missing values
###
# print the summary of the train set
summary(titanic.train)
# print the summary of the test set
summary(titanic.test)
# checking the presence of empty strings or other irregular values in character variables in the train set
# first for one variable, and then for all
sum(titanic.train$Cabin=="" | titanic.train$Cabin==" " | titanic.train$Cabin=="-" | is.na(titanic.train$Cabin))
char_vars = c("Name", "Sex", "Cabin", "Embarked", "Ticket")
apply(X = titanic.train[,char_vars],
MARGIN = 2,
FUN = function(x) sum(x=="" | x==" " | x=="-" | is.na(x)))
# do the same in the test set
apply(X = titanic.test[,char_vars],
MARGIN = 2,
FUN = function(x) sum(x=="" | x==" " | x=="-" | is.na(x)))
# check the irregular values present in the Embarked variable
head(sort(unique(titanic.train$Embarked)))
# set the empty Embarked values to NA in the train set
titanic.train$Embarked[titanic.train$Embarked==""] <- NA
# get indices of observations with no Cabin value from the first class, in the train set
train.class1.no.cabin <- which(titanic.train$Pclass==1 & titanic.train$Cabin=="")
length(train.class1.no.cabin)
# get indices of observations with no Cabin value from the first class, in the test set
test.class1.no.cabin <- which(titanic.test$Pclass==1 & titanic.test$Cabin=="")
length(test.class1.no.cabin)
# set the Cabin value for identified passengers to NA in the train and test sets
titanic.train$Cabin[train.class1.no.cabin] <- NA
titanic.test$Cabin[test.class1.no.cabin] <- NA
# print the number of missing Cabin values in the train and test sets
sum(is.na(titanic.train$Cabin))
sum(is.na(titanic.test$Cabin))
###
# Handling missing values
###
###
## Categorical variables with a small number of missing values
###
# create the contingency table for the values of the Embarked variable
xtabs(~Embarked, data = titanic.train)
# replace all NA values for the Embarked variable with 'S' in the train set
titanic.train$Embarked[is.na(titanic.train$Embarked)] <- 'S'
# print the contingency table for the values of the Embarked variable
xtabs(~Embarked, data = titanic.train)
# transform the Embarked variable into a factor in both sets
titanic.train$Embarked <- factor(titanic.train$Embarked,
levels = c("S","C","Q"),
labels = c("Southampton", "Cherbourg", "Queenstown"))
titanic.test$Embarked <- factor(titanic.test$Embarked,
levels = c("S","C","Q"),
labels = c("Southampton", "Cherbourg", "Queenstown"))
###
# Numerical variables with a small number of missing values
###
# test the Fare variable for normality
shapiro.test(titanic.test$Fare)
# get the class of the observation with missing Fare variable
missing.fare.pclass <- titanic.test$Pclass[is.na(titanic.test$Fare)]
# calculate the median value for the Fare variable of all passengers from the 3rd class
median.fare <- median(x = titanic.test$Fare[titanic.test$Pclass == missing.fare.pclass],
na.rm = T) # we have to set this to true as Fare has one NA value
# set the median value to the Fare variable of the passenger with a missing Fare
titanic.test$Fare[is.na(titanic.test$Fare)] <- median.fare
# print the summary of the test set
summary(titanic.test$Fare)
###
# Feature selection
###
###
## Examining the predictive power of variables from the data set
###
# transform the Sex variable into factor
titanic.train$Sex <- factor(titanic.train$Sex)
# get the summary of the Sex variable
summary(titanic.train$Sex)
# compute the proportions table of the Sex variable
prop.table(summary( titanic.train$Sex ))
# create a contingency table for Sex vs. Survived
sex.survived.counts <- xtabs(~Sex + Survived, data = titanic.train)
sex.survived.counts
# transform the Survived variable into factor
titanic.train$Survived <- factor(titanic.train$Survived,
levels = c(0,1), labels = c('No','Yes'))
# create the table again, now labels for Survived will be available
sex.survived.counts <- xtabs(~Sex + Survived, data = titanic.train)
sex.survived.counts
# compute the proportions for Sex vs. Survived
sex.surv.tbl <- prop.table(sex.survived.counts,
margin = 1) # proportions are computed at the row level (each row sums to 1)
sex.surv.tbl
# transform the Pclass variable into factor
titanic.train$Pclass <- factor(titanic.train$Pclass,
levels = c(1,2,3),
labels = c("1st", "2nd", "3rd"))
# plot the number of passengers for different classes and Survived values
library(ggplot2)
gp1 <- ggplot(titanic.train, aes(x = Pclass, fill=Survived)) +
geom_bar(position = "dodge", width = 0.4) +
labs(x="Passenger class", y="Number of passengers") +
theme_bw()
gp1
# add the Sex facet to the plot
gp2 <- gp1 + facet_wrap(~Sex, nrow=2)
gp2
# Instead of counts, plot the proportions
gp3 <- ggplot(titanic.train, aes(x = Pclass, fill=Survived)) +
geom_bar(position = "fill", width = 0.4) + # note that 'fill' is used instead of 'dodge'
labs(x="Passenger class", y="Proportion of passengers") +
theme_bw()
gp3
gp4 <- gp3 + facet_wrap(~Sex, nrow=2)
gp4
# plot the number of passengers for different ports and Survived values
ggplot(titanic.train, aes(x = Embarked, fill = Survived)) +
geom_bar(position = "dodge", width = 0.5) +
labs(x="Place of embarkment", y="Number of passengers") +
theme_bw()
# examine the relation between Embarked and Survived, but with proportions
ggplot(data = titanic.train,
mapping = aes(x = Embarked, fill=Survived)) +
geom_bar(position = "fill", width = 0.5) +
labs(x="\nPlace of embarkment", y="Proportion of passengers\n") +
theme_minimal()
# examine the relation between Fare and Survived
ggplot(data = titanic.train,
mapping = aes(x = Fare, fill=Survived)) +
geom_density(alpha = 0.5) +
theme_minimal()
###
# Feature engineering
###
# add the Survived variable to the test set
titanic.test$Survived <- NA
# transform the Pclass variable into factor (in the test set)
titanic.test$Pclass <- factor(x = titanic.test$Pclass,
levels = c(1,2,3),
labels = c("1st", "2nd", "3rd"))
# transform the Sex variable into factor (in the test set)
titanic.test$Sex <- factor(titanic.test$Sex)
# merge train and test sets
titanic.all <- rbind(titanic.train, titanic.test)
###
## Creating the FamilySize variable
###
# examine the values of the SibSp variable
summary(titanic.all$SibSp)
table(titanic.all$SibSp)
# examine the values of the Parch variable
summary(titanic.all$Parch)
table(titanic.all$Parch)
# create a new variable FamilySize based on the SibSp and Parch values
titanic.all$FamilySize <- titanic.all$SibSp + titanic.all$Parch
summary(titanic.all$FamilySize)
# print the contingency table for the FamilySize
table(titanic.all$FamilySize)
# compute the proportion of FamilySize >= 3 in all passangers
sum(titanic.all$FamilySize>=3)/length(titanic.all$FamilySize)
# set the FamilySize to 3 to all observations where FamilySize > 3
titanic.all$FamilySize[titanic.all$FamilySize > 3] <- 3
# transform FamilySize into factor
titanic.all$FamilySize <- factor(titanic.all$FamilySize,
levels = 0:3,
labels = c(0:2, "3+"))
table(titanic.all$FamilySize)
# plot the FamilySize vs. Survived
ggplot(titanic.all[is.na(titanic.all$Survived) == FALSE,],
aes(x = FamilySize, fill = Survived)) +
geom_bar(position = "dodge", width = 0.5) +
theme_light()
#####
# kraj 8. i 9. casa
#####
#####
# DT primer
#####
## Napomena: ova skripta sadrzi originalno resenje jednog od studenata.
## Na nekoliko mesta u skripti su dodate napomene kako bi se ukazalo na neka malo bolja resenja
## ili na blagu izmenu nekih formulacija kako bi se iste ucinile jasnijim
# ucitavanje dataseta
data <- read.csv('bank.csv', stringsAsFactors = F)
str(data)
all(!is.na(data$Exited))
# pravljenje izlazne varijable
data$Stayed <- ifelse(data$Exited == 0, 'Yes', 'No')
data$Stayed <- as.factor(data$Stayed)
data$Exited <- NULL
any(is.na(data$Age))
# dopuna nedostajucih vrednosti Age varijable
unique(data$Age)
data$Age[data$Age == "-"] <- NA
data$Age <- as.numeric(data$Age)
shapiro.test(data$Age)
data$Age[is.na(data$Age)] <- median(data$Age, na.rm = T)
# pravljenje subseta koji cemo nadalje koristiti
dataSub <- subset(data, data$Age <= 87)
rm(data)
str(dataSub)
#provera nedostajucih vrednosti
apply(dataSub, 2, function(x) sum(is.na(x)))
apply(dataSub, 2, function(x) sum(x == "", na.rm = T))
apply(dataSub, 2, function(x) sum(x == "-", na.rm = T))
apply(dataSub, 2, function(x) sum(x == " ", na.rm = T))
# varijablu Surname necemo ukljuciti u model iz razloga sto je u pitanju karakter varijabla sa previse
# razlicitih vrednosti
length(unique(dataSub$Surname))
dataSub$Surname <- NULL
# varijabla Geography NA i - vrednosti tako da cemo ih zameniti dominantnom klasom i pretvoriti u faktorsku varijablu
length(unique(dataSub$Geography))
sort(table(dataSub$Geography))
dataSub$Geography[dataSub$Geography == "-" | is.na(dataSub$Geography)] <- 'France'
dataSub$Geography <- as.factor(dataSub$Geography)
str(dataSub)
# varijabla Gender ima samo dve moguce vrednosti tako da je pretvaramo u faktorsku
length(unique(dataSub$Gender))
dataSub$Gender <- as.factor(dataSub$Gender)
# varijabla Card Type slicno Gender ima 4 moguce vrednosti tako da je pretvaramo u faktorsku
length(unique(dataSub$Card.Type))
dataSub$Card.Type <- as.factor(dataSub$Card.Type)
# varijable HasCrCard i IsActiveMember imaju samo dve moguce vrednosti tako da cemo ih takodje pretvoriti u faktorske
dataSub$HasCrCard <- factor(dataSub$HasCrCard, levels = 0:1, labels = c('No', 'Yes'))
dataSub$IsActiveMember <- factor(dataSub$IsActiveMember, levels = 0:1, labels = c('No', 'Yes'))
library(ggplot2)
# na grafiku vidimo da se vrednosti za CreditScore ne razlikuju mnogo u pogledu izlazne varijable tako da ovu
# varijablu necemo koristiti za model
ggplot(dataSub, aes(x = CreditScore, fill = Stayed)) + geom_density(alpha = 0.55) + theme_minimal()
dataSub$CreditScore <- NULL
# na grafiku vidimo da se vrednosti za Age razlikuju u pogledu izlazne varijable tako da cemo ovu
# varijablu koristiti za model
ggplot(dataSub, aes(x = Age, fill = Stayed)) + geom_density(alpha = 0.55) + theme_minimal()
# na grafiku vidimo da se vrednosti za Tenure ne razlikuju mnogo u pogledu izlazne varijable tako da ovu
# varijablu necemo koristiti za model
ggplot(dataSub, aes(x = Tenure, fill = Stayed)) + geom_density(alpha = 0.55) + theme_minimal()
dataSub$Tenure <- NULL
# na grafiku vidimo da se vrednosti za Balace ne razlikuju mnogo, ali dovoljno u pogledu izlazne varijable tako da cemo ovu
# varijablu koristiti za model
ggplot(dataSub, aes(x = Balance, fill = Stayed)) + geom_density(alpha = 0.55) + theme_minimal()
# na grafiku vidimo da se vrednosti za NumOfProducts razlikuju u pogledu izlazne varijable tako da cemo ovu
# varijablu koristiti za model
ggplot(dataSub, aes(x = NumOfProducts, fill = Stayed)) + geom_density(alpha = 0.55) + theme_minimal()
# na grafiku vidimo da se vrednosti za EstimatedSalary ne razlikuju mnogo u pogledu izlazne varijable tako da ovu
# varijablu necemo koristiti za model
ggplot(dataSub, aes(x = EstimatedSalary, fill = Stayed)) + geom_density(alpha = 0.55) + theme_minimal()
dataSub$EstimatedSalary <- NULL
# na grafiku vidimo da se vrednosti za Satisfaction.Score ne razlikuju znacajno u pogledu izlazne varijable tako da necemo ovu
# varijablu koristiti za model
ggplot(dataSub, aes(x = Satisfaction.Score, fill = Stayed)) + geom_density(alpha = 0.55) + theme_minimal()
## Napomena: Bolja opcija za vizuelizaciju je koriscenje bar plot-a, s obzirom da omogucuje bolje sagledavanje
## distribucije izlazne varijable za razlicite vrednosti ulazne varijable:
## ggplot(dataSub, aes(x = Satisfaction.Score, fill = Stayed)) + geom_bar(position = "fill") + theme_minimal()
dataSub$Satisfaction.Score <- NULL
# na grafiku vidimo da se vrednosti za Point.Earned ne razlikuju mnogo u pogledu izlazne varijable tako da ovu
# varijablu necemo koristiti za model
ggplot(dataSub, aes(x = Point.Earned, fill = Stayed)) + geom_density(alpha = 0.55) + theme_minimal()
dataSub$Point.Earned <- NULL
# na grafiku vidimo da se vrednosti za Geography razlikuju u pogledu izlazne varijable tako da cemo ovu
# varijablu koristiti za model
ggplot(dataSub, aes(x = Geography, fill = Stayed)) + geom_bar(position = 'dodge', width = 0.5) + theme_minimal()
## Napomena: Bolja opcija je koriscenje vrednosti "fill" za position parametar jer se time dobija prikaz u
## relativnim vrednostima (proportions), a ne apsolutnim, pa se lakse moze sagledati znacajnost varijable
## Ne odnosi se samo na ovu varijablu, vec vazi generalno u slucaju koriscenja bar plota za proveru znacajnosti
## varijabli
# na grafiku vidimo da se vrednosti za Gender razlikuju u pogledu izlazne varijable tako da cemo ovu
# varijablu koristiti za model
ggplot(dataSub, aes(x = Gender, fill = Stayed)) + geom_bar(position = 'dodge', width = 0.5) + theme_minimal()
# na grafiku vidimo da se vrednosti za HasCrCard ne razlikuju u pogledu izlazne varijable tako da necemo ovu
# varijablu koristiti za model
ggplot(dataSub, aes(x = HasCrCard, fill = Stayed)) + geom_bar(position = 'fill', width = 0.5) + theme_minimal()
dataSub$HasCrCard <- NULL
# na grafiku vidimo da se vrednosti za IsActiveMember razlikuju u pogledu izlazne varijable tako da cemo ovu
# varijablu koristiti za model
ggplot(dataSub, aes(x = IsActiveMember, fill = Stayed)) + geom_bar(position = 'dodge', width = 0.5) + theme_minimal()
# na grafiku vidimo da se proporcije za Card.Type ne razlikuju mnogo u pogledu izlazne varijable tako da ovu
# varijablu necemo koristiti za model
ggplot(dataSub, aes(x = Card.Type, fill = Stayed)) + geom_bar(position = 'fill', width = 0.5) + theme_minimal()
dataSub$Card.Type <- NULL
library(caret)
# delimo dataset na trening(80% pocetnog dataseta) i test(20% pocetnog dataseta)
set.seed(1)
train_indices <- createDataPartition(dataSub$Stayed, p = 0.8, list = FALSE)
train_data <- dataSub[train_indices,]
test_data <- dataSub[-train_indices,]
library(rpart)
library(e1071)
# definisemo da cemo za odredjivanje cp vrednosti koristiti 10-fold cross-validation
# definisemo u kom rasponu vredsnoti trazimo najbolju vrednost cp parametra
tr_ctrl <- trainControl(method = 'cv', number = 10)
cp_grid <- expand.grid(.cp = seq(0.01, 0.05, 0.001))
# pokrecemo kros-validaciju, a kako je u pitanju probabilisticki proces moramo podesiti seed vrednost
set.seed(1)
cv <- train(x = train_data[, -7], y = train_data$Stayed, method = 'rpart', trControl = tr_ctrl, tuneGrid = cp_grid)
cv
# vidimo da nam je najbolja vrednost za cp parametar 0.023
best_cp <- cv$bestTune$cp
# kreiramo model sa nadjenom cp vrednoscu
set.seed(1)
model <- rpart(Stayed ~ ., data = train_data, control = rpart.control(cp = best_cp))
model
# pravimo predikcije
model_pred <- predict(model, newdata = test_data, type = 'class')
#kreiramo matricu konfuzije
cm <- table(actual = test_data$Stayed, predicted = model_pred)
cm
#iz evaluaciona matrice gledamo 4 vrednosti TP, TN, FP i FN:
# True-Positive(TP) - koliko klijenata smo predvideli da su ostali u banci i oni su zaista ostali, u nasem slucaju 749
# True-Negative(TN) - koliko klijenata smo predvideli da su napustili banku i oni su je zaista napustili, u nasem slucaju 91
# False-Positive(FP) - koliko klijenata smo predvideli da su ostali u banci a oni su je ipak napustili, u nasem slucaju 113
# False-Negative(FN) - koliko klijenata smo predvideli da su napustili banku, a oni su ipak ostali, u nasem slucaju 25
# definisemo funkciju za izracunavanje evaluacionih metrika
compute_eval_metrics <- function(cm){
TP <- cm[2,2]
TN <- cm[1,1]
FP <- cm[1,2]
FN <- cm[2,1]
a <- (TP + TN) / sum(cm)
p <- TP / (TP + FP)
r <- TP / (TP + FN)
f1 <- 2*p*r/(p+r)
c(accuracy = a, precision = p, recall = r, F1 = f1)
}
# evaluacione metrike
eval <- compute_eval_metrics(cm)
eval
# racunamo 4 evaluacione matrike za nas model:
#
# accuracy - koliko tacno predvidjamo izlaznu varijablu, u nasem slucaju od 978 klijenata mi smo tacno predvideli 807,
# odnosno, za 85.9% klijenata smo tacno predvideli da li ce ostati u banci ili ne
#
# precision - od svih klijenata koje smo predvideli da ce ostati banci, koji procenta klijenata je zaista ostao,
# kako nam je FP 126, takodje vidimo da nam je ova metrika visoka i iznosi 0.87; odnosno, od svih klijenata za koje smo
# predvideli da ce ostati u banci, 87% njih je stvarno ostalo
#
# recall - od svih klijenata koji su zaista ostali u banci, koji procenat klijenata smo mi predvideli da ce ostati u banci
# kako nam je FN samo 45, mozemo videti nam je i ova metrika visoka odnosno 0.968; odnosno, medju svim klijentima koji su
# zaista ostali u banci, mi smo tacno predvideli 96.8% njih
#
# F1 - predstavlja meru performanski modela za balansirane vrednosti precision i recall, i za njom postoji potreba jer su
# precision i recall u antagonistickom odnosu, povecanje jedne, druga se smanjuje i obrnuto
# kako su i precision i recall visoki, i vrednost ove metrike je visokih 0.92. U kontekstu zadatka, ova metrika pokazuje
# da smo predvidjanjima ne samo u najvecoj meri obuhvatili sve klijente koji su ostali u banci, vec i da nismo pravili
# mnogo gresaka predvidjajuci ostanak u banci onih koji su je zapravo napustili.
## Napomena: neke od prethodnih formulacija su blago izmenjene, bez promene znacenja, kako bi se ucinile jasnijim i/ili
## kompletnijim. VAZNO: date interpretacije tj nacin na koji su iskazane nije jedini korektan nacin, vec pokazuje samo
## jedan od nacina da se formulisu interpretacije metrika kako bi se odgovorilo na zahtev zadatka.
##
#####
# Kraj DT primera
#####
#####
# K Means primer
#####
## Napomena: ova skripta sadrzi originalno resenje jednog od studenata.
## Na nekoliko mesta u skripti su dodate napomene kako bi se ukazalo na neka malo bolja resenja
## ili na blagu izmenu nekih formulacija kako bi se iste ucinile jasnijim
# ucitavamo pocetni dataset
df <- read.csv("spotify.csv", stringsAsFactors = F)
# izbacujemo iz skupa podataka one obzervacije koje imaju vrednost atributa speechiness > 60
indeksi <- which(df$speechiness > 60)
df <- df[-indeksi,]
sum(df$speechiness > 60)
# vidimo da sada imamo 0 obzervacija koje imaju vrednost atributa speechiness > 60
# pre nego sto odredimo znacaj atributa, malo cemo da sredimo skup podataka
# kako kmeans algoritam radi sa numerickim podacima na ulazu, kolonu track_name cemo izbaciti iz daljeg
# razmatranja jer ne mozemo da je pretvorimo u numericku varijablu a da sa logicke strane donese neku vaznu informaciju
df$track_name <- NULL
table(df$mode)
# kako je raspodela promenljive mode slicna za obe vrednosti, a nemamo informaciju o tome kakav je odnos izmedju
# tih vrednosti kako bi ih pretvorili u numericku varijablu, izbacicemo i atribut mode
#Dodatna napomena: varijablu mode smo mogli i da prepoznamo kao ordinalnu pa da je na taj nacin posmatramo.
df$mode <- NULL
# proveravamo postojanje nedostajucih vrednosti
apply(df, 2, FUN = function(x) sum(is.na(x) | x == "" | x == " " | x == "-"))
# kolona in_spotify_charts se pretvara u numericku kolonu prvo pa ce se onda zameniti njene vrednosti odgovarajucom
# vrednoscu
df$in_spotify_charts <- as.numeric(df$in_spotify_charts)
summary(df$in_spotify_charts)
# da bi zamenili na vrednosti prvo proveravamo raspodelu promenljive
shapiro.test(df$in_spotify_charts)
# kako promenljiva dolazi iz raspodele cija distribucija nije normalna, na vrednosti cemo zamentiti medijanom
df$in_spotify_charts[is.na(df$in_spotify_charts)] <- median(df$in_spotify_charts , na.rm = T)
summary(df$in_spotify_charts)
# vidimo da vise nemamo na vrednosti
# isti postupak radimo i za promemnljivu in_apple_charts
shapiro.test(df$in_apple_charts)
# na vrednosti kod in_apple_charts kolone menjamo takodje sa medijanom
df$in_apple_charts[is.na(df$in_apple_charts)] <- median(df$in_apple_charts, na.rm = T)
summary(df$in_spotify_charts)
# izbacili smo na vrednosti iz kolone in_apple_charts
apply(df, 2, FUN = function(x) sum(is.na(x) | x == "" | x == " " | x == "-"))
# jos jednom proveravamo prisustvo na vrednosti i vidimo da nemamo vise nijednu
# sada cemo odraditi znacajnost varijabli kreiranjem matrice korelacije
library(ggplot2)
cm <- cor(df)
cm
library(corrplot)
corrplot.mixed(cm, tl.cex = 0.7, number.cex = 0.7)
# vidimo da atributi kao sto su streams i in_apple_playlist su visoko korelisani sa sa nekim drugim atributima
# tako da cemo njih izbaciti iz daljeg razmatranja
df$streams <- NULL
df$in_apple_playlists <- NULL
# ostale atribute za sada zadrzavamo u daljoj analizi
# pre daljem rada, potrebno je proveriti postojanje autlajera (ekstremnih vrdnosti) u skupu podataka
# zbog nacina rada kmeans algoritma, autlajeri mogu znacajno uticati na dobijene klastere, te je njih potrebno
# obraditi na odgovoarajuci nacin
# proveravamo postojanje autlajera
apply(df, 2, FUN = function(x) length(boxplot.stats(x)$out))
# za kolone koje imaju ukupan broj autlajera veci od 10% od ukupnog broja obzervacija smo se odlucili
# da izbacimo iz daljeg razmatranja jer primenom tehnika za otklanjanje autlajera cemo znacajno uticati
# na raspodelu tih atributa te smatramo da takve kolone treba izbaciti iz daljeg razmatranja
df$released_year <- NULL
df$in_spotify_playlists <- NULL
df$speechiness <- NULL
# za ostale kolone, koje imaju manje od 10% autlajera, cemo primeniti tehniku Winsorize
library(DescTools)
# za svaki od preostalih atributa cemo iscrtati boxplot, videti sa koje strane se nalaze autlajeri
# i zameniti ih odgovarajucim vrednostima (najceste 5 i 95 percentilom u zavisnosti od vrednosti autlajera)
boxplot(df$in_spotify_charts)
tmp <- Winsorize(df$in_spotify_charts, probs = c(0, 0.95))
boxplot(tmp)
# pomericemo jos nize granicu jer idalje postoje autlajeri
tmp <- Winsorize(df$in_spotify_charts, probs = c(0, 0.92))
boxplot(tmp)
# sada cemo ovakve vrednosti dodeliti koloni in_spotify_charts i ovaj postupak primeniti za sve ostale kolone
df$in_spotify_charts <- tmp
boxplot(df$in_spotify_charts)
boxplot(df$in_apple_charts)
tmp <- Winsorize(df$in_apple_charts, probs = c(0, 0.95))
boxplot(tmp)
df$in_apple_charts <- tmp
boxplot(df$in_apple_charts)
boxplot(df$bpm)
tmp <- Winsorize(df$in_apple_charts, probs = c(0, 0.95))
boxplot(tmp)
df$bpm <- tmp
boxplot(df$bpm)
boxplot(df$energy)
tmp <- Winsorize(df$energy , probs = c(0.05, 1))
boxplot(tmp)
df$energy <- tmp
boxplot(df$energy)
boxplot(df$instrumentalness)
table(df$instrumentalness)
tmp <- Winsorize(df$instrumentalness , probs = c(0, 0.95))
boxplot(tmp)
# uzimamo jos manju vrednost jer idalje ima autlajera
tmp <- Winsorize(df$instrumentalness , probs = c(0, 0.9))
boxplot(tmp)
df$instrumentalness <- tmp
boxplot(df$instrumentalness)
table(df$instrumentalness)
# kod atributa instrumentalness vidimo jedno vrlo zanimljivu stvar a to je da nakon uklanjanja autlajera
# sve obezrvacije imaju vrednost 0 za dati atribut
# to znaci da taj atribut ne nosi nikakvu informaciju i zato cemo ga izbaciti iz daljeg razmatranja
df$instrumentalness <- NULL
boxplot(df$liveness)
tmp <- Winsorize(df$liveness , probs = c(0, 0.95))
boxplot(tmp)
df$liveness <- tmp
boxplot(df$liveness)
apply(df, 2, FUN = function(x) length(boxplot.stats(x)$out))
# ponovo izvrsavamo pocetnu funkciju i vidimo da vise nemamo autlajere
summary(df)
# sada cemo odraditi proces normalizacije kako bi sveli sve vrednosti na isti opseg
# a u cilju da zbog razlicitih velicina atributa neki atributi ne bi postali znacajniji samo
# zato sto imaju veci opseg
normalizuj <- function(x) {
if(all(is.na(x) | x == 0)){
return (x)
}
if (max(x) == min(x)) {
return (x)
}
( (x - min(x)) / (max(x) - min(x)) )
}
df2 <- as.data.frame(
apply(df, 2, FUN = function(x) normalizuj(x))
)
summary(df2)
# sada ovako normalizovane vrednosti (vrednosti svih atributa svedeni na opseg od 0 do 1
# gde je 0 minimalna a 1 maksimalna vrednost) vracamo u pocetni skup podataka
df <- df2
summary(df)
# sada sprovodimo elbow metodu kako bi nasli najbolju vrednost za broj klastera
# prvo kreiramo pomocni dataframe u koje cemo smestati razlicite razultate za razlicite vrednosti klastera
df.tmp <- data.frame()
# odradicemo proveru za opseg od 2 do 10 klastera i videti koja vrednost od tih daje najbolje rezultate
for (k in 2:10) {
model <- kmeans(x = df,
centers = k,
iter.max = 20,
nstart = 1000)
df.tmp <- rbind(df.tmp,
c(k, model$tot.withinss, model$betweenss / model$totss))
}
df.tmp
# radi bolje preglednosti promenicemo nazive kolona
colnames(df.tmp)<- c("k", "tot.withinss", "ratio")
df.tmp
# iscrtacemo grafik da zagledamo ponasanje modela za razlicite vrednosti parametra k
library(ggplot2)
ggplot(data = df.tmp, aes (x = k, y = tot.withinss)) + geom_line() + geom_point()
# sagledavajuci grafik kao i vrednosti koje se nalaze u df.tmp setu vidimo da su nam ponajbolje opcije za broj klastera 2 ili 3,
# s tim da za nijansu su bolje vrednosti za broj klastera = 2 te cemo tu vrednost uzeti za kreiranje model
# izvrsavamo klasterizaciju za broj klastera = 2
model <- kmeans(x = df,
centers = 2,
iter.max = 20,
nstart = 1000)
model
# vidimo da sve obzervacije se nalaze ili u klasteru 1 ili u klasteru 2
source("Utility.R")
stats <- summary.stats(df, model$cluster, 2)
stats
# ucitana je skripta Utility i iz nje pozvana funkcija summary.stats koja nam daje odredjene statistike za model koji smo
# kreirali a na osnovu kojih cemo izvrsiti interpretaciju dobijenih klastera
# 1. Broj pesama po klasteru
# U klasteru 1 imamo 378 pesama, dok u klasteru 2 imamo 574, odnosno klaster 2 sadrzi vise pesama
# 2. Centri klastera
# U klasteru 1 se uglavnom nalaze pesme koje su visoko kotirane na applerang listi (in_apple_charts) i imaju veci bpm,
# odnosno brze su pesme. Sa druge strane, klaster 2 uglavnom sadrzi pesme nisko kotirane na applerang listi, koje
# imaju nesto sporiji ritam (bpm) i takodje nisu dobro kotirane na spotify rang listi
# 3. Disperzija od centara
# Za klaster 2 mozemo da kazemo da, iako je veci po svom obimu, je homogeniji jer je disperzija od centra klastera
# za skoro sve atribute manja u odnosu na model 1 a noricto atribute in_apple_charts i bpm koji su se pokazali kao
# znacajni za ovako kreiran model. Kod klastera 1, kojeg karakterise manja homogenost u pogledu na klaster 2,
# vidimo da za atribut in_spotify_charts imamo najvecu disperziju obzervacija u odnosu na centar klastera, sto bi znacilo
# da se u tom klasteru mogu naci i pesme koje imaju visoku poziciju na sporitfy rang listi ali i nisku poziciju
## Napomena: neke od prethodnih formulacija su blago izmenjene, bez promene znacenja, kako bi se ucinile jasnijim i/ili
## kompletnijim. VAZNO: date interpretacije tj nacin na koji su iskazane nije jedini korektan nacin, vec pokazuje samo
## jedan od nacina da se formulisu interpretacije metrika kako bi se odgovorilo na zahtev zadatka.
##
#####
# Kraj K Means primera
#####
#####
# KNN primer
#####
## Napomena: ova skripta sadrzi originalno resenje jednog od studenata.
## Na nekoliko mesta u skripti su dodate napomene kako bi se ukazalo na neka malo bolja resenja
## ili na blagu izmenu nekih formulacija kako bi se iste ucinile jasnijim
#ucitavamo dataset
data <- read.csv("spotify.csv", stringsAsFactors = F)
str(data)
#proveravamo da li varijabla streams ima nedostajuce vrednosti jer na osnovu nje kreiramo izlaznu promenljivu
sum(is.na(data$streams))
sum(data$streams == "", na.rm = T)
sum(data$streams == " ", na.rm = T)
sum(data$streams == "-", na.rm = T)
#nema nedostajucih vrednosti, racunamo treci kvartil
treciKvartil <- quantile(data$streams, 0.75)
treciKvartil
data$High_Streams <- ifelse(data$streams > treciKvartil, yes = "yes", no = "no")
data$High_Streams <- as.factor(data$High_Streams)
#odmah izbacujemo varijablu streams jer smo je koristili za kreiranje izlazne varijable, pa nam nece biti vise potrebna za model
data$streams <- NULL
str(data)
#IZBOR VARIJABLI - KNN koristi samo numericke podatke!
#proveravamo koje cemo varijable ukljuciti u model, proveravamo koliko imaju razlicitih vrednosti
length(unique(data$track_name))
#svaka pesma ima jedinstveno ime i nece uticati na to da li pesma ima veliki broj pustanja pesme
data$track_name <- NULL
length(unique(data$in_spotify_charts))
table(data$in_spotify_charts)
#mozemo je prebaciti u numericku varijablu (jer je to sustinski numericka vrednost), necemo je izbacivati
length(unique(data$mode))
#takodje se moze pretvoriti u numercku s obzirom da je ordinalna
#proveravamo da li promenljive imaju nedostajuce vrednosti
apply(data, MARGIN = 2, FUN = function(x) sum(is.na(x)))
apply(data, MARGIN = 2, FUN = function(x) sum(x == "", na.rm = T))
apply(data, MARGIN = 2, FUN = function(x) sum(x == " ", na.rm = T))
apply(data, MARGIN = 2, FUN = function(x) sum(x == "-", na.rm = T))
#in_spotify_charts i in_apple_charts imaju 1 i 2 NA vrednosti respektivno
#in_spotify_charts ima "-" vrednost
data$in_spotify_charts[data$in_spotify_charts == "-"] <- NA
data$in_spotify_charts <- as.numeric(data$in_spotify_charts, na.rm = T)
#proveravamo da li ima normalnu raspodelu kako bismo znali sa kojim vrednostima menjamo NA vrednosti i za in_spotify_charts i za in_apple_charts
apply(X = data[,c(3,5)], 2, FUN = function(x) shapiro.test(x))
#obe imaju vrednost p < 0.05 sto znaci da nemaju normalnu raspodelu i NA vrednosti menjamo medijanom
medijanaSpotify <- median(data$in_spotify_charts, na.rm = T)
medijanaSpotify
medijanaApple <- median(data$in_apple_charts, na.rm = T)
medijanaApple
#dobili smo da su medijane 3 i 38 za spotify i apple respektivno
data$in_spotify_charts[is.na(data$in_spotify_charts)] <- medijanaSpotify
data$in_apple_charts[is.na(data$in_apple_charts)] <- medijanaApple
#nema vise nedostajucih vrednosti
str(data)
data$mode <- as.factor(data$mode)
#radimo plotovanje kako bismo jos jednom proverili koje cemo varijable ukljuciti u nas model
library(ggplot2)
ggplot(data, mapping = aes(x = released_year, fill = High_Streams)) + geom_density(alpha = 0.5)
#sto je skorija godina izdanja to je manja verovatnoca da pesma ima visok broj pustanja
ggplot(data, mapping = aes(x = in_spotify_charts, fill = High_Streams)) + geom_density(alpha = 0.5)
#sto je veca pozicija pesme na spotify rang listi to je veca verovatnoca da ima veci broj pustanja, ta verovatnoca opada sa smanjenjem pozicije
ggplot(data, mapping = aes(x = in_spotify_playlists, fill = High_Streams)) + geom_density(alpha = 0.5)
#ukoliko se pesma nalazi u malom broju plejlista na spotify platformi to je manja verovatnoca da ima veci broj pustanja
ggplot(data, mapping = aes(x = in_apple_playlists, fill = High_Streams)) + geom_density(alpha = 0.5)
#ukoliko se pesma nalazi u malom broju plejlista na apple platformi to je manja verovatnoca da ima veci broj pustanja
ggplot(data, mapping = aes(x = in_apple_charts, fill = High_Streams)) + geom_density(alpha = 0.5)
#ukoliko pesma nema poziciju na apple rang listi to je veca verovatnoca da nema veliki broj pustanja
ggplot(data, mapping = aes(x = bpm, fill = High_Streams)) + geom_density(alpha = 0.5)
#nema neke prevelike razlike u broju pustanja pesme kada je u pitanju mera brzine pesme, nije nam relevantna za model
data$bpm <- NULL
ggplot(data, mapping = aes(x = liveness, fill = High_Streams)) + geom_density(alpha = 0.5)
#isto vazi i za nivo prisustva elemenata uzivo izvodjenja, ni ona nije relevantna za nas model
data$liveness <- NULL
ggplot(data, mapping = aes(x = energy, fill = High_Streams)) + geom_density(alpha = 0.5)
#ni energetski nivo pesme ne utice na broj pustanja pesme, mozemo je izbaciti
data$energy <- NULL
ggplot(data, mapping = aes(x = instrumentalness, fill = High_Streams)) + geom_density(alpha = 0.5)
#ni instrumentalness nam nije relevantna i nju izbacujemo
data$instrumentalness <- NULL
ggplot(data, mapping = aes(x = speechiness, fill = High_Streams)) + geom_density(alpha = 0.5)
#ni ova varijabla nam ne utice na model tako da cemo je izbaciti
data$speechiness <- NULL
ggplot(data, mapping = aes(x = mode, fill = High_Streams)) + geom_bar(position = 'dodge')
#veca je verovatnoca da pesma ima veci broj pustanja ako je karakter pesme major
str(data)
data$mode <- as.numeric(data$mode)
#prebacili smo varijablu mode u numericku jer knn model radi samo sa numerickim varijablama
#gotovo je sredjivanje podataka
#proveravamo da li neke variujable imaju autlajere
apply(X = data[,c(1,2,3,4,5,6)], 2, FUN = function(x) length(boxplot.stats(x)$out))
#vidimo da sve varijable imaju autlajere, potrebno je da ih standardizujemo
#proveravamo da li imaju normalnu raspodelu kako bismo znali na koji nacin cemo izvrsiti standardizaciju
apply(X = data[,c(1,2,3,4,5,6)], 2, FUN = function(x) shapiro.test(x))
#sve varijable imaju p < 0.05 sto znaci da nemaju normalnu raspodelu
data.std <- apply(X = data[,c(1,2,3,4,5,6)], 2, FUN = function(x) scale(x, center = median(x), scale = IQR(x)))
data.std <- as.data.frame(data.std)
#kreirali smo standardizovani dataframe
#za scale smo stavili IQR - interquartile range sto predstavlja razliku izmedju 3. i 1. kvartila i kao centar smo stavili medijanu
#ubacujemo u standardizovani dataset preostalu varijablu mode i izlaznu varijablu
data.std$mode <- data$mode
data.std$High_Streams <- data$High_Streams
#sa tako standardizovanim dataset-om kreiramo nas model
#kreiramo train i test setove
library(caret)
set.seed(1010) # allows for repeating the randomization process exactly
indexes <- createDataPartition(data.std$High_Streams, p = 0.8, list = FALSE)
train.data <- data.std[indexes, ]
test.data <- data.std[-indexes, ]
#radimo krosvalidaciju sa 10 iteracija
library(e1071)
library(caret)
numFolds = trainControl(method = "cv", number = 10)
kGrid = expand.grid(.k = seq(from = 3, to = 25, by = 2))
#stavili smo neparne brojeve od 3 do 25 jer su nam oni potrebni kako bi jedna klasa bila dominantnija
knn.cv <- train(
x = train.data[, -7],
y = train.data$High_Streams,
method = "knn", # use rpart() to build multiple
trControl = numFolds, # <folds> from above
tuneGrid = kGrid)
knn.cv
#dobili smo da nam je optimalna vrednost za parametar k = 7;
#mozemo je izuci direktno odatle
kValue <- knn.cv$bestTune$k
kValue
library(class)
knn.pred <- knn(train = train.data[,-7], # training data without the output (class) variable
test = test.data[,-7], # test data without the output (class) variable
cl = train.data$High_Streams, # output (class) variable is specified here
k = kValue)
#kreirali smo model knn i sada mozemo kreirati matricu konfuzije
knn.cm <- table(true = test.data$High_Streams, predicted = knn.pred)
knn.cm
# predicted
# true no yes
# no 135 8
# yes 12 35
#na glavnoj dijagonali matirce se nalaze tacne predikcije, a na sporednoj pogresne predikcije
#od ukupno 190 observacija tacno smo predvideli 170 sto znaci da ce nam accuracy biti veoma visok
#od ukupno 143 pesme koje stvarno nece imati veliki broj pustanja, tacno smo predvideli 135
#a od ukupno 47 pesama koje ce imati veliki broj pustanja, tacno smo predvideli 35
#TP - true positive - predstavlja sve pesme koje ce stvarno imati veliki broj pustanja a i mi smo predvideli da hoce
#TN - true negative - predstavlja sve pesme koje nece imati veliki broj pustanja i mi smo predvideli da nece
#FP - false positive - predstavlja sve pesme za koje smo rekli dace imati veliki broj pustanja a zapravo nemaju
#FN - false negative - predstavlja sve pesme za koje smo rekli da nece imati veliki broj pustanja a zapravo imaju
#kreiramo funkciju koja ce nam sluziti za evaluaciju modela
getEvaluationMetrics <- function(cm) {
TP <- cm[2,2]
TN <- cm[1,1]
FP <- cm[1,2]
FN <- cm[2,1]
accuracy <- sum(diag(cm))/sum(cm)
precision <- TP/(TP+FP)
recall <- TP/(TP+FN)
F1 <- (2*precision*recall)/(precision + recall)
c(Accuracy = accuracy,
Precision = precision,
Recall = recall,
F1 = F1)
}
knn.eval <- getEvaluationMetrics(knn.cm)
knn.eval
#accuracy - predstavlja udeo tacno predvidjenih u odnosu na sve observacije.
#nas model je od ukupno 190 observacija tacno predvideo 170 i zbog toga nam je accuracy 89.47%
#za 170 pesama je tacno predvideo da li ce imati veliki broj pustanja pesme
#precision - predstavlja udeo observacija za koje smo predvideli da su pozitivne i stvarno jesu
#za 43 pesme smo predvideli da ce imati veliki broj pustanja, od toga 35 pesama stvarno ima veliki broj pustanja
#i zbog toga nam je precision 81.39%
#recall - predstavlja udeo observacija koje su stvarno pozitivne a mi smo predvideli da jesu
#od ukupno 47 stvarno pozitivnih observacija mi smo tacno predvideli 35
#za 35 pesama smo rekli da ce imati veliki broj pustanja od ukupno 47 koje stvarno imaju veliki broj pustanja
#i zbog toga nam je recall 74.46%
#F1 - predstavlja meru performanski modela za balansirane vrednosti precision i recall, i za njom postoji potreba jer su
# precision i recall u antagonistickom odnosu, povecanje jedne, druga se smanjuje i obrnuto. Kada su i recall i precision visoke,
# visok je i F1.
# U kontekstu ovog zadatka, moze se reci da nas model dobro balansira između sposobnosti da prepozna većinu stvarno popularnih
#pesama i toga da ne označava previše pesama kao veoma popularne kada nisu.
## Napomena: neke od prethodnih formulacija su blago izmenjene, bez promene znacenja, kako bi se ucinile jasnijim i/ili
## kompletnijim. VAZNO: date interpretacije tj nacin na koji su iskazane nije jedini korektan nacin, vec pokazuje samo
## jedan od nacina da se formulisu interpretacije metrika kako bi se odgovorilo na zahtev zadatka.
#####
# Kraj KNN primera
#####
#####
# LR primer
#####
#Učitavanje dataseta
data <- read.csv("fastfood.csv", stringsAsFactors = FALSE)
summary(data)
str(data)
#KREIRANJE PODSKUPA PODATAKA koji sadrži proizvode čije ukupne masti su <= 125
data.subset <- subset(data, data$total_fat <= 125)
#PROVERA NEDOSTAJUĆIH VREDNOSTI (NA,“-”, ” ”, ””)
apply(data.subset, MARGIN = 2, FUN = function(x) length(which(is.na(x))))
#Atributi cal_fat(1), fiber(12), protein(2) i calcium(209) imaju nedostajuće(NA) vrednosti.
apply(data.subset, MARGIN = 2, FUN = function(x) length(which(x == "-" | x == " " | x == "")))
#Atribut restaurant(3) ima nedostajuće vrednosti.
#ZAMENA NEDOSTAJUĆIH VREDNOSTI
column.with.na<-c(4,11,13,14)
apply(data.subset[,column.with.na],MARGIN = 2, FUN = function(x) shapiro.test(x))
#Svi atributi koji imaju NA vrednosti nemaju normalnu raspodelu.
#S obzirom da atributi nemaju normalnu raspodelu (p-value<0.05), menjamo nedostajuće vrednosti medijanom.
cal_fat_median <- median(data.subset$cal_fat, TRUE)
data.subset$cal_fat[is.na(data.subset$cal_fat)] <- cal_fat_median
fiber_median <- median(data.subset$fiber, TRUE)
data.subset$fiber[is.na(data.subset$fiber)] <- fiber_median
protein_median <- median(data.subset$protein, TRUE)
data.subset$protein[is.na(data.subset$protein)] <- protein_median
#Atribut calcium ima dosta nedostajućih vrednosti, pa neće biti ni zamenjene.
#Atribut calcium se izbacuje iz dataseta, zbog prevelikog broja nedostajućih vrednosti.
data.subset$calcium <- NULL
#Zamena - i " " sa NA
data.subset$restaurant[data.subset$restaurant == "-" | data.subset$restaurant == " " | data.subset$restaurant == ""] <- NA
table(data.subset$restaurant)
#S obzirom da je restaurant character varijabla, gledamo koja vrednost je najzastupljenija među podacima.
#Za atribut restaurant najučestalija je vrednost "Taco Bell", pa sa njom manja NA vrednosti.
data.subset$restaurant[is.na(data.subset$restaurant)] <- "Taco Bell"
#IZBOR ATRIBUTA ZA MODEL LINEARNE REGRESIJE
#Model LR podržava samo numeričke atribute.
#Naziv restorana (atribut restaurant) i naziv proizvoda (atribut item) se neće koristiti za kreiranje modela,
# jer su one character varijable koje se ne mogu pretvoriti u numeričke.
data.subset$restaurant <- NULL
data.subset$item <- NULL
#Značajnost ostalih atributa se ispituje pomoću korelacione matrice
cor.mat <- cor(data.subset)
library(corrplot)
corrplot.mixed(cor.mat, tl.cex = 0.75, number.cex = 0.75)
#Na korelacionoj matrici se može primetiti da su za predviđaje atributa protein najznačajni(najveća korelisanost) atributi:
# 1)cholesterol - 0.87 (jedinično povećanje atributa cholesterol dovešće do povećanja atributa protein za 0.87)
# 2)calories - 0.8 (jedinično povećanje atributa calories dovešće do povećanja atributa protein za 0.8)
# 3)sodium - 0.76 (jedinično povećanje atributa sodium dovešće do povećanja atributa protein za 0.76)
#Zaključak: Za predviđanje atributa protein (prvi model LR) koristiće se atributi cholesterol,calories i sodium.
#KREIRANJE MODELA LR
#Podela dataseta na train i test
#Train set sadrzi 80% opservacija,a test ostalih 20%
library(caret)
set.seed(123)
i <- createDataPartition(data.subset$protein, p = 0.8, list = FALSE)
train <- data.subset[i, ]
test <- data.subset[-i, ]
#Kreiranje modela LR
lm1 <- lm(protein ~ cholesterol + calories + sodium, data = train)
summary(lm1)
comment("-Residuali predstavljaju razliku između stvarnih vrednosti zavisne varijable(protein) i
predviđenih vrednosti koje model generiše. Vrednost reziduala bi trebala da bude jednaka nuli.
Medijana reziduala je -0.702. To znači da je polovina naših podataka imala rezidualne vrednosti
manje od -0.702, a druga polovina veće.
-Intercept je 2.4256456, što znači da kada su svi ostali atributi 0, očekujemo da je prosečna
vrednost atributa protein 2.4256456.
-Koeficijenti za varijable cholesterol i sodium su statistički značajni jer su p-vrednosti daleko manje od 0.05,
dok koeficijent za calories nije statistički značajan.
Dakle cholesterol i sodium su relevantne za predikciju nivoa proteina u modelu.
-Koeficijenti: Kada se ostale varijable drže konstantnim:
Jedinično povećanje cholesterol za 1 dovešće do povećanja protein za 0.182 (ovo znači da postoji pozitivna linearna veza između nivoa holesterola i proteina, tj. što je veći nivo holesterola, veći je i nivo proteina).
Jedinično povećanje calories za 1 dovešće do povećanja protein za 0.0034,
međutim s obzirom da calories nije statistički značajan, ne možemo da kažemo da će značajano uticati na nivo proteina).
Jedinično povećanje sodium za 1 dovešće do povećanja protein za 0.0083 (postoji pozitivna linearna veza između nivoa sodium-a i proteina).
-Standardna greška reziduala je 6.492, što nam daje informaciju o prosečnoj udaljenosti između stvarnih vrednosti i predviđenih vrednosti.
-Koeficijent determinacije (R-kvadrat):Model objašnjava oko 83.52% varijabiliteta protein varijable.
-F-statistika: Nulta hipoteza F-statistike glasi da svi koeficijenti nezavisnih varijabli su nula,
što znači da nezavisne varijable ne utiču na zavisnu varijablu.
S obzirom da je p-vrednost manja od 0.05, može se zaključiti da model ima značajnu prediktivnu moć i da ga treba razmatrati.")
lm1.pred <- predict(lm1, test)
lm1.pred
# DIJAGNOSTIČKI GRAFIKONI
par(mfrow = c(2,2))
plot(lm1)
par(mfrow = c(1,1))
comment("-Reziduals vs Fitted- prikazuje raspodelu reziduala (razlike između stvarnih i predviđenih vrednosti) u odnosu na predviđene vrednosti (da li je zadovoljenja pretpostavka lineranosti).
Residuali bi trebali da budu raspoređeni oko horizontalne linije (njihova vrednost treba da bude 0).
Sa grafikona može da se zaključi da postoje manja odstupanja, tj. nelineranost u podacima, što bi moglo dovesti do neadekvatnosti modela.
-Normal Q-Q- govori o normalnoj raspodeli reziduala. Tačke bi trebale da se prate duž dijagonale.
S obzirom da tačke odstupaju od dijagonale na krajevima, ne možemo se u potpunosti osloniti na model.
-Scale-Location- prikazuje da li varijabilitet(greške,odstupanja) reziduala variraju na ujednačen način.
Crvena linija treba da bude horizontalna. Na grafikonu varijansa reziduala nije konstantna, što može dovesti do neadekvatnosti modela.
-Residuals vs. Leverage-prikazuje ekstremno visoke i/ili ekstremno niske vrednosti.
Većina tačaka trebala bi biti unutar granica, a ekstremne vrednosti bi trebale biti retke (nije dobro ukoliko postoje).
Na osnovu Kukova distance, možemo zaključiti da model sadrži outlayer.")
#MULTIKOLINERANOST
library(car)
vif(lm1)
sqrt(vif(lm1))
#S obzirom da je za atribut calories sqrt(vif) vrednost >2, može se zaključiti da postoji multikorelisanost, koja može dovesti do nepouzdanih procena.
#Novi model za predikciju proteina treba kreirati sa atributima cholesterol i sodium.
lm2<-lm(protein ~ cholesterol+ sodium, data = train)
sqrt(vif(lm2))
#Multikolienearnost ne postoji više.
summary(lm2)
#####
# Kraj LR primera
#####
#####
# NB primer
#####
# ucitavanje dataseta
data <- read.csv("apples.csv", stringsAsFactors = F)
# upoznajemo se sa strukturom dataseta
str(data)
# provera da li varijabla od koje cemo kreirati izlaznu sadrzi NA vrednosti
all(!is.na(data$Quality))
# kreiranje izlazne varijable
data$IsGood <- ifelse(data$Quality == "good", "yes", "no")
data$IsGood <- as.factor(data$IsGood)
data$Quality <- NULL
# izbacujemo varijable koje nema smisla iskoristiti za predvidjanje kvaliteta jabuke
data$A_id <- NULL
data$Code <- NULL
# Varijabla Continent ima 6 mogucih vrednosti, tako da je pretvaramo u faktorsku
length(unique(data$Continent))
data$Continent <- as.factor(data$Continent)
# naredne varijable su u osnovi numericki podaci tako da ih pretvaramo
data$Weight <- as.numeric(data$Weight)
data$Sweetness <- as.numeric(data$Sweetness)
data$Acidity <- as.numeric(data$Acidity)
# provera nedostajucih vrednosti
apply(data, MARGIN = 2, FUN = function(x) sum(is.na(x)))
# uocavamo nedostajuce vrednosti kod varijabli Weight, Sweetness, Acidity
apply(data, MARGIN = 2, FUN = function(x) sum(x == "", na.rm=T))
apply(data, MARGIN = 2, FUN = function(x) sum(x == "-", na.rm=T))
apply(data, MARGIN = 2, FUN = function(x) sum(x == " ", na.rm=T))
# nema nedostajucih vrednosti u nekom drugom obliku
# nije neophodno proveravati da li varijable imaju normalnu raspodelu jer ako imaju svakako je medijana jednaka mean-u
# tako da sve NA vrednosti menjamo medijanom
medianW<-median(data$Weight, na.rm=T)
data$Weight[is.na(data$Weight)] <- medianW
medianS<-median(data$Sweetness, na.rm=T)
data$Sweetness[is.na(data$Sweetness)] <- medianS
medianA <- median(data$Acidity, na.rm=T)
data$Acidity[is.na(data$Acidity)] <- medianA
# vise nemamo nedostajuce vrednosti
apply(data, MARGIN = 2, FUN = function(x) sum(is.na(x)))
# kako bismo izvrsili odabir atributa vizualizovacemo podatke
library(ggplot2)
# na grafiku vidimo male razlike u pogledu izlazne varijable u odnosu na Continent
# tako da cemo ovu varijablu koristiti za model
ggplot(data, aes(x = Continent, fill=IsGood)) +
geom_bar(position = "dodge", width = 0.4)
# na grafiku vidimo da se vrednosti za Size razlikuju u pogledu izlazne varijable
# tako da cemo ovu varijablu koristiti za model
ggplot(data, aes(x = Size, fill=IsGood)) +
geom_density(alpha = 0.5)
# na grafiku vidimo da se vrednosti za Weight razlikuju u pogledu izlazne varijable
# tako da cemo ovu varijablu koristiti za model
ggplot(data, aes(x = Weight, fill=IsGood)) +
geom_density(alpha = 0.5)
# na grafiku vidimo da se vrednosti za Sweetness razlikuju u pogledu izlazne varijable
# tako da cemo ovu varijablu koristiti za model
ggplot(data, aes(x = Sweetness, fill=IsGood)) +
geom_density(alpha = 0.5)
# na grafiku vidimo da se vrednosti za Crunchiness razlikuju u pogledu izlazne varijable
# tako da cemo ovu varijablu koristiti za mode
ggplot(data, aes(x = Crunchiness, fill=IsGood)) +
geom_density(alpha = 0.5)
# na grafiku vidimo da se vrednosti za Juiciness razlikuju u pogledu izlazne varijable
# tako da cemo ovu varijablu koristiti za mode
ggplot(data, aes(x = Juiciness, fill=IsGood)) +
geom_density(alpha = 0.5)
# na grafiku vidimo da se vrednosti za Ripeness razlikuju u pogledu izlazne varijable
# tako da cemo ovu varijablu koristiti za mode
ggplot(data, aes(x = Ripeness, fill=IsGood)) +
geom_density(alpha = 0.5)
# na grafiku vidimo da se vrednosti za Acidity ne razlikuju puno u pogledu izlazne varijable
# tako da ovu varijablu necemo koristiti za kreiranje naseg modela
ggplot(data, aes(x = Acidity, fill=IsGood)) +
geom_density(alpha = 0.5)
data$Acidity <- NULL
# obzirom na to da algoritam Naive Bayes radi sa faktorskim i numerickim varijablama ako one imaju normalnu raspodelu
# a ako nemaju potrebno je da se izvrsi diskretizacija
# moramo prvo da izvrsimo proveru normalnosti raspodele za numericke varijable
num.vars <- c(2:7)
apply(data[,num.vars],
MARGIN = 2,
FUN = function(x) shapiro.test(x))
# Na osnovu shapiro testa zakljucujemo da:
# Varijable Size, Juiciness, Ripeness imaju normalnu raspodelu
# Varijable Weight, Sweetness, Crunchiness nemaju normalnu raspodelu
# Varijable koje imaju normalnu raspodelu mozemo da iskoristimo za kreiranje modela
# Varijable koje nemaju normalnu raspodelu moramo prvo da diskretizujemo
to_discretize <- c(3,4,5)
# pokusacemo da podelimo ukupan opseg vrednosti svake varijable na 5 intervala jednake ucestalosti (kvantilna diskretizacija)
library(bnlearn)
discretized <- discretize(data[,to_discretize],
method = 'quantile',
breaks = c(5,5,5))
# vidimo da je diskretizacija uspesno izvrsena
summary(discretized)
data.new <- cbind(data[,c(1,2)], discretized, data[,c(6,7,8)])
# delimo dataset na trening (80% pocetnog dataseta) i test (20% pocetnog dataseta)
library(caret)
set.seed(10)
train_indices <- createDataPartition(data.new$IsGood, p = 0.8, list = FALSE)
train.data <- data.new[train_indices,]
test.data <- data.new[-train_indices,]
library(e1071)
# kreiramo model
nb1 <- naiveBayes(IsGood ~ ., data = train.data)
# vrsimo predvidjanje
nb1.pred <- predict(nb1, newdata = test.data, type = "class")
#kreiramo matricu konfuzije
nb1.cm <- table(true = test.data$IsGood, predicted = nb1.pred)
nb1.cm
# interpretacija matrice konfuzije
# True negative (TN) - za 207 jabuka smo predvideli da nece biti dobrog kvaliteta i one zaista nisu dobrog kvaliteta
# True positive (TP) - za 232 jabuke smo predvideli da ce biti dobrog kvaliteta i one zaista jesu dobrog kvaliteta
# False positive (FP) - za 92 jabuke smo predvideli da ce biti dobrog kvaliteta a one zapravo nisu dobrog kvaliteta
# False negative (FN) - za 68 jabuka smo predvideli da nece biti biti dobrog kvaliteta a one zapravo jesu dobrog kvaliteta
# definisemo funkciju za izracunavanje evalucionih metrika
compute.eval.metrics <- function(cmatrix) {
TP <- cmatrix[2,2]
TN <- cmatrix[1,1]
FP <- cmatrix[1,2]
FN <- cmatrix[2,1]
acc <- sum(diag(cmatrix)) / sum(cmatrix)
precision <- TP / (TP + FP)
recall <- TP / (TP + FN)
F1 <- 2*precision*recall / (precision + recall)
c(accuracy = acc, precision = precision, recall = recall, F1 = F1)
}
#4 metrike koje se koristimo za evaluaciju performansi klasifikacionog modela su:
#Accuracy (tacnost) - procenat tacnih predikcija (i pozitivnih i negativnih).
#Precision (preciznost) - procenat tacno predvidjenih pozitivnih predikcija u odnosu na sve pozitivne predikcije koje je model napravio
#Recall (odziv) - procenat pozitivnih predikcija koje je model tacno identifikovao kao pozitivne u odnosu na ukupan broj stvarno pozitivnih
#F1 - mera za balansiranje vrednosti precision-a i recall-a
# racunamo ove metrike za nas model
nb1.eval <- compute.eval.metrics(nb1.cm)
nb1.eval
# Tumacimo metrike u nasem domenu problema:
# Accuracy - za 73,2% jabuka smo tacno predvideli da li ce ostvariti dobar kvalitet ili ne
# Precision - od svih jabuka za koje smo predvideli da ce biti dobrog kvaliteta, 71,6% je stvarno dobrog kvaliteta
# Recall - od svih jabuka koje su zaista dobrog kvaliteta, mi smo tacno predvideli 77,3% njih
# F1 - mera koja balansira preciznost i odziv i kod nas iznosi 0.74
# ROC kriva
# kreiramo predikcije u formi verovatnoca
nb2.pred.prob <- predict(nb1, newdata = test.data, type = "raw")
nb2.pred.prob
library(pROC)
# kreiramo roc krivu
nb2.roc <- roc(response = as.numeric(test.data$IsGood),
predictor = nb2.pred.prob[,2])
nb2.roc$auc
# AUC je povrsina ispod ROC krive i ona nam govori koliko dobro model razdvaja pozitivnu od negativne klase, kod nas iznosi 0.81
nb2.coords <- coords(nb2.roc,
ret = c("spec", "sens", "thr"),
x = "local maximas", transpose = FALSE)
nb2.coords
# koristimo youden metodu koja ce nam odrediti prag verovatnoce kojim se maksimizuje suma sensitivity i specificity metrika
plot.roc(nb2.roc,
print.thres = TRUE,
print.thres.best.method = "youden")
prob.threshold <- 0.451
# vrsimo ponovo predvidjanje sa novim pragom verovatnoce
nb2.pred2 <- ifelse(test = nb2.pred.prob[,2] >= prob.threshold,
yes = "Yes",
no = "No")
nb2.pred2 <- as.factor(nb2.pred2)
# kreiramo novu matricu konfuzije
nb2.cm2 <- table(actual = test.data$IsGood, predicted = nb2.pred2)
nb2.cm2
# ponovo racunamo evalucione metrike
nb2.eval2 <- compute.eval.metrics(nb2.cm2)
data.frame(rbind(nb1.eval, nb2.eval2), row.names = c("nb1", "nb2"))
# kada uporedimo metrike, mozemo da zakljucimo da je tacnost ostala ista
# preciznost se vrlo malo smanjila
# odziv se povecao, samim tim i metrika F1
#####
# Kraj NB primera
#####