---
# System prepended metadata

title: Multivariate hw4
tags: [multivariate]

---

# 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