基于R软件的车险欺诈识别.docx

企业资源管理与信息化大作业第2组基于R软件的车险欺诈识别1. 保险欺诈简介1.1保险欺诈定义保险欺诈国际上一般也称保险犯罪。严格意义上说，保险欺诈较保险犯罪含义更广。保险当事人双方都可能构成保险欺诈。凡保险关系投保人一方不遵守诚信原则，故意隐瞒有关保险标的的真实情况，诱使保险人承保，或者利用保险合同内容，故意制造或捏造保险事故造成保险公司损害，以谋取保险赔付金的，均属投保方欺诈。凡保险人在缺乏必要偿付能力或未经批准擅自经营业务，并利用拟订保险条款和保险费率的机会，或夸大保险责任范围诱导、欺骗投保人和被保险人的，均属保险人欺诈。保险欺诈一经实施，必然造成危害结果，有必要严加防范。1.2保险欺诈产生原因1.2.1 主体本身问题保险主体有投保人、保险人和被保险人，由于之间的信息不对称，所以存在道德风险，这为欺诈提供了可能。投保人有可能会隐瞒标的真实信息；保险人因为更专业，所以可能在保险条款和宣传中暗藏玄机；被保险人则有可能为了获得更多的理赔费尽心机。在市场交易和合同关系中，这种信息不对称普遍存在，所以欺诈也普遍存在。主体多元化也是保险欺诈的一个因素。由于保险分成产险和人身险，人身险包括寿险、健康险、人身意外伤害险等，而产险有包括企业财产险、工程险、责任险、货运险、保证险、船舶险、家财险等险种。险种的多样化，参与主体就变得复杂，尤其是每一个险种有着巨大的市场受众，从而让风险变得很不好控制。比如车险，目前中国的汽车保有量超过1.37 亿，每辆车肯定都买了保险，这么多的车主道德分数是否都合格？保险公司有没有不诚信？这些都很难说清楚。1.2.2 保险的特性射幸性除了保险主体的因素，保险合同的射幸性也是促成保险欺诈的一个因素。如果保险期间没有出险，被保险人会觉得自己的保费“白交了”，因此他会考虑减少保费，或者得到理赔。这样就会导致一些人机会性欺诈、夸大损失，或者无中生有“制造出险”，甚至职业性欺诈。1.2.3 管理问题保险公司方面，对保险中的核保、核赔过程不严谨，导致没有有效识别保险欺诈。另外，由于公司对欺诈的识别建设不够，加之保险欺诈识别成果行业不共享，将使欺诈不能有效识别，即使识别出来，悲剧可能会在另外的保险公司上演。对业务人员的培训不够也是一个因素；另一方面业务人员的职业道德约束不够，保险欺诈有可能是“内鬼”所为。行业规范也有待加强，其中最重要的是对保险欺诈的惩罚。目前虽然有刑法、保险法等相关法律法规加以约束，但是对欺诈的惩罚较轻，欺诈成本低，使得有些人愿意铤而走险。1.3目前保险欺诈现状2014 年5 月，中国保险行业协会发布了2003 年以来十起反保险欺诈典型案例，并宣布将于9 月份成立反保险欺诈专业委员会，以预防和打击保险欺诈犯罪。据中保协消息，随着我国保险业的发展，保险欺诈案件呈逐年上升之势，保险欺诈呈现团伙化、专业化和职业化等特征，反欺诈形势日益严峻。1近年来，我国保险欺诈活动频繁。据保监会数据，2013年各级稽查部门共累计查实违法违规资金23 亿元、违法违规行为118 项，对639 家机构和820 人实施1764 项次行政处罚，指导协调保险公司完成责任追究172 起案件，组织行业向公安机关移送涉嫌车险欺诈案件2375 起，涉案金额1.37 亿元，挽回经济损失7580 万元；对6.7 亿元股本资金来源、4531 家新设保险公司及其分支机构反洗钱制度进行反洗钱审查，对1 万多名高管进行反洗钱培训测试。2保险欺诈涉案金额占保险市场金额的比例高，对保险市场影响深远。2014 年5 月，据中保协工作人员表示，近年来我国保险欺诈案件呈逐年上升趋势，跨境案例增多。在保险中车险欺诈占比最高，约占50%至70%。在寿险中，高额意外险是重点领域。3在国际上，保险诈骗金额约占赔付总额的10%30%。某些险种的欺诈金额占比甚至高达50%。保险欺诈已成为世界各国保险业不得不面对的共同难题。一直以来，我国重视对保险欺诈的法制规范，随着目前形势的严峻有加强之势。对于保险活动的规范，除了保险法、刑法、民法通则、合同法、证券法、保险公司管理规定、道路交通安全法及其实施条例、交强险条例、关于规范人身保险经营行为有关问题的通知等一系列相关文件，近年来也根据形势，推出了一系列相关规定或解释。2保险欺诈识别针对保险欺诈向全球蔓延，研究方法开始由定性研究向实证研究过渡。早期的识别模型主要是以Logit模型为主。Artis、Ayuso和Guillen建立了AAG欺诈识别模型。Caudill建立了多项分对数模型。Ridit模型通过建立标准组，将其他组的数据与之做对比来计算R值。我国学者叶明华运用Logit回归分析对识别因子进行精炼，将BP神经网络应用在车险欺诈识别中，得到了更高的准确率3。但在实际操作中，有两种原因造成准确率发生偏差。一种是数据样本小，缺乏对于现实中包含微弱信息的大型数据库的数据分析。另一种是以理想的数据平衡状态作为研究背景，且现有分类器的设计都是基于类分布大致平衡这一假设，脱离实际4。针对现实世界中广泛存在不平衡数据这一事实，本研究利用随机森林对欺诈进行识别预测，并同Logistic算法进行对比。3 算法简介3.1 随机森林随机森林，顾名思义，是用随机的方式建立一个森林，森林里面有很多的决策树组成，随机森林的每一棵决策树之间是没有关联的。在得到森林之后，当有一个新的输入样本进入的时候，就让森林中的每一棵决策树分别进行一下判断，看看这个样本应该属于哪一类（对于分类算法），然后看看哪一类被选择最多，就预测这个样本为那一类。由上不难看出，决策树是随机森林的基础，所以首先介绍决策树。3.1.1 决策树决策树实际上是将空间用超平面进行划分的一种方法，每次分割的时候，都将当前的空间一分为二，比如说下面的决策树：就是将空间划分成下面的样子：这样使得每一个叶子节点都是在空间中的一个不相交的区域，在进行决策的时候，会根据输入样本每一维feature的值，一步一步往下，最后使得样本落入N个区域中的一个（假设有N个叶子节点）。决策树生长的核心是确定决策树的分枝准则：n 如何从众多的属性变量中选择一个当前的最佳分支变量；也就是选择能使异质性下降最快的变量。异质性的度量：GINI、TWOING、least squared deviation。前两种主要针对分类型变量，LSD针对连续性变量。n 如何从分支变量的众多取值中找到一个当前的最佳分割点（分割阈值）。分割阈值：A、数值型变量对记录的值从小到大排序，计算每个值作为临界点产生的子节点的异质性统计量。能够使异质性减小程度最大的临界值便是最佳的划分点。B、分类型变量列出划分为两个子集的所有可能组合，计算每种组合下生成子节点的异质性。同样，找到使异质性减小程度最大的组合作为最佳划分点。n 决策树停止生长的条件满足以下一个即停止生长。（1）节点达到完全纯性；（2）数树的深度达到用户指定的深度；（3）节点中样本的个数少于用户指定的个数；（4）异质性指标下降的最大幅度小于用户指定的幅度。3.1.2 随机森林简介随机森林利用随机的方式将许多决策树组合成一个森林，每个决策树在分类的时候投票决定测试样本的最终类别。下面我们再详细说一下随机森林是如何构建的。随机森林主要包括4个部分：随机选择样本；随机选择特征；构建决策树；随机森林投票分类。 随机选择样本给定一个训练样本集，数量为N，我们使用有放回采样到N个样本，构成一个新的训练集。注意这里是有放回的采样，所以会采样到重复的样本。详细来说，就是采样N次，每次采样一个，放回，继续采样。即得到了N个样本。然后我们把这个样本集作为训练集，进入下面的一步。 随机选择特征 在构建决策树的时候，我们前面已经讲过如何在一个节点上，计算所有特征的Information Gain（ID3） 或者 Gain Ratio（C4.5），然后选择一个最大增益的特征作为划分下一个子节点的走向。但是，在随机森林中，我们不计算所有特征的增益，而是从总量为M的特征向量中，随机选择m个特征，其中m可以等于sqrt(M)，然后计算m个特征的增益，选择最优特征（属性）。注意，这里的随机选择特征是无放回的选择！所以，随机森林中包含两个随机的过程：随机选择样本，随机选择特征。 构建决策树有了上面随机产生的样本集，我们就可以使用一般决策树的构建方法，得到一棵分类（或者预测）的决策树。需要注意的是，在计算节点最优分类特征的时候，我们要使用上面的随机选择特征方法。而选择特征的标准可以是我们常见的Information Gain（ID3） 或者 Gain Ratio（C4.5）。 随机森林投票分类 通过上面的三步走，我们可以得到一棵决策树，我们可以重复这样的过程H次，就得到了H棵决策树。然后来了一个测试样本，我们就可以用每一棵决策树都对它分类一遍，得到了H个分类结果。这时，我们可以使用简单的投票机制，获得该测试样本的最终分类结果。随机森林优缺点分析 优点：n 它能够处理很高维度（feature很多）的数据，并且不用做特征选择；n 由于随机选择样本导致的每次学习决策树使用不同训练集，所以可以一定程度上避免过拟合；过拟合（overfittingt）是这样一种现象：一个假设在训练数据上能够获得比其他假设更好的拟合，但是在训练数据外的数据集上却不能很好的拟合数据。 缺点：n 随机森林已经被证明在某些噪音较大的分类或回归问题上会过拟合；n 对于有不同级别的属性的数据，级别划分较多的属性会对随机森林产生更大的影响，所以随机森林在这种数据上产出的属性权值是不可信的3.2 逻辑回归Logistic回归与多重线性回归实际上有很多相同之处，最大的区别就在于它们的因变量不同，其他的基本都差不多。正是因为如此，这两种回归可以归于同一个家族，即广义线性模型（generalizedlinear model）。这一家族中的模型形式基本上都差不多，不同的就是因变量不同。n 如果是连续的，就是多重线性回归；n 如果是二项分布，就是Logistic回归；n 如果是Poisson分布，就是Poisson回归；n 如果是负二项分布，就是负二项回归。Logistic回归的因变量可以是二分类的，也可以是多分类的，但是二分类的更为常用，也更加容易解释。所以实际中最常用的就是二分类的Logistic回归。3.2.1 Logistic回归的主要用途：寻找危险因素：寻找某一疾病的危险因素等；预测：根据模型，预测在不同的自变量情况下，发生某病或某种情况的概率有多大；判别：实际上跟预测有些类似，也是根据模型，判断某人属于某病或属于某种情况的概率有多大，也就是看一下这个人有多大的可能性是属于某病。3.2.2常规步骤Regression问题的常规步骤为：1 寻找h函数（即hypothesis）；2 构造J函数（损失函数）；3 想办法使得J函数最小并求得回归参数（）n 构造预测函数hLogistic回归虽然名字里带“回归”，但是它实际上是一种分类方法，主要用于两分类问题（即输出只有两种，分别代表两个类别），所以利用了Logistic函数（或称为Sigmoid函数），函数形式为：对于线性边界的情况，边界形式如下：构造预测函数为：函数的值有特殊的含义，它表示结果取1的概率，因此对于输入x分类结果为类别1和类别0的概率分别为：n 构造损失函数JCost函数和J函数如下，它们是基于最大似然估计推导得到的。n 梯度下降法求的最小值更新过程：更新过程可以写成：4 算法性能评价指标4.1 ROC曲线 ROC曲线，即受试者工作特征曲线（receiver operating characteristic curve），曲线及其分析以统计决策理论为基础,起源于电子信号观测理论、用于雷达信号接收能力的评价,目前已应用于许多医学、非医学领域,如人类感知和决策研究、工业质量控制、军事监控等。ROC曲线上每个点反映着对同一信号刺激的感受性。 在诊断试验的评价研究中,它是以每一个检测结果作为可能的诊断界值，计算得到相应的真阳性率（true positive fraction，简称TPF）和假阳性率（false positive fraction，简称FPF）。是以假阳性率（即1-特异度）为横坐标、以真阳性率（即灵敏度）为纵坐标绘制而成的曲线，如图1所示，其能够直观的表明诊断试验的准确度,图1 诊断实验ROC曲线范例 对某个分类器而言，我们可以根据其在测试样本上的表现得到一个TPR和FPR点对。这样，这个分类器就可以映射成ROC平面上的一个点。调整这个分类器分类时候使用的阈值，我们就可以得到一个经过(0, 0)，(1, 1)的曲线，这就是此分类器的ROC曲线。一般情况下，这个曲线都应该处于(0, 0)和(1, 1)连线的上方。 针对一个二分类问题，将实例分成正类(postive)或者负类(negative)。实际中分类时，会出现四种情况：(1)若一个实例是正类，并且被预测为正类，即为真正类(True Postive TP)；(2)若一个实例是正类，但是被预测成为负类，即为假负类(False Negative FN)；(3)若一个实例是负类，但是被预测成为正类，即为假正类(False Postive FP)；(4)若一个实例是负类，但是被预测成为负类，即为真负类(True Negative TN)。由此，我们可得到表1。表1 ROC曲线分类表 这个表格中，1代表正数，0代表负数，TP为正确的肯定数目；FN为漏报，即没有找到正确匹配的数目，FP为误报，及没有匹配不正确，TN为正确拒绝的非匹配数目，由此得到计算公式如下：(1)真正类率(True Postive Rate)TPR:TP/(TP+FN),代表分类器预测的正类中实际正实例占所有正实例的比例。Sensitivity(2)负正类率(False Postive Rate)FPR:FP/(FP+TN)，代表分类器预测的正类中实际负实例占所有负实例的比例。1-Specificity(3)真负类率(True Negative Rate)TNR:TN/(FP+TN),代表分类器预测的负类中实际负实例占所有负实例的比例，TNR=1-FPR。Specificity 由计算公式，我们可求出各分类器预测的指，并画出曲线。ROC曲线中，较为重要的四个点： 第一个点，(0,1)，FPR=0，TPR=1，即误纳率为0，误拒率为0，就是负类样本中被判断为正类的比例为0，说明负类样本都被判断为负类，正类样本都被判断正确，所以这是一个完美的分类器，它将所有的样本都正确分类。 第二个点，(1,0)，FPR=1，TPR=0，把该判断为正类的判断为负类，把该判断为负类的判断为正类，是一个最糟糕的分类器，因为他成功的避开了所有的正确分类。 第三个点，(0,0)，FPR=TPR=0，该分类器预测所有的样本都为负样本。 第四个点（1,1），分类器实际上预测所有的样本都为正样本。4.2 AUC曲线 虽然ROC曲线较为直观，但人们仍希望有一个数值来标志分类器的好坏，于是引入AUC曲线。他是用来度量分类模型好坏的一个标准。AUC的值就是处于ROC curve下方的那部分面积的大小。 理论上，当诊断实验完全无诊断价值时，ROC曲线为一条从原点到右上角的对角线，即线段（0，0）（1,1），这条线称为机会线，若获得的ROC曲线落在机会线上，则曲线下面积为0.5；此时诊断无价值。理想的诊断实验ROC曲线应是从原点垂直上升至左上角，然后水平到达右上角，其曲线下面积为1，是最好的诊断结果，也是理想状态下的结果。 通常，AUC的值介于0.5到1.0之间，AUC值越大，越接近于1，说明诊断效果越好，AUC在0.70.9时有一定准确性，AUC在 0.50.7时有较低准确性，AUC=0.5时，说明诊断方法完全不起作用，无诊断价值。AUC0.5不符合真实情况，在实际中极少出现。 AUC值的计算分为参数法和非参数法，均适用于结果为连续性资料或等级资料的诊断试验准确度的评价，但计算较为复杂，大多需要借助统计软件来实现。本文主要介绍非参数法： 非参数法是根据诊断试验的检测结果直接计算绘制曲线所需的工作点（真阳性率,假阳性率）,由此绘制的ROC曲线称为ROC经验曲线,其曲线下面积可由梯形规则计算得到；Bamber于1975年发现：ROC经验曲线下面积等价于患者组和非患者组实验结果秩和检验的Wilcoxon Mann-Whitney检验统计量，因而可由Wilcoxon Mann-Whitney统计量估计曲线下面积的大小。假设患者组有nj个检测结果,记为xi(i=1,2,.,nj)，非患者组有n0个检测结果，记为yj（j=1,2,.,n0）。如果检测结果越大为患者的可能性越大,则根据Wilcoxon Mann-Whitney统计量，ROC曲线下面积就是患者组每个检测结果大于非患者组每个检测结果的概率，非参数法面积（A）估计的公式可表示为：5. 实验和结果5.1原始数据研究样本为美国车辆保险机构提供的1994年的车险索赔数据，一共6142条记录，33个变量。其中欺诈样本有409条记录，占总样本的6.66%，所以是非平衡数据。自变量31个，1个应变量（0是诚实索赔，1是欺诈索赔），1个控制变量年份。5.2数据处理过程：首先在EXCEL中把分类数据变成连续性数值（即分水平）Month事故发生的月份Month2WeekOfMonth第几周DayOfWeek周几DayOfWeek2R里默认周日为1Make车的品牌Make2Honda=1;Toyota=2;Ford=3;Mazda=4;Chevrolet=5;Pontiac=6;Accura=7;Dodge=8;Mercury=9;Jaguar=10;BMW=11;Ferrari=12;Lexus=13;Mecedes=14;Nisson=15;Porche=16;Saab=17;Saturn=18;VW=19AccidentArea事故发生地点AccidentArea2Urban=1；Rural=2DayOfWeekClaimed提出索赔是周几DayOfWeekClaimed2R里默认周日为1MonthClaimed提出索赔的月份MonthClaimed2WeekOfMonthClaimed提出索赔是第几周Sex性别Sex2Male=1；Female=2MaritalStatus婚姻状态MaritalStatus2Single=1；Married=2；Divorced=3;Widow=4Age年龄0-80把0变成17，未成年人Fault过失Fault2Policy Holder=1；Third Party=2PolicyType保险类型PolicyType2Sedan -All Perils=1；Sedan -Collision=2；Sedan - Liability=3；Sport-All Perils=4；Sport - Collision=5；Sport - Liability=6；Utility - All Perils=7; Utility - Collision=8;Utility - Liability=9.VehicleCategory汽车的类型VehicleCategory2Sedan=1；Sport=2；Utility=3VehiclePrice汽车价格区间VehiclePrice2less than 20000=0；20000 to 29000=1；30000 to 39000=2；40000 to 59000=3；60000 to 69000=4；more than 69000=5FraudFound是否为欺诈欺诈=1非欺诈=0PolicyNumber保险单号，相当于ID1-15420RepNumber推销员工号1-16Deductible免赔额300,400,500,700DriverRating驾驶员等级1-4Days:Policy-Accident下保单到事故发生的间隔Days:Policy-Accident2None=0; 1 to 7=1; 8 to 15 =2; 15 to 30 =3; More than 30=4Days:Policy-Claim下保单到索赔发生的间隔Days:Policy-Claim2None=0; 8 to 15 =1; 15 to 30 =2;More than 30=3PastNumberOfClaims过去索赔的次数PastNumberOfClaims2None=0; 1=1; 2 to 4 =2; More than 4=3AgeOfVehicle汽车使用年数AgeOfVehicle2New = 0; 2 years =1; 3 years =2;4 years =3; 5 years =4; 6 years =5;7 years =6; More than 7 = 7AgeOfPolicyHolder投保人的年纪区间AgeOfPolicyHolder216 to 17 =1；18 to 20 =2; 21 to 25 =3;26 to 30 =4; 31 to 35 =5; 36 to 40 =6;41 to 50 =7; 51 to 65 =8; over 65 =9PoliceReportFiled警方报道是否存档PoliceReportFiled2No=1; Yes=2WitnessPresent是否有目击证人WitnessPresent2No=1; Yes=2AgentType代理商类型AgentType2Internal=1; External=2NumberOfSuppliments2补充的代号NumberOfSuppliments2none=0; 1 to 2=1; 3 to 5=2;more than 5=3.AddressChange-Claim索赔和地点变更的间隔年数AddressChange-Claim2no change =0; under 6 months=11 year =2; 2 to 3 years =3;4 to 8 years=4NumberOfCars拥有汽车的数量NumberOfCars21 vehicle=1； 2 vehicles=2； 3 to 4=3；5 to 8=4；more than 8=5Year年1994-19961994:1-61421995:6143-113371996:11338-15420BasePolicy基础保险政策BasePolicy2All Perils=1； Collision=2；Liability=3删去了BasePolicy2,因为是PolicyType的基础.剩下32个变量，即黄色标注的变量。利用自编语句，加上DayOfAcci和DayOfClaim, Claim-accident（事故发生到索赔的时间间隔）和Isdriver（是否为投保人）=36个变量（请见附录）减去以下变量=22个变量，进行实验1Month2事故日期2WeekOfMonth3DayOfWeek24DayOfWeekClaimed2索赔日期5MonthClaimed26WeekOfMonthClaimed7PolicyNumber保单的顺序8DayOfAcci使用索赔-事故的时间间隔即可9DayOfClaim10Make2车的品牌影响不大11PolicyType2保单类型12AddressChange-Claim2地点是否变更13X序号14Year1994年的自检验不需要年份5.3实验过程首先利用随机森林进行实验。具体语句如下：#先将数据随机分成训练集和测试集mydata1994-read.csv(Pyle1994.csv,header=TRUE)b-c(which(names(mydata1994)=X),which(names(mydata1994)=Year)Pyle1994-mydata1994,-bnames(Pyle1994)fix(Pyle1994)ind-sample(c(1,2),nrow(Pyle1994),replace=T,prob=c(0.7,0.3)traindata-Pyle1994ind=1,testdata-Pyle1994ind=2,library(pROC) #安装roc包#使用随机森林library(randomForest)rf-randomForest(traindata$FraudFound.P.,data=traindata,ntree=200)rftestpre-predict(rf,newdata=testdata,type=class)table(testpre,testdata$FraudFound.P.)fity-as.numeric(testpre)rf.roc-roc(testdata$FraudFound.P.,fity)rf.roc结果如图所示：利用Logistic算法进行运算。具体语句为：#使用逻辑回归train.log-glm(FraudFound.P.,data=traindata,family=binomial(link=logit)testpredlog0.5)fitylog-as.numeric(testpredlog)log.roc-roc(testdata$FraudFound.P.,fitylog)Log.roc结果如图：5.4 试验结果将两种算法的ROC曲线进行对比：由以上对比可以清楚地看出来：在非平衡大数据的分类问题中，随机森林表现比Logistic好。参考文献1中国经济网.车险成保险欺诈易发领域中保协加强预防打击(2014-05-21)./14/0521/07/9SOKLFJK00253B0H.html#from=relevant#xwwzy_35_bottomnewskwd.2北京日报.去年全国车险欺诈案2375 起（2014-04-17). /14/0417/ 03/9Q0LMGLN00014AED.html.3叶明华.基于BP神经网络的保险欺诈识别研究J.保险研究.2011（3）：79-864范广哲，刘瑾雯.非平衡大数据下的保险欺诈识别研究J.中国新通信，2015（01）：98-995 廖晓阳, 金成勇. 如何在诊断性试验研究中正确应用有关测量指标J. 华西医学. 2000, 15(2):1466 Metz CE, Herman BA, Shen J-H, Maximum-likelihood estimation of receiver operating characteristic(ROC) curves from continuously-distributed data. Statistics in Medicine, 19987 宋华玲, ROC曲线的评价研究及应用D. 第二军医大学，2006附录：mydata-read.csv(C:/Users/duancaiwei/Desktop/car claim/DATA/Pyle_numori.csv,header=T)dim(mydata)names(mydata)#先对变量进行转换和增加#将三个日期相关变量转变成日期，求出事故日期和索赔日期，并求二者之差for (i in 1:nrow(mydata) if (mydatai,1 %in% c(1,3,5,7,8,10,12) start1-paste(as.character(mydatai,32),as.character(mydatafi,1),as.character(01),sep = -) start1day-as.Date(start1) end1-paste(as.character(mydatafi,32),as.character(mydatafi,1),as.character(31),sep = -) end1day-as.Date(end1) if (mydatai,1 %in% c(4,6,9,11) start1-paste(as.character(mydatai,32),as.character(mydatai,1),as.character(01),sep = -) start1day-as.Date(start1) end1-paste(as.character(mydatai,32),as.character(mydatai,1),as.character(30),sep = -) end1day-as.Date(end1) if (mydatai,1 =2) start1-paste(as.character(mydatai,32),as.character(mydatai,1),as.character(01),sep = -) start1day-as.Date(start1) end1-paste(as.character(mydatai,32),as.character(mydatai,1),as.character(28),sep = -) end1day-as.Date(end1) library(lubridate) seqDate-seq.Date(ymd(start1day),ymd(end1day),by=day) a-seqDatewhich(wday(seqDate)=as.character(mydatai,3)mydatai,2#输出NA，因为第五周，所以得减去1 b-seqDatewhich(wday(seqDate)=as.character(mydatai,3)mydatai,2-1 mydatai,33-ifelse(a %in% seqDate,a,b) if(i=nrow(mydata) mydata$DayOfAcci - as.Date(mydata$DayOfAcci,origin = 1970-01-01) #到此求出了事故日期for (i in 1:nrow(mydata) if (mydatai,7 %in% c(1,3,5,7,8,10,12) start1-paste(as.character(mydatai,32),as.character(mydatafi,7),as.character(01),sep = -) start1day-as.Date(start1) end1-paste(as.character(mydatafi,32),as.character(mydatafi,7),as.character(31),sep = -) end1day-as.Date(end1) if (mydatai,7 %in% c(4,6,9,11) start1-paste(as.character(mydatai,32),as.character(mydatai,7),as.character(01),sep = -) start1day-as.Date(start1) end1-paste(as.character(mydatai,32),as.character(mydatai,7),as.character(30),sep = -) end1day-as.Date(end1) if (mydatai,7 =2) start1-paste(as.character(mydatai,32),as.character(mydatai,7),as.character(01),sep = -) start1day-as.Date(start1) end1-paste(as.character(mydatai,32),as.character(mydatai,7),as.character(28),sep = -) end1day-as.Date(end1) library(lubridate) seqDate-seq.Date(ymd(start1day),ymd(end1day),by=day) a-seqDatewhich(wday(seqDate)=as.character(mydatai,6)mydatai,8#输出NA，因为第五周，所以得减去1 b-seqDatewhich(wday(seqDate)=as.character(mydatai,6)mydatai,8-1 mydatai,34-ifelse(