#title Linear Regression
[[TableOfContents]]

==== 개요 ====
linear regression은 lm() 함수를 이용하면 된다. 이 문서의 목적은 최적화를 알아보기 위함이므로 optim() 함수를 이용하여 cost function()을 만들어 최적화하는 방법을 알아보겠다. 

==== data 준비 ====
trees 데이터셋을 이용할 것인데, 컬럼명이 대문자라 코딩하기 귀찮으니 소문자로 변경하자.
{{{
df <- trees
colnames(df) <- tolower(colnames(df))
head(df)
}}}

데이터가 준비되었다. 

==== simple regression ====
독립 변수가 1개일 경우을 simple regression이라고 한다. lm()을 사용하면 간단히 할 수 있다. 
{{{
lm.out <- lm(volume ~ girth, data=df)
summary(lm.out)
}}}

결과
{{{
> summary(lm.out)

Call:
lm(formula = volume ~ girth, data = df)

Residuals:
   Min     1Q Median     3Q    Max 
-8.065 -3.107  0.152  3.495  9.587 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -36.9435     3.3651  -10.98 7.62e-12 ***
girth         5.0659     0.2474   20.48  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.252 on 29 degrees of freedom
Multiple R-squared:  0.9353,	Adjusted R-squared:  0.9331 
F-statistic: 419.4 on 1 and 29 DF,  p-value: < 2.2e-16
}}}

위 결과는 최소제곱법으로 최적화된 결과로 아래와 같은 회귀식을 만들 수 있다. 

 volume = 5.0659 * girth - 36.9435

5.0659은 가중치고, -36.9435는 바이어스다. 

이것을 수동으로 해보자. 가설은 다음과 같다.  

 H(x)  = wx + b

최적화할 함수 즉, cost function은
 
 cost(w,b) = avg((H(x) - y)^2)

가 되겠다. 추정값(H(x))과 실제값(y)의 차이가 가장 적은 것을 찾는 것이 목적이다. 즉, cost(w,b)가 최소화되는 것이 목적이다. 제곱을 한 이유는 마이너스값이 나오지 않게 하기 위함과 추정값과 실제값의 차이가 큰 값에 패널티를 주기 위함이다. 다음과 같이 cost function을 만들고 optim() 함수를 이용해서 최적화를 해보자. 참고로 optim() 함수는 최소값을 찾는 것이 default다.

{{{
cost.f <- function(par, x, y){
  w <- par[1]
  b <- par[2]
  H <- w*x + b
  mean((H - y)^2)
}

result <- optim(par = c(0, 0), cost.f, x = df$girth, y = df$volume)
result
}}}

결과
{{{
> result
$par
[1]   5.066008 -36.945325

$value
[1] 16.91299

$counts
function gradient 
     103       NA 

$convergence
[1] 0

$message
NULL
}}}

위 결과를 보면

 * w = 5.066008
 * b = -36.945325

로 최적화했다. 위 lm() 함수의 결과인 volume = 5.0659 * girth - 36.9435와 거의 같다. 

==== multi regression ====
독립 변수가 2개 이상인 경우를 multi regression이라고 한다. multi regression도 마찬가지로 lm() 함수를 이용하는 것과 cost function을 만들어 최적하는 방법을 써보자. 먼저 lm() 함수를 이용하면 다음과 같은 결과를 얻을 수 있다. 
{{{
lm.out <- lm(volume ~ girth + height, data=df)
summary(lm.out)
}}}

결과
{{{
> summary(lm.out)

Call:
lm(formula = volume ~ girth + height, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.4065 -2.6493 -0.2876  2.2003  8.4847 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -57.9877     8.6382  -6.713 2.75e-07 ***
girth         4.7082     0.2643  17.816  < 2e-16 ***
height        0.3393     0.1302   2.607   0.0145 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.882 on 28 degrees of freedom
Multiple R-squared:  0.948,	Adjusted R-squared:  0.9442 
F-statistic:   255 on 2 and 28 DF,  p-value: < 2.2e-16
}}}

이제 cost function을 만들어야하는데, 문제는 독립변수가 n개 생겨 다음과 같이 hypothesis가 복잡해졌다는 것이다. 

 H(x)  = w1x1 + w2x2 + ... + b

이것을 단순화시키려면 대수적으로 접근해야하는데, 이것을 행렬로 취급하면 쉬워진다. 즉, 

 H(X) = WX + b 

와 같이 대문자 W와 대문자 X를 행렬로 취급하여 H()를 만들면 simple regression과 같이 생각할 수 있게 된다. b도 행렬 X에 포함시켜 더욱 단순하게

 H(X) = WX

와 같이 만들 수 있다. R로 cost function을 만들기 전에 우선 b를 데이터셋에 포함시키는데, 1로 세팅하면 된다. 

{{{
b <- 1
tmp <- cbind(b, df)
}}}

optim() 함수의 par 파라미터는 초기값을 주는 것인데 행렬의 곱(inner product)을 하려면 수학적 규칙을 맞춰야 한다. 예를 들면,

 * M1 = 3행 2열
 * M2 = 3행 2열

M1•M2를 하려면 M1(3x'''2''')와 M2('''3'''x2) 볼드체 부분(M1의 행수와 M2의 열수가 같아야 함)이 같아야만 되는 수학적 규칙이 있다. 그러므로 위 M1•M2는 할 수 없다. 위 예제의 데이터는 31행 3열(girth, height, b)이고 초기값은 3개가 입력되므로 3행 1열로 만들면 행렬의 곱을 할 수 있게 된다. R에서 행렬의 곱은 %*% 연산자를 쓰면 되므로 다음과 같이 cost function을 작성하면 된다. 
{{{
cost.f <- function(par, X, y){
  W <- as.matrix(par)
  X <- as.matrix(X)
  mean((X%*%W - y)^2)
}
}}}

cost function을 작성했으니 이제 최적화 해보자. 
{{{
result <- optim(par = c(0, 0, 0), cost.f, X = tmp[1:3], y = tmp[4])
result
}}}

결과
{{{
> result
$par
[1] -57.9996505   4.7082120   0.3394336

$value
[1] 13.61037

$counts
function gradient 
     238       NA 

$convergence
[1] 0

$message
NULL
}}}

결과를 보면 lm() 함수와 거의 유사하다. 

==== 응용 ====
다음과 같이 이상치가 섞인 데이터가 있다. 
{{{
set.seed(123)
mydata <- within(data.frame(x=1:10), y <- rnorm(x, mean=x))
mydata$y[2] <- 20
plot(mydata)
}}}
--출처: http://www.alastairsanderson.com/R/tutorials/robust-regression-in-R/
attachment:LinearRegression/outlier.png

lm() 함수는 최소제곱법으로 이상치에 취약하다. lm()함수를 쓰면 다음과 같은 결과를 얻는다. 
{{{
lm.out <- lm(y~x, data=mydata)
plot(mydata)
abline(lm.out, col="red")
}}}
attachment:LinearRegression/lm.png

최소제곱법 말고 다른 방법을 사용해보도록 하자. 실제값과 추정치의 차이에 제곱을하여 차이가 크면 패널티를 줬지만, 여기에서는 평균 절대 편차(MAD, mean absolute deviation)로 패널티를 주지 말아보자. 
{{{
cost.f <- function(par, x, y){
  w <- par[1]
  b <- par[2]
  H <- w*x + b
  mean(abs(H - y))
}

result <- optim(par = c(0, 0), cost.f, x = mydata$x, y = mydata$y)
result
}}}

결과
{{{
> result
$par
[1] 0.8849192 0.7047651

$value
[1] 2.37167

$counts
function gradient 
     249       NA 

$convergence
[1] 0

$message
NULL
}}}

그림으로 비교해보자. 
{{{
plot(mydata)
abline(lm.out, col="red")
abline(a = result$par[1], b = result$par[2], col = "blue")
}}}
attachment:LinearRegression/mad.png

뭔가 lm() 함수를 쓴 것보다 괜찮아 보이긴 한다. 하지만, 만족스럽지는 않다. 
참고로 robust regression의 결과는 다음과 같다. black line이 rlm()의 결과다. 
{{{
library(MASS)
rlm.out <- rlm(y~x, data=mydata)
abline(rlm.out)
}}}
attachment:LinearRegression/rlm.png

음.. 좋다. 잘 fitting됬다.