毕设论文数据分析记录-part3:各变量因子的相对贡献程度

20220311—了解随机森林回归模型、多元非线性回归

Matlab TreeBagger随机森林回归实例_wokaowokaowokao12345的专栏-CSDN博客_matlab treebagger

20220313—RF工具箱的加载

参考以下两篇，下载了编译器；

matlab MinGW-w64 C/C++ Compiler 的配置(附百度云下载资源）_学术小渣渣的博客-CSDN博客

如何在Matlab 2016a中配置MinGW-w64 C/C++ 编译器_辗迟大虾的博客-CSDN博客

下载RF工具箱；

成功！（还有一点小问题）——后续计算中再看吧

（1）tree bagger------各特性说明

Bag of decision trees - MATLAB- MathWorks 中国

详例：Matlab TreeBagger随机森林回归实例_wokaowokaowokao12345的专栏-CSDN博客_matlab treebagger

官方详例（以上博文应该也有参考官方详例） ————超级详细

使用 TreeBagger 对回归树进行 Bootstrap Aggregation （Bagging） - MATLAB & Simulink - MathWorks 中国

load total6_tmax_r.mat;
data=total6;
Y = data(:,1);
X = data(:,2:end);
rng(1945,'twister');
leaf = [5 10 20 50 100];
col = 'rbcmy';
figure
for i=1:length(leaf)b = TreeBagger(100,X,Y,'Method','regression','OOBPredictorImportance','On','MinLeafSize',leaf(i));plot(oobError(b),col(i))hold on
end
xlabel('Number of Grown Trees')
ylabel('Mean Squared Error')
legend({'5' '10' '20' '50' '100'},'Location','NorthEast')
hold offfigure
plot(oobError(b))
xlabel('Number of Grown Trees')
ylabel('Out-of-Bag Mean Squared Error')figure
bar(b.OOBPermutedPredictorDeltaError)
xlabel('Feature Number')
ylabel('Out-of-Bag Feature Importance')

运行我的数据后，得到误差曲线以及四个自变量重要性分布；

参考其他博文理解袋外误差含义：

随机森林进行特征重要性度量的详细说明_bbbeoy的专栏-CSDN博客

随机森林变量重要性评分及其研究进展 - 道客巴巴

变量的重要性；

参考文章：

论文在线阅读—中国知网

（2）RF工具箱的使用

matlab下利用随机森林包做回归拟合_认真努力，做一只会飞的毛毛虫。-CSDN博客_matlab随机森林回归

自己尝试输代码：随机森林工具箱（调用前需要setenv('MW_MINGW64_LOC','C:\TDM-GCC-64')；mex -setup;）：

load total6_tmax_r.mat;a=randperm(2050);
data=total6;
Train=data(a(1:round(2050*0.75)),:);
Test=data(a(round(2050*0.75)+1:2050),:);P_train=Train(:,2:end);
T_train=Train(:,1);P_test=Test(:,2:end);
T_test=Test(:,1);setenv('MW_MINGW64_LOC','C:\TDM-GCC-64')；
mex -setup;model=regRF_train(P_train,T_train);
Y_hat=regRF_predict(P_test,model);err1=((Y_hat-T_test)./T_test;前期数据准备（6月份各空气污染物对最大温度影响程度）
t=xlsread('tmax_1989.10-2014.09.xlsx','maize_');
a=t(:,2:end);
julyt=[];
for j=9:12:300jt=a(:,j);julyt=[julyt;jt];
end
o3=xlsread('to3_1989.10-2014.09.xlsx','maize.1');
b=o3(2:end,:);
julyo=[];
for j=1:4:100jo=b(:,j);julyo=[julyo;jo];
end
aod=xlsread('aod_1989.10-2014.09.xlsx','maize.1');
c=aod(2:end,:);
julya=[];
for j=1:4:100ja=c(:,j);julya=[julya;ja];
end
so2=xlsread('so2_1989.10-2014.09.xlsx','maize.1');
d=so2(2:end,:);
julys=[];
for j=1:4:100js=d(:,j);julys=[julys;js];
end
co=xlsread('co_1989.10-2014.09.xlsx','maize.1');
e=co(2:end,:);
julyc=[];
for j=1:4:100jc=e(:,j);julyc=[julyc;jc];
end
total6=[zscore(julyt) zscore(julya) zscore(julyo) zscore(julys) zscore(julyc)];

目前不晓得RF回归后，各项参数的意义；

为什么model工作区没有这项呢？

为什么运行后工作区importance只有一列呢？

（偏最小二乘回归了解）

matlab偏最小二乘法及其检验_陆嵩的博客-CSDN博客_偏最小二乘法matlab

偏最小二乘法PLS（matlab自带代码）_billy145533的博客-CSDN博客_pls算法matlab源码

20220316——偏最小二乘回归（集合主成分分析、回归分析等的优点）

（1）主成分分析

特征值和特征向量的理解：https://jingyan.baidu.com/article/27fa7326afb4c146f8271ff3.html

（2）偏最小二乘原理以及matlab实现

pls回归分析，腰围影响了你的体能水平。（matlab与python两种语言编程）_哔哩哔哩_bilibili

（3）代码实战：运用matlab自带plsregression函数运行；

具体参考：博主的代码非常详细；基本复制运行，没有发现有什么问题；傻瓜攻略（十二）——MATLAB实现偏最小二乘回归PLS_素观江湖真的博客-CSDN博客_matlab偏最小二乘回归

最终结果的检验方式，参考：

MATLAB回归预测模型的结果展示和效果检验_素观江湖真的博客-CSDN博客_回归预测模型matlab

（完全搬运以上博主代码，换成了自己的数据）
load('total_t_r_p.mat');
ab0=total66;
mu=mean(ab0);sig=std(ab0); %求均值和标准差
rr=corrcoef(ab0);   %求相关系数矩阵
ab=zscore(ab0); %数据标准化
a=ab(:,[4:end]);    %提出标准化后的自变量数据
b=ab(:,[1:3]);  %提出标准化后的因变量数据
%% 判断提出成分对的个数
[XL,YL,XS,YS,BETA,PCTVAR,MSE,stats] =plsregress(a,b);
xw=a\XS;  %求自变量提出成分的系数，每列对应一个成分，这里xw等于stats.W
yw=b\YS;  %求因变量提出成分的系数
a_0=PCTVAR(1,:);b_0=PCTVAR(2,:);% 自变量和因变量提出成分的贡献率
a_1=cumsum(a_0);b_1=cumsum(b_0);% 计算累计贡献率
i=1;%赋初始值while ((a_1(i)<0.9)&(a_0(i)>0.05)&(b_1(i)<0.9)&(b_0(i)>0.05)) % 提取主成分的条件i=i+1;
end
ncomp=i;% 选取的主成分对的个数
fprintf('%d对成分分别为：\n',ncomp);% 打印主成分的信息
for i=1:ncompfprintf('第%d对成分：\n',i);fprintf('u%d=',i);for k=1:size(a,2)%此处为变量x的个数fprintf('+(%f*x_%d)',xw(k,i),k);endfprintf('\n');fprintf('v%d=',i);for k=1:size(b,2)%此处为变量y的个数fprintf('+(%f*y_%d)',yw(k,i),k);endfprintf('\n');
end
%% 确定主成分后的回归分析
[XL2,YL2,XS2,YS2,BETA2,PCTVAR2,MSE2,stats2] =plsregress(a,b,ncomp);
n=size(a,2); m=size(b,2);%n是自变量的个数,m是因变量的个数
beta3(1,:)=mu(n+1:end)-mu(1:n)./sig(1:n)*BETA2([2:end],:).*sig(n+1:end); %原始数据回归方程的常数项
beta3([2:n+1],:)=(1./sig(1:n))'*sig(n+1:end).*BETA2([2:end],:); %计算原始变量x1,...,xn的系数，每一列是一个回归方程
fprintf('最后得出如下回归方程：\n')
for i=1:size(b,2)%此处为变量y的个数fprintf('y%d=%f',i,beta3(1,i));for j=1:size(a,2)%此处为变量x的个数fprintf('+(%f*x%d)',beta3(j+1,i),j);endfprintf('\n');
end
%% 求预测值
y1 = repmat(beta3(1,:),[size(a,1),1])+ab0(:,[1:n])*beta3([2:end],:);  %求y1,..,ym的预测值
y0 = ab0(:,1:3);  % 真实值%% 画回归系数的直方图
bar(BETA2','k')
%% 贡献率画图
figure
percent_explained1 = 100 * PCTVAR(1,:) / sum(PCTVAR(1,:));
pareto(percent_explained1);
xlabel('主成分')
ylabel('贡献率(%)')
title('主成分对自变量的贡献率')figure
percent_explained = 100 * PCTVAR(2,:) / sum(PCTVAR(2,:));
pareto(percent_explained);
xlabel('主成分')
ylabel('贡献率(%)')
title('主成分对因变量的贡献率')%% 绘制预测结果和真实值的对比
N = size(a,1);% 样本个数
for i =1:size(b,2)yz = y0(:,i);% 真实值yc = y1(:,i);% 预测值    R2 = (N*sum(yc.*yz)-sum(yc)*sum(yz))^2/((N*sum((yc).^2)-(sum(yc))^2)*(N*sum((yz).^2)-(sum(yz))^2)); %计算R方figureplot(1:N,yz,'b:*',1:N,yc,'r-o')legend('真实值','预测值','location','best')xlabel('预测样本')ylabel('值')string = {['第',num2str(i),'个因变量预测结果对比'];['R^2=' num2str(R2)]};title(string)
end%% 三种方法检验网络性能
for i =1:size(b,2)yz = y0(:,i);% 真实值yc = y1(:,i);% 预测值   % 第一种方法，均方误差perf = mse(yz,yc)% 第二种方法，回归图figure;plotregression(yz,yc,['第',num2str(i),'个回归图'])% 第三种方法，误差直方图e = yz-yc;figure;ploterrhist(e,['第',num2str(i),'个误差直方图'])
end

（4）详细复现pls算法原理，实现偏最小二乘回归：
参考：matlab偏最小二乘法及其检验_陆嵩的博客-CSDN博客_偏最小二乘法matlab

（完全搬运以上博主代码，数据换为自己的）
clc % 清屏
clear all; % 删除workplace变量
close all; % 关掉显示图形窗口
load('total_t_r_p.mat');
pz=total66;
mu=mean(pz); %求均值
sig=std(pz); %求标准差
rr=corrcoef(pz); %求相关系数矩阵
data=zscore(pz); %数据标准化
n=4; % n 是自变量的个数
m=3; % m 是因变量的个数
x0=pz(:,m+1:end);y0=pz(:,1:m);%定义自变量为前n列，因变量为n+1到m列
e0=data(:,m+1:end);f0=data(:,1:m);%e0为自变量标准化值，f0为因变量标准化值
num=size(e0,1);%求样本点的个数
chg=eye(n); % w 到 w* 变换矩阵的初始化
for i=1:n%计算 w，w* 和t 的得分向量，matrix=e0'*f0*f0'*e0;[vec,val]=eig(matrix); %求特征值和特征向量val=diag(val); %提出对角线元素，即特征值[val,ind]=sort(val,'descend');%降序排列，ind为排序后原下标序号w(:,i)=vec(:,ind(1)); %提出最大特征值对应的特征向量w_star(:,i)=chg*w(:,i); %计算w*的取值t(:,i)=e0*w(:,i); %计算成分ti 的得分alpha=e0'*t(:,i)/(t(:,i)'*t(:,i)); %计算alpha_ichg=chg*(eye(n)-w(:,i)*alpha'); %计算w 到w*的变换矩阵e=e0-t(:,i)*alpha'; %计算残差矩阵e0=e;%将残差矩阵带进下次循环%计算ss(i)的值beta=[t(:,1:i),ones(num,1)]\f0; %求回归方程的系数beta(end,:)=[]; %删除回归分析的常数项cancha=f0-t(:,1:i)*beta; %求残差矩阵ss(i)=sum(sum(cancha.^2)); %求误差平方和%计算p(i)for j=1:numt1=t(:,1:i);f1=f0;she_t=t1(j,:);she_f=f1(j,:); %把舍去的第j 个样本点保存起来t1(j,:)=[];f1(j,:)=[]; %删除第j 个观测值beta1=[t1,ones(num-1,1)]\f1; %求回归分析的系数beta1(end,:)=[]; %删除回归分析的常数项cancha=she_f-she_t*beta1; %求残差向量p_i(j)=sum(cancha.^2);endp(i)=sum(p_i);if i>1Q_h2(i)=1-p(i)/ss(i-1);elseQ_h2(1)=1;endif Q_h2(i)<0.0975fprintf('提出的成分个数r=%d',i);fprintf('   ');fprintf('交叉的有效性=%f',Q_h2(i));r=i;breakend
end
beta_z=[t(:,1:r),ones(num,1)]\f0; %求Y 关于t 的回归系数
beta_z(end,:)=[]; %删除常数项
xishu=w_star(:,1:r)*beta_z; %求Y 关于X 的回归系数，且是针对标准数据的回归系数，每一列是一个回归方程
mu_x=mu(1:n);mu_y=mu(n+1:end);
sig_x=sig(1:n);sig_y=sig(n+1:end);
for i=1:mch0(i)=mu_y(i)-mu_x./sig_x*sig_y(i)*xishu(:,i); %计算原始数据的回归方程的常数项
end
for i=1:mxish(:,i)=xishu(:,i)./sig_x'*sig_y(i); %计算原始数据的回归方程的系数，每一列是一个回归方程
end
sol=[ch0;xish] %显示回归方程的系数，每一列是一个方程，每一列的第一个数是常数项,每一列为一个因变量与自变量们的回归方程
%此为还原为原始变量后的方程
%% 感觉用途不大，用到的时候再查询怎么使用
save mydata x0 y0 num xishu ch0 xish
w1=w(:,1)
w2=w(:,2)
wx1=w_star(:,1)
wx2=w_star(:,2)
tx1=t(:,1)'
tx2=t(:,2)'
beta_z %回归系数
xishu%系数矩阵,即未还原原始变量的系数,每一列为一个因变量与自变量的回归方程%% 用法：分别计算出第四列第五列第六列和前三列的线性回归关系，给出系数，系数以列的方式给出，
%%sol分别为常数项系数，x1 x2 x3的系数clc % 清屏
clear all; % 删除workplace变量
close all; % 关掉显示图形窗口
format short
load('mydata.mat');%mydata为计算偏最小二乘保存的数据集，可以用于检验
%% 更直观的解释各个自变量的作用
figure
bar(xishu')%分别画出三个自变量对三个因变量标准化后回归方程的系数的的长度图
axis tight
hold on
annotation('textbox',[0.26 0.14 0.086 0.07],'String',{'温度'},'FitBoxToText','off');
annotation('textbox',[0.56 0.14 0.086 0.07],'String',{'辐射量'},'FitBoxToText','off');
annotation('textbox',[0.76 0.14 0.086 0.07],'String',{'降雨'},'FitBoxToText','off');%在指定位置加注释
%% 拟合效果的确定
%所有点都在对角线附近均匀分布，则效果较好
ch0=repmat(ch0,num,1);%repmat起复制矩阵组合为新矩阵的作用
yhat=ch0+x0*xish; %计算y 的预测值
y1max=max(yhat);
y2max=max(y0);
ymax=max([y1max;y2max]);
cancha=yhat-y0; %计算残差
figure
subplot(2,2,1)
plot(0:ymax(1),0:ymax(1),yhat(:,1),y0(:,1),'*')
title('温度成绩预测')
subplot(2,2,2)
plot(0:ymax(2),0:ymax(2),yhat(:,2),y0(:,2),'O')
title('辐射量成绩预测')
subplot(2,1,2)
plot(0:ymax(3),0:ymax(3),yhat(:,3),y0(:,3),'H')
title('降雨量成绩预测')
%% 绘制拟合效果图和权重比重图

偏最小二乘路径模型（plspm）

偏最小二乘法路径模型分析plspm；

科学网—R语言统计：偏最小二乘路径模型（plspm） - 涂波的博文

这篇也比较具体：进一步解释潜在变量；

R数据分析：PLS结构方程模型介绍，论文报告方法和实际操作_公众号Codewar原创作者的博客-CSDN博客

下载R、RSTUDIO

在r中运行该程序包；

读取数据时遇到问题：R语言中xlsx加载时出现报错的解决办法 - 哔哩哔哩

plspm包运行具体步骤：

偏最小二乘法路径模型分析plspm

实际操作有点问题~~~

20220317—继续plspm分析

参考以上博文了解更基本知识：

PLS-PM with R 偏最小二乘路径建模 (R语言)_小天使甲的博客-CSDN博客_plspm包

具体的实现过程：参考以下博文（超级详细，附带返回值的一些例子）

如何用R实现同时拥有形成性、反映型结构的PLS-SEM? - 知乎

R语言，plspm片最小二乘法做结构方程模型 - LISREL、AMOS等结构方程模型分析软件 - 经管之家(原人大经济论坛) 操作中，一些小问题比如工作路径的设置以及数据的加载参考

【R语言】随机森林模型构建_哔哩哔哩_bilibili

setwd("E:/不过尔尔/毕业设计-中期/data-sorting")
getwd()
install.packages("pacman")
pacman::p_load(randomForest,caret,pROC)
library(xlsx)
data<-read.xlsx("./aosc_trp_y_6.xlsx",1)
aod<-c(0,0,0,0,0,0,0,0)
o3<-c(0,0,0,0,0,0,0,0)
so2<-c(0,0,0,0,0,0,0,0)
co<-c(0,0,0,0,0,0,0,0)
tmax<-c(1,1,1,1,0,0,0,0)
rad<-c(1,1,1,1,0,0,0,0)
pre<-c(1,1,1,1,0,0,0,0)
myield<-c(0,0,0,0,1,1,1,0)
sat_path<- rbind(aod,o3,so2,co,tmax,rad,pre,myield)
innerplot(sat_path)
sat_blocks <- list(1,2,3,4,5,6,7,8)
sat_mod<-rep("A",8)
satpls <- plspm(data, sat_path,sat_blocks, modes = sat_mod,scaled = FALSE)
summary(satpls)
innerplot(satpls)前期数据准备：（循环6-9四个月aosc-trp-my数据对应）
list=dir('aosc*.xlsx');
for z=1:4
t=xlsread('tmax_1989.10-2014.09.xlsx','maize_');
a=t(:,2:end);
julyt=[];
for j=8+z:12:300jt=a(:,j);julyt=[julyt jt];
end
juct=[];
for k=1:length(julyt)
yt=julyt(k,:);
avaluet=[];
for i=1:length(yt)-1
valuea = yt(i+1) - yt(i);
avaluet=[avaluet;valuea];
end
juct=[juct;avaluet];
end
o3=xlsread('to3_1989.10-2014.09.xlsx','maize.1');
b=o3(2:end,:);
julyo=[];
for j=z:4:100jo=b(:,j);julyo=[julyo jo];
end
juco=[];
for k=1:length(julyo)
xo=julyo(k,:);
avalueo=[];
for i=1:length(xo)-1
valueb = xo(i+1) - xo(i);
avalueo=[avalueo;valueb];
end
juco=[juco;avalueo];
end
aod=xlsread('aod_1989.10-2014.09.xlsx','maize.1');
c=aod(2:end,:);
julya=[];
for j=z:4:100ja=c(:,j);julya=[julya ja];
end
juca=[];
for k=1:length(julya)
xa=julya(k,:);
avaluea=[];
for i=1:length(xa)-1
valuec = xa(i+1) - xa(i);
avaluea=[avaluea;valuec];
end
juca=[juca;avaluea];
end
so2=xlsread('so2_1989.10-2014.09.xlsx','maize.1');
d=so2(2:end,:);
julys=[];
for j=z:4:100js=d(:,j);julys=[julys js];
end
jucs=[];
for k=1:length(julys)
xs=julys(k,:);
avalues=[];
for i=1:length(xs)-1
valued = xs(i+1) - xs(i);
avalues=[avalues;valued];
end
jucs=[jucs;avalues];
end
co=xlsread('co_1989.10-2014.09.xlsx','maize.1');
e=co(2:end,:);
julyc=[];
for j=z:4:100jc=e(:,j);julyc=[julyc jc];
end
jucc=[];
for k=1:length(julyc)
xc=julyc(k,:);
avaluec=[];
for i=1:length(xc)-1
valuee = xc(i+1) - xc(i);
avaluec=[avaluec;valuee];
end
jucc=[jucc;avaluec];
end
p=xlsread('pre_1989.10-2014.09.xlsx','maize_');
f=p(:,2:end);
julyp=[];
for j=8+z:12:300jp=f(:,j);julyp=[julyp jp];
end
jucp=[];
for k=1:length(julyp)
yp=julyp(k,:);
avaluep=[];
for i=1:length(yp)-1
valuef= yp(i+1) - yp(i);
avaluep=[avaluep;valuef];
end
jucp=[jucp;avaluep];
end
rad=xlsread('rad_1989.10-2014.09.xlsx','maize.1');
g=rad;
julyr=[];
for j=z:4:100jr=g(:,j);julyr=[julyr jr];
end
jucr=[];
for k=1:length(julyr)
yr=julyr(k,:);
avaluer=[];
for i=1:length(yr)-1
valueg = yr(i+1) - yr(i);
avaluer=[avaluer;valueg];
end
jucr=[jucr;avaluer];
end
maize=xlsread('gai-maizey.xls','gai-maizey');
h=maize(:,5:end-1);
jucm=[];
for k=1:length(h)
ym=h(k,:);
avaluem=[];
for i=1:length(ym)-1
valueh = ym(i+1) - ym(i);
avaluem=[avaluem;valueh];
end
jucm=[jucm;avaluem];
end
total66=[juca juco jucs jucc juct jucr jucp jucm ];
total666=zscore(total66);
xlswrite(list(z).name,total66,1);
xlswrite(list(z).name,total666,2);
end

带入我的数据后，首先是外部、内部模型的一个反馈：

参考：【[R] 【教程】教你如何读懂线性回归lm的结果summary（判断显著性）[转]】_lq497028254的博客-CSDN博客_r语言summary结果解读

各加载分数情况：

第一次试运行完，在观测变量和潜变量这块，我们这里的观测变量就是温度降雨等这些潜变量实际的数据，相当于每一个潜变量只对应一种观测变量，然后在模型模式选择的时候，是应该选reflective还是formative呢？自己A和B模式都试过之后，发现结果没有大的差别；数据做标准化处理前后也没有大的差别；

个人理解，PLSPM模型，有PLS回归方法的原理在其中，最后的加载分数，从回归的角度反应两个变量之间的相关性，可以显示出多对一时某个变量较大的影响；

R数据分析：PLS结构方程模型介绍，论文报告方法和实际操作_公众号Codewar原创作者的博客-CSDN博客

20220318——R语言里随机森林，PLSPM，整合数据

R中数据的写入

R语言--数据框导出excel文件 - 知乎

正则：匹配以某字符串开头或不以某字符串开头的字符串 - 百度文库

R语言：提取路径中的文件名字符串（basename函数）_weixin_33858336的博客-CSDN博客

各种声明工作蒲等参考下文

【原创推荐】使用R写入Excel方法总结

R语言心得说：R语言之xlsx包读写Excel数据_狼の牙的博客-CSDN博客_r语言读取excel数据

自己运行时遇到问题：

wb<-loadWorkbook("./aosc_pre.xlsx")
sheets <- getSheets(wb)
r<-as.data.frame(defaultSummary(data.frame(obs=testdata$pre,pred=testpred)))
addDataFrame(r, sheet = "Sheet2" ,row.names=FALSE,col.names = F,startColumn = 6)
saveWorkbook(wb, file="./aosc_pre.xlsx")

$ operator is invalid for atomic vectors???

参考以下博文也没有解决

R语言收捲时出错: $ operator is invalid for atomic vectors 求大神们指教-CSDN论坛

R中随机森林学习

R语言机器学习模型-随机森林模型解决回归、二分类、多分类问题（已完结）_哔哩哔哩_bilibili

！！！for循环语句+if（if语句中的格式!!!!）！！！

R语言中if判断语句 - 小鲨鱼2018 - 博客园

以下代码参考R语言机器学习模型-随机森林模型解决回归、二分类、多分类问题（已完结）_哔哩哔哩_bilibili

自己加了循环（太难了）

随机森林：

library(randomForest)
library(tidyverse)
library(skimr)
library(DataExplorer)
library(caret)
library(pROC)
library(xlsx)
setwd("E:/不过尔尔/毕业设计-中期/data-sorting/RF")
filelist <- list.files(pattern="^6.*")
num<-c(1,2,3)
for (i in num)
{
mydata<-read.xlsx(filelist[i],sheetName = "Sheet1")
if(i==1)
{y1=mydata$pre
}else
if(i==2)
{y1=mydata$rad
}else
{y1=mydata$tma
}
set.seed(42)
trains<-createDataPartition(y1<-y1,p=0.75,list=F
)
traindata<-mydata[trains,]
testdata<-mydata[-trains,]
# hist(traindata$pre,breaks = 50)
# hist(testdata$pre,breaks = 50)colnames(mydata)s=substr(basename(filelist[i]), start = 7, stop = 9)st=paste(s,"~")form_reg<-as.formula(paste0(st,paste(colnames(traindata)[1:4],collapse="+")))form_reg
set.seed(42)
fit_rf_reg<-randomForest(form_reg,data=traindata,ntree=500,mtry=4,importance=T
)
fit_rf_reg
# plot(fit_rf_reg,main="ERROR&TREES")
importance(fit_rf_reg)
#
# varImpPlot(fit_rf_reg,main="Variable Importance Plot")
# varImpPlot(fit_rf_reg,main="Variable Importance Plot",type=1)
# varImpPlot(fit_rf_reg,main="Variable Importance Plot",type=2)# partialPlot(x=fit_rf_reg,
#             pred.data = traindata,
#             x.var=aod)
# plot(pre~aod,data=traindata)testpred<-predict(fit_rf_reg,newdata = testdata)
if(i==1)
{t=testdata$pre
}else
if(i==2)
{t=testdata$rad
}else
{t=testdata$tma
}
defaultSummary(data.frame(obs=t,pred=testpred))
# plot(x=testdata$pre,
#      y=testpred,
#      xlab = "Actual",
#      ylab = "Prediction",
#      main="随机森林-实际值与预测值比较",
#      sub="测试集")
# testlinmod<-lm(testpred~testdata$pre)
# abline(testlinmod,col="blue",lwd=2.5,lty="solid")
# abline(a=0,b=1,col="red",lwd=2.5,lty="dashed")
# legend("topleft",
#        legend=c("Model","Base"),
#        col=c("bule","red"),
#        lwd=2.5,
#        lty=c("solid","dashed"))
#
# trainpred<-predict(fit_rf_reg,newdata = traindata)
# defaultSummary(data.frame(obs=traindata$pre,pred=trainpred))
# plot(x=traindata$pre,
#      y=trainpred,
#      xlab = "Actual",
#      ylab = "Prediction",
#      main="随机森林-实际值与预测值比较",
#      sub="训练集")
# trainlinmod<-lm(trainpred-traindata$pre)file.path(filelist[i])
path<-file.path(filelist[i])
r<-as.data.frame(defaultSummary(data.frame(obs=t,pred=testpred)))
write.xlsx(fit_rf_reg[["importance"]],file=path,sheetName = "Sheet2",append = TRUE)
write.xlsx(r,file=path,sheetName = "Sheet3",append = TRUE)
}

plspm路径分析

setwd("E:/不过尔尔/毕业设计-中期/data-sorting")
getwd()
library(plspm)
library(xlsx)
filelist <- list.files(pattern="^aosc.*")
num<-c(1,2,3,4)
for (i in num)
{
data<-read.xlsx(filelist[i],sheetName = "Sheet2")
aod<-c(0,0,0,0,0,0,0,0)
o3<-c(0,0,0,0,0,0,0,0)
so2<-c(0,0,0,0,0,0,0,0)
co<-c(0,0,0,0,0,0,0,0)
tmax<-c(1,1,1,1,0,0,0,0)
rad<-c(1,1,1,1,0,0,0,0)
pre<-c(1,1,1,1,0,0,0,0)
myield<-c(0,0,0,0,1,1,1,0)
sat_path<- rbind(aod,o3,so2,co,tmax,rad,pre,myield)
sat_blocks <- list(1,2,3,4,5,6,7,8)
sat_mod<-rep("A",8)
satpls<- plspm(data, sat_path,sat_blocks, modes = sat_mod,scaled = FALSE)
summary(satpls)
if(i==1)
{a<-satpls[["inner_model"]]
}else
if(i==2)
{b<-satpls[["inner_model"]]
}else
if(i==3)
{c<-satpls[["inner_model"]]
}else
{d<-satpls[["inner_model"]]
}file.path(filelist[i])
path<-file.path(filelist[i])
write.xlsx(satpls[["effects"]],file=path,sheetName = "Sheet3",append = TRUE)
}

试着出图：：：

20220319——出图分析

关于Plspm详细的一些系数说明

R数据分析：PLS结构方程模型介绍，论文报告方法和实际操作_公众号Codewar原创作者的博客-CSDN博客_pls结构方程模型

Is there any goodness of fit index for PLS?

出图：

首先是PLSPM路径图（在VISIO中画的）；

还有随机森林各变量相对贡献率（尝试用双层饼图来表示）

20220320——结合最小二乘路径以及随机森林得到的变量贡献率分析

20220322——前述方法分析冬小麦

完整的流程：生长期月份数据分别整理进EXCEL表格（matlab），为plspm数据处理做准备（R），然后进行RF分析（R）

list=dir('waosc*.xlsx')（该步骤没有用在小麦分析中，选取的月份较为分散）t=xlsread('tmax_1989.10-2014.09.xlsx','wheat_');
a=t(:,2:end);
julyt=[];
for j=5:12:300jt=a(:,j:j+1);jt=mean(jt,2);julyt=[julyt jt];
end
juct=[];
for k=1:length(julyt)
yt=julyt(k,:);
avaluet=[];
for i=1:length(yt)-1
valuea = yt(i+1) - yt(i);
avaluet=[avaluet;valuea];
end
juct=[juct;avaluet];
end
o3=xlsread('to3_1989.10-2014.09.xlsx','wheat.1');
b=o3(2:end,:);
julyo=[];
for j=5:9:225jo=b(:,j:j+1);jo=mean(jo,2);julyo=[julyo jo];
end
juco=[];
for k=1:length(julyo)
xo=julyo(k,:);
avalueo=[];
for i=1:length(xo)-1
valueb = xo(i+1) - xo(i);
avalueo=[avalueo;valueb];
end
juco=[juco;avalueo];
end
aod=xlsread('aod_1989.10-2014.09.xlsx','wheat.1');
c=aod(2:end,:);
julya=[];
for j=5:9:225ja=c(:,j:j+1);ja=mean(ja,2);julya=[julya ja];
end
juca=[];
for k=1:length(julya)
xa=julya(k,:);
avaluea=[];
for i=1:length(xa)-1
valuec = xa(i+1) - xa(i);
avaluea=[avaluea;valuec];
end
juca=[juca;avaluea];
end
so2=xlsread('so2_1989.10-2014.09.xlsx','wheat.1');
d=so2(2:end,:);
julys=[];
for j=5:9:225js=d(:,j:j+1);js=mean(js,2);julys=[julys js];
end
jucs=[];
for k=1:length(julys)
xs=julys(k,:);
avalues=[];
for i=1:length(xs)-1
valued = xs(i+1) - xs(i);
avalues=[avalues;valued];
end
jucs=[jucs;avalues];
end
co=xlsread('co_1989.10-2014.09.xlsx','wheat.1');
e=co(2:end,:);
julyc=[];
for j=5:9:225jc=e(:,j:j+1);jc=mean(jc,2);julyc=[julyc jc];
end
jucc=[];
for k=1:length(julyc)
xc=julyc(k,:);
avaluec=[];
for i=1:length(xc)-1
valuee = xc(i+1) - xc(i);
avaluec=[avaluec;valuee];
end
jucc=[jucc;avaluec];
end
p=xlsread('pre_1989.10-2014.09.xlsx','wheat_');
f=p(:,2:end);
julyp=[];
for j=5:12:300jp=f(:,j:j+1);jp=mean(jp,2);julyp=[julyp jp];
end
jucp=[];
for k=1:length(julyp)
yp=julyp(k,:);
avaluep=[];
for i=1:length(yp)-1
valuef= yp(i+1) - yp(i);
avaluep=[avaluep;valuef];
end
jucp=[jucp;avaluep];
end
rad=xlsread('rad_1989.10-2014.09.xlsx','wheat_');
g=rad(:,2:end);
julyr=[];
for j=5:12:300jr=g(:,j:j+1);jr=mean(jr,2);julyr=[julyr jr];
end
jucr=[];
for k=1:length(julyr)
yr=julyr(k,:);
avaluer=[];
for i=1:length(yr)-1
valueg = yr(i+1) - yr(i);
avaluer=[avaluer;valueg];
end
jucr=[jucr;avaluer];
endwheat=xlsread('gai-wheaty.xls','gai-wheaty');
h=wheat(:,4:end);
jucw=[];
for k=1:length(h)
yw=h(k,:);
avaluew=[];
for i=1:length(yw)-1
valueh = yw(i+1) - yw(i);
avaluew=[avaluew;valueh];
end
jucw=[jucw;avaluew];
endtotal66=[juca juco jucs jucc juct jucr jucp jucw ];
total666=zscore(total66);
xlswrite('waosc_trp_wy_2_3.xlsx',total66,1);
xlswrite('waosc_trp_wy_2_3.xlsx',total666,2);

（plspm、rf分析都如前所述）

20220323——主成分分析、偏最小二乘原理

整理完了夏玉米、冬小麦数据，认真想了一下自己用的方法，为什么用？一些原理好像懂了，又好像没懂；决定今天理解其中的原因！！！

1.相关系数与回归系数；【数学式有联系（协方差、各自方差），但是含义不同；回归系数会受到量纲数量级的影响】

2.回归分析；（最常用的是最小二乘法）（回归系数的显著性检验，最终模型的R2解释度不高不影响，主要是各变量能否通过显著性检验，认为其在概率的角度会对因变量产生影响）；

回归分析|笔记整理（2）——一元线性回归（下） - 知乎

3.主成分分析（多重共线性，自变量之间存在关系使得因变量与自变量的表达式出现很多种从而导致y与x1、、、xn之间的关系不清楚）

一看就懂的多重共线性 - 知乎

消除共线性的影响，降维，并且构造的新变量后，使其对y的方差贡献最大，且各新变量之间线性无关，在原理的过程中，方差、协方差之间的计算，特征值与特征变量；

4.结构方程：内在模型，潜变量（不容易直接测量的，比如：、、能力）；外在模型，表示潜变量的观测数据；内容来自于：结构方程模型原理及应用_哔哩哔哩_bilibili

（）原生矩阵，计算机根据模型构建的矩阵相差不大（协方差矩阵）；拟合优度指数（卡方越小越好，其他指数）

（）探索性因子分析——验证性因子分析，分类依据不确定；（主成分分析提取类别，各个因子在类别中的关系）

（）SEM；【例如：学科相关性】（各因子的拟合指数；自由度（大）、拟合程度；输出：原矩阵与模型矩阵差异最小（协方差）、各路径参数，因子负荷、相关系数，总的拟合指数）

（）组成：外部模型，内部模型；优点：多个因子关系共同表示；传统自变量与因变量假设没有误差（观察得分、真分数，误差），可做的分析多样（回归分析、因子分析、检验、交互作用）

（）多质多法；

（）结构方程分析步骤

最大似然回归

也参考另外一个课程：01-结构方程模型（SEM）基础介绍_哔哩哔哩_bilibili

标准化和非标准化的路径系数（看P值即可；两者代表的含义有一丝丝差别）

总结：偏最小二乘路径分析，本质上是回归，但可以同时进行多对多的回归，偏最小二乘路径结合偏最小而成回归原理，结合主成分分析和传统回归，对相应的自变量因变量进行回归；

20220324—结合显著路径和随机森林，选择各月份典型污染物与最终的单产作多元线性拟合，看其影响程度大小

vif+fitlm(经过因子膨胀分析后，不存在多重共线性，可用一般的多元线性拟合)

（在分析小麦时，4月份co和5月份的co两个自变量vif>10，尝试去除其中的一个进行拟合）

% maize
load('pz6-pls.mat');
total=zscore(total);
o311=reshape(total(:,1),[82,24]);
aod23=reshape(total(:,2),[82,24]);
o34=reshape(total(:,3),[82,24]);
so28=reshape(total(:,4),[82,24]);
aod5=reshape(total(:,5),[82,24]);
aod9=reshape(total(:,6),[82,24]);
wy=reshape(total(:,7),[82,24]);bm11=[];bm23=[];bm4=[];bm5=[];bm55=[];bm555=[];
v=[];
xxi=[];p=[];r2=[];
for k=1:size(o311,1)
bm11=o311(k,:);bm23=aod23(k,:);bm4=o34(k,:);bm5=so28(k,:);bm55=aod5(k,:);bm555=aod9(k,:);
x=[bm11(:) bm23(:) bm4(:) bm5(:) bm55(:) bm555(:)];
rx=corrcoef(x);
vif=diag(inv(rx));vif=vif';
v=[v;vif];y=wy(k,:);y=y(:);
md1=fitlm(x,y);
xi=md1.Coefficients.Estimate;xi=xi';
xxi=[xxi;xi];
pi=md1.Coefficients.pValue;pi=pi';
p=[p;pi];
r22=md1.Rsquared.Ordinary;
r2=[r2;r22];
end% wheat
load('pw8-pls.mat');
total=zscore(total);
o311=reshape(total(:,1),[65,24]);
aod23=reshape(total(:,2),[65,24]);
o34=reshape(total(:,3),[65,24]);
% co4=reshape(pw8(:,4),[65,24]);
aod5=reshape(total(:,5),[65,24]);
so25=reshape(total(:,6),[65,24]);
co5=reshape(total(:,7),[65,24]);
wy=reshape(total(:,8),[65,24]);bm11=[];bm23=[];bm4=[];bm5=[];bm55=[];bm555=[];
% bm44=[];
v=[];
xxi=[];p=[];r2=[];
for k=1:size(o311,1)
bm11=o311(k,:);bm23=aod23(k,:);bm4=o34(k,:);bm5=aod5(k,:);bm55=so25(k,:);bm555=co5(k,:);
% bm44=co4(k,:);
x=[bm11(:) bm23(:) bm4(:) bm5(:) bm55(:) bm555(:)];
rx=corrcoef(x);
vif=diag(inv(rx));vif=vif';
v=[v;vif];y=wy(k,:);y=y(:);
md1=fitlm(x,y);
xi=md1.Coefficients.Estimate;xi=xi';
xxi=[xxi;xi];
pi=md1.Coefficients.pValue;pi=pi';
p=[p;pi];
r22=md1.Rsquared.Ordinary;
r2=[r2;r22];
end

毕设论文数据分析记录-part3:各变量因子的相对贡献程度相关推荐

毕设论文数据分析记录-part2：相关性分析
20220225-相关性分析 polyfit线性拟合后进行t检验判断线性关系显著水平-一阶差分去趋势 a=xlsread('gai-maizey.xls','gai-maizey'); a1=a(:, ...
毕设论文数据分析记录-part1：长时间序列的趋势分析
20220210-MK非参数检验 (刘烁楠. 变化环境下月水量平衡模型参数时变特征研究[D]. 2019.) MK检验_哔哩哔哩_bilibili 通过UP主视频里的详细EXCEL过程,理解了MK具体 ...
【计算机视觉】Mip-nerf 论文精读记录
[计算机视觉]Mip-nerf 论文精读记录本人是刚入门的计算机视觉小白,此系列为nerf论文精读系列笔记记录,感兴趣的朋友可以关注一下,共同成长! Mip-NeRF: A Multiscale R ...
t检验的p值对照表_论文数据分析实战 | 如何对汇总数据进行t检验
在SPSS统计分析交流群中有学员在阅读论文的过程中看到下面的这张表格: 这张表中记录了第16届世界男子篮球锦标赛中国队与前8名球队进攻指标比较的结果,其中这份表格并没有给出详细的P值,而只是告诉我们P ...
【秃头系列】-【本科生毕设论文格式Word】自动生成目录并调整目录
文章目录 01 - 论文目录 02 - 如何自动生成 03 - 调整目录 04 - 使用效果 05 - 总结上一文:[秃头系列]-[本科生毕设论文格式Word]自动生成论文多级标题并排版正文 ...
论文数据分析-1(论文数据统计)
这是在学习数据分析的一个实例,论文数据分析,这是第一部分,笔者刚学习此项内容,有问题大家提出来,不喜勿喷. 任务1:论文数据统计1 1.1 任务说明任务主题:论文数量统计,即统计2019年全年计算机 ...
ICCV2017 论文浏览记录(转)
mark一下,感谢作者分享! 作者将ICCV2017上的论文进行了汇总,在此记录下来,平时多注意阅读积累. 之前很早就想试着做一下试着把顶会的论文浏览一遍看一下自己感兴趣的,顺便统计一下国内高校或者研 ...
ICCV2017 论文浏览记录
之前很早就想试着做一下试着把顶会的论文浏览一遍看一下自己感兴趣的,顺便统计一下国内高校或者研究机构的研究方向,下面是作为一个图像处理初学者在浏览完论文后的觉得有趣的文章: ICCV2017 论文浏览 ...
数据分析记录（六）--多元线性回归在SPSS中的实现（步骤及指标含义）
数据分析记录(六)–多元线性回归在SPSS中的实现(步骤及指标含义) 本文仅作为自己的学习记录以备以后复习查阅在回归分析中,如果有两个或两个以上的自变量,就称为多元回归.事实上,一种现象常常是与多个 ...

毕设论文数据分析记录-part3:各变量因子的相对贡献程度

毕设论文数据分析记录-part3:各变量因子的相对贡献程度相关推荐

最新文章

热门文章