R语言数学建模（１）：Regression analysis

某建模萌新半年来学习建模的一些记录加总结

一，linear regression

对于这样一个有明显线性趋势的数据,我们要找到一条直线使其具备预测该数据走向的能力，采用最小二乘准则（least square estimates）即，拟合出的直线符合残差平方和最小的准则。将上述翻译成数学语言即可得到如下模型

$min\sum _i{}(y_{i}-\hat{y}_{i})^{2}\\ s.t.\quad\hat{y}_{i}=ax_{i}+b$

这是一个凸优化问题，我们可以用高等数学的方法将约束条件代入目标函数得到一个无约束优化，再对a，b分别求偏导数，即可得到a与b的估计值，在统计学上可以证明LSE得到的是参数a与b的无偏估计即 $E(\hat{a}) = a \quad E(\hat{b}) = b$

此外，我们对于参数a与b的估计，我们还可以通过极大似然估计（Maxmum Likelihood Estimate）求得，其结果与LSE是一样的。

总离差平方和：SST

$\sum_{i}(y_{i}-\hat{y}_{i})^2$

回归平方和： SSR

$\sum_{i}(y_{i}-\bar{y}_{i})^2$

残差平方和： SSE

$\sum_{i}(y_{i}-\hat{y}_{i})^2$

可以证明三个平方和满足 $SST=SSE+SSR$

根据最小二乘准则，对于不同的模型，可以比较不同模型的残差平方和来选择最优拟合能力的模型

定义拟合优度 $R^2=1-SSE/SST$ SSE越小其值越大，但不会超过1。

supplyment: 在机器学习中，模型的误差，分为训练误差和泛化误差，这里的R方只能反应模型的训练误差，而训练误差小未必该模型就拥有优良的预测能力，即训练误差小不等于泛化误差小。当模型过拟合时，我们就认为其不具备预测数据的能力。具体参考李航的统计机器学习，本文不做详述。在对二维数据进行多项式回归时，使用的多项式次数过高往往会造成过拟合现象（可以想象成插值那种极端情形）。

二，R语言实现+ggplot可视化

lm函数不仅可以进行一元线性拟合，还可以进行多元线性和非线性拟合。

##输入数据集
x <- c(0.1, 0.11 ,0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.20, 0.21, 0.23)
y <- c(42 ,43 ,45 ,45 ,45 ,47.5 ,49 ,53 ,50 ,55 ,55 ,60)
##建立数据框
linear <- data.frame(x,y)
##拟合
model <- lm(y~x,data=linear)
##查看
summary(model)##输出结果Call:
lm(formula = y ~ x, data = linear)Residuals:Min       1Q   Median       3Q      Max 
-2.00449 -0.63600 -0.02401  0.71297  2.32451 Coefficients:Estimate Std. Error t value Pr(>|t|)    
(Intercept)   28.083      1.567   17.92 6.27e-09 ***
x            132.899      9.606   13.84 7.59e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 1.309 on 10 degrees of freedom
Multiple R-squared:  0.9503,    Adjusted R-squared:  0.9454 
F-statistic: 191.4 on 1 and 10 DF,  p-value: 7.585e-08

模型给出了各参数的估计值 Estimate y = 28.083+132.899x 可以看到p-value很小，可以认为我们的回归方程是显著的

我们可以用函数 confint(model)得到该模型在2.5%-97.5%范围的置信区间

confint(model)2.5 %    97.5 %
(Intercept)  24.59062  31.57455
x           111.49556 154.30337

可以用ggplot中的geom_smooth函数进行模型的可视化

geom_smooth(data,formula,method,se=T,colour,size)

data~数据集

formula~拟合的法则，可参照lm函数

method~拟合方法 loess：局部加权回归 lm:线性回归 glm:广义线性回归 gam：广义加性回归

colour~线条颜色

size~线条宽度

se~是否添加置信区间，默认为True

ggplot(linear,aes(x=x,y=y))+geom_point(shape=1,colour="black")+##geom_abline(intercept = 28.083,slope=132.899)geom_smooth(method = lm,se = F,colour = "steelblue",size=0.5)+theme(plot.title = element_text(hjust=0.5,size=15))+labs(title="Linear Regression")+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.border = element_blank(),axis.line =element_line(colour = "black"))+annotate("text",x=0.2,y=58,label="y=132.9x+28.08")+annotate("text",x=0.2,y=59,parse=TRUE,  label = "atop(R^2==0.9503)",size=4)

若要拟合的式子为 y=ax^2+b 可用如下代码

x <- c(0.1, 0.11 ,0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.20, 0.21, 0.23)
y <- c(42 ,43 ,45 ,45 ,45 ,47.5 ,49 ,53 ,50 ,55 ,55 ,60)
##建立数据框
unlinear <- data.frame(x,y)
model <- lm(y~I(x^2),unlinear)
summary(model)Call:
lm(formula = y ~ I(x^2), data = linear)Residuals:Min       1Q   Median       3Q      Max 
-1.46660 -0.59878 -0.07453  0.32904  2.95051 Coefficients:Estimate Std. Error t value Pr(>|t|)    
(Intercept)  38.3482     0.8144   47.09  4.5e-13 ***
I(x^2)      404.8879    27.5137   14.72  4.2e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 1.234 on 10 degrees of freedom
Multiple R-squared:  0.9559,	Adjusted R-squared:  0.9514 
F-statistic: 216.6 on 1 and 10 DF,  p-value: 4.201e-08

可以看到此时的 multiple R-squared = 0.9559 竟然比线性模型的0.9503还要高，此时

y = 38.3482+404.8879x^2