肺活量(LungCapData)是衡量一个人身体健康程度的重要指标。一般而言，影响肺活量的因素有很多，例如性别、年龄、身材、是否经常运动、是否吸烟与酗酒等。本文欲探究年龄(Age)、身高(Height)、是否吸烟(Smoke)、性别(Gender)与是否是剖腹产(Caesarean)对肺活量是否存在影响，若有影响，影响程度如何。根据分析结果，本文提出了相应的建议。

二、数据来源和相关说明

本文利用LungCapData数据，对肺活量的影响因素进行探究。数据共计725条，涉及的变量6个，分别为：

表 2-1 变量说明

变量	符号	解释
肺活量	Y	反应体质的指标
对数肺活量	ln⁡(Y)	肺活量取对数
年龄	X1	连续变量
身高	X2	连续变量
是否吸烟	X3	否=0；是=1
性别	X4	女=0；男=1
是否剖腹产	X5	否=0；是=1

三、描述性分析

为了获得对数据的整体了解，本文先对数据进行了描述性统计分析。

3.1 样本描述

表 3-1 样本描述

变量名	最小值	中位数	最大值	平均值
LungCap(Y)	0.507	8.000	14.675	7.863
Age(X1)	3.000	13.000	19.000	12.330
Height(X2)	45.300	65.400	81.800	64.840
Smoke(X3)	0.000	0.000	1.000	0.106
Gender(X4)	0.000	1.000	1.000	0.506
Caesarean(X5)	0.000	0.000	1.000	0.226

从表3-1可以得出：肺活量（LungCapData）介于0.507-14.675之前，其平均水平约为7.863（平均值）和8.000（中位数）；年龄介于3-19岁之间，其平均水平约为12.330（平均值）和13.000（中位数）；身高介于45.300-81.800之间，其平均水平为64.840（平均值）和65.400（中位数）；是否吸烟取值为0&1，其平均水平为0.106（平均值）和0.000（中位数），由此可见不吸烟的人数更多；性别取值为0&1，其平均水平为0.506（平均值）和1（中位数），由此可见男性人数多一些；是否是剖腹产取值为0&1，其平均水平为0.226（平均值）和0（中位数），由此可见非剖腹产人数较多。

3.2 数据可视化

为了更直观地看数据的分布情况，本文用R软件，分别绘制了肺活量、年龄和体重的直方图，是否吸烟、性别和是否是剖腹产的饼图。为了分析各个变量之间的关系，本文又绘制了各个变量的相关性热力图，是否吸烟、性别和是否是剖腹产与肺活量的箱线图，具体结果见图3-1到图3-4。

图 3-1 直方图

从图3-1直方图可以得出：肺活量、年龄和身高均大致服从正态分布。其中，肺活量数值多集中在6~9之间，年龄多集中在10~16岁，身高多集中在60~70之间。

图 3-2 饼图

从图3-2饼图可以得出：样本中，绝大多数人群不吸烟，占比89.4%，少数人吸烟，占比10.9%；男性与女性人数基本相同，男性略多余女性，分别占比50.6%和49.4%；人群中大部分人不是剖腹产，占比77.4%，少数人是剖腹产，占比22.6%。

图 3-3 相关性热力图

从图3-3热力图可以得出：肺活量与年龄和身高都有较强的相关性，相关系数分别为0.82和0.91。年龄和身高之间也存在较强的相关性，相关系数为0.84。

图 3-4 箱线图

从图3-4箱线图可以得出：从直观上看，是否吸烟与性别在对数肺活量上存在一定的差异性，而是否是剖腹产在对数量上基本不存在差异性。此外，不吸烟人群肺活量的最大值要高于吸烟人群；男性肺活量要高于女性，这与常识相符；非剖腹产人群的肺活量也高于剖腹产人群。

四、数据建模

4.1 回归模型A

为了分析肺活量的影响因素，本文首先建立了回归模型A （不考虑交互效应），即以ln⁡(Y) 作为因变量，X1~X5 作为自变量，建立回归模型。

lnY=β0+β1X1+β2X2+β3X3+β4X4+β5X5+ε

模型结果如下表4-1所示。

表 4-1 回归模型A结果

变量	Estimate	Std.Error	t value	Pr(>\|t\|)
β0	-0.9746	0.097	-10.099	<2e-16***
X1	0.0251	0.004	6.797	2.25e-11***
X2	0.0407	0.002	19.743	<2e-16
X3	-0.0678	0.026	-2.626	0.009***
X4	0.0526	0.016	3.224	0.001**
X5	-0.0345	0.019	-1.854	0.064
R-squared	0.7697
Adjusted R-squared p-value	0.7681 <2.2e-16

从表4-1可以得出：X1~X5 中，只有X5 的P 值>0.05 ，即在显著性水平α=0.05 的条件下，是否是剖腹产对肺活量不存在显著性影响。因而，本文考虑去除变量X5 后再建立回归模型B 。

4.2 回归模型B

以ln⁡(Y) 作为因变量，X1~X4 作为自变量，建立回归模型。

lnY=β0+β1X1+β2X2+β3X3+β4X4+ε

模型结果如下表4-2所示。同时，对模型A和模型B进行方差分析，以进一步判断是否舍弃变量X5 ，分析结果如表4-2所示。

表 4-2 回归模型B结果

变量	Estimate	Std.Error	t value	Pr(>\|t\|)
	-0.9763	0.097	-10.099	<2e-16***
X1	0.0250	0.004	6.768	2.71e-11***
X2	0.0406	0.002	19.676	<2e-16
X3	-0.0691	0.026	-2.674	0.008***
X4	0.0524	0.016	3.202	0.001**
R-squared	0.7686
Adjusted R-squared p-value	0.7673 <2.2e-16
P-value	2.2e-6

由表4-2可以得出：变量X1~X4 的P 值均<0.05 ，即在显著性水平α=0.05 的条件下，这四个变量对肺活量均存在显著性影响。

表 4-3 回归模型A和回归模型B方差分析

模型	Res.Df	RSS	Df	Sum of Sq	F	Pr(>f)
A	719	31.404
B	720	31.554	-1	-0.150	3.436	0.064

由表4-3可以得出：模型A 和模型B 方差分析P 值为0.064>0.05 ，因而在显著性水平α=0.05 的条件下接受原假设，即认为模型A 和模型B 不存在显著性差异，即变量X5 可以剔除。因而，本文采用模型B 的分析结果对肺活量的影响因素进行定量分析。模型B 对应的回归模型如下。

lnY=-0.9763+0.0250X1+0.0406X2-0.0691X3+0.0524X4

4.3 结果分析

根据表4-2可以得出：β1=0.0250 ，意味着其他变量保持不变时，X1 (年龄)每增加1岁，肺活量将相较于之前增加2.50%；β2=0.0406 ，意味着其他变量保持不变时，X2 (身高)高度每增加1，肺活量将相较于之前增加4.06%；β3=-0.0691 ，意味着其他变量保持不变时，X3 (是否吸烟)每增加1，即吸烟的人相较于不吸烟的人，肺活量将降低6.91%；β4=0.0524 ，意味着其他变量保持不变时，X4 (性别)每增加1，即男性相较于女性，肺活量将增加5.24%。

综上所述，年龄、身高和性别均对肺活量有正向影响，其中性别的影响最大，身高次之，年龄影响最小；是否吸烟对肺活量存在负向影响，即从平均意义上而言，吸烟会降低人的肺活量，且影响较大。

五、结论及建议

5.1 结论

样本中，男性略多于女性，吸烟者和剖腹产仅占一小部分，分别为10.6%和22.6%。肺活量、年龄和身高都大致服从正态分布。从图形上看，年龄和身高与肺活量存在强相关性，是否吸烟与性别在肺活量上有较大差异性，是否是剖腹产差异性较小。

根据回归分析可以得出：在显著性水平α=0.05 下，五个因素中，四个因素均对肺活量存在显著性影响，分别是年龄、身高、是否吸烟和性别。随着年龄和身高的增长，肺活量均会增加；男性肺活量的均值要高于女性；吸烟者的肺活量均值要低于不吸烟者。

5.2 建议

综上所述，本文提出以下两点建议：

吸烟者应尽量戒烟：根据本文分析，吸烟会明显降低肺活量。此外，吸烟对身体的危害不言而喻，因而为了自身身体健康考虑，应考虑尽早戒烟。
肺活量是测量心肺功能的一项重要指标，肺活量的大小反映了人体的心肺功能和健康程度。据网上资料显示，可以通过做有氧运动，如登山、慢跑、游泳等；或者做增强扩胸的运动等都可以增加自身肺活量。为了自身的身体健康，可以考虑经常做运动来增加自己的肺活量，改善自身体质。

六、代码

a=read.csv("D:/个人成长/学业/课程/大三下课程/统计模型/线性回归模型和方差分析模型/LungCapData.csv",header=1)
a[1:5,]
b=a[,1:3]
a$Y=log(a$LungCap)
dim(a)
par(mfrow=c(1,3))
hist(a$LungCap,xlab="肺活量",ylab = "Frequency",main=NULL)
#hist(log(a$LungCap),xlab="肺活量",ylab = "Frequency",main=NULL)
hist(a$Age,xlab="年龄",ylab = "Frequency",main=NULL)
hist(a$Height,xlab="身高",ylab = "Frequency",main=NULL)
#hist(log(a$Age),xlab="年龄",ylab = "Frequency",main=NULL)
#hist(log(a$Height),xlab="身高",ylab = "Frequency",main=NULL)
x=c(648,77)
piepercent1=round(100*x/sum(x),1)
y=c(358,367)
piepercent2=round(100*y/sum(x),1)
z=c(561,164)
piepercent3=round(100*z/sum(x),1)
par(mfrow=c(1,3))
pie(x,labels=piepercent1,main="是否吸烟",col=rainbow(length(x)))
legend("topright",c("不吸烟","吸烟"),cex=0.8,fill=rainbow(length(x)))
pie(y,labels=piepercent2,main="性别",col=rainbow(length(y)))
legend("topright",c("女性","男性"),cex=0.8,fill=rainbow(length(y)))
pie(z,labels=piepercent3,main="是否是剖腹产",col=rainbow(length(z)))
legend("topright",c("否","是"),cex=0.8,fill=rainbow(length(z)))
library(corrplot)
##默认参数
k=cor(b,use='everything',method='pearson')
par(mfrow=c(1,1))
corrplot(k,addCoef.col = "black")
summary(a)
descrb = function(var)
  Z=a[,var]
  N=tapply(a$Y,Z,length)
  MU=tapply(a$Y,Z,mean)
  SD=tapply(a$Y,Z,sd)
  MIN=tapply(a$Y,Z,min)
  MED=tapply(a$Y,Z,median)
  MAX=tapply(a$Y,Z,max)
  out=cbind(N,MU,SD,MIN,MED,MAX)
  out

descrb("Smoke") 
descrb("Gender")
descrb("Caesarean")

par(mfrow=c(1,3))
boxplot(Y~Smoke,a,ylab="对数肺活量",xlab="是否吸烟",names=c("否","是"))
boxplot(Y~Gender,a,ylab="对数肺活量",xlab="性别",names=c("女","男"))
boxplot(Y~Caesarean,a,ylab="对数肺活量",xlab="是否是剖腹产",names=c("否","是"))

Model.A=lm(Y~Age+Height+Smoke+Gender+Caesarean,data=a)
summary(Model.A)

Model.B=lm(Y~Age+Height+Smoke+Gender,data=a)
summary(Model.B)

anova(Model.A,Model.B)

个人见解，请各位读者批评指正！

统计模型产品磨损老化单因素方差分析

产品磨损老化单因素方差分析

摘要本文通过描述性统计分析与单因素方差分析，分析国内外工厂某产品磨损老化实验数据，得出如下结论：四个工厂之间、国内外之间、国内各工厂之间均存在显著性差异；工厂A3与A4、A2与A4、A1与A3、A1与A2之间也存在显著性差异。即国内与国外，国内之间产品质量参差不齐。四个工厂之间，变化率均值由低到高分别是A4、A1、A3和A2，即可以认为国外产品质量最好，国内工厂乙质量最差。基于此，本文提出了加强质量监控与检测标准以提高国内工厂产品质量整体水平的建议。

一、研究目的

有四种产品，分别来自国内甲、乙、丙三个工厂与国外同类产品，现从各厂分别取10，6，6和2个产品做300小时连续磨损老化检验，得到变化率。通过分析国内外产品质量，对国内外产品优化提出建议。

二、数据来源和相关说明

数据来源于300小时连续磨损老化试验数据，具体数据如表2-1所示。

表 2-1 300小时连续磨损老化实验数据

产品	老化率
A1	20 18 19 17 15 16 13 18 22 17
A2	26 19 26 28 23 25
A3	24 25 18 22 27 24
A4	12 14

【注】A1、A2、A3、A4分别代表国内甲、乙、丙三个国内工厂与国外同类产品

三、描述性分析

为了对数据有更直观地了解，本文首先进行了描述性统计分析，绘制了各个工厂产品磨损老化变化率的直方图与散点图，结果如图3-1所示。

由图3-1可以得出：四个工厂产品变化率均值由低到高分别是A4、A1、A3和A2，即可以认为国外产品质量最好，国内工厂甲次之，其次是工厂丙，工厂乙质量最差。有箱线图可以看出，各个工厂生产的产品变化率差别较大，可能存在显著性差异。

图 3-1 直方图与散点图

四、数学建模

为了分析国内外产品、国内各产品与各个工厂之间的产品是否存在显著性差异，本文用R语言进行了单因素方差分析。

进行单因素方差分析前，首先对数据的正态性与方差齐性进行检验，本文采用Shapiro检验和bartlett检验，具体检验结果如4-1所示。两个检验的P值均>0.05，即在显著性水平α=0.05 下都通过了显著性检验。

表 4-1 正态性与方差齐性检验

	正态性检验	方差齐性检验
p-value	0.4163	0.8326

其次，对四个工厂的产品进行单因素方差分析。

Yik=μi+εik,i=1,2,3,4;k=1,2,⋯,ni

即检验μ1=μ2=μ3=μ4 是否成立，检验结果如表4-2所示。由P值<0.05可以得出，在显著性水平α=0.05 下拒绝原假设，即认为四个工厂的产品存在显著性差异。

表 4-2 四个工厂间的单因素方差分析

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
factory	3	346	115.33	14.66	2.79e-05***
Residuals	20	157.3	7.87

第三，对国内外产品的差异性进行检验。合并甲、乙、丙三个工厂的产品磨损数据，进而与国外工厂A4的数据进行单因素方差分析，分析结果如表4-3所示。由P值=0.0169<0.05可以得出，在显著性水平α=0.05 下拒绝原假设，即认为国内外的产品也存在显著性差异。

表 4-3 国内外产品单因素方差分析

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
factory	1	117.3	117.33	6.687	0.0169*
Residuals	22	386	17.55

第四，对国内甲、乙、丙工场的产品进行单因素方差分析，分析结果如表4-4所示。由P值=0.000184<0.05可以得出，在显著性水平α=0.05 下拒绝原假设，即认为国内产品之间也存在显著性差异。

表 4-4 国内产品单因素方差分析

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
factory	2	228.7	114.33	13.98	0.000184***
Residuals	19	155.3	8.18

第五，对四个工厂的产品进行两两比较，进行进一步分析。分析结果如图4-1所示。六组中，只有A4-A1与A3-A2两组的95%置信区间包含0，因而认为只有这两组之间不存在显著性差异，其他四组之间均存在显著性差异。

图 4-1 两两比较

五、结论与建议

5.1 结论

综上所述，四个工厂之间、国内外之间、国内各工厂之间均存在显著性差异；工厂A3与A4、A2与A4、A1与A3、A1与A2之间存在显著性差异。即国内与国外，国内之间产品质量参差不齐。四个工厂之间，变化率均值由低到高分别是A4、A1、A3和A2，即可以认为国外产品质量最好，国内工厂乙质量最差。

5.2建议

国外产品质量要高于国内质量，因而国内工厂要考虑加强质量监控，提高产品质量。同时政府要加强对各企业的监控，提高检测标准，提高国内企业产品质量整体水平。

六、代码

c=read.csv("D:/个人成长/学业/课程/大三下课程/统计模型/作业/第二次作业/two.csv",header=1)
c

#检验正态性和方差齐性
library(car)
shapiro.test(c$rate)
bartlett.test(rate~factory, data = c)

par(mfrow=c(1,2))
boxplot(rate~factory,c,ylab="变化率",xlab="工厂")
stripchart(rate ~ factory, c, vertical=TRUE, method="stack",xlab="工厂",ylab="变化率")
fit3=aov(rate~factory,data=c)
summary(fit3)

par(mfrow=c(1,1))
d=read.csv("D:/个人成长/学业/课程/大三下课程/统计模型/作业/第二次作业/three.csv",header=1)
d
boxplot(rate~factory,d,ylab="变化率",xlab="工厂",names=c("国内","国外"))
fit4=aov(rate~factory,data=d)
summary(fit4)

e=read.csv("D:/个人成长/学业/课程/大三下课程/统计模型/作业/第二次作业/four.csv",header=1)
e
fit5=aov(rate~factory,data=e)
summary(fit5)

##方差分析后，多重比较，继续探寻两两分组间的差异
#Tukey HSD 检验
tuk=TukeyHSD(fit3, conf.level = 0.95)
plot(tuk)
tuk

以上是关于统计模型肺活量影响因素分析报告的主要内容，如果未能解决你的问题，请参考以下文章

R语言Logistic回归模型案例:低出生婴儿体重的影响因素分析（列线图校准曲线）

在回归分析中，怎样用SPSS求解完全二次模型