在线时间:8:00-16:00
迪恩网络APP
随时随地掌握行业动态
扫描二维码
关注迪恩网络微信公众号
用R语言进行简单线性回归分析。数据出自何晓群--应用回归分析,语言例如以下所看到的: x y 3.4 26.2 1.8 17.8 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 1.1 17.3 6.1 43.2 4.8 36.4 3.8 26.1 #-------------------------------------------------------------#数据准备 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) #学生化残差 plot(fire.sre) abline(h = 0) 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) #取消连接 -------------------------------------------------------------------------------------------------- #附自编的过程程序:(R最大的优点是能够自己编想要的程序和函数,尤其没有内置函数的时候) fire <- read.table('D:/fire.txt', head = T) attach(fire) -------------------------------------------- lxy <- function(x){ sum <- 0 sum0 <- 0 for(i in 1:length(x)){ sum0 <- (x[i] - mean(x)) * (y[i]-mean(y)) sum <- sum + sum0} sum} --------------------------------------------------------------------------------- #用这个就不须要循环了 lxy <- function(x){ mid <- (x - mean(x)) * (y-mean(y)) sum <- sum(mid) sum} #对于数据框、列表等数据对象要善用apply()函数。 --------------------------------------------------------------------------------- lxx <- function(x){ sum <- 0 sum0 <- 0 for(i in 1:length(x)){ sum0 <- (x[i] - mean(x))^2 sum <- sum + sum0} 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) #调整后的决定系数 ---------------------------------------------------------------------------------- esrequre <- function(x){ #求标准差平方预计值 sum <- 0 sum0 <- 0 for(i in 1:length(x)){ sum0 <- residu[i]^2 sum <- sum + sum0} 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 sum0 <- 0 for(i in 1:length(x)){ sum0 <- residu[i]^2 sum <- sum + sum0} sum} SSE <- SSe(x); SSE #残差平方和 MSE <- SSE/(length(x)-2); MSE #残差均方和 SSr <- function(x){ sum <- 0 sum0 <- 0 for(i in 1:length(x)){ sum0 <- ((b0 + b1*x[i]) - mean(y))^2 sum <- sum + sum0} sum} s-s-r <- SSr(x); s-s-r #回归平方和 MSR <- s-s-r/1; MSR #回归均方和 val_F <- s-s-r / MSE; val_F #检验回归方程F值 hi <- 1/length(x) + (x-mean(x))^2/Lxx #杠杆值 ZRE <- residu / ester; ZRE #标准化残差 SRE <- residu/(ester*sqrt(1-hi)); SRE #学生化残差 Y <- function(x){b0 + b1 * x} #点预计 Y(3.5) |
请发表评论