1963-1992年美国香烟消费数据研究
李雅茜 218025200067 应用统计
摘要
本文通过对1963-1992年美国香烟消费数据进行研究，进行探索性分析，寻找影响香烟销售量的重要因素，以及这些因素对香烟销售量所造成影响。然后通过对寻找到的时间和地区这两个重要影响因素建立线性回归模型，来对香烟销售量数据进行拟合，并对拟合结果进行评估与可视化。最后，对研究结果进行总结与反思，依据美国控烟成功的经验，针对我国的控烟问题提出了相应建议。
关键词：美国；香烟消费；香烟税；线性模型；
1. 绪论
1.1 研究背景
总所周知，吸烟有害健康。吸烟对肺的伤害最大，原因在于它是烟草吸入体内后较先到达的器官。其实，吸烟不仅伤肺，它对人体的伤害几乎是“全方位”的：呼吸道系统、心脑血管系统、内分泌系统等，只要有烟草的“侵入”，就存在损伤的可能。有统计数字表明，全球每年死于烟草的人数高达500万，其中30%的癌症死亡，包括90%的肺癌死亡，都和吸烟有关。除了肺癌之外，还有16种癌症和吸烟直接有关。
2018年美国癌症统计年报显示，美国的癌症死亡率创历史新低，1991年到2015年14年间，癌症的死亡率整体下降了26%。美国癌症协会把这算成控烟成功的功劳。从1990年到2015年间，癌症之王肺癌的死亡率在男性中下降了45%。而美国的吸烟者比率从50年代近乎一半人口都吸烟的42.2%，降到了2014年的17.8%，而这个数字还在持续下降。从吸烟作为一种全民追逐的潮流，到全民抵制烟草，美国用了五十年时间。
但在仅仅半个世纪前，美国许多吸烟者甚至声称吸烟对他们的健康有益。二战之后的美国，吸烟是一种流行文化和社会风尚，40%的美国人经常吸烟，有三分之一的女性是烟民，香烟也非常便宜，一包不到0.5美元，谁都能买得起。当时美国香烟消费量达3820亿支一年，平均每个消费者每年要抽2262支之多。吸烟这种习惯跨越社会经济、性别、种族和种族界限。烟民在可以在医院、餐馆、飞机等各种各样的公共场所中吸，你每走几步路就可以看到赶时髦的夹着烟头的人。
从20世纪30年代起，美国学者和社会活动家已经意识到癌症死亡率的增加。50年代，科学界开始出现对吸烟和健康之间关系的真正的科学理性的研究。1964年美国联邦卫生局发布了第一份《吸烟与健康》报告，第一次以政府名义给出了“吸烟危害健康”的结论。1965年联邦香烟标签和广告法开始规定香烟包装上必须出现：“警告：吸烟会危害你的健康”。1969年的《公共卫生吸烟法》禁止烟草公司在1971年后在广播频道播出广告。但是，在整个20世纪60年代和70年代，人们对吸烟行为和尼古丁成瘾的普遍了解依旧非常有限。70年代吸烟率仍为40%左右。但政府开始逐渐意识到吸烟的后果不仅仅是私人的问题，它可能会影响每一个纳税人。在经济上，政府也逐渐“觉醒”，烟草行业的繁荣给政府带来的巨大收入，很可能会在未来给政府带来巨额账单。1976年《吸烟与健康》报告有了升级版，这一报告的重大意义在于，揭露了二手烟对非自愿吸烟者的健康伤害。80年代中期，几乎所有的州都颁布了限制人们在公共场所吸烟的条例。1988年，国会对美国所有2小时以内的国内航班实施禁烟令。两年后，这项禁令延长到6小时以内的航班，实际上禁止在国内所有航班上吸烟。90年代，FDA开展了一项大规模调查，调查烟草企业操纵香烟中的尼古丁含量，和它们向年轻人推销烟草的行为。
美国自豪于控烟成功的同时，全球的70%烟民已经转移到发展中国家，而全球近一半因吸烟致死的死亡人数集中在中国。在中国，2015年的统计数据显示，15岁以上人群吸烟率为28.1%，其中男性吸烟率高达52.9%。中国的20多个城市虽然发布了全面禁烟的规定，但所覆盖的人口范围不足10%，因为对经济利益的追求，烟草企业及相关主管部门很难积极主动控烟。因此，美国的控烟史，对于我国极具借鉴意义，研究其控烟历程是较为有价值的。
1.2 研究意义
通过对1963-1992年美国各州香烟消费数据进行分析，寻找影响香烟销售量的重要因素，建立线性回归模型，分析依据美国控烟成功的经验，对我国的控烟问题提出有效建议。
1.3 研究内容和方法
第一步，从本文的数据背景入手，阐明本文研究的意义，确定本文的研究内容、方法、思路。
第二步，对1963-1992年美国各州香烟消费数据进行预处理和描述性统计。
第三步，探索性分析，寻找影响香烟销售量的重要因素，建立线性回归模型并进行比较，并对香烟销售量数据进行拟合。
第四步，对研究结果进行总结与反思，依据美国控烟成功的经验，针对我国的控烟问题提出相应建议。
2. 数据预处理与描述
2.1 数据预处理
数据为Ecdat包中的1963-1992年美国各州香烟消费的数据(Cigar)，具体变量包括：州（state）、年份（year）、每包香烟价格（price）、各州人口（pop）、大于16岁的人口（pop16）、消费者价格指数（cpi）、人均可支配收入（ndi）、人均香烟销售量（sales）、每包香烟在相邻州的最低价格（pimin），下面将原始数据的前几行进行展示：
Cigar=Ecdat::Cigar
head(Cigar) state year price pop pop16 cpi ndi sales pimin
1 1 63 28.6 3383 2236.5 30.6 1558.305 93.9 26.1
2 1 64 29.8 3431 2276.7 31.0 1684.073 95.4 27.5
3 1 65 29.8 3486 2327.5 31.5 1809.842 98.5 28.9
4 1 66 31.5 3524 2369.7 32.4 1915.160 96.4 29.5
5 1 67 31.6 3533 2393.7 33.4 2023.546 95.5 29.6
6 1 68 35.6 3522 2405.2 34.8 2202.486 88.4 32.02.1.1 缺失值
由于数据较多，数据时间跨度较长，可能存在隐含缺失，即某些州存在某一年份数据的缺失，因此，通过complete()来查看数据中是否有隐含缺失，从结果来看数据中并没有出现此种缺失，说明该数据是比较整齐的数据。
check=complete(Cigar,state,year)从前面展示的原始数据可以看到，数据中的年份只包含后两位数，如“1963”年缩写为“63”，不便于后面进行图表展示等工作，因此将其进行补全。
Cigar=Cigar%>%mutate(century=19)%>%
unite(year,century,year,sep = '')2.1.2 转换变量类型
由于原始数据中各州用数字代表，R导入数据时将其识别为数值型变量，因此需要将这一变量转换为因子变量。同时，上一步生成的年份由于通过unite()进行合并，新生成的年份变量为字符型变量，因此，将其变量类型转换为数值型变量。
Cigar$state=as.factor(Cigar$state)
Cigar$year=as.numeric(Cigar$year)2.1.3 生成新变量
后面经过分析，发现人均香烟销售量可能与人口比例和本州与相邻州香烟价格的差异有较大关系，为了方便后面的数据分析，在这里先根据原始数据生成以下两个新变量：成年人口比例proportion，本州与相邻州的最低价格的差价pricegap。下面将处理后的数据进行展示：
Cigar=Cigar%>%mutate(proportion=pop16/pop,pricegap=price-pimin)
head(Cigar) state year price pop pop16 cpi ndi sales pimin proportion
1 1 1963 28.6 3383 2236.5 30.6 1558.305 93.9 26.1 0.6610996
2 1 1964 29.8 3431 2276.7 31.0 1684.073 95.4 27.5 0.6635675
3 1 1965 29.8 3486 2327.5 31.5 1809.842 98.5 28.9 0.6676707
4 1 1966 31.5 3524 2369.7 32.4 1915.160 96.4 29.5 0.6724461
5 1 1967 31.6 3533 2393.7 33.4 2023.546 95.5 29.6 0.6775262
6 1 1968 35.6 3522 2405.2 34.8 2202.486 88.4 32.0 0.6829074
pricegap
1 2.5
2 2.3
3 0.9
4 2.0
5 2.0
6 3.62.1.4 生成新表
由于年份和州为本次分析的两个重要影响因素，为方便后续分析，通过总汇函数summarise()得到时间层面的1963-1992美国全国的人均香烟销售量、每包香烟价格等信息和地区层面的46个州30年平均的人口数、本州与相邻州的最低价格的差价等信息，生成以下两个表：
year=Cigar%>%group_by(year)%>%
summarise(sales=mean(sales),price=mean(price),ndi=mean(ndi),cpi=mean(cpi))
head(year)# A tibble: 6 x 5
year sales price ndi cpi
<dbl> <dbl> <dbl> <dbl> <dbl>
1 1963 127. 27.1 2100. 30.6
2 1964 124. 27.9 2230. 31
3 1965 126. 28.1 2395. 31.5
4 1966 126. 29.6 2554. 32.4
5 1967 126. 30.3 2697. 33.4
6 1968 124. 32.4 2889. 34.8state=Cigar%>%group_by(state)%>%
summarise(sales=mean(sales),price=mean(price),pimin=mean(pimin),pop=mean(pop),proportion=mean(proportion))%>%
mutate(pricegap=price-pimin)
head(state)# A tibble: 6 x 7
state sales price pimin pop proportion pricegap
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 107. 69.1 65.0 3768. 0.725 4.18
2 3 113. 67.5 65.1 2499. 0.719 2.33
3 4 115. 69.4 63.1 2172. 0.733 6.36
4 5 115. 71.5 67.0 23149. 0.745 4.48
5 7 118. 78.5 70.4 3081. 0.754 8.10
6 8 151. 69.3 63.6 590. 0.732 5.692.2 描述性统计
(1)年份
从时间上来看，该数据集时间跨度为1963-1992共30年，覆盖了美国香烟消费从兴盛到成功控烟的主要时期，每年包含46个州的香烟消费数据，包括了美国的绝大部分地区，因此，该数据还是十分具有分析价值的。
n_distinct(Cigar$year)[1] 30Cigar%>%group_by(year)%>%count()# A tibble: 30 x 2
# Groups: year [30]
year n
<dbl> <int>
1 1963 46
2 1964 46
3 1965 46
4 1966 46
5 1967 46
6 1968 46
7 1969 46
8 1970 46
9 1971 46
10 1972 46
# ... with 20 more rows由于时间因素对美国每包香烟价格、人均可支配收入、消费者价格指数、差价有较大影响，这些变量都有明显的随时间变化的趋势，因此将年份与这几个变量结合进行可视化。从下面绘制的几张图可以看到，随着时间的变化，美国每包香烟价格、人均可支配收入、消费者价格指数都呈现快速上升的趋势。
(2)价格
可以看到，随着时间的增加，每包香烟的价格大体呈现类似二次曲线的上升趋势。随着各州政府对控烟问题的重视，各州开始在香烟销售上施加不同比例的税收，导致从1970年开始，各州的每包香烟价格的差异渐渐扩大，1990年之后，随着控烟力度的加强，各州的每包香烟价格的差异进一步拉大。
ggplot(Cigar,aes(x=year,y=price,group=year))+geom_boxplot()+
labs(x='年份(year)',y='每包香烟价格(price)',title='1963-1992年美国香烟价格')
(3)人均可支配收入
随着时间增加，美国经济迅速发展，人均可支配收入迅速增加，但与此同时，由于各州经济发展的不均衡，各州的发展水平不同，从1970年开始，各州的收入差距渐渐拉大，1985年之后，各州的收入差距进一步扩大，这可能也是各州的每包香烟价格存在较大差异除不同税收比例之外的另一个原因。
ggplot(Cigar,aes(x=year,y=ndi,group=year))+geom_boxplot()+
labs(x='年份(year)',y='人均可支配收入(ndi)',title='1963-1992年人均可支配收入')
(4)消费者价格指数
随着时间增加和经济的发展，消费者价格指数也迅速攀升，但由于消费者价格指数为全国性的指标，各州的数值都一样，因此无法对各州数据进行分析。
ggplot(Cigar,aes(x=year,y=cpi))+geom_point()+
labs(x='年份(year)',y='消费者价格指数(cpi)',title='1963-1992年消费者价格指数')
(5)差价
随时间增加，如上文提到的，由于各州征收香烟税的幅度不同，各州与相邻州最低价格的差异逐渐增大，由此也带来了香烟走私的问题。从下图可以看到1965年之后，差价以0为中心向两边扩大，整体呈喇叭形。
ggplot(Cigar,aes(x=year,y=pricegap,group=year))+geom_boxplot()+
labs(x='年份(year)',y='差价(pricegap)',title='1963-1992年美国香烟差价')
(6)销量
30年中各州香烟销量近似正态分布，主要集中在100-150的区间内，均值为124，方差为31。
ggplot(Cigar)+geom_histogram(aes(sales))`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
sd(Cigar$sales)[1] 30.99105summary(Cigar$sales) Min. 1st Qu. Median Mean 3rd Qu. Max.
53.4 107.9 121.2 124.0 133.2 297.9 (7)人口
由于各州人口相对稳定，因此画出各州人口平均值的分布，可以看到大多数州人口在5000以下，有少数几个超过10000的人口大省。
ggplot(state,aes(x=pop))+geom_histogram(binwidth =1000)+
labs(x='人口(pop)',title='各州人口分布')
利用生成的成年人口比例绘制箱线图，可以看到，成年人口比例分布区间在0.6-0.9，可见各州成年人口比例还是存在较大差异。
ggplot(Cigar,aes(y=proportion))+geom_boxplot(group=1)+
labs(x='成年人口比例(proportion)',title='各州成年人口比例')
(8)州
数据中共包含46个州，各州经济发展水平、风俗习惯各不相同，例如，犹他州拥有高比例的摩门教徒人口（一种禁止吸烟的宗教），1988年人均香烟销量为55包，略低于全国平均113包的一半；内华达州是一个高度旅游的州，1988年人均香烟量为142包，比全国平均水平多出29包。同时，各州征收香烟税的比例也不相同，例如，由于印第安保留地出售免税香烟。蒙大拿州，新墨西哥州和亚利桑那州等印第安保留地州是非印第安人从保留地购买免税香烟的税收收入最大的输家；由于很多人从佛罗里达州，德克萨斯州，华盛顿州和佐治亚州的免税军事基地购买香烟，因此他们也是税收收入最大的输家。
各种因素综合导致美国各州人均香烟销售量存在较大差异，因此，建立模型时需考虑地区的影响，但可以假设这些州对香烟销量的特殊影响是相对固定的。
ggplot(Cigar,aes(state,sales))+geom_boxplot()+
labs(x='各州(state)',y='人均香烟销售量(sales)',title='1963-1992年美国各州人均香烟销售量')+scale_x_discrete(labels=NULL)
3. 探索性数据分析
3.1 变量分析
(1)时间
总体来看，1963-1992年美国人均香烟销售量随时间先增加后降低，符合美国自80年代开始加强控烟力度的总体背景，其中，1963-1976年人均香烟销售量整体呈上升趋势，而1976-1992年销售量呈下降趋势，因此考虑在模型中加入平方项进行拟合。
ggplot(year,aes(x=year,y=sales,group=1))+geom_line()+
labs(x='年份(year)',y='人均香烟销售量(sales)',title='1963-1992年美国人均香烟销售量')
(2)各州
比较各州人均香烟销售量与价格，发现整体上各州人均香烟销售量随价格升高而降低，但有几个异常值点需要进行进一步探究。
ggplot(state,aes(x=price,y=sales))+geom_point()+
labs(x='价格(price)',y='人均香烟销售量(sales)',title='美国各州人均香烟销售量与价格')
将各州人均香烟销售量与价格同时进行比较。通过观察发现人均香烟销售量较高的几个州价格相对较低，但优势并不十分明显，因此考虑引入每包香烟与相邻州的最低价格的差价这一变量。
ggplot(state,aes(x=state))+geom_point(aes(y=sales))+geom_point(aes(y=price),color='red')+scale_x_discrete(labels = NULL)+
labs(x='各州(state)',y='人均香烟销售量(sales)与价格(price)',title='美国各州人均香烟销售量与价格')
下面绘制各州人均香烟销售量与相邻州差价的关系图，可以发现当本州与相邻州的最低价格的差价小于0时，即本州价格更为优惠时，周边州的消费者可能更倾向于从该州购买香烟。价格差异引发香烟走私的问题。这可以对部分异常值点进行解释。
ggplot(state,aes(x=pricegap,y=sales))+geom_point()+
labs(x='相邻州差价(pricegap)',y='人均香烟销售量(sales)',title='美国各州人均香烟销售量与相邻州差价')
3.2 建立模型
3.2.1 数据处理
为方便后面的计算，将时间进行标准化。同时，根据州进行分组，通过嵌套生成列表变量，便于后面进行多模型同步的操作。
Cigar_=Cigar%>%mutate(year_s=scale(year))
staten=Cigar_%>%group_by(state)%>%nest()3.2.2 建立模型
(1)模型1
把时间和时间的平方项作为变量建立线性模型，通过glance()计算可决系数，各州可决系数在0.36-0.93区间内，将其进行可视化。
state_model1=function(df){
lm(sales~I(year_s^2)+year_s,data=df)
}
statem1=staten%>%mutate(model=map(data,state_model1))%>%
mutate(resids=map2(data,model,add_residuals))
glance1=statem1%>%mutate(glance=map(model,broom::glance))%>%unnest(glance)
glance1%>%ggplot(aes(reorder(state,r.squared),r.squared))+geom_point()+scale_x_discrete(labels = NULL)+
labs(x='州(state)',y='可决系数(r.squared)',title='美国各州模型可决系数')
(2)模型2
在前一个模型的基础上加入差价这一变量，可以发现可决系数有了明显提升，各州可决系数较之前缩小到了0.53-0.94区间内，将其进行可视化。
state_model2=function(df){
lm(sales~I(year_s^2)+year_s+pricegap,data=df)
}
statem2=staten%>%mutate(model=map(data,state_model2))%>%
mutate(resids=map2(data,model,add_residuals))
glance2=statem2%>%mutate(glance=map(model,broom::glance))%>%unnest(glance)
glance2%>%ggplot(aes(reorder(state,r.squared),r.squared))+geom_point()+scale_x_discrete(labels = NULL)+
labs(x='州(state)',y='可决系数(r.squared)',title='美国各州模型可决系数')
3.2.3 拟合效果
将可决系数小于0.6的两个州筛选出来，可以看到这两个州在1970年有非常明显的下降，导致了误差的产生，对于这两个州来说，可能在现有模型基础上还需进行其他调整。经分析，从1968年人均香烟销售量突然出现较为明显的下降趋势，可能由于60年代美国联邦卫生局发布了第一份《吸烟与健康》报告，让部分烟民认识到了“吸烟危害健康”，然而1970年后又开始迅速回升，直到80年代政府加大控烟力度，人均香烟销售量才开始迅速下降。
predict=glance2%>%mutate(pred=map2(data,model,add_predictions))%>%unnest(pred)
predict%>%filter(state==c('26','39'))%>%ggplot(aes(x=year,y=sales,group=state,color=state))+geom_point()+geom_line(aes(y=pred))+
labs(x='年份(year)',y='人均香烟销售量(sales)',color='各州(state)',title='1963-1992年美国各州人均香烟销售量')
下面将各州销售量数据散点图与拟合曲线进行展示，大部分州的模型拟合程度还是比较好的，与该州人均香烟销售量的变化趋势相符。
predict%>%ggplot(aes(x=year,y=sales,group=state,color=state))+geom_point()+geom_line(aes(y=pred))+guides(color=FALSE)+
labs(x='年份(year)',y='人均香烟销售量(sales)',color='各州(state)',title='1963-1992年美国各州人均香烟销售量')
4. 结论与反思
1.从数据本身来看，该数据集受时间因素和地区因素的影响明显，反映了不同年代，政府总体上对控烟问题的重视程度和采取征收香烟税这一措施的有效性，与美国控烟的具体时代背景相符。同时，也反映了不同地区文化的差异性，以及各州政府控烟力度不同而导致的香烟走私问题。此外，还反映出随着时间的变化，美国经济迅速发展，物价和人均可支配收入迅速上升，但与此同时，也存在地区发展不平衡的问题。
2.通过探索性分析，发现了时间和地区这两个重要影响因素，在考虑了时间因素和地区因素的影响后，进行建模，通过可决系数可以看出，对于大部分州来说，建立的模型还是比较合适的，但仍存在改进的空间。在后面的研究中，可加入其它影响变量对模型进行进一步的调整，使模型能对香烟销售量进行更好的拟合。
3.同时，由于部分影响因素难易量化或通过已有变量体现，例如各州政府的控烟力度和对香烟税的征收力度，以及人们对吸烟危害的认识程度，都对香烟销售量有较大影响，而这些因素较难通过已有变量进行体现，因此，也需要进一步的研究，可能需要增加各州具体征收香烟税的数值和香烟税与其他税收的比较来体现控烟力度的大小，或根据各州政府出台的法规政策来增加虚拟变量，从而对政府控烟力度进行量化。
4.美国控烟的成功对于我国极具借鉴意义，其征收香烟税等控烟方式事实证明是非常有效的。香烟作为一种证明对健康有害又可以合法销售的商品，世界各国都对烟草制品课以重税，通过此种方式，抬高香烟售价抑制烟民购买。欧盟国家这一比例普遍达到75%以上，我国的卷烟税负水平只有56%，而世卫组织建议的税负水平为75%。中国烟草行业每年给财政贡献的税收占财政总收入的6%，这个贡献背后是残酷而高昂的成本。数据显示，2005年中国人群中归因于烟草使用的死亡已达120万人，其中有33.8%的人在40-69岁之间死去；2030年烟草归因死亡估计占40岁以上人口死亡的25%，超过300万人。中国的控烟问题刻不容缓。