```
#bibliography: [../References/BILab.bib, packages.bib]
```
# 回归分析
## 一元回归
自己构建一套回归的方法
```{r}
set.seed(1)
x=1:100
y=x+rnorm(100,sd=5)
target_fun=function(a0=1,b0=1){
y_est=a0+b0*x
target=sum(abs(y-y_est))
return(target)
}
#暴力破解法
re=c()
as=c()
bs=c()
i=1
for(ai in 0:10){
for(bi in 0:10){
re[i]=target_fun(ai,bi)
as=c(as,ai)
bs=c(bs,bi)
i=i+1
}
}
mini=which(re==min(re))
a_best=as[mini]
b_best=bs[mini]
plot(x=x,y=y)
abline(a=a_best, b=b_best)
```
```{r}
x=1:100
y=1.2345*x+rnorm(100,mean=3.1415,sd=1)
#反向梯度求解
target_fun2=function(par){
target_fun(par[1],par[2])
}
#这里默认是采用牛顿梯度法
re2=optim(c(1,1),target_fun2)
plot(x=x,y=y)
abline(a=re2$par[1], b=re2$par[2])
```
内在的统计框架
$$y=a+bx+\epsilon$$
1. x测量没有任何误差
2. $\epsilon$是一个符合正态分布随机变量,均值为0,标准差为$\delta$。
2. y一个符合正态分布的随机变量,其中均值为$a+bx$,标准差为$\delta$
R中还内置了一个非常方便的函数来求解回归系数。
```{r}
md1=lm(y~x)
md1
```
```{r}
plot(x=x,y=y)
abline(md1)
```
```{r}
plot(x=x,y=y)
abline(md1)
y_pred=predict(md1)
points(x=x,y=y_pred,col=2,pch=19)
```
```{r}
delta=sd(y-y_pred)
delta
```
参数的置信区间
```{r}
confint(md1)
```
参数与零是否有显著差异
```{r}
summary(md1)
```
### 模型的选择
针对最简单的回归模型,我们可以有三个备选模型
$$y=a+bx+\epsilon$$
$$y=bx+\epsilon$$
$$y=a+\epsilon$$
到底选哪个模型?
```{r}
y=x+rnorm(100,sd=10)
md1=lm(y~x)
md1
```
```{r}
summary(md1)
```
```{r}
md2=lm(y~x-1)
md2
```
```{r}
md3=lm(y~1)
md3
```
```{r}
plot(x=x,y=y)
abline(md1)
abline(md2,col=2)
abline(md3,col=3)
```
### 模型评价
解释度
```{r}
x=1:100
y=5.5*x+rnorm(100,sd=20,mean=10)
md1=lm(y~x)
md2=lm(y~1)
plot(x=x,y=y)
abline(md1)
abline(md2,col=2)
```
残差和
```{r}
var(predict(md1)-y)
var(predict(md2)-y)
```
```{r}
sm1=summary(md1)
sm2=summary(md2)
```
诊断模型
1. 残差是否符合正态
```{r}
shapiro.test(resid(md1))
```
2. 方差是否齐性
```{r}
y2=x*1.3+rnorm(100,sd=2*x)
#par(mfrow=c(2,1))
plot(x=x,y=y)
plot(x=x,y=y2)
```
```{r}
plot(md1)
```
```{r}
plot(lm(y2~x))
```
```{r}
y3=y
y3[3]=1000
plot(lm(y3~x))
```
## 多元回归
## 结构方程
## 作业
‘treegrowth.txt’文件中包含了一个样地内所有树木,两次的胸径测量值。请你建立一个模型,用第一次的胸径来预测解释第二次的胸径大小。