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=