本节书摘来自华章计算机《数学建模：基于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 从计算结果来看,2010年与第2012年,2011年与2012年微生物方面的食品安全情况有显著差异.> prop.trend.test(x = rt$microorganism, n = rt$sampling)
Chi-squared Test for Trend in Proportions
data: rt$microorganism out of rt$sampling,
using scores: 1 2 3
X-squared = 70.1004, df = 1, p-value < 2.2e-16从趋势检验来看,拒绝原假设,说明微生物方面的食品不合格率在逐年减少.
(2) 关于重金属方面的食品安全情况的变化趋势> prop.test(x = rt$metal, n = rt$sampling)
3-sample test for equality of proportions without
continuity correction
data: rt$metal out of rt$sampling
X-squared = 310.7125, df = 2, p-value < 2.2e-16
alternative hypothesis: two.sided
sample estimates:
prop 1prop 2prop 3
0.0030959750.0337855130.003299708 拒绝原假设,说明关于重金属方面的食品安全情况这三年是有差异的.> pairwise.prop.test(x = rt$metal, n = rt$sampling)
Pairwise comparisons using Pairwise comparison of
proportions
data: rt$metal out of rt$sampling
12
27.1e-15 -
31< 2e-16
P value adjustment method: holm 从计算结果来看,2010年与第2011年,2011年与2012年重金属方面的食品安全情况有显著差异.> prop.trend.test(x = rt$metal, n = rt$sampling)
Chi-squared Test for Trend in Proportions
data: rt$metal out of rt$sampling,
using scores: 1 2 3
X-squared = 65.8371, df = 1, p-value = 4.898e-16从趋势检验来看,拒绝原假设,说明重金属方面的食品不合格率在逐年减少.
(3) 关于添加剂含量方面的食品安全情况的变化趋势> prop.test(x = rt$additive, n = rt$sampling)
3-sample test for equality of proportions without
continuity correction
data: rt$additive out of rt$sampling
X-squared = 245.0698, df = 2, p-value < 2.2e-16
alternative hypothesis: two.sided
sample estimates:
prop 1prop 2prop 3
0.0137107470.0305861270.002734044拒绝原假设,说明添加剂含量方面的食品安全情况这三年是有差异的.> pairwise.prop.test(x = rt$additive, n = rt$sampling)
Pairwise comparisons using Pairwise comparison of
proportions
data: rt$additive out of rt$sampling
12
21.7e-05-
32.3e-11< 2e-16
P value adjustment method: holm 从计算结果来看,2010年与第2011年,2010年与第2012年,2011年与2012年添加剂含量方面的食品安全情况有显著差异.> prop.trend.test(x = rt$additive, n = rt$sampling)
Chi-squared Test for Trend in Proportions
data: rt$additive out of rt$sampling,
using scores: 1 2 3
X-squared = 111.1509, df = 1, p-value < 2.2e-16从趋势检验来看,拒绝原假设,说明添加剂含量方面的食品不合格率在逐年减少.
1.7.3 问题2：找出某些规律性的东西
1.食品生产到销售各环节与食品质量的关系
食品在生产到销售的过程中,大致可分为:餐饮、储存、生产、使用、运输和销售等环节,将2010～2012年食品在各环节的抽检情况汇总(见表1.19).
2.食品产地与食品质量的关系
不同食品产地对食品质量的影响主要体现在三个方面：运输距离的长短,不同食品产地本身的质量和深圳市对不同产地产品流入市场时的监管力度.
在数据处理中,由于部分产地的检测数目较小,该省份或城市的合格率会产生较大的偶然性,因此统一剔除检测数目小于50的产地.筛选后所得的36个食品产地与合格率信息,如表1.20所示.
考查食品产地与食品质量之间的关系,可从近距离(800km以内)和远距离(800km以上)两方面考虑.考查的指标是食品产地到销售地的距离和食品抽检的不合格率.由于这两项指标的数据不一定满足正态分布,所以这里采用秩相关检验——Spearman相关检验.
将数据按表格形式存入数据文件(文件名:summer_22.data),在表头中,distance表示食品产地距深圳市的距离,sampling表示抽检总数,unqualified表示食品不合格数.以下是R程序(文件名:summer_22.R)及计算结果.
先计算近距离(800km以内)的食品安全情况：rt <- read.table("summer_22.data")
rt$ratio <- rt$unqualified / rt$sampling
cor.test(~ distance + ratio, data = rt, subset = 1:18,
method = "spearman")
Spearman's rank correlation rho
data: distance and ratio
S = 804.9153, p-value = 0.5018
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.169334接受原假设,说明食品产地与食品质量是独立的.
再计算远距离(800km以上)的食品安全情况：cor.test(~ distance + ratio, data = rt, subset = 19:36,
method = "spearman")
Spearman's rank correlation rho
data: distance and ratio
S = 1256.987, p-value = 0.231
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.2972007接受原假设,说明食品产地与食品质量是独立的.
将近距离和远距离的食品放在一起讨论：cor.test(~ distance + ratio, data = rt,
method = "spearman")
Spearman's rank correlation rho
data: distance and ratio
S = 11417.17, p-value = 0.003876
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.4693908拒绝原假设,说明食品产地与食品质量是有关的.从相关系数来看,相关系数的符号为“-”,也就是说,食品产地越远,食品抽检的不合格率越小,食品质量越高.这一结论与我们的生活常识似乎是矛盾的,应该是产地越远,食品抽检的不合格率越高.
仔细分析,上述计算得到的结论是合理的,因为产地越远,食品的包装就会越好,运输条件也相对近距离要好,否则到了销售地后,食品的包装损坏严重,会给生产厂家带来不必要的损失.
3.食品销售地点（即抽检地点）与食品质量的关系
为了研究这个问题,将深圳市按照城区划分,各个城区的合格数与不合格数如表1.21所示.
表1.21 2010—2012年各地区食品安全情况汇总表
抽检地区
宝安 福田 光明 龙岗 罗湖 南山 坪山 盐田抽检总数 5773 2682 690 4899 2754 2893 734 986
不合格数 293 128 67 256 113 107 32 53
作列联表检验,检验食品质量与抽检地是否独立.以下是程序(文件名:summer_23.R)和计算结果：x <- scan()
5773 2682 690 4899 2754 2893 734 986
293 128 67 256 113 107 32 53
X <- matrix(c(x［9:16］, x［1:8］-x［9:16］), nr = 2, byrow = T)
chisq.test(X)
Pearson's Chi-squared test
data: X
X-squared = 49.5067, df = 7, p-value = 1.805e-08拒绝原假设,说明食品的抽检地点与食品质量是有关系的.这个结论或许与每个地区的食品来源有关.
也可以作比率检验(prop.test()函数),其结果是相同的,这一检验留给读者完成.
4.季节与食品质量的关系
为了研究季节与食品质量的关系,整理出饮品、果蔬、粮食、焙烤、肉蛋和水产品共6类主要食品的抽检情况,如表1.22所示.
表1.22 六类主要食品分季度汇总表
第一季度 第二季度 第三季度 第四季度
合格 不合格 合格 不合格 合格 不合格 合格 不合格饮品 777 15 254 7 1225 29 988 25
果蔬 861 33 275 29 518 25 1049 50
粮食 809 47 721 32 1171 16 1156 36
焙烤 476 130 205 17 989 18 491 31
肉蛋 426 62 391 36 280 16 766 76
水产 98 2 222 10 373 22 330 18
将数据以表格形式存放在文本文件(summer_24.data)中,使用列联表作独立性检验.
(1) 关于饮品类食品的检验情况rt <- read.table("summer_24.data")
X <- matrix(as.numeric(rt［1, ］), nc = 4)
chisq.test(X)
Pearson's Chi-squared test
data: X
X-squared = 0.8809, df = 3, p-value = 0.83程序中的as.numeric()函数是将其他类型数据强制转换成数值型.计算结果是,接受原假设,即饮品类食品的安全情况与季节无关.
(2) 关于果蔬类食品的检验情况X <- matrix(as.numeric(rt［2, ］), nc = 4)
chisq.test(X)
Pearson's Chi-squared test
data: X
X-squared = 17.4588, df = 3, p-value = 0.0005687拒绝原假设,即果蔬类食品的安全情况与季节有关.
(3) 关于粮食类食品的检验情况X <- matrix(as.numeric(rt［3, ］), nc = 4)
chisq.test(X)
Pearson's Chi-squared test
data: X
X-squared = 29.5963, df = 3, p-value = 1.678e-06拒绝原假设,即粮食类食品的安全情况与季节有关.
(4) 关于焙烤类食品的检验情况X <- matrix(as.numeric(rt［4, ］), nc = 4)
chisq.test(X)
Pearson's Chi-squared test
data: X
X-squared = 197.4468, df = 3, p-value < 2.2e-16拒绝原假设,即焙烤类食品的安全情况与季节有关.
(5) 关于肉蛋类食品的检验情况X <- matrix(as.numeric(rt［5, ］), nc = 4)
chisq.test(X)
Pearson's Chi-squared test
data: X
X-squared = 12.5369, df = 3, p-value = 0.005753拒绝原假设,即肉蛋类食品的安全情况与季节有关.
(6)关于水产类食品的检验情况X <- matrix(as.numeric(rt［6, ］), nc = 4)
fisher.test(X)
Fisher's Exact Test for Count Data
data: X
p-value = 0.5258
alternative hypothesis: two.sided接受原假设,即水产类食品的安全情况与季节无关.注:由于第一季度的不合格数为2,小于5,因此,使用Fisher精确检验.
1.7.4 问题3：如何改进食品的抽检办法
这个问题实际上相当的复杂,改进抽检办法的目的就是在保证抽检质量的前提下,减少抽检次数,从而降低抽检成本.这里涉及抽检方法、抽检成本,以及优化等方面的知识.
本小节以简单随机抽样为基础,简单计算出抽样的样本容量.
设y1,y2,…,yn为随机抽样的样本,y=∑ni=1yi为样本均值,N为总体的总数,f=n/N为抽样比.由抽样理论可知,var(y)=1－fnS2=1n－1NS2(1.79)其中S2为总体的方差.由式(1.79)得到1n=1N+var(y)S2(1.80)式(1.80)虽然明确地给出了样本容量n的计算公式,但式中的总体方差S2是未知的,而抽样误差var(y)也与总体方差S2有关.
通常人们用置信度1－α和绝对误差限y－Y≤d来计算抽样的样本容量,即考虑P{y－Y≤d}=1－α（1.81）其中Y为总体的均值.由式(1.81)和双侧分位点uα/2的定义分别得到Py－Yvar(y)≤dvar(y)=1－α（1.82）
Py－Yvar(y)≤uα/2=1－α（1.83）对比式(1.82)和式(1.83)得到,uα/2=d/var(y),即var(y)=d2/u2α/2（1.84）将式(1.84)代入式(1.85),得到1n=1N+d2u2α/2S2（1.85）如果考虑置信度和相对误差Py－YY≤r=1－α（1.86）类似前面的推导,可以得到1n=1N+r2Y2u2α/2S2（1.87）对于食品的安全检查,人们最关心的是食品的不合格率p,因此,用p替换Y,并令S2=pq,其中q=1－p.代入式(1.87)得到1n=1N+r2pu2α/2q（1.88）在实际应用中,通常取α=0.05,用正态分布的分位点近似uα/2.如果某类食品的总数N=10000,食品的不合格率p=0.05,检验的相对误差r=0.1,即10%,按式(1.88)计算得到n=4220约为总数的一半.
这个抽检次数远远大于目前的抽检次数,因此,只使用简单随机抽检的方法是不能减少抽检次数的.关于其他抽检方法的讨论已超出本书的内容,所以不再进行讨论.
1.7.5 结论
(1) 深圳市2010年至2012年三年来的食品安全状况逐年变好;
(2) 食品的餐饮、储存等各环节与食品质量有关;
(3) 食品产地与食品质量有关,而且是负相关;
(4) 食品的抽检地点与食品质量有关;
(5) 饮品和水产类食品的质量与季节无关,果蔬、粮食、焙烤和肉蛋类食品的质量与季节有关.
习题1
求学生身高的平均值、中位数、方差、极差和四分位极差.
2.假定某地区成年男子的身高(单位:cm)服从正态分布X~N(170,7.692),求该地区成年男子身高超过175cm的概率.如果希望有95%以上的人在通过门时不低头,门的高度至少应为多少?
3.绘出习题1中学生身高数据的直方图、核密度估计曲线.(要求:直方图高度为频率/间距,并将核密度估计曲线与直方图画在一张图上).
4.绘出习题1中学生身高数据的正态Q-Q图,说明该数据是否接近正态分布.
5.正常男子血小板计数均值为225×109/L,今测得20名男性油漆作业工人的血小板计数值(单位:109/L)如下：
220 188 162 230 145 160 238 188 247 113
126 245 164 231 256 183 190 158 224 175
假定数据服从正态分布,问油漆工人的血小板计数与正常成年男子有无显著差异?
6.已知某种灯泡寿命服从正态分布,在某星期所生产的该灯泡中随机抽取10只,测得其寿命(单位:h)为
1067 919 1196 785 1126 936 918 1156 920 948
求这个星期生产出的灯泡能使用1000h以上的概率.
7.为研究国产四类新药阿卡波糖胶囊效果,某医院用40名Ⅱ型糖尿病病人进行同期随机对照实验.试验者将这些病人随机等分到试验组(阿卡波糖胶囊组)和对照组(拜唐苹胶囊组),分别测得试验开始前和8周后空腹血糖,计算得到空腹血糖下降值如表1.23所示.假设数据服从正态分布,试用t检验(讨论方差相同和方差不同两种情况)和成对数据的t检验来判断:国产四类新药阿卡波糖胶囊与拜唐苹胶囊对空腹血糖的降糖效果是否相同.分析三种检验方法哪一种最好,并说明理由.
表1.23 实验组与对照组空腹血糖下降值
试验组-0.70-5.602.002.800.703.504.005.807.10-0.50
2.50-1.601.703.000.404.504.602.506.00-1.40
对照组3.706.505.005.200.800.200.603.406.60-1.10
6.003.802.001.602.002.201.203.101.70-2.00
8.一项调查显示某城市老年人口比重为14.7％.该市老年研究协会为了检验该项调查是否可靠,随机抽选了400名居民,发现其中有57位是老年人.问调查结果是否支持该市老年人口比重为14.7％的看法(α=0.05).
9.研究者发现,妇女患乳腺癌可能与初次分娩时的年龄有关,表1.24给出国际卫生组织在1970年的报告.试分析初次分娩妇女各年龄段患乳腺癌的比例是否相同.
表1.24 初次分娩年龄与乳腺癌患病人数
初次分娩年龄<2020~2425~2930~34≥3
5乳腺癌数32012061011463220
调查总数1742563839041555626
10.某项调查询问了2000名青年人.问题是:“你认为我们的生活环境比过去更好、更差,还是没有变化?”有800人觉得“越来越好”,有720人感觉“一天不如一天”,有400人表示“没有变化,一直如此”,还有80人说不知道.试分析:在总体中是否可以认为“我们的生活环境比过去更好”的人比“我们的生活环境比过去更差”的人多?
11.在某养鱼塘中,根据过去经验,鱼的长度的中位数为14.6cm,现对鱼塘中鱼的长度进行一次估测,随机地从鱼塘中取出10条鱼，长度如下：
13.32 13.06 14.02 11.86 13.58
13.77 13.51 14.42 14.44 15.43
将它们作为一个样本进行检验.试分析,该鱼塘中鱼的长度是在中位数之上,还是在中位数之下？(1)用符号检验;(2)用Wilcoxon符号秩检验.
12.考虑例1.19若新方法与原方法得到排序结果改为表1.25所示的情形,能否说明新方法比原方法显著提高了教学效果?
表1.25 学生数学能力排序结果
新方法467910
原方法12358
13.为了比较一种新疗法对某种疾病的治疗效果,将40名患者随机地分为两组,每组20人,一组采用新疗法,另一组采用原标准疗法.经过一段时间的治疗后,对每个患者的疗效作仔细地评估,并划分为差、较差、一般、较好和好五个等级.两组中处于不同等级的患者人数如表1.26所示.试分析,由此结果能否认为新方法的疗效显著地优于原疗法(α=0.05)?
表1.26 不同方法治疗后的结果
等级差较差一般较好好
新疗法组01973
原疗法组221141
14.据经验知,机床发生故障的概率服从均匀分布,某车间在一周内统计所有车床发生故障的频数资料如表1.27所示,试分析,这组数据是否服从均匀分布?
表1.27 某车间在一周内所有车床发生故障的频数资料
星期一二三四五六
频数78391617
- Mendel用豌豆的两对相对性状进行杂交实验,黄色圆滑种子与绿色皱缩种子的豌豆杂交后,第二代根据自由组合规律,理论分离比为黄圆∶黄皱∶绿圆∶绿皱=916∶316∶316∶116实际实验值为:黄圆315粒,黄皱101粒,绿圆108粒,绿皱32粒,共556粒,问此结果是否符合自由组合规律?
16.在一次实验中,每隔一定的时间观察一次由某种铀所放射到达计数器的α粒子数X,共观察了100次,其数据如表1.28所示.试分析,能否认为α粒子数X服从Poisson分布.
表1.28 观察到的α粒子数
粒子数01234567891011
频数 1516172611992121
17.一般认为长途电话通过电话总机的过程是一个随机过程,其间打进电话的时间间隔服从指数分布.某个星期下午1:00以后最先打进的10个电话的时间为
1:06 1:08 1:16 1:22 1:23 1:34 1:44 1:47 1:51 1:57
试用Kolmogorov-Smirnov检验分析打进电话的时间间隔是否服从指数分布.
18.试用Shapiro-Wilk检验判断习题1中的身高数据是否服从正态分布.
19.在高中一年级男生中抽取300名考察其两个属性:B是1500米长跑,C是每天平均锻炼时间,得到4×3列联表,如表1.29所示.试对α=0.05,检验B与C是否独立.
表1.29 300名高中学生体育锻炼的考察结果1500米长跑记录锻炼时间2小时以上1~2小时1小时以下合计5″01′~5″30′451210675″31′~6″00′462028946″01′~6″30′282330816″31′~7″00′11123558合计13067103300
20.为研究分娩过程中使用胎儿电子监测仪对剖腹产率有无影响,对5824例分娩的经产妇进行回顾性调查,结果如表1.30所示,试进行分析.
表1.30 5824例经产妇回顾性调查结果剖腹产胎儿电子监测仪使用未使用合计是358229587否249227455237合计285029745824
21.为研究人脑的左右半球恶性肿瘤的发病率是否有显著差异,对人脑恶性肿瘤和良性肿瘤的发病情况作了调查,调查结果如表1.31所示.试进行分析.
表1.31 人脑左右半球恶性肿瘤和良性肿瘤的发病情况
良性 恶性 合计左半球9312
右半球134
合计10616
22.表1.32列出某高中18名学生某门课程的高考成绩和模拟考试成绩,这组数据能否说明:高考成绩与模拟考试成绩是相关的？
表1.32 高考成绩和模拟考试成绩
学号 高考成绩 模拟成绩 学号 高考成绩 模拟成绩 学号 高考成绩 模拟成绩1 87 90 7 78 65 13 90 100
2 76 98 8 91 90 14 92 97
3 77 92 9 76 84 15 100 97
4 85 87 10 100 92 16 100 95
5 89 87 11 96 100 17 90 94
6 83 62 12 96 98 18 99 100
23.调查某大学学生每周学习时间(单位:h)与得分的平均等级之间关系,现抽查10名学生的资料如表1.33所示.表中等级10表示最好,1表示最差.试用秩相关检验(Spearman检验和Kendall检验)分析学习时间与学习等级有无关系.
表1.33 每周学习时间与学习等级(单位:h)
学习时间24172041522346181529
学习等级81479510326
