数据科学24 | 回归模型-基本概念与最小二乘法

2020-07-03 16:49:19 浏览数 (1)

回归分析在统计学中非常重要，目的在于了解两个或多个变量间是否相关、相关方向与强度，并建立数学模型以便观察特定变量来预测研究者感兴趣的变量。回归分析可以帮助人们了解在只有一个自变量变化时因变量的变化量。

用一个简单的例子介绍最小二乘回归法拟合线性模型：

例：UsingR包的galton数据集，包括配对的父母和孩子的身高。

查看父母身高和孩子身高的边缘分布，父母性别的身高差异通过母亲身高乘1.8校正：

代码语言：javascript复制
library(UsingR)
data(galton)
library(reshape)
long <- melt(galton)
g <- ggplot(long, aes(x = value, fill = variable))
g <- g   geom_histogram(color = "black", binwidth =1)
g <- g   facet_grid(.~variable)
g

图1.孩子和父母身高的边缘分布

用父母的身高预测孩子的身高，不考虑父母的身高时，利用最小二乘法求孩子身高的最佳预测?：令

Y_i

为第

i

个孩子的身高，

i=1,…,n

，当

sumlimits_{i=1}^n(Y_i-?)^2

最小时，孩子身高的真实值与预测值的差值最小，即残差平方和最小，此时?即为孩子身高的最佳预测，等于孩子身高分布估计的均值。

用manipulate()函数查看不同?值下残差平方的平均值变化：

代码语言：javascript复制
library(manipulate)
myHist <- function(mu){
  mse <- mean((galton$child - mu)^2) #对残差平方取均值而不是求和
  g <- ggplot(galton, aes(x = child))   geom_histogram(fill = "salmon", color = "black", binwidth = 1)
  g <- g   geom_vline(xintercept = mu, size = 3)
  g <- g   ggtitle(paste("mu = ", mu, "MSE = ", round(mse, 2), sep = ""))
  g
}
manipulate(myHist(mu), mu = slider(62, 74, step = 0.5))

图2.不同?值下残差平方的平均值变化

可以看到?值变大向分布中心靠近时，残差平方的均值变小；?值从分布中心继续变大时，残差平方的均值重新变大。当?等于孩子身高均值时，残差平方的均值最小，即孩子身高的最小二乘估计是孩子身高的均值。

代码语言：javascript复制
g <- ggplot(galton, aes(x = child))   geom_histogram(fill = "salmon", color = "black", binwidth = 1)
g <- g   geom_vline(xintercept = mean(galton$child), size = 3)
g

图3.孩子身高的均值

证明孩子身高的均值

bar{Y}

是使公式

sumlimits_{i=1}^n(Y_i-?)^2

最小的?值：

begin{split} sumlimits_{i=1}^n(Y_i-?)^2 &=sumlimits_{i=1}^n(Y_i-bar{Y} bar{Y}-?)^2 \ &=sumlimits_{i=1}^n(Y_i-bar{Y})^2 2sumlimits_{i=1}^n(Y_i-bar{Y})(bar{Y}-?) sumlimits_{i=1}^n(bar{Y}-?)^2\ &=sumlimits_{i=1}^n(Y_i-bar{Y})^2 2(bar{Y}-?)sumlimits_{i=1}^n(Y_i-bar{Y}) sumlimits_{i=1}^n(bar{Y}-?)^2\ &=sumlimits_{i=1}^n(Y_i-bar{Y})^2 2(bar{Y}-?)sumlimits_{i=1}^n(Y_i-nbar{Y}) sumlimits_{i=1}^n(bar{Y}-?)^2\ &=sumlimits_{i=1}^n(Y_i-bar{Y})^2 sumlimits_{i=1}^n(bar{Y}-?)^2\ &≥sumlimits_{i=1}^n(Y_i-bar{Y})^2 end{split}

即?等于孩子身高均值

bar{Y}

时，残差平方和最小。

比较配对的父母身高和孩子身高：

代码语言：javascript复制
ggplot(galton, aes(x=parent, y=child))  geom_point()

图4.父母身高及相应的孩子身高的散点图

这个图中有许多点被重复绘制，数据的频数信息没有被展示出来。

代码语言：javascript复制
library(UsingR)
data(galton)
library(dplyr)
library(ggplot2)
freqData <- as.data.frame(table(galton$parent,galton$child)) 
names(freqData) <- c("child","parent","freq")
freqData$parent<- as.numeric(as.character(freqData$parent))
freqData$child <- as.numeric(as.character(freqData$child))
g <- ggplot(filter(freqData, freq > 0), aes(x = parent, y = child))
g <- g   scale_size(range = c(2, 20), guide = "none")
g <- g   geom_point(color = "grey50", aes(size = freq))
g <- g   geom_point(aes(color=freq,size = freq))
g <- g   scale_color_gradient(low = "lightblue", high = "white")
g

图5.父母身高与孩子身高关系的气泡图

气泡大小及颜色深浅表示在特定父母身高与相应孩子身高的配对组合的数量。

最小二乘法拟合线性模型解释父母身高与孩子身高的关系，令回归线经过原点，即截距为0，这条线可用

y=?x

表示。令

X_i

为父母身高，最适合的线性模型的斜率?使实际观测值与预测值之间的残差平方和

sumlimits_{i=1}^n(Y_i-X_i?)^2

最小。

使用manipulate()函数查看不同?值的残差平方和变化：

代码语言：javascript复制
y <- galton$child - mean(galton$child) 
x <- galton$parent - mean(galton$parent) #通过减去均值使数据回归线经过原点
freqData <- as.data.frame(table(x,y)) 
names(freqData) <- c("child","parent","freq")
freqData$parent<- as.numeric(as.character(freqData$parent))
freqData$child <- as.numeric(as.character(freqData$child))
myPlot <- function(beta){
  g <- ggplot(filter(freqData, freq > 0), aes(x = parent, y = child))
  g <- g   scale_size(range = c(2, 20), guide = "none")
  g <- g   geom_point(color = "grey50", aes(size = freq 20))
  g <- g   geom_point(aes(color=freq,size = freq))
  g <- g   scale_color_gradient(low = "lightblue", high = "white")
  g <- g   geom_abline(intercept = 0, slope = beta, size = 3)
  mse <- mean((y - beta * x)^2)
  g <- g   ggtitle(paste("beta = ", beta, "mse = ", round(mse, 3)))
  g
}
manipulate(myPlot(beta), beta = slider(0.6, 1.2, step = 0.02))

图6.不同?值的残差平方和变化

可以看到，斜率?=0.64时，残差平方和最小。可以用

y=0.64x

预测孩子的身高。

在R中可以用lm()函数快速拟合线性模型。

代码语言：javascript复制
lm(I(child-mean(child))~I(parent-mean(parent))-1,data=galton)

Call:
lm(formula = I(child - mean(child)) ~ I(parent - mean(parent)) - 
    1, data = galton)

Coefficients:
I(parent - mean(parent))  
                  0.6463

可以在图5基础上重新绘制线性回归线：

代码语言：javascript复制
freqData <- as.data.frame(table(galton$parent,galton$child)) 
names(freqData) <- c("child","parent","freq")
freqData$parent<- as.numeric(as.character(freqData$parent))
freqData$child <- as.numeric(as.character(freqData$child))
g <- ggplot(filter(freqData, freq > 0), aes(x = parent, y = child))
g <- g   scale_size(range = c(2, 20), guide = "none")
g <- g   geom_point(color = "grey50", aes(size = freq))
g <- g   geom_point(aes(color=freq,size = freq))
g <- g   scale_color_gradient(low = "lightblue", high = "white")
g <- g   geom_smooth(method = lm, formula=y~x, se=FALSE, size = 3, color="red")
g

图7.添加回归线

基本概念

1. 经验均值
bar{X}

  1. 定义经验均值为
bar{X}=frac{1}{n}sumlimits_{i=1}^nX_i
  1. 样本数据点减去平均值会得到均值为0的数据，定义
tilde{X_i}=X_i-bar{X}

，则

tilde{X_i}

的均值为0。这个过程称为"居中"随机变量。

  1. 均值是使
sumlimits_{i=1}^n(X_i-?)^2

最小的最小二乘解

2. 经验标准差和方差

  1. 定义经验方差为
  2. 定义经验标准差为
S= sqrt{S^2}

，注意标准差与数据有相同单位

X_i/s

的经验标准差为1，这个过程称为"缩放"数据。

  1. 有时选择以分母
n

代替分母

n-1

，后者为无偏估计

3. 标准正态分布

  1. 定义标准正态分布
Z_i=frac{X_i-bar{X}}{s}

，经验均值为0，经验标准差为1。

  1. 将数据“居中”并“缩放”的过程称为“标准化“
4. 经验协方差

  1. 对于成对的数据
(X_i, Y_i)

，定义经验协方差为

  1. 同样，有时选择以分母
n

代替分母

n-1

，后者为无偏估计

5. 相关系数

  1. 定义相关系数，其中
S_x

和

S_y

分别是

X

观测值和

Y

观测值的标准差的估计值

  1. 相关系数
Cor(X,Y)=Cor(Y,X)
-1≤Cor(X,Y)≤1
  1. 当且仅当
X

或

Y

观测值分别恰好落在正斜率线或负斜率线时，

Cor(X，Y)= 1

和

Cor(X，Y)= -1
Cor(X,Y)

度量

X

和

Y

数据之间线性关系的强度

Cor(X，Y)= 0

表示没有线性关系

线性方程的普通最小二乘法（OLS）

回顾前面galton数据集中父母与孩子身高的例子

令

Y_i

为第

i

个孩子的身高，

X_i

为父母身高，线性回归

Y=?_0 X?_1

，最小二乘法要求

sumlimits_{i=1}^n(Y_i-(?_0 ?_1X_i))^2

最小。

  • 最优解为，
hat{?_0}=bar{Y}-hat{?_1}bar{X}

，回归线为

Y=hat{?_0} hat{?_1}X

，经过点

(bar{X},bar{Y})

。

  • 若已知
Y

预测

X

，此时回归线斜率为

  • 如果将数据居中，
(X_i-bar{X},Y_i-bar{Y})

，回归线斜率相同，并经过原点

  • 如果标准化数据，，斜率为
Cor(Y,X)
代码语言：javascript复制
y<-galton$child
x<-galton$parent 
beta1<-cor(y,x)* sd(y)/sd(x) 
beta0<-mean(y)-beta1*mean(x)
rbind(c(beta0,beta1),coef(lm(y~x)))
     (Intercept)      x
[1,]       23.94 0.6463
[2,]       23.94 0.6463

在R中检查计算，根据公式计算的斜率和截距与lm()函数拟合回归线得到的结果一样。

代码语言：javascript复制
beta1<-cor(y,x)* sd(x)/sd(y) 
beta0<-mean(x)-beta1*mean(y) 
rbind(c(beta0,beta1),coef(lm(x~y)))
     (Intercept)      y
[1,]       46.14 0.3256
[2,]       46.14 0.3256

用变量Y预测X，计算结果一致。

代码语言：javascript复制
yc<-y-mean(y) 
xc<-x-mean(x) 
beta1<-sum(yc*xc)/sum(xc^2) 
c(beta1,coef(lm(y~x))[2])
            x 
0.6463 0.6463

将数据居中，

(X_i-bar{X},Y_i-bar{Y})

，回归线斜率相同。

代码语言：javascript复制
yn<-(y-mean(y))/sd(y) 
xn<-(x-mean(x))/sd(x) 
c(cor(y,x),cor(yn,xn),coef(lm(yn~xn))[2])
                  xn 
0.4588 0.4588 0.4588

标准化数据，回归线斜率为相关系数

Cor(Y,X)

。

线性回归 编程算法

0 人点赞