用R语言进行简单线性回归分析

数模天下

用R语言进行简单线性回归分析，数据出自何晓群--应用回归分析，语言如下所示：

4 I! G; i1 T; y0 K" p

x y

) E0 v/ U0 z7 k5 w

3.4 26.2

1.8 17.8

! s$ Z( @* v" U6 I. o$ H  m0 s

4.6 31.3

2.3 23.1

3.1 27.5

5.5 36

0.7 14.1

3 22.3

2.6 19.6

4.3 31.3

2.1 24

  }( x# [4 K  U" W% g

1.1 17.3

6.1 43.2

4.8 36.4

+ Q) A4 V  W* `& ?. ]4 W$ f3 b

3.8 26.1

9 U( f6 z* s( Q) U1 [& u3 _/ R

#-------------------------------------------------------------#数据准备

fire <- read.table('D:/fire.txt', head = T)

#-------------------------------------------------------------#回归分析

plot(fire$y ~ fire$x)

fire.reg <- lm(fire$y ~ fire$x, data = fire) #回归拟合

summary(fire.reg) #回归分析表

anova(fire.reg) #方差分析表

abline(fire.reg, col = 2, lty = 2) #拟合直线

#-------------------------------------------------------------#残差分析

fire.res <- residuals(fire.reg) #残差

fire.sre <- rstandard(fire.reg) #学生化残差

* d2 K+ ?& r* n- \; `

plot(fire.sre)

abline(h = 0)

5 ~$ @2 }$ z: o

text(11, fire.sre[11], label = 11, adj = (-0.3), col = 2) #标注点

#-------------------------------------------------------------#预测与控制

attach(fire) #连接

fire.reg <- lm(y ~ x) #这种回归拟合简单

fire.points <- data.frame(x = c(3.5, 4))

fire.pred <- predict(fire.reg, fire.points, interval = 'prediction', level = 0.95) #预测：置信区间

fire.pred

detach(fire) #取消连接

& e9 o2 i: F" N! O9 ]

--------------------------------------------------------------------------------------------------

#附自编的过程程序：（R最大的好处是可以自己编想要的程序和函数，尤其没有内置函数的时候）

8 e- @+ J% ^' Q6 ~, S1 V9 H; ^

fire <- read.table('D:/fire.txt', head = T)

3 U5 t* I# i$ H

attach(fire)

--------------------------------------------

1 O# |/ Y9 w- |; K3 @2 v8 I: f

lxy <- function(x){

sum <- 0

; V. E; c, P. B, W# G7 o0 e

sum0 <- 0

for(i in 1:length(x)){

, Q( D4 O" t( n7 R- U0 y( V

sum0 <- (x - mean(x)) * (y-mean(y))

sum <- sum + sum0}

sum}

---------------------------------------------------------------------------------

#用这个就不需要循环了

3 L9 F" N4 R3 l) ]' j

lxy <- function(x){

mid <- (x - mean(x)) * (y-mean(y))

* V9 y+ E0 x  G5 [1 j8 d  U  a

sum <- sum(mid)

5 L9 l' {' [9 ~5 ?* X6 Y; R8 N, d

sum}

#对于数据框、列表等数据对象要善用apply()函数。

---------------------------------------------------------------------------------

lxx <- function(x){

sum <- 0

sum0 <- 0

; V  d% C2 E. F# H6 p5 q, _2 h0 g- P

for(i in 1:length(x)){

sum0 <- (x - mean(x))^2

sum <- sum + sum0}

1 X8 j. m8 I. {5 a( M

sum}

Lxx <- lxx(x)

Lyy <- lxx(y)

Lxy <- lxy(x)

b1 <- Lxy / Lxx; b1 #回归系数斜率

b0 <- mean(y) - b1 * mean(x); b0 #回归系数截距

residu <- y - (b0 + b1*x); residu #残差

r <- Lxy / sqrt(Lxx * Lyy); r #相关系数

rsqure <- r^2; rsqure #决定系数

adrsqure <- 1 - ((length(x)-1)/(length(x)-2))*(1-r^2) #调整后的决定系数

----------------------------------------------------------------------------------

6 [0 T* h2 U6 Z0 u, f- a! b

esrequre <- function(x){ #求标准差平方估计值

4 h- q2 j8 r1 q7 b7 l

sum <- 0

sum0 <- 0

for(i in 1:length(x)){

sum0 <- residu^2

sum <- sum + sum0}

4 j+ Y! f7 R( l4 y  b

residusqure <- sum/(length(x)-2)

residusqure}

esterreq <- esrequre(x); esterreq #标准差平方估计值(MSE)

ester <- sqrt(esrequre(x)); ester #标准差估计值(回归分析表给出的标准误差)

val_t <- b1*sqrt(Lxx) / ester; val_t #检验回归系数斜率b1的t值

SSe <- function(x){ #求残差平方和

sum <- 0

2 n# j* m, V* F/ a7 g# \/ T

sum0 <- 0

for(i in 1:length(x)){

5 i6 m5 j* u' n6 `8 _2 t! N

sum0 <- residu^2

sum <- sum + sum0}

* h; N7 u: l. G" ]

sum}

SSE <- SSe(x); SSE #残差平方和

. y/ g2 h$ Z/ W" b- j1 X3 u: n

MSE <- SSE/(length(x)-2); MSE #残差均方和

SSr <- function(x){

sum <- 0

9 t  }/ ?2 |, v- [- b# ~9 H/ C

sum0 <- 0

for(i in 1:length(x)){

sum0 <- ((b0 + b1*x) - mean(y))^2

sum <- sum + sum0}

sum}

SSR <- SSr(x); SSR #回归平方和

MSR <- SSR/1; MSR #回归均方和

val_F <- SSR / MSE; val_F #检验回归方程F值

hi <- 1/length(x) + (x-mean(x))^2/Lxx #杠杆值

ZRE <- residu / ester; ZRE #标准化残差

SRE <- residu/(ester*sqrt(1-hi)); SRE #学生化残差

  b+ @0 B0 ]' @( v

Y <- function(x){b0 + b1 * x} #点估计

Y(3.5)

& i  q$ c9 K& n9 U
0 I# F/ ~5 y2 V+ _