1.1介绍
将数据以若干等分方式进行现象回归建模，如此可独立观察每一等分模型的趋势，并且与一般线性回归模型做比较。
假设y依赖于x的一个连续性因变量，建模后的标准线性回归模型表示为yi=b0+b1xi+ei
ei服从均值为0的正态分布。
1.2回归分析
回归分析是确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法。
1.3步骤
- 数据用excel整理成“.csv”文件，在计算机的“C:\”中新建目录，命名为“data”， 将整理的“.csv”复制到该目录中，并改名为“testdata.csv”
- 登录http://eplrm.byethost31.com
【注意，这里用的潘文超老师的做法噢！我仅仅是个搬运工，然后对代码进行了一些整理而已！】
1.4代码（这部分是不用改的，直接用就好了）
#使用前先运行前面粗体字行代码（从标注start开始到end结束）
#=====================start====================
install.packages("SparseM")
install.packages("quantreg")
install.packages("boot")
install.packages("ggplot2")
install.packages("grid")
install.packages("Rmisc")
library(SparseM)
library(quantreg)
library(boot)
library(ggplot2)
library(grid)
library("Rmisc")
#定义一个函数equal_lm来处理回归分析并解析数据
#四个参数，第一个formulation是方程式，第二个datak是数据，
#第三个是modelname模型名字，第四个C是置信度,K表示变量个数（y和x一起）
equal_lm<-function(formulation,datak,modelname,c,k){
variable<-c("constant","x1","x2","x3","x4","x5","x6","x7","x8","x9","x10")
names<-rep(0,k)
for (i in 1:k){
names[i]<-variable[i]
}
LRM=lm(formulation,data=datak)
output=summary(LRM)
confident=confint(LRM,level=c) #置信度区间
beta<-output$coefficients #参数估计
R_saure<-output$r.squared #R方的值
#将系数存储到变量中
conf<-rep(0,length(names))
std.error<-rep(0,length(names))
t_value<-rep(0,length(names))
sig<-rep(0,length(names))
low_confident<-rep(0,length(names))
up_confidient<-rep(0,length(names))
r_squre<-c(modelname,rep("-",(length(names)-1)))
for(i in 1:length(beta[,1])){
conf[i]<-beta[i,1]
std.error[i]<-beta[i,2]
t_value[i]<-beta[i,3]
sig[i]<-beta[i,4]
low_confident[i]<-confident[i,1]
up_confidient[i]<-confident[i,2]
r_squre[2]<-R_saure
}
result<-data.frame(names,r_squre,conf,std.error,t_value,sig,low_confident,up_confidient)
return(result)
}
#=====================end====================
#2.DF_model将数据等分划分,并计算最终的结果
#formulation 为回归模型方程，alldata为所有数据，c为置信度，k为模型变量个数
#alpha为分成了几等分
DF_model<-function(formulation,alldata,c,k,alpha){
n=length(alldata[,1]) #数据的样本个数
m=round(n/alpha) #表示每一等分数据的个数
#m0为原始模型（未进行等分处理时的模型拟合结果）
result<-equal_lm(formulation,alldata,"model0",c,k)
for(i in 1:alpha){
if(i<alpha){
temp1=m*(i-1)+1
temp2=m*i
}
if(i==alpha){
temp1=m*(i-1)+1
temp2=n
}
#拿到每一次等分的数据
datak=alldata[temp1:temp2,]
modelname<-paste("model",i)
result_temp<-equal_lm(formulation,datak,modelname,c,k)
result<-rbind(result,result_temp)
}
return(result)
}
#3,定义一个函数PY_data将数据平移划分,并计算模型结果
#formulation 为回归模型方程，alldata为所有数据，c为置信度，k为模型变量个数
#alpha为移动数据样本的个数，比如一次为5个数据，移动一步，alpha=5
py_model<-function(formulation,alldata,c,k,alpha){
n=length(alldata[,1]) #数据的样本个数, alpha为每一份样本的数据个数
move<-n-alpha+1 #表示样本为alpha，一共需要移动的步数
#m0为最开始的所有数据的模型
m0<-equal_lm(formulation,data,"model0",c,k)
for (i in 1:move){
temp1=i
temp2=i+alpha-1
datak=alldata[temp1:temp2,]
modelname<-paste("model",i)
mk<-equal_lm(formulation,datak,modelname,c,k)
m0<-rbind(m0,mk)
}
return(m0)
}
#===============end================
#4,画图置信区间
plot_conf<-function(k,mn,variable){
names<-rep(0,k)
for (i in 1:k){
names[i]<-variable[i]
}
n=length(mn[,1])
delta<-round(n/k)
up_beta_x=rep(0,delta)
low_beta_x=rep(0,delta)
estima_beta_x=rep(0,delta)
step=seq(1,delta,1)
p<-NULL
for (i in 1:k){
for(j in 1:delta){
temp=(j-1)*k+i
up_beta_x[j]=mn[temp,8]
low_beta_x[j]=mn[temp,7]
estima_beta_x[j]=mn[temp,3]
}
name=names[i]
plotdata0<-data.frame(up_beta_x,low_beta_x,estima_beta_x,step)
len<-length(plotdata0[,1])
plotdata<-plotdata0[2:len,]
layer1=ggplot(plotdata, aes(x = step, y=estima_beta_x)) +
geom_ribbon(aes(ymin = low_beta_x, ymax=up_beta_x), alpha = 0.2) +
geom_line(linetype = 1, size =1.5, colour = "black")+geom_point(size = 1.8, shape = 21,fill = "white") # 折线图函数
# 带状图函数：ymin设置下界，ymax设置上界；
layer2=geom_hline(size=1.0,aes(yintercept=mn[i,3]),linetype=4,colour="red")
layer3=geom_hline(size=0.7,aes(yintercept=mn[i,7]),linetype=2,colour="black")
layer4=geom_hline(size=0.7,aes(yintercept=mn[i,8]),linetype=2,colour="black")
plots=layer1+layer2+layer3+layer4+ xlab("step") + ylab(name)
p[[i]]=plots
}
multiplot( plotlist=p ,cols=2)
}
#========end=========
#==样本的方差检验函数F_test==
#与DF_model过程很像
#variable<-c("GDP","X1","X2")，data为所有数据，alpha为几等分
F_test<-function(variable,data,alpha){
k<-length(variable)
data_compare<-rep("-",k)
f_value<-rep(0,k)
p_value<-rep(0,k)
result<-data.frame(data_compare,variable,f_value,p_value)
#alpha=3
alldata=data
n=length(alldata[,1]) #数据的样本个数
m=round(n/alpha) #表示每一等分数据的个数
#m0为原始模型（未进行等分处理时的模型拟合结果）
temp=alpha-1
for(i in 1:temp){
temp11=m*(i-1)+1
temp12=m*i
#拿到每一次等分的数据
datak1=alldata[temp11:temp12,]
temp2=i+1
for(j in temp2:alpha){
if(i<alpha){
temp21=m*(j-1)+1
temp22=m*j
}
if(i==alpha){
temp21=m*(j-1)+1
temp22=n
}
datak2=alldata[temp21:temp22,]
modelname<-paste("model",i,"and",j)
data_compare[1]<-modelname
for(line in 1:k){
linedata1<-datak1[,line]
linedata2<-datak2[,line]
varout<-var.test(linedata1,linedata2)
f_value[line]<-varout$statistic
p_value[line]<-varout$p.value
}
result_temp<-data.frame(data_compare,variable,f_value,p_value)
result<-rbind(result,result_temp)
}
}
len_res<-length(result[,1])
len_str<-k+1
Result<-result[len_str:len_res,]
return(Result)
}
#============平移模型系数F检验==========
#三个参数，第一个平移结果，alpah为几等分，varivable为变量
py_xishu_f_test<-function(mode_pingyi_result,alpha,variavle){
data<-mode_pingyi_result
k=length(variavle) #变量的个数
n=length(mode_pingyi_result[,1]) #数据的长度
m<-round(n/k)-1 #一共有多少个模型
data_compare<-rep("-",k)
f_value<-rep(0,k)
p_value<-rep(0,k)
delta<-round(m/alpha) #每一等分的数量
result<-data.frame(data_compare,variavle,f_value,p_value)
temp1=alpha-1
for(i in 1:temp1){
temp2=i+1
line11<-((i-1)*delta+1)*k+1
line12<-(i*delta+1)*k
for(j in temp2:alpha){
line21<-((j-1)*delta+1)*k+1
line22<-(j*delta+1)*k
modename<-paste("model",i,"and",j)
datak1=data[line11:line12,3]
datak2=data[line21:line22,3]
data_compare[2]=modename
for(beta in 1:k){
index<-seq(1,m,k)+beta-1
varout<-var.test(datak1[index],datak2[index])
f_value[beta]<-varout$statistic
p_value[beta]<-varout$p.value
}
result_temp<-data.frame(data_compare,variavle,f_value,p_value)
result<-rbind(result,result_temp)
}
}
len_res<-length(result[,1])
len_str<-k+1
Result<-result[len_str:len_res,]
return(Result)}
#=======end======
【注意：以上代码直接用！！！！】
1.5代码（这部分是要根据自己的数据进行修改的）
#从这里开始依序执行等分线性回归分析
#====從這裡開始依序執行等分線性回歸分析================
#==1.读取数据==
#设置工作目录括号内参数为存放数据的目录，并读取数据
setwd("C:/data")
data=read.csv("testdata.csv")
#====2.简单线性回归====
#对所有数据
#c=0.95表示置信度为0.95，k=10表示一共10个变量（9个X变量+1个常数项）
data=read.csv("testdata.csv")
model0_result1<equal_lm(Y~X1+X2+X3+X4+X5+X6+X7+X8+X9,datak=data,modelname="model0",c=0.95,k=10)
model0_result1
#======3.等分回归======
#等分回归，分3等分
#c=0.90表示置信度，k=10表示（9个X变量+1个常数项）alpha=3表示3等分
data=read.csv("testdata.csv")
model_dengfen_result<-DF_model(Y~X1+X2+X3+X4+X5+X6+X7+X8+X9,alldata=data,c=0.90,k=10,alpha = 3)
model_dengfen_result
#保存结果到csv文件，mode_pingyi_result为变量名称，后面参数为存为文件名称
write.csv(model_dengfen_result,"dengfeng.csv")
variable<-c("Y","X1","X2","X3","X4","X5","X6","X7","X8","X9")
plot_conf(k=10, mn=model_dengfen_result,variable)
#=========4.平移模型======
# c=0.95表示置信度为0.95，k=10表示9个变量和1个常数项
data=read.csv("testdata.csv")
mode_pingyi_result<-py_model(Y~X1+X2+X3+X4+X5+X6+X7+X8+X9,data,c=0.90,k=10,alpha = nrow(data)/3)
mode_pingyi_result #保存结果到csv文件
write.csv(mode_pingyi_result,"pingyimodel.csv")
variable<-c("Y","X1","X2","X3","X4","X5","X6","X7","X8","X9")
plot_conf(k=10, mn=mode_pingyi_result,variable)
#========5.使用F检验测试样本数据=====
#两个参数，一个是data，一个是变量variable，一个是等分数量
data=read.csv("testdata.csv")
variable<-c("Y","X1","X2","X3","X4","X5","X6","X7","X8","X9")
F_test_result<-F_test(variable,data,3)
F_test_result
write.csv(F_test_result,"data_F_test.csv")
#===6.平移模型回归系数的f检验 ===
data=read.csv("testdata.csv")
mode_pingyi_result<-py_model(Y~X1+X2+X3+X4+X5+X6+X7+X8+X9,data,c=0.95,k=10, alpha= nrow(data)/3)
mode_pingyi_result
data_beta_result<-py_xishu_f_test(mode_pingyi_result,3,c("constant","x1","x2","X3","X4","X5","X6","X7","X8","X9"))
data_beta_result
write.csv(data_beta_result,"data_F_test.csv")
扫一扫
921
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
