Rozwiązanie raportu 5

Wgranie danych, sprawdzenie braków, podział zbioru

data(advertising)

data <- advertising
nan_values <- sum(is.na(data))
nan_values
## [1] 0
x <- model.matrix(sales~.,data)[, -1]
y <- data$sales
set.seed(2147)

trainset <- data %>%
  sample_frac(0.7)

testset <- data %>%
  setdiff(trainset)

x_train <- model.matrix(sales~., trainset)[,-1]
x_test <- model.matrix(sales~., testset)[,-1]

y_train <- trainset$sales
y_test <- testset$sales

Regresja grzbietowa - Ridge (L2)

cv_out <- cv.glmnet(x_train, y_train, alpha=0, folds=5)
cv_out
## 
## Call:  cv.glmnet(x = x_train, y = y_train, alpha = 0, folds = 5) 
## 
## Measure: Mean-Squared Error 
## 
##     Lambda Index Measure     SE Nonzero
## min 0.4048   100   2.526 0.3606       3
## 1se 0.8521    92   2.851 0.4303       3
best_lambda <- cv_out$lambda.min
cat(paste('Wartość lambdy, która najmocniej minimalizuje MSE na zbiorze treningowym to:', best_lambda))
## Wartość lambdy, która najmocniej minimalizuje MSE na zbiorze treningowym to: 0.404811764760053
plot(cv_out)

predict_ridge <- predict(cv_out, s=best_lambda, newx=x_test)
mse_test <- mean((predict_ridge - y_test)^2)

cat(paste('Metryka MSE na zbiorze testowym w regresji grzbietowej wyniosła:', mse_test))
## Metryka MSE na zbiorze testowym w regresji grzbietowej wyniosła: 4.425334281384
ridge_gen <- glmnet(x, y, alpha=0)
ridge_gen_coef <- predict(ridge_gen, type='coefficients', s=best_lambda)

cat(paste('Współczynniki przy użyciu najlepszej lambdy: \n'))
## Współczynniki przy użyciu najlepszej lambdy:
ridge_gen_coef
## 4 x 1 sparse Matrix of class "dgCMatrix"
##                      s1
## (Intercept) 3.634052677
## TV          0.042519738
## radio       0.173864503
## newspaper   0.002993224

Żaden ze współczynników nie jest zerowy, zatem wyniki pokazują, że regresja grzbietowa wykonała swoją pracę w poprawny sposób z wynikiem MSE na zbiorze testowym równym 4.425334281384.

Regresja LASSO (L1)

model_lasso <- glmnet(x_train, y_train, alpha=1, lamba=(10^seq(10,-2,length=100)))

plot(model_lasso)

cv_out_lasso <- cv.glmnet(x_train, y_train, alpha=1, folds=5)
cv_out_lasso
## 
## Call:  cv.glmnet(x = x_train, y = y_train, alpha = 1, folds = 5) 
## 
## Measure: Mean-Squared Error 
## 
##     Lambda Index Measure     SE Nonzero
## min 0.0741    44   2.465 0.3758       2
## 1se 0.4341    25   2.838 0.5049       2
best_lambda_lasso <- cv_out_lasso$lambda.min
cat(paste('Wartość lambdy, która najmocniej minimalizuje MSE na zbiorze treningowym to:', best_lambda_lasso))
## Wartość lambdy, która najmocniej minimalizuje MSE na zbiorze treningowym to: 0.0741104394060153
plot(cv_out_lasso)

predict_lasso <- predict(cv_out_lasso, s=best_lambda_lasso, newx=x_test)
mse_test_lasso <- mean((predict_lasso - y_test)^2)

cat(paste('Metryka MSE na zbiorze testowym w regresji LASSO wyniosła:', mse_test_lasso))
## Metryka MSE na zbiorze testowym w regresji LASSO wyniosła: 4.23978569962593
lasso_gen <- glmnet(x, y, alpha=1)
lasso_gen_coef <- predict(lasso_gen, type='coefficients', s=best_lambda)

cat(paste('Współczynniki przy użyciu najlepszej lambdy: \n'))
## Współczynniki przy użyciu najlepszej lambdy:
lasso_gen_coef
## 4 x 1 sparse Matrix of class "dgCMatrix"
##                     s1
## (Intercept) 4.18290846
## TV          0.04127349
## radio       0.16208021
## newspaper   .

Regresja Lasso wykorzystała swoją zdolność do redukcji zmiennych niezależnych i usunęła zmienną newspaper (współczynnik stojący przy tej zmiennej okazał sie 0).

model_ols <- lm(sales~., data=trainset)
predict_ols <- predict(model_ols, testset)
mse_ols <- mean((predict_ols - y_test)^2)

cat(paste('Metryka MSE na zbiorze testowym wyniosła:', mse_ols))
## Metryka MSE na zbiorze testowym wyniosła: 4.22464823572012
summary(model_ols)
## 
## Call:
## lm(formula = sales ~ ., data = trainset)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.5583 -0.5352  0.1459  1.1234  2.6577 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.172866   0.341101   9.302 3.12e-16 ***
## TV           0.044804   0.001524  29.391  < 2e-16 ***
## radio        0.192747   0.009526  20.233  < 2e-16 ***
## newspaper   -0.003063   0.006040  -0.507    0.613    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.507 on 136 degrees of freedom
## Multiple R-squared:  0.9156, Adjusted R-squared:  0.9138 
## F-statistic: 491.9 on 3 and 136 DF,  p-value: < 2.2e-16

Wnioski

Ad.1 Wybrałem zbiór danych Advertising ze strony ISLR.

Ad.2 Modelowałem sprzedaż z różnych źródeł reklamowania.

Ad.3 Oczekiwałem, że regresja grzbietowa okaże się lepsza od Lasso, gdyż charakterystyka wybranego zbioru danych sugerowała, że każda ze zmiennych niezależnych ma istotny wpływ na model. Niemniej jednak to regresja Lasso okazała się lepsza pod względem metryki błędu średniokwadratowego (MSE). Co ciekawe stało się tak, za sprawą przypisania zmiennej newspaper współczynnika równego 0. Co ciekawe na tle tych modeli i tak lepszy wynik osiągnęła metoda najmniejszych kwadratów, choć nieznacznie niższy od regresji Lasso.

Ad.4 W modelu regresji grzbietowej każda zmienna niezależna okazała się istotna, a kolejność ich wpływu na wynik to:

  • radio (1),
  • TV (2),
  • newspaper (3).

W przypadku regresjii Lasso zmienna newspaper została wyeliminowana z modelu. Warto dodać, że zmienna radio w tym modelu również okazała się najważniejsza. W modelu regresji liniowej (OLS) zmienna newspaper także okazała się statystycznie nieistotna.

LS0tDQp0aXRsZTogIk5pZWtsYXN5Y3puZSBtZXRvZHkgc3RhdHlzdHlraSINCnN1YnRpdGxlOiAiUklER0UgaSBMQVNTTyINCmF1dGhvcjogIlBpb3RyIFN6eXB1bGnFhHNraSINCmRhdGU6ICIyMDI0LTEyLTE0Ig0Kb3V0cHV0Og0KICBodG1sX2RvY3VtZW50Og0KICAgIHRoZW1lOiBjZXJ1bGVhbg0KICAgIGhpZ2hsaWdodDogdGV4dG1hdGUNCiAgICBmb250c2l6ZTogOHB0DQogICAgdG9jOiB0cnVlDQogICAgY29kZV9kb3dubG9hZDogdHJ1ZQ0KICAgIHRvY19mbG9hdDoNCiAgICAgIGNvbGxhcHNlZDogZmFsc2UNCmVkaXRvcl9vcHRpb25zOiANCiAgbWFya2Rvd246IA0KICAgIHdyYXA6IDcyDQotLS0NCg0KIyMgUm96d2nEhXphbmllIHJhcG9ydHUgNQ0KDQoNCmBgYHtyIHNldHVwLCBpbmNsdWRlPUZBTFNFfQ0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFKQ0KYGBgDQoNCmBgYHtyLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFLCBlY2hvPUZBTFNFfQ0KI2luc3RhbGwucGFja2FnZXMoJ2dsbW5ldCcpDQpsaWJyYXJ5KGdsbW5ldCkNCiNpbnN0YWxsLnBhY2thZ2VzKCdnbG10b29sYm94JykNCmxpYnJhcnkoZ2xtdG9vbGJveCkNCiNpbnN0YWxsLnBhY2thZ2VzKCdkcGx5cicpDQpsaWJyYXJ5KGRwbHlyKQ0KbGlicmFyeShnZ3Bsb3QyKQ0KYGBgDQoNCiMjIyBXZ3JhbmllIGRhbnljaCwgc3ByYXdkemVuaWUgYnJha8OzdywgcG9kemlhxYIgemJpb3J1DQoNCmBgYHtyfQ0KZGF0YShhZHZlcnRpc2luZykNCg0KZGF0YSA8LSBhZHZlcnRpc2luZw0KYGBgDQoNCmBgYHtyfQ0KbmFuX3ZhbHVlcyA8LSBzdW0oaXMubmEoZGF0YSkpDQpuYW5fdmFsdWVzDQpgYGANCg0KYGBge3J9DQp4IDwtIG1vZGVsLm1hdHJpeChzYWxlc34uLGRhdGEpWywgLTFdDQp5IDwtIGRhdGEkc2FsZXMNCmBgYA0KDQpgYGB7cn0NCnNldC5zZWVkKDIxNDcpDQoNCnRyYWluc2V0IDwtIGRhdGEgJT4lDQogIHNhbXBsZV9mcmFjKDAuNykNCg0KdGVzdHNldCA8LSBkYXRhICU+JQ0KICBzZXRkaWZmKHRyYWluc2V0KQ0KDQp4X3RyYWluIDwtIG1vZGVsLm1hdHJpeChzYWxlc34uLCB0cmFpbnNldClbLC0xXQ0KeF90ZXN0IDwtIG1vZGVsLm1hdHJpeChzYWxlc34uLCB0ZXN0c2V0KVssLTFdDQoNCnlfdHJhaW4gPC0gdHJhaW5zZXQkc2FsZXMNCnlfdGVzdCA8LSB0ZXN0c2V0JHNhbGVzDQpgYGANCg0KIyMjIFJlZ3Jlc2phIGdyemJpZXRvd2EgLSBSaWRnZSAoTDIpDQoNCmBgYHtyfQ0KY3Zfb3V0IDwtIGN2LmdsbW5ldCh4X3RyYWluLCB5X3RyYWluLCBhbHBoYT0wLCBmb2xkcz01KQ0KY3Zfb3V0DQoNCmJlc3RfbGFtYmRhIDwtIGN2X291dCRsYW1iZGEubWluDQpjYXQocGFzdGUoJ1dhcnRvxZvEhyBsYW1iZHksIGt0w7NyYSBuYWptb2NuaWVqIG1pbmltYWxpenVqZSBNU0UgbmEgemJpb3J6ZSB0cmVuaW5nb3d5bSB0bzonLCBiZXN0X2xhbWJkYSkpDQpgYGANCg0KYGBge3J9DQpwbG90KGN2X291dCkNCmBgYA0KDQpgYGB7cn0NCnByZWRpY3RfcmlkZ2UgPC0gcHJlZGljdChjdl9vdXQsIHM9YmVzdF9sYW1iZGEsIG5ld3g9eF90ZXN0KQ0KbXNlX3Rlc3QgPC0gbWVhbigocHJlZGljdF9yaWRnZSAtIHlfdGVzdCleMikNCg0KY2F0KHBhc3RlKCdNZXRyeWthIE1TRSBuYSB6YmlvcnplIHRlc3Rvd3ltIHcgcmVncmVzamkgZ3J6YmlldG93ZWogd3luaW9zxYJhOicsIG1zZV90ZXN0KSkNCmBgYA0KDQpgYGB7cn0NCnJpZGdlX2dlbiA8LSBnbG1uZXQoeCwgeSwgYWxwaGE9MCkNCnJpZGdlX2dlbl9jb2VmIDwtIHByZWRpY3QocmlkZ2VfZ2VuLCB0eXBlPSdjb2VmZmljaWVudHMnLCBzPWJlc3RfbGFtYmRhKQ0KDQpjYXQocGFzdGUoJ1dzcMOzxYJjenlubmlraSBwcnp5IHXFvHljaXUgbmFqbGVwc3plaiBsYW1iZHk6IFxuJykpDQpyaWRnZV9nZW5fY29lZg0KYGBgDQrFu2FkZW4gemUgd3Nww7PFgmN6eW5uaWvDs3cgbmllIGplc3QgemVyb3d5LCB6YXRlbSB3eW5pa2kgcG9rYXp1asSFLCDFvGUgcmVncmVzamEgZ3J6YmlldG93YSB3eWtvbmHFgmEgc3dvasSFIHByYWPEmSB3IHBvcHJhd255IHNwb3PDs2IgeiB3eW5pa2llbSBNU0UgbmEgemJpb3J6ZSB0ZXN0b3d5bSByw7N3bnltIDQuNDI1MzM0MjgxMzg0LiANCg0KIyMjIFJlZ3Jlc2phIExBU1NPIChMMSkNCg0KYGBge3J9DQptb2RlbF9sYXNzbyA8LSBnbG1uZXQoeF90cmFpbiwgeV90cmFpbiwgYWxwaGE9MSwgbGFtYmE9KDEwXnNlcSgxMCwtMixsZW5ndGg9MTAwKSkpDQoNCnBsb3QobW9kZWxfbGFzc28pDQpgYGANCg0KYGBge3J9DQpjdl9vdXRfbGFzc28gPC0gY3YuZ2xtbmV0KHhfdHJhaW4sIHlfdHJhaW4sIGFscGhhPTEsIGZvbGRzPTUpDQpjdl9vdXRfbGFzc28NCg0KYmVzdF9sYW1iZGFfbGFzc28gPC0gY3Zfb3V0X2xhc3NvJGxhbWJkYS5taW4NCmNhdChwYXN0ZSgnV2FydG/Fm8SHIGxhbWJkeSwga3TDs3JhIG5ham1vY25pZWogbWluaW1hbGl6dWplIE1TRSBuYSB6YmlvcnplIHRyZW5pbmdvd3ltIHRvOicsIGJlc3RfbGFtYmRhX2xhc3NvKSkNCmBgYA0KDQpgYGB7cn0NCnBsb3QoY3Zfb3V0X2xhc3NvKQ0KYGBgDQoNCg0KDQpgYGB7cn0NCnByZWRpY3RfbGFzc28gPC0gcHJlZGljdChjdl9vdXRfbGFzc28sIHM9YmVzdF9sYW1iZGFfbGFzc28sIG5ld3g9eF90ZXN0KQ0KbXNlX3Rlc3RfbGFzc28gPC0gbWVhbigocHJlZGljdF9sYXNzbyAtIHlfdGVzdCleMikNCg0KY2F0KHBhc3RlKCdNZXRyeWthIE1TRSBuYSB6YmlvcnplIHRlc3Rvd3ltIHcgcmVncmVzamkgTEFTU08gd3luaW9zxYJhOicsIG1zZV90ZXN0X2xhc3NvKSkNCmBgYA0KDQpgYGB7cn0NCmxhc3NvX2dlbiA8LSBnbG1uZXQoeCwgeSwgYWxwaGE9MSkNCmxhc3NvX2dlbl9jb2VmIDwtIHByZWRpY3QobGFzc29fZ2VuLCB0eXBlPSdjb2VmZmljaWVudHMnLCBzPWJlc3RfbGFtYmRhKQ0KDQpjYXQocGFzdGUoJ1dzcMOzxYJjenlubmlraSBwcnp5IHXFvHljaXUgbmFqbGVwc3plaiBsYW1iZHk6IFxuJykpDQpsYXNzb19nZW5fY29lZg0KYGBgDQoNClJlZ3Jlc2phIExhc3NvIHd5a29yenlzdGHFgmEgc3dvasSFIHpkb2xub8WbxIcgZG8gcmVkdWtjamkgem1pZW5ueWNoIG5pZXphbGXFvG55Y2ggaSB1c3VuxJnFgmEgem1pZW5uxIUgbmV3c3BhcGVyICh3c3DDs8WCY3p5bm5payBzdG9qxIVjeSBwcnp5IHRlaiB6bWllbm5laiBva2F6YcWCIHNpZSAwKS4gDQoNCmBgYHtyfQ0KbW9kZWxfb2xzIDwtIGxtKHNhbGVzfi4sIGRhdGE9dHJhaW5zZXQpDQpwcmVkaWN0X29scyA8LSBwcmVkaWN0KG1vZGVsX29scywgdGVzdHNldCkNCm1zZV9vbHMgPC0gbWVhbigocHJlZGljdF9vbHMgLSB5X3Rlc3QpXjIpDQoNCmNhdChwYXN0ZSgnTWV0cnlrYSBNU0UgbmEgemJpb3J6ZSB0ZXN0b3d5bSB3eW5pb3PFgmE6JywgbXNlX29scykpDQpgYGANCg0KYGBge3J9DQpzdW1tYXJ5KG1vZGVsX29scykNCmBgYA0KDQojIyBXbmlvc2tpIA0KDQpBZC4xIFd5YnJhxYJlbSB6YmnDs3IgZGFueWNoIEFkdmVydGlzaW5nIHplIHN0cm9ueSBJU0xSLg0KDQpBZC4yIE1vZGVsb3dhxYJlbSBzcHJ6ZWRhxbwgeiByw7PFvG55Y2ggxbpyw7NkZcWCIHJla2xhbW93YW5pYS4gDQoNCkFkLjMgT2N6ZWtpd2HFgmVtLCDFvGUgcmVncmVzamEgZ3J6YmlldG93YSBva2HFvGUgc2nEmSBsZXBzemEgb2QgTGFzc28sIGdkecW8IGNoYXJha3RlcnlzdHlrYSB3eWJyYW5lZ28gemJpb3J1IGRhbnljaCBzdWdlcm93YcWCYSwgxbxlIGthxbxkYSB6ZSB6bWllbm55Y2ggbmllemFsZcW8bnljaCBtYSBpc3RvdG55IHdwxYJ5dyBuYSBtb2RlbC4gTmllbW5pZWogamVkbmFrIHRvIHJlZ3Jlc2phIExhc3NvIG9rYXphxYJhIHNpxJkgbGVwc3phIHBvZCB3emdsxJlkZW0gbWV0cnlraSBixYLEmWR1IMWbcmVkbmlva3dhZHJhdG93ZWdvIChNU0UpLiBDbyBjaWVrYXdlIHN0YcWCbyBzacSZIHRhaywgemEgc3ByYXfEhSBwcnp5cGlzYW5pYSB6bWllbm5laiAqKm5ld3NwYXBlcioqIHdzcMOzxYJjenlubmlrYSByw7N3bmVnbyAwLiBDbyBjaWVrYXdlIG5hIHRsZSB0eWNoIG1vZGVsaSBpIHRhayBsZXBzenkgd3luaWsgb3NpxIVnbsSZxYJhIG1ldG9kYSBuYWptbmllanN6eWNoIGt3YWRyYXTDs3csIGNob8SHIG5pZXpuYWN6bmllIG5pxbxzenkgb2QgcmVncmVzamkgTGFzc28uIA0KDQpBZC40IFcgbW9kZWx1IHJlZ3Jlc2ppIGdyemJpZXRvd2VqIGthxbxkYSB6bWllbm5hIG5pZXphbGXFvG5hIG9rYXphxYJhIHNpxJkgaXN0b3RuYSwgYSBrb2xlam5vxZvEhyBpY2ggd3DFgnl3dSBuYSB3eW5payB0bzogDQoNCi0gcmFkaW8gKDEpLCANCi0gVFYgKDIpLCANCi0gbmV3c3BhcGVyICgzKS4NCg0KVyBwcnp5cGFka3UgcmVncmVzamlpIExhc3NvIHptaWVubmEgbmV3c3BhcGVyIHpvc3RhxYJhIHd5ZWxpbWlub3dhbmEgeiBtb2RlbHUuIFdhcnRvIGRvZGHEhywgxbxlIHptaWVubmEgKipyYWRpbyoqIHcgdHltIG1vZGVsdSByw7N3bmllxbwgb2themHFgmEgc2nEmSBuYWp3YcW8bmllanN6YS4gVyBtb2RlbHUgcmVncmVzamkgbGluaW93ZWogKE9MUykgem1pZW5uYSAqKm5ld3NwYXBlcioqIHRha8W8ZSBva2F6YcWCYSBzacSZIHN0YXR5c3R5Y3puaWUgbmllaXN0b3RuYS4gDQoNCg0KDQo=