《数学建模：基于R》一一1.7 数学建模案例分析——食品质量安全抽检数据分析...

本节书摘来自华章计算机《数学建模：基于R》一书中的第1章，第1.7节,作者：薛　毅 更多章节内容可以访问云栖社区“华章计算机”公众号查看。

1.7　数学建模案例分析——食品质量安全抽检数据分析

1.7.1　问题的提出　该题选自2013年“深圳杯”数学建模夏令营A题.
“民以食为天”,食品安全关系到千家万户的生活与健康.随着人们对生活质量的追求和安全意识的提高,食品安全已成为社会关注的热点,也是政府民生工程的一个主题.城市食品的来源越来越广泛,人们消费加工好的食品的比例也越来越高,因此除食材的生产收获外,食品的运输、加工、包装、贮存、销售以及餐饮等每一个环节皆可能影响食品的质量与安全.另一方面,食品质量与安全又是一个专业性很强的问题,其标准的制定和抽样检测及评价都需要科学有效的方法.深圳是食品抽检、监督最统一、最规范、最公开的城市之一.请下载2010年、2011年和2012年深圳市的食品抽检数据（注意蔬菜、鱼类、鸡鸭等抽检数据的获取）　数据下载网站：www.szaic.gov.cn（深圳市市场监督管理局网站）.并根据这些资料来讨论：
1） 如何评价深圳市这三年各主要食品领域微生物、重金属、添加剂含量等安全情况的变化趋势.
2） 从这些数据中能否找出某些规律性的东西：如食品产地与食品质量的关系；食品销售地点（即抽检地点）与食品质量的关系；季节因素等.
3） 能否改进食品抽检的办法,使之更科学更有效地反映食品质量状况且不过分增加监管成本（食品抽检是需要费用的）,例如对于抽检结果稳定且抽检频次过高的食品领域该作怎样的调整？
解决上述问题的最大困难是数据整理,因为这三年中的数据格式不统一（三年的Excel表的格式不同）,文件不统一（有Excel表和Word文档）,有的表只有不合格食品,没有说明本次抽检共抽检了多少食品.这些都给数据分析带来困难,需要花大量的时间来分析整理　本节的数据是根据2013年“深圳杯”数学建模夏令营的几篇优秀论文整理得到的..
1.7.2　问题1：三年各主要食品领域安全情况的变化趋势
对深圳市2010年至2012年的所有能够得到的不合格数据进行汇总,并将检测指标大致分为以下四类：微生物、重金属、食品添加剂和其他类,汇总结果如表1.18所示.
表1.18　2010—2012年各类主要食品抽检情况汇总表
年度抽检总数抽检不合格种类微生物重金属添加剂其他类总数
2010 2261 60 7 31 40 138
2011 7814 250 264 239 79 832
2012 10607 116 35 29 78 258
食品抽检只有两种结果,一是合格,二是不合格.因此,选择二项分布的方法进行分析,这里采用prop.test()函数,pairwise.prop.test()函数和prop.trend.test()函数完成安全情况分析.
用prop.test()函数检验这三年抽检不合格的比率是否相同,即检验H0:p1=p2=p3,　　H1:p1,p2,p3不全相同如果拒绝原假设,则表明这三年的不合格比率有显著差异.
用pairwise.prop.test()函数作多重比率检验,即检验谁与谁有显著差异,该函数的使用格式为pairwise.prop.test(x, n,
　　p.adjust.method = p.adjust.methods, ...)参数x为整数向量,表示试验成功的次数,或者为一个2列的矩阵,第1列表示成功的次数,第2列表示失败的次数.
n为整数向量,表示试验的次数,当x为矩阵时,该值无效.
p.adjust.method为P值的调整方法(细节见2.2.2节).
...为附加参数.
用prop.trend.test()函数作比率趋势性检验,即检验H0:p1=p2=…=pk,　　H1:p1≤p2≤…≤pk或者检验H0:p1=p2=…=pk,　　H1:p1≥p2≥…≥pk其使用格式为prop.trend.test(x, n, score = seq_along(x))参数x为事件数.n为试验次数.score为分组得分,默认值为自然顺序.
将数据按表格形式存入数据文件(文件名:summer_1.data),在表头中,sampling表示抽检总数,microorganism表示微生物方面的食品不合格数,metal表示重金属方面的食品不合格数,additive表示添加剂含量方面的食品不合格数,others表示其他方面的食品不合格数,unqualified表示食品不合格的总数.用read.table()函数读取表格形式的数据文件.
(1) 关于微生物方面的食品安全情况的变化趋势> rt <- read.table("summer_1.data")

prop.test(x = rt$microorganism, n = rt$sampling)

　　3-sample test for equality of proportions without
　　continuity correction
data: rt$microorganism out of rt$sampling
X-squared = 103.3388, df = 2, p-value < 2.2e-16
alternative hypothesis: two.sided
sample estimates:
prop 1prop 2prop 3
0.026536930.031993860.01093617拒绝原假设,说明关于微生物方面的食品安全情况这三年是有差异的.> pairwise.prop.test(x = rt$microorganism, n = rt$sampling)

　　Pairwise comparisons using Pairwise comparison of
　　proportions
data: rt$microorganism out of rt$sampling
12
20.21-
32.4e-08< 2e-16
P value adjustment method: holm 从计算结果