R语言--线性回归分析（1）--simple linear regression

此博客仅为个人理解和上课笔记梳理，如果发现有误，请各位大佬不吝赐教！

一、回归函数

---

通过不同的X（predictor，independent variable, covariate, input, etc.)去估计Y值(the dependent variable, output, response, etc.)

回归分析： 一种通过分析两个或多个变量间关系，以达到可以通过其他量的值来预测所需量的值的统计学方法。 也就是说，等同于预测一个函数关系：$Y = f(x)$

模型(model)：变量间关系的估计 （函数关系，统计关系）
1. 函数关系： x和Y之间的关系可以通过函数来表达$Y_i = f(x_i)$, 如果我们知道$x_i$的值，就可以直接获得$Y_i$的值。这种关系是具有确定性的(deterministic)。
2. 统计关系： 一般描述为 trend + error。 例如：信号和噪声。 这种关系是具有随机性的(stochastic)。它可以认为是一个函数关系加上噪声。

经验模型(empirical model)： $Y = f(x) + \epsilon(x)$， 其中$\epsilon(x)$是均值为0，彼此不相关$corr(\epsilon_i,\epsilon_j) = 0$的噪声。

$corr(x,y) = 0$ ： 相关关系是一种非确定性的关系，相关系数是研究变量之间线性关系程度的量。

1. 若 X和Y不相关，corr(x,y) = 0。只能说明X和Y之间不存在线性关系，但不排除X和Y之间可能存在其他关系。
2. 若X和Y独立，则一定有corr(x,y) = 0，若X和Y不相关，则仅仅是不存在线性关系。

使用有限的数据去估计函数关系几乎是不可能的，为了使得问题可解，我们定义：
$y_i = \beta_0 + \beta_1x_i + \epsilon_i, i = 1,...,n$

Note:
1. $x_i$不仅仅局限于线性，我们可以定义$x^* = x_i^2, x^* = logx_i$,等等。
2. $E(\epsilon_i)=0$ and $Var(\epsilon_i) = \delta^2$
3. 因为$\epsilon$的存在，即使假设$x_i$是常数，Y依然是一个随机量。

二、对回归参数的解释

---

example：

$Y = 9.5 + 2.1x + \epsilon$ 。其中x的范围为20-80
$\beta_0$ : intercept（截距）： 当x = 0 时， Y的值。 但是当x = 0 在观测值范围以外时，$\beta_0$即没有相关的意义。
$\beta_1$ : slope（斜率），每当x增加1时，Y的变化。

最小均方误差

为了确定回归函数，我们必须在一定程度上尽可能的降低预测误差（$Y_i - \beta_0 - \beta1x_i$)
$Q_1(\beta_0,\beta_1) = \sum_{i=0}^n|y_i-(\beta_0+\beta_1x_i)|$
$Q_2(\beta_0,\beta_1) = \sum_{i=0}^n\{y_i-(\beta_0+\beta_1x_i)\}^2$
$Q_3(\beta_0,\beta_1) = \sum_{i=0}^n\{y_i-(\beta_0+\beta_1x_i)\}$
至于为什么线性回归问题选择了最小均方误差，即$Q_2$，请参考最小均方误差

对$Q_2$求导等于0之后可以获得：
$\hat\beta_1 = \sum_{i=1}^n\frac{(x_i-\bar x)(y_i - \bar y)}{\sum_{i=1}^n(x_i-\bar x)}$
$\hat\beta_0 = \bar y - \hat\beta_1\bar x$

最小均方估计回归曲线： $\hat Y = \hat\beta_0+\hat\beta_1x$

为什么最小均方估计是好的?
1. 他们是无偏估计量：$E(\hat\beta_0) = \beta_0$ , $E(\hat\beta_1) = \beta_1$
2. 在所有的线性无偏估计量中，他们具有最小的方差，他们是最佳线性无偏估计量（best linear unbiased estimators (BLUEs) )

偏差和估计值

residual（偏差）： $e_i = y_i - \hat y_i = y_i - (\hat\beta_0+\beta_1x_i)$
fitted value（估计值）: $\hat y_i = \hat\beta_0 + \hat\beta_1x_i, i = 1,2...,n$
Note:
1. $e_i = y_i - \hat y_i \neq \epsilon_i$ 其中$\hat y_i$是基于$x_i$的，而$\epsilon_i$的值是随机的

Properties of the Fitted Regression Line

1. $\sum_{i=0}^ne_i = 0$

2. $\sum_{i=0}^ny_i = \sum_{i=0}^n\hat y_i$

3. $\sum_{i=0}^nx_ie_i = 0$

4. $\sum_{i=0}^n\hat y_ie_i = 0$
   Note:
   在证明这些公式的过程中，有个小技巧可以用
   $\sum(y_i-\bar y) = 0$

Estimating the Error Variance

sum of square : $\sum_{i=1}^n(y_i-\bar y)^2$

mean square : $s^2 = \frac{1}{n-1}\sum_{i=1}^n(y_i - \bar y)^2$

SSE( sum of squared residuals) : $SSE = \sum_{i=1}^n(y_i-\hat y_i)^2$

MSE( mean squared residuals) : $MSE = \frac{\sum_{i=1}^n(y_i-\hat y_i)^2}{n-2} = \frac{SSE}{n-2} = \frac{\sum_{i=1}^ne_i^2}{n-2}$

Note:
1. MSE 是 $\delta^2$的无偏估计量， $E(MSE) = \delta^2$
2. $\hat\delta = \sqrt {MSE}$ , 很多情况下我们可以知道MSE，即可用此公式求出$\hat\delta$

The Normal Error Regression Model

$\epsilon_i \sim^{iid} N(0,\delta^2) , i = 1,2,...,n$ 独立同分布(iid)

$E(\epsilon) = 0$

$Var(\epsilon) = \delta^2$

$Corr(\epsilon_i,\epsilon_j) = 0, i \neq j$

rm(list=ls())

#####################################################
# Toluca example
#载入x和y的数据，默认为多列一行
x<-c(80,30,50,90,70,60,120,80,100,50,40,70,90,20,110,100,30,50,90,110,30,90,40,80,70)
Y<-c(399,121,221,376,361,224,546,352,353,157,160,252,389,113,435,420,212,268,377,421,273,468,244,342,323)

# 绘制散点图，其中xlab和ylab分别定义了坐标名称 
plot(x,Y,xlab="Lot Size (Parts/lot)",ylab="Work Hours" ) 

# Fit 一个model，该方法会返回beta_0和beta_1的值
# b获得了fit的model的参数，即直接组成了beta的矩阵
# seq函数返回的是20-120范围内，以1为间隔的值，b[1]为beta_0的值，b[2]为beta_1的值，lines（x范围，线性方程）绘制直线。
lm(Y~x)
b<-lm(Y~x)$coefficients
lines(seq(20,120,1),b[1]+b[2]*seq(20,120,1))

# 不使用lm方法，而是通过矩阵的方法获得beta的值
# 其中只需要注意构造x的矩阵，第一列全为1
# rep(1,length(seq(20,120,1))) 构造了第一列值全为1，且行数等于seq(20,120,1)的长度，cbind()将两个矩阵左右相连 
cov.p<-cbind(rep(1,length(seq(20,120,1))),seq(20,120,1))
plot(x,Y,xlab="Lot Size (Parts/lot)",ylab="Work Hours")
# %*%表示矩阵的点乘，col控制颜色
lines(seq(20,120,1), cov.p%*%(lm(Y~x)$coefficients),col="red")

从图中可以很方便的获得$\beta_0 = 62.37, \beta_1 = 3.57$

# 在上面的散点图中绘制出误差值 
for(i in 1:length(x)){
lines(c(x[i],x[i]),c(b[1]+b[2]*x[i],Y[i]),col="green")
}

# 计算MSE，SSE
Y.hat<-b[1]+b[2]*x
n<-length(Y)
SSE<-sum( (Y-Y.hat)^2)
MSE<-SSE/(n-2)
sqrt(MSE)

# We can also get this from the summary of the regression output given by R
summary(lm(Y~x))

对此图的详解将放到后面的博客中。