线性回归in R

回归分析

有的变量间有一定的关系,但由于情况错综复杂,无法精确研究,或由于存在不可避免的误差等原因,以致无法使用函数形式来表示。为研究这类变量之间的关系就需要大量试验或观测获得数据,用统计方法寻找它们之间的关系,这种关系反应了变量之间的统计规律。使用统计方法寻找统计规律便是回归分析

在回归分析中研究的主要问题是:

1)确定 Y Y Y X 1 , X 2 , ⋯   , X p X_1,X_2,\cdots,X_p X1,X2,,Xp的定量关系表达式。这种表达式成为回归方程。

2)对求得的回归方程的可信度进行检验(是否有线性关系)

3)判断自变量 X j ( j = 1 , 2 , ⋯   , p ) 对 X_j(j=1,2,\cdots,p)对 Xj(j=1,2,,p) Y Y Y是否有影响(是否显著)

4)利用所求得的回归方程进行检验和控制。


一元线性回归

Y = β 0 + β 1 X + ε Y=\beta_0+\beta_1X+\varepsilon Y=β0+β1X+ε

其中, β 0 + β 1 X \beta_0+\beta_1X β0+β1X表示 Y Y Y X X X的变化而线性变化的部分, ε \varepsilon ε式随机误差,代表所有误差的总和,其值不可观测。一般假设 ε ∼ N ( 0 , σ 2 ) \varepsilon\sim N(0,\sigma^2) εN(0,σ2) 。称 f ( X ) = β 0 + β 1 X 为 一 元 线 性 回 归 函 数 , f(X)=\beta_0+\beta_1X为一元线性回归函数, f(X)=β0+β1X线 β 0 \beta_0 β0为回归常数, β 1 \beta_1 β1为回归系数。

一元线性回归的模型可以表示为:
y i = β 0 + β 1 x i + ε i i = 1 , 2 , ⋯   , n . y_i=\beta_0+\beta_1x_i+\varepsilon_i \quad i=1,2,\cdots,n. yi=β0+β1xi+εii=1,2,,n.

回归参数的估计

最小二乘法估计:
Q ( β 0 , β 1 ) = ∑ i = 1 n ( y i − β 0 − β 1 x i ) 2 Q(\beta_0,\beta_1)=\sum_{i=1}^{n}(y_i-\beta_0-\beta_1x_i)^2 Q(β0,β1)=i=1n(yiβ0β1xi)2
关于 Q Q Q分别求 β 0 \beta_0 β0 β 1 \beta_1 β1的偏导:
β 1 ^ = ∑ i = 1 n ( x 1 − x ˉ ) ( y i − y ˉ ) ∑ i = 1 n ( x i − x ˉ ) 2 = S x y S x x , β 0 ^ = y ˉ − β 1 ^ x ˉ \hat{\beta_1}=\frac{\sum_{i=1}^{n}(x_1-\bar{x})(y_i-\bar y)} {\sum_{i=1}^n(x_i-\bar x)^2}=\frac{S_{xy}} {S_{xx}},\quad \hat{\beta_0}=\bar{y}-\hat{\beta_1}\bar x β1^=i=1n(xixˉ)2i=1n(x1xˉ)(yiyˉ)=SxxSxy,β0^=yˉβ1^xˉ
通常取
σ ^ 2 = ∑ i = 1 n ( y i − β 0 − β 1 x i ) 2 n − 2 \hat{\sigma}^2=\frac{\sum_{i=1}^{n}(y_i-\beta_0-\beta_1x_i)^2 }{n-2} σ^2=n2i=1n(yiβ0β1xi)2
为参数 σ 2 \sigma^2 σ2的估计量,也称为 σ 2 \sigma^2 σ2的最小二乘估计。

关于 β 0 \beta_0 β0 β 1 \beta_1 β1估计的方差为
V a r ( β 0 ) = σ 2 ( 1 n + x ˉ 2 S x x ) , V a r ( β 1 ) = σ 2 S x x Var(\beta_0)=\sigma^2(\frac{1}{n}+\frac{\bar{x}^2}{S_{xx}}),\quad Var(\beta_1)=\frac{\sigma^2}{S_{xx}} Var(β0)=σ2(n1+Sxxxˉ2),Var(β1)=Sxxσ2
如果 σ 2 \sigma^2 σ2未知,用 σ ^ \hat{\sigma} σ^替换 σ \sigma σ,得到
s d ( β 0 ^ ) = σ ^ 1 n + x ˉ 2 S x x , s d ( β 1 ^ ) = σ ^ S x x sd(\hat{\beta_0}) = \hat{\sigma}\sqrt{\frac{1}{n}+\frac{\bar{x}^2}{S_{xx}}},\quad sd(\hat{\beta_1})=\frac{\hat{\sigma}}{\sqrt{S_{xx}}} sd(β0^)=σ^n1+Sxxxˉ2 ,sd(β1^)=Sxx σ^

回归方程的显著性检验

​ 使用最小二乘法并不一定要知道 Y Y Y X X X是否有线性相关的关系。如果不存在线性相关的关系,得出的模型,就毫无意义。因此,需要对回归方程进行检验。

​ 从统计上讲, β 1 \beta_1 β1 E ( Y ) E(Y) E(Y) X X X线性变化的变化率,若 β 1 = 0 \beta_1=0 β1=0,则 E ( Y ) E(Y) E(Y)不随着 X X X线性变化。

假设检验为:
H 0 : β 1 = 0 , H 1 : β 1 ≠ 0 H_0:\beta_1=0,\quad H_1:\beta_1\neq0 H0:β1=0,H1:β1=0
1. t t t检验法。当 H 0 H_0 H0成立时,统计量
T = β ^ 1 s d ( β ^ 1 ) = β 1 ^ S x x σ ^ ∼ t ( n − 2 ) T=\frac{\hat{\beta}_1}{sd(\hat{\beta}_1)}=\frac{\hat{\beta_1}\sqrt{S_{xx}}}{\hat{\sigma}}\sim t(n-2) T=sd(β^1)β^1=σ^β1^Sxx t(n2)
2. F F F检验。当 H 0 H_0 H0成立时,统计量
F = β 1 ^ 2 S x x σ ^ 2 ∼ F ( 1 , n − 2 ) F=\frac{\hat{\beta_1}^2S_{xx}} {\hat{\sigma}^2}\sim F(1,n-2) F=σ^2β1^2SxxF(1,n2)
3.相关系数检验法。
R = S x y S x x S y y R=\frac{S_{xy}}{\sqrt{S_{xx}S_{yy}}} R=SxxSyy Sxy
称为样本相关系数,对于给定的显著性水平 α \alpha α,拒绝域为
∣ R ∣ > r α ( n − 2 ) |R|>r_{\alpha}(n-2) R>rα(n2)


【例题】利用R语言中的lm()建立一元线性回归方程

image-20210827150439799
> x <- c(0.10,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.0,43.5,45.0,45.5,45.0,47.5,
+        49.0,53.0,50.0,55.0,55.0,60.0)
> linemod <- lm(y ~ x)
> summary(linemod)

模型summary

Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.0431 -0.7056  0.1694  0.6633  2.2653 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   28.493      1.580   18.04 5.88e-09 ***
x            130.835      9.683   13.51 9.50e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.319 on 10 degrees of freedom
Multiple R-squared:  0.9481,	Adjusted R-squared:  0.9429 
F-statistic: 182.6 on 1 and 10 DF,  p-value: 9.505e-08
  • Std. Error为 s d ( β 1 ^ ) sd(\hat{\beta_1}) sd(β1^),是标准差,不是标准误

  • t t t value 为 Estimate/Std. Error, 18.04 = 28.493/1.580

  • F F F value为 Estimate^2/Std. Error^2. 182.6 = 130.835^2 / 9.683^2

  • Multiple R-squared为 R 2 = S x y 2 S x x S y y R^2=\frac{S_{xy}^2}{S_{xx}S_{yy}} R2=SxxSyySxy2,

  • Residual standard error 为残差的标准差,即 σ ^ \hat{\sigma} σ^

参数 β 0 \beta_0 β0 β 1 \beta_1 β1的区间估计

​ 在得到 β 0 ^ 、 β 1 ^ \hat{\beta_0}、\hat{\beta_1} β0^β1^之后,有时还需要它们的区间估计。
T i = β i ^ − β i s d ( β ^ i ) ∼ t ( n − 2 ) T_i = \frac{\hat{\beta_i}-\beta_i}{sd(\hat{\beta}_i)}\sim t(n-2) Ti=sd(β^i)βi^βit(n2)
对给定的置信水平1- α \alpha α,有
p { β i ^ − s d ( β i ^ ) t α / 2 ( n − 2 ) , β i ^ + s d ( β i ^ ) t α / 2 ( n − 2 ) } p\left\{\hat{\beta_i}-sd(\hat{\beta_i})t_{\alpha/2}(n-2) ,\hat{\beta_i}+sd(\hat{\beta_i})t_{\alpha/2}(n-2) \right\} p{βi^sd(βi^)tα/2(n2),βi^+sd(βi^)tα/2(n2)}
【代码实现】

beta.int <- function(fm,alpha=0.05){
  A <- summary(fm)$coefficients
  df <- fm$linemod$df.residual
  left <- A[,1]-A[,2]*qt(1-alpha/2,df)
  right <- A[,1]+A[,2]*qt(1-alpha/2,df)
  rowname <- dimnames(A)[[1]]
  colname <- c('Estimate','Left','Right')
  matrix(c(A[,1],left,right),ncol=3,dimnames = list(rowname,colname))
  
}
> beta.int(linemod)
             Estimate      Left     Right
(Intercept)  28.49282  24.97279  32.01285
x           130.83483 109.25892 152.41074

所以线性模型linemod的置信区间为[24.97279 32.01285]和[109.25892 152.41074]

预测

对与 X = x 0 , Y = y 0 X=x_0,Y=y_0 X=x0,Y=y0的置信度为 1 − α 1-\alpha 1α的预测区间为
[ y ^ 0 − l , y ^ 0 + l ] [\hat{y}_0-l,\hat{y}_0+l] [y^0l,y^0+l]
其中
l = t α / 2 ( n − 2 ) σ ^ 1 + 1 n + ( x ˉ − x 0 ) 2 S x x l = t_{\alpha/2}(n-2)\hat{\sigma}\sqrt{1+\frac{1}{n}+\frac{(\bar{x}-x_0 )^2}{S_{xx}}} l=tα/2(n2)σ^1+n1+Sxx(xˉx0)2
当样本容量很大时,根式中近似等于1,因此,置信区间近似等于[ y 0 − σ ^ Z α / 2 , y 0 + σ ^ Z α / 2 y_0-\hat{\sigma}Z_{\alpha/2},y_0+\hat{\sigma}Z_{\alpha/2} y0σ^Zα/2,y0+σ^Zα/2]


【代码】

predict(model,newx,interval='prediction')

回归诊断(Regression diagnostics)

​ 回归诊断的主要内容有:

  • 关于误差项是否满足独立性、等方差性、正态性

  • 选择线性模型是否合适

  • 是否存在异常样本

  • 回归分析的结果是否对某些样本的依赖过重?也就是说,回归模型是否具备稳定性

  • 自变量之间是否存在高度相关?是否存在多重共线性问题

运用R语言进行贝叶斯多元线性回归可以使用runjags软件包来实现。runjags是一个R语言包,它使用JAGS(Just Another Gibbs Sampler)引擎进行贝叶斯分析。 以下是一个简单的示例代码,展示如何使用runjags软件包进行贝叶斯多元线性回归: ``` library(runjags) # 设定数据 x1 <- rnorm(100) x2 <- rnorm(100) y <- 2*x1 + 3*x2 + rnorm(100) # 构建模型 model_string <- " model { # 数据 for (i in 1:n) { y[i] ~ dnorm(mu[i], tau) mu[i] <- b[1] + b[2]*x1[i] + b[3]*x2[i] } # 先验 for (j in 1:3) { b[j] ~ dnorm(0, 0.001) } tau ~ dgamma(0.001, 0.001) }" # 运行模型 n <- length(y) data_list <- list("x1"=x1, "x2"=x2, "y"=y, "n"=n) params <- c("b", "tau") model <- jags.model(textConnection(model_string), data=data_list, n.chains=3) update(model, 10000) samples <- coda.samples(model, params, n.iter=50000) # 结果展示 summary(samples) plot(samples) ``` 在上面的代码中,我们首先生成一些随机数据。然后,我们构建了一个模型字符串,该模型在数据的基础上定义了一个多元线性回归模型。随后,我们使用runjags包中的函数来运行这个模型,得到了参数的后验分布。最后,我们展示了参数的摘要和一些可视化结果。 需要注意的是,这只是一个简单的示例,实际应用中可能需要更复杂的模型和更多的参数。同时,贝叶斯分析需要对先验分布进行合理的选择,这对于初学者来说可能是一个挑战。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值