Projekt zaliczeniowy
Założenia
Projekt polega na wykonaniu analizy jednego z ogólnodostępnych
zbiorów danych metodami studiowanymi na zajęciach. Zbiór danych powinien
zawierać co najmniej około 200 obserwacji dotyczących co najmniej
kilkunastu zmiennych. Analiza powinna dotyczyć odpowiednio określonego
problemu regresji lub klasyfikacji. Analiza powinna obejmować:
- ocenę jakości predykcji zbudowanych modeli, czyli oszacowanie błędu
testowego predykcji metodą zbioru walidacyjnego lub metodą walidacji
krzyżowej;
- ocenę istotności każdego z predyktorów;
- ranking ważności predyktorów - o ile model na to pozwala;
- ocenę charakteru (dodatni/ujemny) i wielkości (o ile jednostek
zmienia się zmienna odpowiedzi przy jednostkowej zmianie predyktora)
wpływu istotnych predyktorów na odpowiedź - o ile model na to
pozwala.
Rozważany problem powinien być przeanalizowany przy pomocy
następujących modeli:
- model liniowy lub logistyczny,
- lasso,
- drzewo decyzyjne (odpowiednio przycięte, jeśli się da),
- las losowy.
Końcowy dokument może mieć format R Markdown lub podobny.
Wczytanie danych i przygotowanie danych
df <- read.csv("dataset/cars-simple.csv")
head(df)
#summary(df)
Wizualizacja danych
# Histogram ceny\
ggplot(df, aes(x = price)) +
geom_histogram(binwidth = 1000, fill = "skyblue", color = "black") +
labs(title = "Rozkład cen samochodów")

# Wykres zależności roku produkcji od przebiegu
ggplot(df, aes(x = year, y = odometer)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, color = "red") +
labs(title = "Zależność roku produkcji od przebiegu", x = "Rok produkcji", y = "Przebieg")

# Wykres zależności roku produkcji od ceny
ggplot(df, aes(x = year, y = price)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, color = "red") +
labs(title = "Zależność roku produkcji od ceny", x = "Rok produkcji", y = "Cena")

Model liniowy
subset <- df[c("price", "year", "odometer")]
# Fit linear regression model with price as the response variable and year and odometer as predictors
lm_compact <- lm(price ~ year + odometer, data = subset)
# Summary of the linear regression model
plot(lm_compact)




# Istotność predyktorów w modelu liniowym
summary(lm_compact)
Call:
lm(formula = price ~ year + odometer, data = subset)
Residuals:
Min 1Q Median 3Q Max
-3221.5 -1128.3 -44.9 969.4 3278.5
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.328e+06 5.060e+05 -6.577 1.86e-08 ***
year 1.662e+03 2.504e+02 6.637 1.48e-08 ***
odometer -5.760e-02 2.649e-02 -2.174 0.034 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1600 on 55 degrees of freedom
Multiple R-squared: 0.9451, Adjusted R-squared: 0.9431
F-statistic: 473.7 on 2 and 55 DF, p-value: < 2.2e-16
Ocena jakości predykcji
# Przykładowa walidacja krzyżowa dla modelu liniowego
library(caret)
# Przygotowanie danych do walidacji krzyżowej
set.seed(123)
indexes <- createDataPartition(df$price, p = 0.7, list = FALSE)
train_data <- df[indexes, ]
test_data <- df[-indexes, ]
# Budowa modelu liniowego na danych treningowych
linear_model <- lm(price ~ year + odometer, data = train_data)
# Predykcja na danych testowych
predictions <- predict(linear_model, newdata = test_data)
# Obliczenie błędu testowego
test_error <- RMSE(predictions, test_data$price)
test_error
[1] 1224.051
# wynik 1224.051 czy jest to akceptowalne???
Lasso
dfSmall <- df[c("price", "year", "odometer")]
x <- model.matrix(price ~ ., data = dfSmall)[,-1]
y <- dfSmall$price
lasso_model <- cv.glmnet(x, y, alpha = 1)
plot(lasso_model)

summary(lasso_model)
Length Class Mode
lambda 55 -none- numeric
cvm 55 -none- numeric
cvsd 55 -none- numeric
cvup 55 -none- numeric
cvlo 55 -none- numeric
nzero 55 -none- numeric
call 4 -none- call
name 1 -none- character
glmnet.fit 12 elnet list
lambda.min 1 -none- numeric
lambda.1se 1 -none- numeric
index 2 -none- numeric
Ocena jakości predykcji
# Przygotowanie danych
set.seed(123)
indexes <- createDataPartition(df$price, p = 0.7, list = FALSE)
train_data <- df[indexes, ]
test_data <- df[-indexes, ]
# Tworzenie macierzy modelowej dla danych treningowych
x_train <- model.matrix(price ~ ., data = train_data)[,-1]
Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
contrasts can be applied only to factors with 2 or more levels
Drzewo decyzyjne
tree_model <- rpart(price ~ year + manufacturer + odometer, data = df, method = "anova")
plot(tree_model)

Las losowy
library(randomForest)
random_forest_model <- randomForest(price ~ year + odometer + condition + state, data = df)
print(random_forest_model)
Call:
randomForest(formula = price ~ year + odometer + condition + state, data = df)
Type of random forest: regression
Number of trees: 500
No. of variables tried at each split: 1
Mean of squared residuals: 2122828
% Var explained: 95.2
varImpPlot(random_forest_model)

LS0tDQp0aXRsZTogIkFuYWxpemEgZGFueWNoIHNhbW9jaG9kb3d5Y2giDQphdXRob3I6ICJFbG1pcmEgUmVzaGV0aXVrIg0KZGF0ZTogImByIGZvcm1hdChTeXMudGltZSgpLCAnJVktJW0tJWQnKWAiDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KDQojIFByb2pla3QgemFsaWN6ZW5pb3d5DQoNCiMjIFphxYJvxbxlbmlhDQoNClByb2pla3QgcG9sZWdhIG5hIHd5a29uYW5pdSBhbmFsaXp5IGplZG5lZ28geiBvZ8OzbG5vZG9zdMSZcG55Y2ggemJpb3LDs3cgZGFueWNoIG1ldG9kYW1pIHN0dWRpb3dhbnltaSBuYSB6YWrEmWNpYWNoLiBaYmnDs3IgZGFueWNoIHBvd2luaWVuIHphd2llcmHEhyBjbyBuYWptbmllaiBva2/Fgm8gMjAwIG9ic2Vyd2FjamkgZG90eWN6xIVjeWNoIGNvIG5ham1uaWVqIGtpbGt1bmFzdHUgem1pZW5ueWNoLiBBbmFsaXphIHBvd2lubmEgZG90eWN6ecSHIG9kcG93aWVkbmlvIG9rcmXFm2xvbmVnbyBwcm9ibGVtdSByZWdyZXNqaSBsdWIga2xhc3lmaWthY2ppLiBBbmFsaXphIHBvd2lubmEgb2Jlam1vd2HEhzoNCg0KLSBvY2VuxJkgamFrb8WbY2kgcHJlZHlrY2ppIHpidWRvd2FueWNoIG1vZGVsaSwgY3p5bGkgb3N6YWNvd2FuaWUgYsWCxJlkdSB0ZXN0b3dlZ28gcHJlZHlrY2ppIG1ldG9kxIUgemJpb3J1IHdhbGlkYWN5am5lZ28gbHViIG1ldG9kxIUgd2FsaWRhY2ppIGtyennFvG93ZWo7DQotIG9jZW7EmSBpc3RvdG5vxZtjaSBrYcW8ZGVnbyB6IHByZWR5a3RvcsOzdzsNCi0gcmFua2luZyB3YcW8bm/Fm2NpIHByZWR5a3RvcsOzdyAtIG8gaWxlIG1vZGVsIG5hIHRvIHBvendhbGE7DQotIG9jZW7EmSBjaGFyYWt0ZXJ1IChkb2RhdG5pL3VqZW1ueSkgaSB3aWVsa2/Fm2NpIChvIGlsZSBqZWRub3N0ZWsgem1pZW5pYSBzacSZIHptaWVubmEgb2Rwb3dpZWR6aSBwcnp5IGplZG5vc3Rrb3dlaiB6bWlhbmllIHByZWR5a3RvcmEpIHdwxYJ5d3UgaXN0b3RueWNoIHByZWR5a3RvcsOzdyBuYSBvZHBvd2llZMW6IC0gbyBpbGUgbW9kZWwgbmEgdG8gcG96d2FsYS4NCg0KUm96d2HFvGFueSBwcm9ibGVtIHBvd2luaWVuIGJ5xIcgcHJ6ZWFuYWxpem93YW55IHByenkgcG9tb2N5IG5hc3TEmXB1asSFY3ljaCBtb2RlbGk6DQoNCi0gbW9kZWwgbGluaW93eSBsdWIgbG9naXN0eWN6bnksDQotIGxhc3NvLA0KLSBkcnpld28gZGVjeXp5am5lIChvZHBvd2llZG5pbyBwcnp5Y2nEmXRlLCBqZcWbbGkgc2nEmSBkYSksDQotIGxhcyBsb3Nvd3kuDQoNCktvxYRjb3d5IGRva3VtZW50IG1vxbxlIG1pZcSHIGZvcm1hdCBSIE1hcmtkb3duIGx1YiBwb2RvYm55Lg0KDQoNCmBgYHtyIHNldHVwLCBpbmNsdWRlPUZBTFNFfSANCmtuaXRyOjpvcHRzX2NodW5rJHNldChlY2hvID0gVFJVRSkgDQpsaWJyYXJ5KHJtYXJrZG93bikNCmxpYnJhcnkodGlkeXZlcnNlKSANCmxpYnJhcnkoZ2xtbmV0KQ0KbGlicmFyeShycGFydCkNCmBgYA0KDQojIyBXY3p5dGFuaWUgZGFueWNoIGkgcHJ6eWdvdG93YW5pZSBkYW55Y2gNCmBgYHtyfQ0KZGYgPC0gcmVhZC5jc3YoImRhdGFzZXQvY2Fycy1zaW1wbGUuY3N2IikNCmhlYWQoZGYpDQojc3VtbWFyeShkZikNCmBgYA0KIyMgV2l6dWFsaXphY2phIGRhbnljaA0KYGBgIHtyfQ0KIyBIaXN0b2dyYW0gY2VueVwNCmdncGxvdChkZiwgYWVzKHggPSBwcmljZSkpICsNCiAgZ2VvbV9oaXN0b2dyYW0oYmlud2lkdGggPSAxMDAwLCBmaWxsID0gInNreWJsdWUiLCBjb2xvciA9ICJibGFjayIpICsNCiAgbGFicyh0aXRsZSA9ICJSb3prxYJhZCBjZW4gc2Ftb2Nob2TDs3ciKQ0KDQojIFd5a3JlcyB6YWxlxbxub8WbY2kgcm9rdSBwcm9kdWtjamkgb2QgcHJ6ZWJpZWd1DQpnZ3Bsb3QoZGYsIGFlcyh4ID0geWVhciwgeSA9IG9kb21ldGVyKSkgKw0KICBnZW9tX3BvaW50KCkgKw0KICBnZW9tX3Ntb290aChtZXRob2QgPSAibG0iLCBzZSA9IEZBTFNFLCBjb2xvciA9ICJyZWQiKSArDQogIGxhYnModGl0bGUgPSAiWmFsZcW8bm/Fm8SHIHJva3UgcHJvZHVrY2ppIG9kIHByemViaWVndSIsIHggPSAiUm9rIHByb2R1a2NqaSIsIHkgPSAiUHJ6ZWJpZWciKQ0KDQoNCiMgV3lrcmVzIHphbGXFvG5vxZtjaSByb2t1IHByb2R1a2NqaSBvZCBjZW55DQpnZ3Bsb3QoZGYsIGFlcyh4ID0geWVhciwgeSA9IHByaWNlKSkgKw0KICBnZW9tX3BvaW50KCkgKw0KICBnZW9tX3Ntb290aChtZXRob2QgPSAibG0iLCBzZSA9IEZBTFNFLCBjb2xvciA9ICJyZWQiKSArDQogIGxhYnModGl0bGUgPSAiWmFsZcW8bm/Fm8SHIHJva3UgcHJvZHVrY2ppIG9kIGNlbnkiLCB4ID0gIlJvayBwcm9kdWtjamkiLCB5ID0gIkNlbmEiKQ0KYGBgDQojIyBNb2RlbCBsaW5pb3d5DQpgYGB7cn0NCnN1YnNldCA8LSBkZltjKCJwcmljZSIsICJ5ZWFyIiwgIm9kb21ldGVyIildDQojIEZpdCBsaW5lYXIgcmVncmVzc2lvbiBtb2RlbCB3aXRoIHByaWNlIGFzIHRoZSByZXNwb25zZSB2YXJpYWJsZSBhbmQgeWVhciBhbmQgb2RvbWV0ZXIgYXMgcHJlZGljdG9ycw0KbG1fY29tcGFjdCA8LSBsbShwcmljZSB+IHllYXIgKyBvZG9tZXRlciwgZGF0YSA9IHN1YnNldCkNCg0KIyBTdW1tYXJ5IG9mIHRoZSBsaW5lYXIgcmVncmVzc2lvbiBtb2RlbA0KcGxvdChsbV9jb21wYWN0KQ0KDQojIElzdG90bm/Fm8SHIHByZWR5a3RvcsOzdyB3IG1vZGVsdSBsaW5pb3d5bQ0Kc3VtbWFyeShsbV9jb21wYWN0KQ0KYGBgDQoNCiMjIyBPY2VuYSBqYWtvxZtjaSBwcmVkeWtjamkNCg0KYGBge3J9DQojIFByenlrxYJhZG93YSB3YWxpZGFjamEga3J6ecW8b3dhIGRsYSBtb2RlbHUgbGluaW93ZWdvDQpsaWJyYXJ5KGNhcmV0KQ0KDQojIFByenlnb3Rvd2FuaWUgZGFueWNoIGRvIHdhbGlkYWNqaSBrcnp5xbxvd2VqDQpzZXQuc2VlZCgxMjMpDQppbmRleGVzIDwtIGNyZWF0ZURhdGFQYXJ0aXRpb24oZGYkcHJpY2UsIHAgPSAwLjcsIGxpc3QgPSBGQUxTRSkNCnRyYWluX2RhdGEgPC0gZGZbaW5kZXhlcywgXQ0KdGVzdF9kYXRhIDwtIGRmWy1pbmRleGVzLCBdDQoNCiMgQnVkb3dhIG1vZGVsdSBsaW5pb3dlZ28gbmEgZGFueWNoIHRyZW5pbmdvd3ljaA0KbGluZWFyX21vZGVsIDwtIGxtKHByaWNlIH4geWVhciArIG9kb21ldGVyLCBkYXRhID0gdHJhaW5fZGF0YSkNCg0KIyBQcmVkeWtjamEgbmEgZGFueWNoIHRlc3Rvd3ljaA0KcHJlZGljdGlvbnMgPC0gcHJlZGljdChsaW5lYXJfbW9kZWwsIG5ld2RhdGEgPSB0ZXN0X2RhdGEpDQoNCiMgT2JsaWN6ZW5pZSBixYLEmWR1IHRlc3Rvd2Vnbw0KdGVzdF9lcnJvciA8LSBSTVNFKHByZWRpY3Rpb25zLCB0ZXN0X2RhdGEkcHJpY2UpDQp0ZXN0X2Vycm9yDQojIHd5bmlrIDEyMjQuMDUxIGN6eSBqZXN0IHRvIGFrY2VwdG93YWxuZT8/Pw0KYGBgDQojIyBMYXNzbw0KYGBge3J9DQpkZlNtYWxsIDwtIGRmW2MoInByaWNlIiwgInllYXIiLCAib2RvbWV0ZXIiKV0NCg0KeCA8LSBtb2RlbC5tYXRyaXgocHJpY2UgfiAuLCBkYXRhID0gZGZTbWFsbClbLC0xXQ0KeSA8LSBkZlNtYWxsJHByaWNlDQpsYXNzb19tb2RlbCA8LSBjdi5nbG1uZXQoeCwgeSwgYWxwaGEgPSAxKQ0KcGxvdChsYXNzb19tb2RlbCkNCnN1bW1hcnkobGFzc29fbW9kZWwpDQpgYGANCiMjIyBPY2VuYSBqYWtvxZtjaSBwcmVkeWtjamkNCg0KYGBge3J9DQojVE9ETw0KYGBgDQoNCiMjIERyemV3byBkZWN5enlqbmUNCg0KYGBge3J9DQp0cmVlX21vZGVsIDwtIHJwYXJ0KHByaWNlIH4geWVhciArIG1hbnVmYWN0dXJlciArIG9kb21ldGVyLCBkYXRhID0gZGYsIG1ldGhvZCA9ICJhbm92YSIpDQpwbG90KHRyZWVfbW9kZWwpDQpgYGANCiMjIExhcyBsb3Nvd3kNCmBgYHtyfQ0KbGlicmFyeShyYW5kb21Gb3Jlc3QpDQpyYW5kb21fZm9yZXN0X21vZGVsIDwtIHJhbmRvbUZvcmVzdChwcmljZSB+IHllYXIgKyBvZG9tZXRlciArIGNvbmRpdGlvbiArIHN0YXRlLCBkYXRhID0gZGYpDQpwcmludChyYW5kb21fZm9yZXN0X21vZGVsKQ0KdmFySW1wUGxvdChyYW5kb21fZm9yZXN0X21vZGVsKQ0KYGBgDQo=