# Multivariate hw4 ###### tags: `multivariate` ## Basic LM code ``` R oridata = read.csv("cars.dat", sep="") data = oridata[, c("M", "W", "D", "C")] data$C = factor(data$C) lm(formula = M ~ W + D + C, data = data) ``` result ``` Multiple R-squared: 0.6958, Adjusted R-squared: 0.6782 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 42.060469 2.223302 18.918 < 2e-16 *** W -0.007329 0.001138 -6.443 1.34e-08 *** D 0.008365 0.009928 0.843 0.402 C2 -0.326517 1.278132 -0.255 0.799 C3 -2.017549 1.259087 -1.602 0.114 ``` We found that Weight is the most significant coefficient in this regression. ![](https://i.imgur.com/Y0dsmvb.png) ## Feature selection Foward with AIC `foward = step(lm(M ~ 1, data), ~ W + D + C, data=data, direction="forward")` ``` Step: AIC=179.12 Call: lm(formula = M ~ W, data = data) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 39.5829916 1.5369186 25.75 <2e-16 *** W -0.0060733 0.0004942 -12.29 <2e-16 *** Multiple R-squared: 0.6772, Adjusted R-squared: 0.6727 ``` Foward with BIC `foward = step(lm(M ~ 1, data), ~ W + D + C, data=data, direction="forward", k=log(n))` `k` is the value to calculate the loss, default AIC is `k=2` The result seletion is same as above ``` Step: AIC=183.73 M ~ W Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 39.5829916 1.5369186 25.75 <2e-16 *** W -0.0060733 0.0004942 -12.29 <2e-16 *** Multiple R-squared: 0.6772, Adjusted R-squared: 0.6727 ``` Backward ``` R backward_aic = step(fit, data=data, direction="backward") backward_bic = step(fit, data=data, direction="backward", k=log(n)) ``` Backward AIC share same result with forward AIC Backward BIC share same result with forward BIC By permutation ``` R library(leaps) fit = regsubsets(M ~ ., data=data, intercept=TRUE, method="exhaustive") summary(fit) ``` ``` 1 subsets of each size up to 3 Selection Algorithm: exhaustive W D C2 C3 1 ( 1 ) "*" " " " " " " 2 ( 1 ) "*" " " " " "*" 3 ( 1 ) "*" "*" " " "*" 4 ( 1 ) "*" "*" "*" "*" ``` By R_adj ``` R library("SignifReg") fit_adj = add1SignifReg(lm(M ~ 1, data), criterion="r-adj") ``` result ``` lm(formula = M ~ W, data = data) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 39.5829916 1.5369186 25.75 <2e-16 *** W -0.0060733 0.0004942 -12.29 <2e-16 *** Residual standard error: 3.31 on 72 degrees of freedom Multiple R-squared: 0.6772, Adjusted R-squared: 0.6727 F-statistic: 151 on 1 and 72 DF, p-value: < 2.2e-16 ``` This result is same as above. ## Regulation L1 loss(LASSO) and l2 loss(Ridge) ``` R library(glmnet) fit_lasso <- cv.glmnet(as.matrix(norm_data[, 2:6]), as.matrix(norm_data[, 1]), alpha=1) fit_ridge <- cv.glmnet(as.matrix(norm_data[, 2:6]), as.matrix(norm_data[, 1]), alpha=0) ``` ![](https://i.imgur.com/rZfdaVr.png) ![](https://i.imgur.com/wi7swnf.png) The lambda of ridge is ``` $lambda.min [1] 0.01258925 $lambda.1se [1] 0.7943282 ``` And of lasso is ``` $lambda.min [1] 0.1258925 $lambda.1se [1] 0.3162278 ``` Then calculate the R2 manually? ``` R test_x = as.matrix(norm_data[,2:6]) test_y = as.matrix(norm_data[,1]) pred_y = predict(fit_lasso, s=fit_lasso$lambda.min, newx=test_x) R2 = 1 - sum((test_y - pred_y) ** 2) / sum(var(test_y) * (n - 1)) ``` The R2 with lambda.min and lambda.1se ``` 0.6889514 0.6610924 ``` which is larger than feature selection method, because the model parameters is bigger. Reference * https://www.rdocumentation.org/packages/glmnet/versions/2.0-18/topics/cv.glmnet * https://www.rdocumentation.org/packages/glmnet/versions/2.0-18/topics/glmnet # Code https://gist.github.com/linnil1/570fd6ef5424313fff8b8151256b9864