实战5 - 抗胰腺癌候选药物的优化建模

1 题目简介

题目来源于2021年中国研究生数学建模竞赛D题——抗乳腺癌候选药物的优化建模。

2 涉及内容

在本次实战的数据分析过程中，涉及以下技术内容：
（1）数据预处理
（2）相关系数
（3）NCA算法
（4）重采样
（5）遗传算法
（6）求AUC值

3 实战步骤

3.1 数据预处理

这里包括数据读取、删除相关性大的列、删除值全都是0的列、将超出mean±4σ的异常值替换为mean±4σ。

global Mdl_Caco_2 Mdl_CYP3A4 Mdl_hERG Mdl_HOB Mdl_MN Mdl_nca_pIC50_ensemble
%读入数据
dir_data1='Molecular_Descriptor.xlsx';
dir_data2='ERα_activity.xlsx';
warning('off');
sheets_1=sheetnames(dir_data1);
data1=readtable(dir_data1,'Sheet',sheets_1(1));
data2=readtable(dir_data2);
data3=readtable(dir_data1,'Sheet',sheets_1(2));

%求出各特征变量之间的相关系数，并删除相关性大的特征变量
table_1=data1(:,2:end);
num_1=data1{:,2:end};
%首先删除全为0的特征变量
del_0=[];
for i=1:size(num_1,2)
    if sum(num_1(:,i))==0
        del_0=[del_0,i];
    end 
end
table_1(:,del_0)=[];
num_1(:,del_0)=[];

%删除相关系数大于0.9的特征变量
corr_data=abs(corrcoef(num_1));
del_x=[];
for i=1:size(corr_data,1)
    if ~ismember(i,del_x) %已经在删除列表里的列不需要再求与其相关性大的列号
        tem1=corr_data(i,:);
        find_x=find(tem1>=0.9);
        find_x(find_x==i)=[];%该列与自己相关系数为1，所以去掉自己的列号
        del_x=[del_x,find_x];
    end
end

u_del_x=unique(del_x);
table_1_ir=table_1;
table_1_ir(:,u_del_x)=[];%table_1_ir为去掉相关性大的列后的表

%将异常值替换为mean±4σ
tem3=table_1_ir{:,:};
tem4=[];
for i=1:size(tem3,2)
    tem1=tem3(:,i);
    tem1((tem1-mean(tem1))>4*std(tem1))=mean(tem1)+4*std(tem1);
    tem1((tem1-mean(tem1))<-4*std(tem1))=mean(tem1)-4*std(tem1);
    tem4=[tem4,tem1];
    
end
tem1=table_1_ir.Properties.VariableNames;
table_1_ir=array2table(tem4);
table_1_ir.Properties.VariableNames=tem1;

3.2 求特征重要度

本来我是使用了NCA算法和树模型两种方法求特征重要度，后来发现用NCA算法求出的特征，训练的模型精度相对较高，所以后来确定的用NCA算法求特征重要的。写论文时大家可以多尝试几种求特征重要度的方法，以使论文内容丰富。

%NCA算法求特征重要度，使用了交叉验证防止过拟合
rng(1) % For reproducibility 
Xtrain=table_1_ir{:,:};
ytrain=data2{:,2};

n = length(ytrain);
cvp = cvpartition(length(ytrain),'kfold',5);
numvalidsets = cvp.NumTestSets;

lambdavals = linspace(0,50,20)*std(ytrain)/n;

lossvals = zeros(length(lambdavals),numvalidsets);

for i = 1:length(lambdavals)
    for k = 1:numvalidsets
        X = Xtrain(cvp.training(k),:);
        y = ytrain(cvp.training(k),:);
        Xvalid = Xtrain(cvp.test(k),:);
        yvalid = ytrain(cvp.test(k),:);

        nca = fsrnca(X,y,'FitMethod','exact', ...
             'Solver','minibatch-lbfgs','Lambda',lambdavals(i), ...
             'GradientTolerance',1e-4,'IterationLimit',30);
        
        lossvals(i,k) = loss(nca,Xvalid,yvalid