数海拾荒数海拾荒记录数据分析学习历程，并适当涉及编程以及AI话题关注专栏更多R（编程语言）{&debug&:false,&apiRoot&:&&,&paySDK&:&https:\u002F\u002Fpay.zhihu.com\u002Fapi\u002Fjs&,&wechatConfigAPI&:&\u002Fapi\u002Fwechat\u002Fjssdkconfig&,&name&:&production&,&instance&:&column&,&tokens&:{&X-XSRF-TOKEN&:null,&X-UDID&:null,&Authorization&:&oauth c3cef7c66aa9e6a1e3160e20&}}{&database&:{&Post&:{&&:{&title&:&配煤动力学特性回归模型的R语言实现&,&author&:&realphy&,&content&:&\u003Cp\u003E之前在研究生毕业论文当中用SPSS采用多元线性及非线性模型用于预测配煤动力学特性，本节主要结合R语言提供的函数和工具重现一下。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E1 几个概念：\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E单煤：就是从一个矿区开采出来的煤炭，一般性质相近；\u003C\u002Fp\u003E\u003Cp\u003E配煤：为了取得特定的特性，由不同种单煤按照不同的比例混合而成，由于煤炭性质的复杂性，通常配煤的动力学特性无法由各种单煤加权平均求得；\u003C\u002Fp\u003E\u003Cp\u003E煤的工业分析（\u003Ca href=\&http:\u002F\u002Fbaike.baidu.com\u002Fitem\u002F%E7%85%A4%E7%9A%84%E5%B7%A5%E4%B8%9A%E5%88%86%E6%9E%90\&\u003E煤的工业分析_百度百科\u003C\u002Fa\u003E）：煤的工业分析是指包括煤的水分（M ）、灰分（A ）、挥发分（V ）和固定碳（Fc )四个分析项目指标的测定的总称。煤的工业分析是了解煤质特性的主要指标，也是评价煤质的基本依据。通常煤的水分、灰分、挥发分是直接测出的，而固定碳是用差减法计算出来的。广义上讲，煤的工业分析还包括煤的全硫分和发热量的测定， 又叫煤的全工业分析。本文的实验数据介于狭义与广义之间，测量了发热量，但是全硫在元素分析中进行。\u003C\u002Fp\u003E\u003Cp\u003E煤的元素分析（\u003Ca href=\&http:\u002F\u002Fbaike.baidu.com\u002Fitem\u002F%E7%85%A4%E7%9A%84%E5%85%83%E7%B4%A0%E5%88%86%E6%9E%90\&\u003E煤的元素分析_百度百科\u003C\u002Fa\u003E）：煤的元素分析是对煤中的元素含量进行检测和分析（一般用质量百分数表示），包括常规的C、H、O、N、S、Al、Si、Fe、Ca等元素含量，还可检测煤中的痕量元素包括Ti、Na、K等。本文的实验数据主要测量了最主要的C、H、O、N、S五种。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E煤的着火温度：煤的着火点是随着热力条件变化而变化的，并不是一个物理常数，只是在一定条件下得出的相对特征值。通常把由缓慢氧化状态转变为高速燃烧状态的瞬间过程称为着火，此时的温度即称为着火温度（Ignition Temperature）。\u003C\u002Fp\u003E\u003Cp\u003E煤的活化能（Activation Energy）：即让某一种煤开始发生\u003Ca href=\&https:\u002F\u002Fzh.wikipedia.org\u002Fwiki\u002F%E5%8C%96%E5%AD%A6%E5%8F%8D%E5%BA%94\&\u003E化学反应\u003C\u002Fa\u003E（通常指燃烧）所需要克服的能量障碍。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E2 数据描述：\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E有16种单煤及48种这些单煤按不同比例混合而成配煤的工业分析以及元素分析数据，并通过热重分析、热重微分分析及差热分析得出的曲线结合Arrhenius方程及Doyle单曲线法求得各个不同煤样的着火温度（）和活化能。\u003C\u002Fp\u003E\u003Cp\u003E这样一共64组数据。\u003C\u002Fp\u003E\u003Cp\u003E考虑到从实验的角度求取着火温度及活化能需要特殊的仪器并耗费较长的时间，希望能基于上述64组数据找到一个用于预测工业分析、元素分析已知煤样的动力学特性，包括着火温度、活化能等，特别是在混配煤动力学特性预测方面进行应用。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E3 数据导入：\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Ccode lang=\&splus\&\u003E& library(xlsx)\n载入需要的程辑包：rJava\n载入需要的程辑包：xlsxjars\n\n载入程辑包：‘xlsx’\n\nThe following objects are masked from ‘package:openxlsx’:\n\n
createWorkbook, loadWorkbook, read.xlsx, saveWorkbook,\n
write.xlsx\n\n& filepath &- \&G:\u002FDataCruiser\u002Fworkspace\u002FCoal Analysis\u002FCoal Analysis data.xlsx\&\n& CoalData &- read.xlsx(filepath, 1, encoding = \&UTF-8\&)\n& class(CoalData)\n[1] \&data.frame\&\n& head(CoalData)\n
煤样编号 煤样名称 Ignition_Temp Activation_Ene Mad.
Vad. FCad.\n1
.85 23.68 26.98 45.49\n2
.00 14.62 10.44 73.94\n3
.50 14.42 26.06 57.02\n4
.47 14.19 30.02 46.32\n5
.21 18.26 29.39 50.14\n6
.62 25.47 23.32 49.59\n
Qnet.ad.J.g.
Cad. Had. Nad. St.add.
.93 3.42 1.12
.90 3.14 1.15
.55 3.92 0.81
.44 2.92 0.83
0.44 13.71\n5
.13 3.74 1.18
0.71 10.77\n6
.26 3.32 1.09
6.72\n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cb\u003E4 数据重命名：\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E着火温度：ignitionTemperature\u003C\u002Fp\u003E\u003Cp\u003E活化能：activationEnergy\u003C\u002Fp\u003E\u003Cp\u003E水分：moisture\u003C\u002Fp\u003E\u003Cp\u003E灰分：ash\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E挥发分：volatileMatter\u003C\u002Fp\u003E\u003Cp\u003E固定碳：fixedCarbon\u003C\u002Fp\u003E\u003Cp\u003E发热量：LHV\u003C\u002Fp\u003E\u003Cp\u003E碳：carbon\u003C\u002Fp\u003E\u003Cp\u003E氢：hydrogen\u003C\u002Fp\u003E\u003Cp\u003E氮：nitrogen\u003C\u002Fp\u003E\u003Cp\u003E硫：sulfur\u003C\u002Fp\u003E\u003Ccode lang=\&splus\&\u003E& names(CoalData) &- c(\&sampleNO\&,\&sampleName\&,\&ignitionTemperature\&,\&activationEnergy\&,\&moistrue\&,\n+
\&ash\&,\&volatileMatter\&,\&fixedCarbon\&,\&LHV\&,\&carbon\&,\&hydrogen\&,\&nitrogen\&,\&sulfur\&,\&oxygen\&)\n& head(CoalData)\n
sampleNO sampleName ignitionTemperature activationEnergy moistrue
3.85 23.68\n2
1.00 14.62\n3
2.50 14.42\n4
9.47 14.19\n5
2.21 18.26\n6
1.62 25.47\n
volatileMatter fixedCarbon
LHV carbon hydrogen nitrogen sulfur\n1
6.72\n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cb\u003E4 数据分组：\u003C\u002Fb\u003E \u003C\u002Fp\u003E\u003Cp\u003E按不同因变量（ignitionTemperature和activationEnergy）进行数据分组。\u003C\u002Fp\u003E\u003Ccode lang=\&splus\&\u003E& IT_Estimate &- as.data.frame(CoalData[,c(\&ignitionTemperature\&,\&moistrue\&,\n+
\&ash\&,\&volatileMatter\&,\&fixedCarbon\&,\&LHV\&,\&carbon\&,\&hydrogen\&,\&nitrogen\&,\&sulfur\&,\&oxygen\&)])\n& head(IT_Estimate)\n
ignitionTemperature moistrue
ash volatileMatter fixedCarbon
3.85 23.68
1.00 14.62
2.50 14.42
9.47 14.19
2.21 18.26
1.62 25.47
49.59 25702.6\n
carbon hydrogen nitrogen sulfur oxygen\n1
6.72\n& AE_Estimate &- as.data.frame(CoalData[,c(\&activationEnergy\&,\&moistrue\&,\n+
\&ash\&,\&volatileMatter\&,\&fixedCarbon\&,\&LHV\&,\&carbon\&,\&hydrogen\&,\&nitrogen\&,\&sulfur\&,\&oxygen\&)])\n& head(AE_Estimate)\n
activationEnergy moistrue
ash volatileMatter fixedCarbon
3.85 23.68
1.00 14.62
2.50 14.42
9.47 14.19
2.21 18.26
1.62 25.47
49.59 25702.6\n
carbon hydrogen nitrogen sulfur oxygen\n1
6.72\n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cb\u003E5 数据相关性检验：\u003C\u002Fb\u003E \u003C\u002Fp\u003E\u003Cp\u003E这里我们采用Pearson相关性分析，需要调用psych包中的corr.test()函数，结果如下：\u003C\u002Fp\u003E\u003Ccode lang=\&splus\&\u003E& corr.test(IT_Estimate, use = \&complete\&)\nCall:corr.test(x = IT_Estimate, use = \&complete\&)\nCorrelation matrix \n
ignitionTemperature moistrue
ash volatileMatter\nignitionTemperature
-0.65\nmoistrue
1.00 -0.42
-0.20\nvolatileMatter
0.31 -0.20
1.00\nfixedCarbon
-0.12 -0.59
-0.65\nLHV
-0.14 -0.79
-0.14\ncarbon
-0.06 -0.82
-0.22\nhydrogen
-0.33 -0.25
0.60\nnitrogen
-0.49 -0.23
0.00\nsulfur
0.30\noxygen
0.61 -0.29
fixedCarbon
LHV carbon hydrogen nitrogen sulfur\nignitionTemperature
-0.07\nmoistrue
-0.12 -0.14
-0.11\nash
-0.59 -0.79
0.29\nvolatileMatter
-0.65 -0.14
0.30\nfixedCarbon
-0.44\nLHV
-0.31\ncarbon
-0.39\nhydrogen
0.26\nnitrogen
-0.06\nsulfur
-0.44 -0.31
1.00\noxygen
-0.40 -0.16
oxygen\nignitionTemperature
-0.57\nmoistrue
-0.29\nvolatileMatter
0.67\nfixedCarbon
-0.40\nLHV
-0.16\ncarbon
-0.27\nhydrogen
0.04\nnitrogen
-0.40\nsulfur
-0.04\noxygen
1.00\nSample Size \n[1] 64\nProbability values (Entries above the diagonal are adjusted for multiple tests.) \n
ignitionTemperature moistrue
ash volatileMatter\nignitionTemperature
0.00\nmoistrue
1.00\nvolatileMatter
0.00\nfixedCarbon
0.27\ncarbon
0.08\nhydrogen
0.00\nnitrogen
0.99\nsulfur
0.02\noxygen
fixedCarbon
LHV carbon hydrogen nitrogen sulfur\nignitionTemperature
1.00\nmoistrue
0.46\nvolatileMatter
0.40\nfixedCarbon
0.33\ncarbon
0.05\nhydrogen
0.80\nnitrogen
1.00\nsulfur
0.00\noxygen
oxygen\nignitionTemperature
0.00\nmoistrue
0.51\nvolatileMatter
0.00\nfixedCarbon
1.00\ncarbon
0.68\nhydrogen
1.00\nnitrogen
0.03\nsulfur
1.00\noxygen
0.00\n\n To see confidence intervals of the correlations, print with the short=FALSE option\n&
corr.test(AE_Estimate, use = \&complete\&)\nCall:corr.test(x = AE_Estimate, use = \&complete\&)\nCorrelation matrix \n
activationEnergy moistrue
ash volatileMatter\nactivationEnergy
-0.40 -0.03
-0.56\nmoistrue
1.00 -0.42
-0.20\nvolatileMatter
0.31 -0.20
1.00\nfixedCarbon
-0.12 -0.59
-0.65\nLHV
-0.14 -0.79
-0.14\ncarbon
-0.06 -0.82
-0.22\nhydrogen
-0.33 -0.25
0.60\nnitrogen
-0.49 -0.23
0.00\nsulfur
0.30\noxygen
0.61 -0.29
fixedCarbon
LHV carbon hydrogen nitrogen sulfur oxygen\nactivationEnergy
-0.46\nmoistrue
-0.12 -0.14
-0.59 -0.79
-0.29\nvolatileMatter
-0.65 -0.14
0.67\nfixedCarbon
-0.40\nLHV
-0.16\ncarbon
-0.27\nhydrogen
0.04\nnitrogen
-0.40\nsulfur
-0.44 -0.31
-0.04\noxygen
-0.40 -0.16
1.00\nSample Size \n[1] 64\nProbability values (Entries above the diagonal are adjusted for multiple tests.) \n
activationEnergy moistrue
ash volatileMatter fixedCarbon\nactivationEnergy
0.00\nmoistrue
0.00\nvolatileMatter
0.00\nfixedCarbon
0.00\ncarbon
0.00\nhydrogen
0.22\nnitrogen
0.01\nsulfur
0.00\noxygen
LHV carbon hydrogen nitrogen sulfur oxygen\nactivationEnergy 0.04
0.01\nmoistrue
0.51\nvolatileMatter
0.00\nfixedCarbon
1.00\ncarbon
0.68\nhydrogen
1.00\nnitrogen
0.03\nsulfur
1.00\noxygen
0.00\n\n To see confidence intervals of the correlations, print with the short=FALSE option\u003C\u002Fcode\u003E\u003Cp\u003E从相关性检验中可以得出无论是着火温度还是活化能，与水分、挥发分、固定碳、发热量以及氧这五个变量相关较大，因此，以下考虑以此五个自变量采用多元线性回归模型分别对着火温度和活化能进行模拟，首先剔除不相关数据，仅包含我们感兴趣的变量。\u003C\u002Fp\u003E\u003Ccode lang=\&splus\&\u003E& AE_Estimate &- as.data.frame(AE_Estimate[,c(\&activationEnergy\&,\&moistrue\&,
\&volatileMatter\&,\&fixedCarbon\&,\&LHV\&,\&oxygen\&)])\n& IT_Estimate &- as.data.frame(IT_Estimate[,c(\&ignitionTemperature\&,\&moistrue\&,
\&volatileMatter\&,\&fixedCarbon\&,\&LHV\&,\&oxygen\&)])\n& \u003C\u002Fcode\u003E\u003Cp\u003E用car包中的scatterplotMatrix()函数生成散点图矩阵。\u003C\u002Fp\u003E\u003Ccode lang=\&splus\&\u003E& library(car)\n\n载入程辑包：‘car’\n\nThe following object is masked from ‘package:psych’:\n\n
logit\n\n& scatterplotMatrix(IT_Estimate, spread = FALSE, smoother.args = list(lty = 2), main = \&Scatter Plot Matrix ignitionTemperature\&)& scatterplotMatrix(AE_Estimate, spread = FALSE, smoother.args = list(lty = 2), main = \&Scatter Plot Matrix of activationEnergy\&)\n\n\u003C\u002Fcode\u003E\u003Cimg src=\&v2-439fdbbc09f0eacb81a62.png\& data-caption=\&\& data-rawwidth=\&916\& data-rawheight=\&900\&\u003E\u003Cimg src=\&v2-dd050ad128.png\& data-caption=\&\& data-rawwidth=\&916\& data-rawheight=\&900\&\u003E\u003Cp\u003E从上图我们可以看出，着火温度是一个单峰曲线，随着水分、挥发分、氧含量的升高而降低，随着固定碳以及发热量升高而升高。活化能也有类似的规律。\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E5 回归分析\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E下面使用lm()函数分别拟合各自的多元回归模型，代码如下：\u003C\u002Fp\u003E\u003Ccode lang=\&splus\&\u003E& fit_AE &- lm(activationEnergy ~ moistrue + volatileMatter + fixedCarbon + LHV + oxygen, data = AE_Estimate)\n& summary(fit_AE)\n\nCall:\nlm(formula = activationEnergy ~ moistrue + volatileMatter + fixedCarbon + \n
LHV + oxygen, data = AE_Estimate)\n\nResiduals:\n
\n\nCoefficients:\n
Estimate Std. Error t value Pr(&|t|)
\n(Intercept)
3.535 0.000808 ***\nmoistrue
-1.090 0.280255
\nvolatileMatter
-2.493 0.015536 *
\nfixedCarbon
-1.099 0.276134
1.875 0.065797 .
0.408 0.684715
\n---\nSignif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n\nResidual standard error: 33440 on 58 degrees of freedom\nMultiple R-squared:
0.465,\tAdjusted R-squared:
0.4188 \nF-statistic: 10.08 on 5 and 58 DF,
p-value: 5.577e-07\n\n& fit_IT &- lm(ignitionTemperature ~ moistrue + volatileMatter + fixedCarbon + LHV + oxygen, data = IT_Estimate)\n& summary(fit_IT)\n\nCall:\nlm(formula = ignitionTemperature ~ moistrue + volatileMatter + \n
fixedCarbon + LHV + oxygen, data = IT_Estimate)\n\nResiduals:\n
Max \n-57.360
57.860 \n\nCoefficients:\n
Estimate Std. Error t value Pr(&|t|)
\n(Intercept)
& 2e-16 ***\nmoistrue
0.00580 ** \nvolatileMatter
5.9e-05 ***\nfixedCarbon
0.00271 ** \noxygen
\n---\nSignif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n\nResidual standard error: 18.61 on 58 degrees of freedom\nMultiple R-squared:
0.6926,\tAdjusted R-squared:
0.6661 \nF-statistic: 26.13 on 5 and 58 DF,
p-value: 1.03e-13\u003C\u002Fcode\u003E\u003Cp\u003E采用标准方法对模型进行诊断：\u003C\u002Fp\u003E\u003Ccode lang=\&splus\&\u003E& confint(fit_AE)\n
97.5 %\n(Intercept)
9. 66\nmoistrue
\nvolatileMatter -1.
-\nfixedCarbon
27.62129\noxygen
\n& confint(fit_IT)\n
97.5 %\n(Intercept)
429.0.\nmoistrue
-1.\nvolatileMatter
-3.\nfixedCarbon
0.\noxygen
5.\n& par(mfrow=c(2,2))\n& plot(fit_AE)\n& plot(fit_AE,main = \&活化能诊断图\&)\n& plot(fit_IT,main = \&着火温度诊断图\&)\n\u003C\u002Fcode\u003E\u003Cimg src=\&v2-959cb3206a36cdb09d4fa1.png\& data-caption=\&\& data-rawwidth=\&916\& data-rawheight=\&900\&\u003E\u003Cp\u003E\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cimg src=\&v2-9e28a354e2166b9bcad1fdd0d7f92fed.png\& data-caption=\&\& data-rawwidth=\&916\& data-rawheight=\&900\&\u003E\u003Cp\u003E我们可以看出总体上数据符合正态性、线性的假设，在活化能的模型当中可能需要适当的加入一些二次项，另外，两个模型都有部分离群点，可能是实验的误差导致。\u003C\u002Fp\u003E&,&updated&:new Date(&T13:52:56.000Z&),&canComment&:false,&commentPermission&:&anyone&,&commentCount&:8,&likeCount&:8,&state&:&published&,&isLiked&:false,&slug&:&&,&isTitleImageFullScreen&:false,&rating&:&none&,&sourceUrl&:&&,&publishedTime&:&T21:52:56+08:00&,&links&:{&comments&:&\u002Fapi\u002Fposts\u002F2Fcomments&},&url&:&\u002Fp\u002F&,&titleImage&:&https:\u002F\u002Fpic1.zhimg.com\u002Fv2-5f310e9e1d3129edafbc76a8f7841c6f_r.jpg&,&summary&:&&,&href&:&\u002Fapi\u002Fposts\u002F&,&meta&:{&previous&:null,&next&:null},&snapshotUrl&:&&,&commentsCount&:8,&likesCount&:8},&&:{&title&:&《R语言实战》第三部分第十章-功效分析学习笔记&,&author&:&realphy&,&content&:&功效分析可以帮助在给定置信度的情况下，判断检测到给定效应值时所需的样本量。在给定置信度水平情况下，计算在某样本量内能检测到给定效应值的概率。在开始学习R语言中功效分析工具包之前首先回顾下统计学中假设检验的概念。\u003Cp\u003E\u003Cb\u003E1 假设检验回顾\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E在统计假设检验中，首先要对总体分布参数设定一个假设（零假设，H0），然后从总体分布中抽样，通过样本计算所得的统计量来对总体参数进行推断。假定零假设为真，如果计算获得观测样本的统计量的概率非常小，便可以拒绝原假设，接受它的对立面（称作备择假设或者研究假设，H1）。\u003C\u002Fp\u003E\u003Cp\u003E在使用取自总体的样本数据来对总体做推断时，通常会有以下四种情况：\u003C\u002Fp\u003E\u003Cp\u003E\u003Cimg src=\&v2-eaf4f0aaf4.png\& data-rawwidth=\&793\& data-rawheight=\&210\&\u003E在研究过程当中，我们一般关注以下四个量：样本大小、显著性水平、功效和效应值。\u003C\u002Fp\u003E\u003Cul\u003E\u003Cli\u003E样本大小指的是实验设计中每种条件\u002F组中观测的数目。\u003C\u002Fli\u003E\u003Cli\u003E显著性水平（也称为alpha）由Ⅰ型错误的概率来定义。也可以把它看作发现效应不发生的概率\u003C\u002Fli\u003E\u003Cli\u003E功效通过1减去Ⅱ型错误的概率来定义。我们可以把它看作真实效应发生的概率。\u003C\u002Fli\u003E\u003Cli\u003E效应值指的是在备择或研究假设下效应的量。效应值的表达式依赖于假设检验中使用的统计方法。\u003C\u002Fli\u003E\u003C\u002Ful\u003E\u003Cbr\u003E其中，样本大小和显著性水平可以直接控制，而功效和效应值的影响时间接的。四个量紧密相关，给定其中任意三个量，便可以推算第四个量。这也是下面学习如何使用R语言中pwr包实现功效分析的基础。\u003Cimg src=\&v2-34eb026decff9d5cb494593.png\& data-rawwidth=\&395\& data-rawheight=\&419\&\u003E\u003Cp\u003E\u003Cbr\u003E\u003Cb\u003E2 pwr包做功效分析\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003Et检验\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E格式：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Epwr.t.test(n=, d=, sig.level=, power=, type=, alternative=)\n\u003C\u002Fcode\u003E\u003Cp\u003E其中：\u003C\u002Fp\u003E\u003Cbr\u003E1, n为样本大小\u003Cbr\u003E\u003Cbr\u003E2, d为效应值，即标准化的均值之差
\u003Cequation\u003Ed=\\frac{\\mu _{1}- \\mu _{2}}{\\sigma } \u003C\u002Fequation\u003E\u003Cbr\u003E\u003Cbr\u003E3, sig.level表示显著性水平（默认为0.05）\u003Cp\u003E4, power为功效水平\u003Cbr\u003E\u003Cbr\u003E5, type指检验类型：双样本t检验（“two.sample”）,单样本t检验（\&one.sample\&）或相依样本t检验（\&paired\&）。默认为双样本t检验\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cbr\u003E6, \&alternative\&指统计检验是双侧检验（\&two.sided\&）还是单侧检验（\&less\&或\&greater\&）。默认为双侧检验\u003Cbr\u003E\u003Cbr\u003E实例：\u003Cp\u003E\u003Cu\u003E求样本大小：\u003C\u002Fu\u003E\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cbr\u003E\u003Ccode lang=\&text\&\u003E& library(pwr)\n& pwr.t.test(d=.8,sig.level = .05,power = .9,type = \&two.sample\&,alternative = \&two.sided\&)\n\n
Two-sample t test power calculation \n\n
n = 33.82554\n
sig.level = 0.05\n
power = 0.9\n
alternative = two.sided\n\nNOTE: n is number in *each* group\n\n& \n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cu\u003E求功效水平：\u003C\u002Fu\u003E\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& pwr.t.test(n=20,d=.5,sig.level = .01)\n\n
Two-sample t test power calculation \n\n
sig.level = 0.01\n
power = 0.1439551\n
alternative = two.sided\n\nNOTE: n is number in *each* group\n\n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cb\u003E方差分析\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E格式：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Epwr.anova.test(k=, n=, f=, sig.level=, power=)\n\u003C\u002Fcode\u003E\u003Cp\u003E其中：\u003C\u002Fp\u003E1, k为组的个数\u003Cp\u003E2, n为各组中样本大小\u003Cbr\u003E\u003Cbr\u003E2, f为效应值，对于单因素方差分析，\u003Cequation\u003Ef=\\sqrt{\\frac{\\sum_{i-1}^{k}{p_{i} x} \\times \\left(\\mu _{i} -\\mu
\\right) ^{2} }{\\sigma ^{2} } } \u003C\u002Fequation\u003E\u003Cbr\u003E\u003Cbr\u003E3, sig.level表示显著性水平（默认为0.05）\u003C\u002Fp\u003E\u003Cp\u003E4, power为功效水平\u003C\u002Fp\u003E\u003Cp\u003E实例：\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cbr\u003E\u003Cp\u003E\u003Cu\u003E求样本大小：\u003C\u002Fu\u003E\u003Cbr\u003E\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& library(pwr)\n& pwr.anova.test(k = 5, f = .25, sig.level = .05,power = .8)\n\n
Balanced one-way analysis of variance power calculation \n\n
n = 39.1534\n
f = 0.25\n
sig.level = 0.05\n
power = 0.8\n\nNOTE: n is number in each group\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cb\u003E相关性\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E格式：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Epwr.r.test(n=, r=, sig.level=, power=, alternative=)\u003C\u002Fcode\u003E\u003Cp\u003E其中：\u003C\u002Fp\u003E\u003Cp\u003E1, n为观测数目\u003C\u002Fp\u003E\u003Cp\u003E2, r为效应值（通过线性相关系数衡量）\u003C\u002Fp\u003E\u003Cp\u003E3, sig.level是显著性水平\u003C\u002Fp\u003E\u003Cp\u003E4, power为功效水平\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E5, \&alternative\&指统计检验是双侧检验（\&two.sided\&）还是单侧检验（\&less\&或\&greater\&）。默认为双侧检验\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E实例：\u003C\u002Fp\u003E\u003Cp\u003E\u003Cu\u003E求样本大小：\u003C\u002Fu\u003E\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& library(pwr)\n& pwr.r.test(r = .25, sig.level = 0.05, power = .90,alternative = \&greater\&)\n\n
approximate correlation power calculation (arctangh transformation) \n\n
n = 133.2803\n
r = 0.25\n
sig.level = 0.05\n
power = 0.9\n
alternative = greater\n\n& \n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cb\u003E线性模型\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E格式：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Epwr.f2.test(u=, v=, f2=, sig.level=, power=)\u003C\u002Fcode\u003E\u003Cp\u003E其中：\u003C\u002Fp\u003E\u003Cp\u003E1, u为分子自由度，v为分母自由度\u003C\u002Fp\u003E\u003Cp\u003E2,\u003Cequation\u003Ef^{2} \u003C\u002Fequation\u003E为效应值\u003C\u002Fp\u003E\u003Cp\u003E当评价目的的不同，\u003Cequation\u003Ef^{2} \u003C\u002Fequation\u003E有两种不同的表达式，如果评价一组预测变量对结果的影响程度时，则：\u003C\u002Fp\u003E\u003Cp\u003E\u003Cequation\u003Ef^{2} =\\frac{R^{2} }{1-R^{2} } \u003C\u002Fequation\u003E，其中\u003Cequation\u003ER^{2} \u003C\u002Fequation\u003E为多重相关性的总体平方值\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E如果要评价的是预测变量对结果的影响超过第二组变量—即协变量—多少时，\u003C\u002Fp\u003E\u003Cp\u003E则：\u003C\u002Fp\u003E\u003Cp\u003E\u003Cequation\u003Ef^{2} =\\frac{R_{AB}^{2}
-R_{A}^{2} }{1-R_{AB}^{2}
} \u003C\u002Fequation\u003E，其中\u003Cequation\u003ER_{A}^{2}\u003C\u002Fequation\u003E为集合A中变量对总体方差的解释率，\u003Cequation\u003ER_{AB}^{2}\u003C\u002Fequation\u003E为集合A和B中变量对总体方差的解释率\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cbr\u003E\u003Cp\u003E3, sig.level是显著性水平\u003C\u002Fp\u003E\u003Cp\u003E4, power为功效水平\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E5, \&alternative\&指统计检验是双侧检验（\&two.sided\&）还是单侧检验（\&less\&或\&greater\&）。默认为双侧检验\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E实例：\u003C\u002Fp\u003E\u003Cp\u003E研究老板的领导风格对员工满意度的影响，是否超过薪水和工作小费对员工满意度的影响。\u003C\u002Fp\u003E\u003Cp\u003E已知领导风格能够解释35%的方差，薪水和消费能够解释约30%的员工满意度的方差。\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E假定显著性水平为0.05，在90%的置信度下。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cu\u003E求样本大小：\u003C\u002Fu\u003E\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& pwr.f2.test(u=3,f2 = 0.0769,sig.level = 0.05,power = 0.9)\n\n
Multiple regression power calculation \n\n
v = 184.2426\n
f2 = 0.0769\n
sig.level = 0.05\n
power = 0.9\n\n&\u003C\u002Fcode\u003E\u003Cp\u003E在多元回归中，分母的自由度等于N-k-1（\u003Ca href=\&https:\u002F\u002Fwww.zhihu.com\u002Fquestion\u002F\& data-editable=\&true\& data-title=\&多元线性回归函数残差平方和的自由度是怎样确定的？ - 多元线性回归 - 知乎\&\u003E多元线性回归函数残差平方和的自由度是怎样确定的？ - 多元线性回归 - 知乎\u003C\u002Fa\u003E），其中N是总观测数，k是预测变量数。\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E本例中，N-7-1 = 185，则样本大小为193。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E比例检验\u003C\u002Fb\u003E\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E格式（各组中n相同）：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Epwr.2p.test(h=, n=, sig.level=, power=)\u003C\u002Fcode\u003E\u003Cp\u003E其中：\u003C\u002Fp\u003E\u003Cp\u003E1，h为效应值\u003C\u002Fp\u003E\u003Cp\u003Eh的定义如下：\u003C\u002Fp\u003E\u003Cp\u003E\u003Cequation\u003E2arcsin(\\sqrt{p_{1} } )-2arcsin(\\sqrt{p_{2} } )\u003C\u002Fequation\u003E，在R语言当中可以用ES.h(p1,p2)函数来计算。\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E2，n为各组相同的样本量\u003C\u002Fp\u003E\u003Cp\u003E格式（各组中n不相同）：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Epwr.2p2n.test(h=, n1=, n2=, sig.level=, power=)\u003C\u002Fcode\u003E\u003Cp\u003E实例：\u003C\u002Fp\u003E\u003Cp\u003E药物差异的药物受试者样本量问题。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cu\u003E求样本大小：\u003C\u002Fu\u003E\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& pwr.2p.test(h=ES.h(.65,.6), sig.level = .05, power = .9,alternative = \&greater\&)\n\n
Difference of proportion power calculation for binomial distribution (arcsine transformation) \n\n
h = 0.1033347\n
sig.level = 0.05\n
power = 0.9\n
alternative = greater\n\nNOTE: same sample sizes\n\n& \u003C\u002Fcode\u003E\u003Cp\u003E\u003Cb\u003E卡方检验\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E卡方检验常用于检验两个\u003Cb\u003E类型变量\u003C\u002Fb\u003E之间的关系，典型的零假设是变量之间独立，备择假设是不独立。\u003C\u002Fp\u003E\u003Cp\u003E格式：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Epwr.chisq.test(w=, N=, df=, sig.level=, power=)\u003C\u002Fcode\u003E其中：\u003Cp\u003E1，N是总样本大小\u003C\u002Fp\u003E\u003Cp\u003E2，df是自由度\u003C\u002Fp\u003E\u003Cp\u003E3，w是效应值，定义如下：\u003C\u002Fp\u003E\u003Cp\u003E\u003Cequation\u003Ew=\\sqrt{\\sum_{i=1}^{m}{\\frac{\\left( p0_{i} -p1_{i} x\\right) ^{b} }{p0_{i} } } } \u003C\u002Fequation\u003E，其中\u003Cequation\u003Ep0_{i} =H_{0} \u003C\u002Fequation\u003E时第i单元格中的概率，\u003Cequation\u003Ep1_{i} =H_{1} \u003C\u002Fequation\u003E时第i单元格中的概率，并且从1到m进行求和，m指的是列联表中单元格的数目。R语言当中ES.w2(P)函数可以计算双因素列联表中的备则假设的效应值。\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E实例：\u003C\u002Fp\u003E\u003Cp\u003E研究对象的双因素列联表如下：\u003C\u002Fp\u003E\u003Cp\u003E\u003Cimg src=\&v2-aa496df4a0f28a2af86694.png\& data-rawwidth=\&777\& data-rawheight=\&157\&\u003E则，效应值计算如下：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& prob &- matrix(c(.42, .28, .03, .07, .10, .10), byrow = TRUE, nrow=3)\n& ES.w2 (prob)\n[1] 0.1853198\n\u003C\u002Fcode\u003E\u003Cp\u003E另外，双因素列联表的自由度为(r-1)(c-1)。\u003Cbr\u003E\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& pwr.chisq.test(w=.1853, df = 2, sig.level = .05, power = .9)\n\n
Chi squared power calculation \n\n
w = 0.1853\n
N = 368.5317\n
sig.level = 0.05\n
power = 0.9\n\nNOTE: N is the number of observations\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cb\u003E效应值的选择\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E功效分析当中，难在预期效应值的选择。在行为科学领域，Cohen曾经提出一个基准。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cimg src=\&v2-70bffc07d8ef0f8ecaaa50d76f502560.png\& data-rawwidth=\&776\& data-rawheight=\&229\&\u003E当然，Cohen的基准值只是根据许多社科类研究提出的一般性建议，特定问题的研究还是需要具体问题具体分析。\u003C\u002Fp\u003E\u003Cp\u003E另外，我们也可以改变研究参数，并记录其对诸如样本大小和功效等方面的影响。现以五分组的单因素方差分析为例，检测一系列的效应值所需要的样本大小。\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Elibrary(pwr)\nes &- seq(.1, .5, .01)\nnes &- length(es)\nsamsize &- NULL\nfor(i in 1:nes){\n
result &- pwr.anova.test(k = 5, f = es[i], sig.level = .05, power = .9)\n
samsize[i] &- ceiling(result$n)\n}\n\nplot(samsize, es, type = \&l\&, lwd = 2, col = \&red\&,\n
ylab = \&Effect Size\&,\n
xlab = \&Sample Size (per cell)\&,\n
main = \&One Way ANOVA with Power = .90 and Alpha = .05\&)\n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cimg src=\&v2-faa34aa63c.png\& data-rawwidth=\&899\& data-rawheight=\&900\&\u003E2 绘制功效分析图形\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E其实跟上一节当中绘制的图形类似，用for循环来计算一系列效应值和功效水平下所需的样本量，并以相关性的功效分析为例。\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E#生成相关系数和功效值系列\nlibrary(pwr)\nr &- seq(.1, .5, .01)\nnr &- length(r)\n\np &- seq(.4, .9, .1)\nnp &- length(p)\n\n#获取样本大小\nsamsize &- array(numeric(nr*np), dim=c(nr,np))\nfor (i in 1:np){\n
for (j in 1:nr){\n
result &- pwr.r.test(n = NULL, r = r[j],\n
sig.level = .05, power = p[i],\n
alternative = \&two.sided\&)\n
samsize[j,i] &- ceiling(result$n)\n
}\n}\n\n#创建图形\nxrange &- range(r)\nyrange &- round(range(samsize))\ncolors &- rainbow(length(p))\nplot(xrange, yrange, type=\&n\&,\n
xlab=\&Correlation Coefficient (r)\&,\n
ylab=\&Sample Size (n)\& )\n\n#添加功效曲线\nfor (i in 1:np){\n
lines(r, samsize[,i], type=\&l\&, lwd=2, col=colors[i])\n}\n\n#添加网格线\nabline(v=0, h=seq(0,yrange[2],50), lty=2, col=\&grey89\&)\nabline(h=0, v=seq(xrange[1],xrange[2],.02), lty=2, col=\&gray89\&)\n\n#添加注释\ntitle(\&Sample Size Estimation for Correlation Studies\\n\nSig=0.05 (Two-tailed)\&)\nlegend(\&topright\&, title=\&Power\&, as.character(p),\n
fill=colors)\n\n\u003C\u002Fcode\u003E\u003Cimg src=\&v2-cb56ebaca456e266c6987.png\& data-rawwidth=\&899\& data-rawheight=\&900\&\u003E\u003Cp\u003E\u003Cb\u003E4 其他软件包\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E除了pwr包，在研究的规划阶段，R还提供了不少其他有用的软件包可供选择，限于篇幅，原文并没有做过多的叙述和举例说明。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cimg src=\&v2-ebcd91d13ea712af074ad8.png\& data-rawwidth=\&715\& data-rawheight=\&551\&\u003E其中，后面五个包主要聚焦于基因研究中的功效分析，识别基因与可观测特征的关联性的研究称为全基因组关联研究（GWAS）。此外，MBESS包也包含了可供各种形式功效分析所用的函数，有兴趣的朋友可以自行help了解。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E5 小结\u003C\u002Fb\u003E\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E在前面三章，主要探索了各种各样的统计假设检验函数，本章更加关注研究的准备阶段，功效分析不仅可以帮助你判断在给定置信度和效应值的前提下所需的样本量，也能说明在给定样本量时检测到要求效应值的概率。对于限定误报效应显著性的可能性（Ⅰ型错误）和正确检测真实效应（功效）的可能性的平衡，也有了初步的了解。\u003C\u002Fp\u003E\u003Cp\u003E本章主要内容是接受pwr包中各个函数的具体用法，这些函数可以对常见的统计方法进行功效分析，特别是计算样本量。最后也给出了一些专业化的功效分析方法。\u003C\u002Fp\u003E\u003Cp\u003E我们必须记住的是，典型的功效分析是一个交互性的过程。研究者会通过改变样本量、效应值、预期显著性水平和预期功效水平等参数，来观测它们对于其他参数的影响。\u003C\u002Fp\u003E\u003Cp\u003EPS：最近连环出差，而且是日常排得满满当当的差，这篇笔记其实上上周就开始提笔了，直到现在才完成。另外，虽然书已经看到第13章了，但是复习笔记才更新到第十章，也是一个回顾总结的过程，今天因为\u003Ca href=\&http:\u002F\u002Fwww.zhihu.com\u002Fpeople\u002F2dcbc480abaa\& data-hash=\&2dcbc480abaa\& class=\&member_mention\& data-hovercard=\&p$b$2dcbc480abaa\&\u003E@不在意\u003C\u002Fa\u003E 同学昨晚在群里的一个问题又回过头去复习了一遍各个检验之间的区别，甚好甚好。\u003C\u002Fp\u003E&,&updated&:new Date(&T13:18:41.000Z&),&canComment&:false,&commentPermission&:&anyone&,&commentCount&:2,&likeCount&:4,&state&:&published&,&isLiked&:false,&slug&:&&,&isTitleImageFullScreen&:false,&rating&:&none&,&sourceUrl&:&&,&publishedTime&:&T21:18:41+08:00&,&links&:{&comments&:&\u002Fapi\u002Fposts\u002F2Fcomments&},&url&:&\u002Fp\u002F&,&titleImage&:&https:\u002F\u002Fpic3.zhimg.com\u002Fv2-e3f3897dbdb0ae83df5b16_r.jpg&,&summary&:&&,&href&:&\u002Fapi\u002Fposts\u002F&,&meta&:{&previous&:null,&next&:null},&snapshotUrl&:&&,&commentsCount&:2,&likesCount&:4},&&:{&title&:&《R语言实战》第三部分第十一章-中级绘图学习笔记&,&author&:&realphy&,&content&:&本章主要在第六章的基础上，在后续几章零星介绍了部分绘图的方法，本章将是之前图形主题的延伸与扩展，并且主要关注用于展示双变量关系（二元关系）和多变量间关系（多元关系）的绘图方法。\u003Cp\u003E\u003Cb\u003E1 散点图\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E散点图比较简单，在之前的内容当中有过介绍，可用来描述两个连续型变量间的关系。\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Eattach(mtcars)\nplot(wt, mpg,\n
main = \&Basic Scatter plot of MPG vs. Weight\&,\n
xlab = \&Car Weight (lbs\u002F1000)\&,\n
ylab = \&Miles Per Gallon\&, pch = 19)\n\nabline(lm(mpg~wt), col = \&red\&, lwd = 2, lty = 1)\nlines(lowess(wt, mpg), col = \&blue\&, lwd = 2, lty = 2)\n\n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cimg src=\&v2-cd0c927fb4b7227c4faa33.png\& data-rawwidth=\&899\& data-rawheight=\&900\&\u003E需要注意的是，R中有两个平滑曲线拟合函数：lowess()函数和loess()函数。loess()函数是基于lowess()函数表达式版本的更新和更强大的拟合函数。由于这两个函数默认值不同，使用时需要加倍小心。\u003C\u002Fp\u003E\u003Cp\u003E此外，car包中的scatterplot()函数增强了散点图的许多功能，以下实例做了一个演示。\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Elibrary(car)\nscatterplot(mpg ~ wt | cyl, data = mtcars, lwd = 2, span = 0.75,\n
main = \&Scatter Plot of MPG vs. Weight by # Cylinders\&,\n
xlab=\&Weight of Car (lbs\u002F1000)\&,\n
ylab=\&Miles Per Gallon\&,\n
legend.plot=TRUE,\n
id.method=\&identify\&,\n
labels=row.names(mtcars),\n
boxplots=\&xy\&\n
)\n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cimg src=\&v2-04fba77fed64371b90acea8ae5705b08.png\& data-rawwidth=\&899\& data-rawheight=\&900\&\u003E\u003Cb\u003E散点图矩阵\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E除了一个散点图，可以同时生成多个散点图，组成散点图矩阵，比如采用R语言中pairs()函数。\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Epairs(~mpg+disp+drat+wt, data=mtcars,main=\&Basic Scatter Plot Matrix\&)\n\u003C\u002Fcode\u003E\u003Cimg src=\&v2-389adea52ac.png\& data-rawwidth=\&899\& data-rawheight=\&900\&\u003E\u003Cp\u003E同样的，car包中有对应的scatterplotMatrix()函数也可以功能更为丰富的散点图矩阵。比如：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Elibrary(car)\nscatterplotMatrix(~ mpg + disp + drat + wt, data=mtcars,\nspread=FALSE, smoother.args=list(lty=2),\nmain=\&Scatter Plot Matrix via car Package\&)\n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cimg src=\&v2-830d422ccb74fa.png\& data-rawwidth=\&899\& data-rawheight=\&900\&\u003E从结果上看，线性和平滑（loess）拟合曲线被默认添加，而且主对角线处还添加了核密度曲线和轴须图。此外，spread = FALSE选项表示不添加展示分散度和对称信息的直线，smoother.args=list(lty=2)则设定平滑（loess）拟合曲线使用虚线而不是实线。\u003C\u002Fp\u003E\u003Cp\u003ER语言中散点图矩阵多种多样，除了上述两种，glus包中的cpars()函数，TeachingDemos包中的pairs2()函数，HH包中的xysplom()函数，ResourceSelection包中的kepairs()函数和SMPracticals包中的pairs.mod()函数都可以绘制各自不同功能的散点图矩阵。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E高密度散点图\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E当数据点重叠非常严重的时候，相关用于绘制高密度散点图的函数便排上用场了。首先是smoothScatter()函数，利用核密度估计生成用颜色密度来表示点分布的散点图。\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& set.seed(1234)\n& n &- 10000\n& c1 &- matrix(rnorm(n,mean = 0,sd = .5), ncol = 2)\n& c2 &- matrix(rnorm(n,mean = 3,sd = 2), ncol = 2)\n& mydata &- rbind(c1, c2)\n& mydata &- as.data.frame(mydata)\n& names(mydata) &- c(\&x\&,\&y\&)\n& with(mydata,\n+
plot(x, y, pch=19, main=\&Scatter Plot with 10,000 Observations\&))\n\n\n\u003C\u002Fcode\u003E\u003Cimg src=\&v2-477ef63efe01.png\& data-rawwidth=\&899\& data-rawheight=\&900\&\u003E\u003Ccode lang=\&text\&\u003E& with(mydata, smoothScatter(x,y,main = \&Scatter Plot Colored by Smoothed Densities\&))\n& \n\u003C\u002Fcode\u003E\u003Cimg src=\&v2-184cf3edbb637faa830acc.png\& data-rawwidth=\&899\& data-rawheight=\&900\&\u003E\u003Cp\u003E为便于比较，以上也给出了未处理时的原图。类似的，hexbin包中的hexbin()函数可以将二元变量的封箱放到六边形单元格中，图形更加直观。示例如下：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Einstall.packages(\&hexbin\&)\nlibrary(hexbin)\nwith(mydata, {\nbin &- hexbin(x, y, xbins=50)\nplot(bin, main=\&Hexagonal Binning with 10,000 Observations\&)\n})\n\u003C\u002Fcode\u003E\u003Cimg src=\&v2-62b004ee3ccb2bfea6b513.png\& data-rawwidth=\&680\& data-rawheight=\&900\&\u003E\u003Cp\u003E\u003Cb\u003E三维散点图\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E三维散点图，顾名思义，就是可以用来对三个定量变量的交互关系进行可视化。\u003C\u002Fp\u003E\u003Cp\u003E格式：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Escatterplot3d(x, y, z)\u003C\u002Fcode\u003E\u003Cp\u003E举个例子，如果对于mtcars包中的汽车英里数、车重和排量间的关系感兴趣。可以用以下代码分别生成以下三幅三维散点图。\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& install.packages(\&scatterplot3d\&)\ntrying URL 'https:\u002F\u002Fmirrors.tuna.tsinghua.edu.cn\u002FCRAN\u002Fbin\u002Fwindows\u002Fcontrib\u002F3.3\u002Fscatterplot3d_0.3-38.zip'\nContent type 'application\u002Fzip' length 440226 bytes (429 KB)\ndownloaded 429 KB\n\npackage ‘scatterplot3d’ successfully unpacked and MD5 sums checked\n\nThe downloaded binary packages are in\n\tC:\\Users\\panhuayin\\AppData\\Local\\Temp\\Rtmp8e4ncd\\downloaded_packages\n& library(scatterplot3d)\n& attach(mtcars)\n& scatterplot3d(wt, disp, mpg, main = \&Basic 3D Scatter Plot\&)\n& \n\u003C\u002Fcode\u003E\u003Cimg src=\&v2-bb9fe7e6fc5c0e2fe5825.png\& data-rawwidth=\&899\& data-rawheight=\&900\&\u003E\u003Cp\u003E利用satterplot3d()函数的一些选项。\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& scatterplot3d(wt, disp, mpg, \n+
pch = 16,\n+
highlight.3d = TRUE,\n+
type = \&h\&,\n+
main = \&3D Scatter Plot with Vertical Lines\&)\n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cbr\u003E\u003Cimg src=\&v2-1fb784cee4c571e9bd7d1.png\& data-rawwidth=\&899\& data-rawheight=\&900\&\u003E此外，还可以借助之前的回归，给图形添加一个回归面，其中平面代表预测值，图中的点是实际值。\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Elibrary(scatterplot3d)\nattach(mtcars)\ns3d &- scatterplot3d(wt, disp, mpg,\n
pch = 16,\n
highlight.3d = TRUE,\n
type = \&h\&,\n
main = \&3D Scatter Plot with Vertical Lines and Regression Plane\&)\nfit &- lm(mpg ~ wt + disp)\ns3d$plane3d(fit)\n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cimg src=\&v2-6b7aa0b1d054c416b697f80e2243b9fd.png\& data-rawwidth=\&899\& data-rawheight=\&900\&\u003E为了更好的进行交互，R还提供了一些旋转图形功能，比如旋转三维散点图。这里我们采用rgl包中的plot3d()函数。\u003C\u002Fp\u003E\u003Cp\u003E格式：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Eplot3d(x, y, z)\u003C\u002Fcode\u003E\u003Cp\u003E其中x，y和z是数值型向量，代表各个点。继续上面的例子，代码和图形如下：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& library(rgl)\n& plot3d(wt, disp, mpg, col = \&red\&, size = 5)\n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cimg src=\&v2-cfdbc4df550104cbf0ae2d4273bca108.png\& data-rawwidth=\&548\& data-rawheight=\&530\&\u003Ecar包中的scatter3d()函数也提供类似的功能，代码和图形如下：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& library(car)\n& with(mtcars, scatter3d(wt, disp, mpg))\n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cimg src=\&v2-36acf16ceaf7fe49578b3c.png\& data-rawwidth=\&667\& data-rawheight=\&676\&\u003E不难发现，scatter3d()函数自带回归面。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E气泡图\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E所谓气泡图（bubble plot），就是三维散点图当中的第三个维度用点的大小来表示，这样便形成了气泡图，因此气泡图是展示三个定量变量间关系的另外一种思路。\u003C\u002Fp\u003E\u003Cp\u003E格式：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Esymbols(x, y, circle=radius)\nsymbols(x, y, circle=sqrt(z\u002Fpi))\u003C\u002Fcode\u003E\u003Cp\u003E其中z为第三个需要绘制的变量，并且，\u003Cb\u003Ecircle\u003C\u002Fb\u003E表示用圆圈来表示气泡，此外还有squares, rectangles, stars, thermometers, boxplots等选项。\u003C\u002Fp\u003E\u003Cp\u003E代码和图形如下：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Eattach(mtcars)\nr &- sqrt(disp\u002Fpi)\n\nsymbols(wt,mpg,circles = r, inches = 0.30,\n
fg = \&white\&, bg = \&lightblue\&,\n
main = \&Bubble Plot with point size proportional to displacement\&,\n
ylab = \&Miles Per Gallon\&,\n
xlab = \&Weight of Car (lbs\u002F1000)\&)\ntext(wt, mpg, rownames(mtcars), cex = 0.6)\ndetach(mtcars)\n\n\u003C\u002Fcode\u003E\u003Cimg src=\&v2-38aef4a902f96ff126adb0e49f8e4eaa.png\& data-rawwidth=\&899\& data-rawheight=\&900\&\u003E\u003Cp\u003E以上图形都可以算是散点图的范畴，虽然散点图比较简单，可能一般的理工科学生在一些实验数据的处理上也会用到过，不过却正是它能够以最为直接的方式来展示数据，因此非常常用，本文也用了大量的篇幅来介绍。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E2 折线图\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E折线图就是将散点图中的所有点用线连接起来而得到，添加了线条以后，便成了一个刻画数据变动的优秀工具。\u003C\u002Fp\u003E\u003Cp\u003E格式：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Eplot(x, y, type = )\nlines(x, y, type = )\u003C\u002Fcode\u003E\u003Cp\u003E其中x和y是要连接的(x, y)点的数值型向量，参数type = 的可选值如下表所示：\u003C\u002Fp\u003E\u003Cimg src=\&v2-d9eaecbf9f79c7e29cf80.png\& data-rawwidth=\&647\& data-rawheight=\&227\&\u003E\u003Cp\u003E下面以基础安装包中的Orange数据集为例，代码和图形如下：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Eopar &- par(no.readonly = TRUE)\npar(mfrow = c(1,2))\nt1 &- subset(Orange, Tree == 1)\nplot(t1$age, t1$circumference,\n
xlab=\&Age (days)\&,\n
ylab=\&Circumference (mm)\&,\n
main=\&Orange Tree 1 Growth\&)\nplot(t1$age, t1$circumference,\n
xlab=\&Age (days)\&,\n
ylab=\&Circumference (mm)\&,\n
main=\&Orange Tree 1 Growth\&,\n
type=\&b\&)\npar(opar)\n\u003C\u002Fcode\u003E\u003Cimg src=\&v2-6b3a26ea010a3fa84a2a7b347af9fd79.png\& data-rawwidth=\&899\& data-rawheight=\&900\&\u003E\u003Cp\u003E下面结合不同的参数（线型、颜色等）在一张图幅上绘制五条折线。\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E#将因子转换为数值型\nOrange$Tree &- as.numeric(Orange$Tree)\nntrees &- max(Orange$Tree)\n\n#构建x，y数据\nxrange &- range(Orange$age)\nyrange &- range(Orange$circumference)\n\n#绘制二维散点图\nplot(xrange, yrange,\n
type = \&n\&,\n
xlab = \&Age (days)\&,\n
ylab = \&Circumference (mm)\&)\n\n#构建颜色、线型、散点的形状的向量，在后续循环当中遍历\ncolors &- rainbow(ntrees)\nlinetype &- c(1:ntrees)\nplotchar &- seq(18, 18+ntrees, 1)\n\n#循环添加设置不同设置参数的线条\nfor (i in 1:ntrees) {\n
tree &- subset(Orange, Tree==i)\n
lines(tree$age, tree$circumference,\n
type=\&b\&,\n
lty=linetype[i],\n
col=colors[i],\n
pch=plotchar[i]\n
)\n}\n\n#设置图名\ntitle(\&Tree Growth\&, \&example of line plot\&)\n\n#添加图例\nlegend(xrange[1], yrange[2],\n
1:ntrees,\n
cex=0.8,\n
col=colors,\n
pch=plotchar,\n
lty=linetype,\n
title=\&Tree\&\n)\n\n\n\n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cimg src=\&v2-c084f28d8d30e51e6edb9d9.png\& data-rawwidth=\&899\& data-rawheight=\&900\&\u003E我们也可以将上述代码中的步骤作为一般绘图的程序，先准备数据，然后创建底图，再添加需要添加的线条（带设置信息的）等其他图形，然后添加图名和图例。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E3 相关图\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E相关图本质上就是对相关系数矩阵进行可视化，以便于更好的发现变量间的相关性、独立性及某种聚集模式，需要用到corrgram包中的corrgram()函数。\u003C\u002Fp\u003E\u003Cp\u003E格式：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Ecorrgram(x, order=, panel=, text.panel=, diag.panel=)\u003C\u002Fcode\u003E\u003Cp\u003E其中x是一行一个的观测数据框，order = TRUE时，相关矩阵将使用主成分分析法（\u003Ca href=\&https:\u002F\u002Fzh.wikipedia.org\u002Fwiki\u002F%E4%B8%BB%E6%88%90%E5%88%86%E5%88%86%E6%9E%90\& data-editable=\&true\& data-title=\&wikipedia.org 的页面\& class=\&\&\u003E主成分分析法\u003C\u002Fa\u003E）对变量进行排序。panel 设定非对角线面板使用的元素类型，而text.panel和diag.panel选项控制着主对角线元素类型。可用的panel选项如下图所示：\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cimg src=\&v2-25e34886eadafe6eddfc7.png\& data-rawwidth=\&652\& data-rawheight=\&238\&\u003E\u003Cp\u003E下面以mtcars包中的数据为例绘制相关图。代码及图形如下：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& install.packages(\&corrgram\&)\n& library(corrgram)\n& corrgram(mtcars, order = TRUE, lower.panel = panel.shade,\n+
upper.panel = panel.pie, text.panel = panel.txt,\n+
main = \&Corrgram of mtcars intercorrelations\&)\n& \n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cimg src=\&v2-9dce76c6fad5fdc2fc58.png\& data-rawwidth=\&899\& data-rawheight=\&900\&\u003E这里有必要对图形进行一些解释，红色表示图形负相关、蓝色表示正相关，颜色的深浅则表示相关性的大小。\u003C\u002Fp\u003E\u003Cp\u003E下面是一个相同数据不同panel设置的相关图，代码和图形如下：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& corrgram(mtcars, order=TRUE, lower.panel=panel.ellipse,\n+
upper.panel=panel.pts, text.panel=panel.txt,\n+
diag.panel=panel.minmax,\n+
main=\&Corrgram of mtcars data using scatter plots\n+
and ellipses\&)\n& \n\u003C\u002Fcode\u003E\u003Cimg src=\&v2-0a147cd286e0e.png\& data-rawwidth=\&899\& data-rawheight=\&900\&\u003E\u003Cp\u003E此处，我们在下三角区域使用平滑拟合曲线和置信椭圆，上三角区域使用散点图。而且，由于档位数必须取3、4或5等离散值，汽缸数也只能取4、6和8等离散值，相应有这两个变量的散点图看起来略显奇怪。\u003Cbr\u003E在结束相关图之前，为了进一步了解corrgram()函数的用法不妨再看几个例子。\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& corrgram(mtcars, lower.panel = panel.shade, \n+
upper.panel = NULL, text.panel = panel.txt,\n+
main = \&Car Mileage Data (unsorted)\&)\n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cimg src=\&v2-8f400cdd7c9a45607b5b.png\& data-rawwidth=\&899\& data-rawheight=\&900\&\u003E在这个例子当中没有对变量进行排序，并且上三角留白，实际上能够给出的信息其实是相似的。\u003C\u002Fp\u003E\u003Cp\u003E下面这个例子是用colorRampPallette()函数自主控制corrgram()函数绘图时的颜色。代码和图形如下：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& cols &- colorRampPalette(c(\&darkgoldenrod4\&,\&burlywood1\&,\n+
\&darkkhaki\&,\&darkgreen\&))\n& corrgram(mtcars, order = TRUE, col.regions = cols,\n+
lower.panel = panel.shade,\n+
upper.panel = panel.conf,\n+
text.panel = panel.txt,\n+
main = \&A Corrgram (or Horse) of a Different Color\&)\n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cimg src=\&v2-e0b6ae7b0dfcaefe20427d.png\& data-rawwidth=\&899\& data-rawheight=\&900\&\u003E\u003Cb\u003E4 马赛克图\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E马赛克图主要用于多变量的类别型变量的可视化。虽然，R基础安装中的mosaicplot()函数也可以绘制马赛克图，考虑到更多的扩展功能，我们更推荐用vcd包中的mosaic()函数。\u003C\u002Fp\u003E\u003Cp\u003E格式：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Emosaic(table)\nmosaic(formula, data = )\u003C\u002Fcode\u003E\u003Cp\u003E其中上面的table是数组形式的列联表，而formula是标准的R表达式，data设定了一个数据框或者表格。\u003C\u002Fp\u003E\u003Cp\u003E下面以基础安装包中的Titanic数据集为例作说明，代码及图形如下：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& library(vcd)\n载入需要的程辑包：grid\n& mosaic(Titanic, shade = TRUE, legend = TRUE)\n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cimg src=\&v2-bdc72fbd6de905bae355ee.png\& data-rawwidth=\&899\& data-rawheight=\&900\&\u003E马赛克图隐含大量的数据信息。每小块方形都可以表示什么性别、什么仓位、是否幸存、年龄分类等数据。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E5 小结\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E本章主要介绍了一些更高级别的数据可视化的工具，包括二维和三维散点图、散点图矩阵、气泡图、折线图、相关系数图和马赛克图。在实践运用当中视情况进行选取。\u003C\u002Fp\u003E&,&updated&:new Date(&T06:01:16.000Z&),&canComment&:false,&commentPermission&:&anyone&,&commentCount&:2,&likeCount&:7,&state&:&published&,&isLiked&:false,&slug&:&&,&isTitleImageFullScreen&:false,&rating&:&none&,&sourceUrl&:&&,&publishedTime&:&T14:01:16+08:00&,&links&:{&comments&:&\u002Fapi\u002Fposts\u002F2Fcomments&},&url&:&\u002Fp\u002F&,&titleImage&:&https:\u002F\u002Fpic3.zhimg.com\u002Fv2-e3f3897dbdb0ae83df5b16_r.jpg&,&summary&:&&,&href&:&\u002Fapi\u002Fposts\u002F&,&meta&:{&previous&:null,&next&:null},&snapshotUrl&:&&,&commentsCount&:2,&likesCount&:7},&&:{&title&:&win 10系统下R语言的升级&,&author&:&realphy&,&content&:&今天在操练《R in Action》第十一章中的实例时惊奇的发现R已经升级到3.3.3版了，作为一个版本控，于是搜了一下R语言的升级步骤，现在综合各方总结给大家。\u003Cp\u003E在ubuntu系统下，非常方便，只要命令行下运行以下代码，如果镜像上面有更新便会自动覆盖安装。\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Esudo apt-get dist-upgrade\u003C\u002Fcode\u003E\u003Cp\u003E在Windows系统下，根据官方说明2.8条（\u003Ca href=\&https:\u002F\u002Fmirrors.tuna.tsinghua.edu.cn\u002FCRAN\u002F\& data-editable=\&true\& data-title=\&The Comprehensive R Archive Network\& class=\&\&\u003EThe Comprehensive R Archive Network\u003C\u002Fa\u003E），步骤如下：\u003C\u002Fp\u003E\u003Cp\u003E1，首先卸载旧版的R，此时会保留安装目录下的library文件夹，这个文件夹里面保存着手动安装的各种包；\u003C\u002Fp\u003E\u003Cp\u003E2，下载最新版的R语言安装文件（\u003Ca href=\&https:\u002F\u002Fcran.r-project.org\u002Fbin\u002Fwindows\u002Fbase\u002FR-3.3.3-win.exe\& class=\&\& data-editable=\&true\& data-title=\&r-project.org 的页面\&\u003Ehttps:\u002F\u002Fcran.r-project.org\u002Fbin\u002Fwindows\u002Fbase\u002FR-3.3.3-win.exe\u003C\u002Fa\u003E），常规安装；\u003C\u002Fp\u003E\u003Cp\u003E3，将旧版R安装目录下的library文件夹拷贝到新版R的安装目录下，部分覆盖此目录下的library文件夹，同文件名的不覆盖，因为同文件的是新版R语言的默认自带安装包；\u003C\u002Fp\u003E\u003Cp\u003E4，此时可以打开Rstudio，顺利的话应该能够识别出新版的R语言，在命令窗口运行\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Eupdate.packages(checkBuilt=TRUE, ask=FALSE)\u003C\u002Fcode\u003E此时会下载更新新版R语言下原旧版R语言所安装的各种包，会有部分更新失败的，此时可以手动用install.packages()函数进行更新。\u003Cp\u003E5，至此，更新就完毕了，可以将旧版的文件夹彻底删除了。\u003C\u002Fp\u003E\u003Cp\u003E除此之外，还有另外一个方法，\u003Cb\u003E使用installr包，\u003C\u002Fb\u003E步骤如下：\u003C\u002Fp\u003E\u003Cp\u003E1，安装installr包，install.packages(\&installr\&)\u003C\u002Fp\u003E\u003Cp\u003E2，载入installr包，library(installr)\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& library(installr)\n\nWelcome to installr version 0.18.0\n\nMore information is available on the installr project website:\nhttps:\u002F\u002Fgithub.com\u002Ftalgalili\u002Finstallr\u002F\n\nContact: &tal.&\nSuggestions and bug-reports can be submitted at: https:\u002F\u002Fgithub.com\u002Ftalgalili\u002Finstallr\u002Fissues\n\n\t\t\tTo suppress this message use:\n\t\t\tsuppressPackageStartupMessages(library(installr))\u003C\u002Fcode\u003E\u003Cp\u003E这个包的作者已经将相关代码在github（\u003Ca href=\&https:\u002F\u002Fgithub.com\u002Ftalgalili\u002Finstallr\u002F\& data-editable=\&true\& data-title=\&talgalili\u002Finstallr\& class=\&\&\u003Einstallr\u003C\u002Fa\u003E）上托管了，在使用上面如果有什么问题可以直接与其联系。\u003C\u002Fp\u003E\u003Cp\u003E3，运行updateR()函数，此时会检查你的版本，如果版本不是最新则会提醒你进行升级，并且建议你在原生的RGui下面进行操作，确认后电脑会自动下载安装文件。\u003C\u002Fp\u003E\u003Cimg src=\&v2-c8023a7fbc30dace4f5f5.png\& data-rawwidth=\&425\& data-rawheight=\&174\&\u003E\u003Cp\u003E不过，可能因为网络的原因，我没有下载成功，Rgui的镜像一直是海外的镜像，我试图改到国内的，可是一直没成功，下载总会半途中断，因此这种方法的更新最终没能成功。另外，上图是我更新成功以后运行updateR()函数的提示。\u003C\u002Fp\u003E&,&updated&:new Date(&T14:28:54.000Z&),&canComment&:false,&commentPermission&:&anyone&,&commentCount&:5,&likeCount&:2,&state&:&published&,&isLiked&:false,&slug&:&&,&isTitleImageFullScreen&:false,&rating&:&none&,&sourceUrl&:&&,&publishedTime&:&T22:28:54+08:00&,&links&:{&comments&:&\u002Fapi\u002Fposts\u002F2Fcomments&},&url&:&\u002Fp\u002F&,&titleImage&:&https:\u002F\u002Fpic1.zhimg.com\u002Fv2-c7e6ab3c31a043e708b27f0df6471919_r.jpg&,&summary&:&&,&href&:&\u002Fapi\u002Fposts\u002F&,&meta&:{&previous&:null,&next&:null},&snapshotUrl&:&&,&commentsCount&:5,&likesCount&:2},&&:{&title&:&《R语言实战》第三部分第十二章-重抽样与自助法学习笔记&,&author&:&realphy&,&content&:&本章主要介绍了两种在统计假设不满足时广泛应用的统计方法：置换检验和自助法。这两种方法依据随机化的思想设计，过去往往只有编程高手或者娴熟的统计专家才有能力使用，如今R语言中已经有对应的软件包。\u003Cp\u003E以下我们将重温一下用传统方法分析过的问题，然后对比看看本章介绍的新方法如何解决它们。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E1 置换检验\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E置换检验（\u003Ca href=\&http:\u002F\u002Fwww.biostat.ucsf.edu\u002Fsen\u002Fstatgen14\u002Fpermutation-tests.htmlpermutation\& class=\&\& data-title=\&permutation test\& data-editable=\&true\&\u003Epermutation test\u003C\u002Fa\u003E），也称随机化检验或重随机化检验。书本上举了一个例子，如果看得不是特别明白可以再看看以下这个例子\u003C\u002Fp\u003E\u003Cp\u003E\u003Ca href=\&https:\u002F\u002Fwww.plob.org\u002Farticle\u002F3176.html\& data-editable=\&true\& data-title=\&Permutation Test 置换检验\& class=\&\&\u003EPermutation Test 置换检验\u003C\u002Fa\u003E\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E首先，它也与之前的参数方法类似，构成并计算一个观测数据的t统计量，并计算出t0，后面就有所不同，也是它的本质所在。它将所有的样本放在一起，重新随机分配的两个不同的组，计算每次随机到的新观测的统计量，如果有N种可能的随机分配组织，则有N个新观测的统计量，我们把每次计算得到的新观测的统计量的分布记为\u003Cb\u003E经验分布。\u003C\u002Fb\u003E最后根据t0落在经验分布中的位置来确定是否拒绝之前的零假设。\u003C\u002Fp\u003E\u003Cp\u003E需要注意的是，虽然置换方法和参数方法都计算了相同的t统计量，但是参数方法是将统计量与理论分布进行比较，置换方法是将统计量与置换观测数据后获得的经验分布进行比较，再根据统计量值得极端性来判断是否有足够的理由拒绝零假设。\u003C\u002Fp\u003E\u003Cp\u003E根据经验分布所依据的数据是否是所有可能的排列组合，置换检验可以分为“精确”检验和非“精确”检验，对于后者可以使用蒙特卡罗模拟（\u003Ca href=\&http:\u002F\u002Fwww.ruanyifeng.com\u002Fblog\u002F\u002Fmonte-carlo-method.html\& class=\&\& data-editable=\&true\& data-title=\&蒙特卡罗方法入门 - 阮一峰的网络日志\&\u003E蒙特卡罗方法入门 - 阮一峰的网络日志\u003C\u002Fa\u003E），从所有可能的排列当中进行抽样，获得一个近似的检验，特别是在样本量很大，减少可能的计算开销的情况下。\u003C\u002Fp\u003E\u003Cp\u003E在介绍完了具体的原理以后，下面结合R中提供的两个软件包来具体实践一下置换检验。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E2 用coin包做置换检验\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E包的安装：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E&install.packages(c(\&coin\&))\n& install.packages(\&lmPerm)\u003C\u002Fcode\u003E\u003Cp\u003E幸运的是，lmPerm包已经可以更新，而且版本也更新到了2.1，目前是Marco Torchiano在维护。在此向原作者Bob Wheeler致意！\u003C\u002Fp\u003E\u003Cp\u003Ecoin包提供了一个针对独立性问题进行置换检验的一般性框架，主要回答以下问题：\u003C\u002Fp\u003E\u003Cp\u003Ea: 响应值与组的分配独立吗？\u003C\u002Fp\u003E\u003Cbr\u003E\u003Cp\u003Eb: 两个数值变量独立吗？\u003C\u002Fp\u003E\u003Cbr\u003E\u003Cp\u003Ec: 两个类型变量独立吗？\u003C\u002Fp\u003E\u003Cbr\u003E\u003Cp\u003E以下为可选置换检验的coin函数：\u003C\u002Fp\u003E\u003Cp\u003E\u003Cimg src=\&v2-7bd95a70a3ba2e7a9d1e16b1c8751038.png\& data-rawwidth=\&714\& data-rawheight=\&343\&\u003E函数格式：\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003Efunction_name(formula, data, distribution = )\u003C\u002Fcode\u003E\u003Cp\u003E其中：\u003C\u002Fp\u003E\u003Cp\u003Ea: formula描述待检验变量间的关系，而且类别型变量和序数变量必须分别转化为因子和有序因子\u003C\u002Fp\u003E\u003Cp\u003Eb: data是一个数据框，而且必须得是数据框\u003C\u002Fp\u003E\u003Cbr\u003E\u003Cp\u003Ec: distribution指定经验分布在零假设条件下的形式，有exact、asymptotic和approximate三个选项。\u003C\u002Fp\u003E\u003Cbr\u003E\u003Cp\u003E若distribution=\&exact\&，那么在零假设条件下，分布的计算是精确的（即依据所有可能的排列组合）。当然，也可以根据它的渐进分布（distribution=\&asymptotic\&）或蒙特卡洛重抽样（distribution=\&approxiamate(B=#)\&）来做近似计算，其中#指所需重复的次数。distribution=\&exact\&当前仅可用于两样本问题。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E2.1 独立两样本和K样本检验\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& library(coin)\n载入需要的程辑包：survival\n& score &- c(40, 57, 45, 55, 58, 57, 64, 55, 62, 65)\n& treatment &- factor(c(rep(\&A\&, 5),rep(\&B\&,5)))\n& treatment\n [1] A A A A A B B B B B\nLevels: A B\n& mydata &- data.frame(treatment, score)\n& t.test(score~treatment,data = mydata, var.equal = TRUE)\n\n\tTwo Sample t-test\n\ndata:
score by treatment\nt = -2.345, df = 8, p-value = 0.04705\nalternative hypothesis: true difference in means is not equal to 0\n95 percent confidence interval:\n -19.0405455
-0.1594545\nsample estimates:\nmean in group A mean in group B \n
60.6 \n\n& oneway_test(score~treatment, data = mydata, distribution = \&exact\&)\n\n\tExact Two-Sample Fisher-Pitman Permutation Test\n\ndata:
score by treatment (A, B)\nZ = -1.9147, p-value = 0.07143\nalternative hypothesis: true mu is not equal to 0\n\n& \u003C\u002Fcode\u003E\u003Cp\u003E\u003Cb\u003E传统t检验表明存在显著性差异（p&0.05），而精确检验却表明差异并不显著（p&0.072）。\u003Cbr\u003E\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& library(MASS)\n& UScrime &- transform(UScrime, So = factor(So))\n& wilcox_test(Prob ~ So, data = UScrime, distribution = \&exact\&)\n\n\tExact Wilcoxon-Mann-Whitney Test\n\ndata:
Prob by So (0, 1)\nZ = -3.7493, p-value = 8.488e-05\nalternative hypothesis: true mu is not equal to 0\n\n& wilcox.test(Prob ~ So, data = UScrime)\n\n\tWilcoxon rank sum test\n\ndata:
Prob by So\nW = 81, p-value = 8.488e-05\nalternative hypothesis: true location shift is not equal to 0\n\u003C\u002Fcode\u003E此处wilcox_test()计算结果与第七章的wilcox.test()函数是一致的，这是因为exact也是wilcox.test()函数的默认算法。\u003Ccode lang=\&text\&\u003E& library(multcomp)\n载入需要的程辑包：mvtnorm\n载入需要的程辑包：survival\n载入需要的程辑包：TH.data\n载入需要的程辑包：MASS\n\n载入程辑包：‘MASS’\n\nThe following object is masked _by_ ‘.GlobalEnv’:\n\n
UScrime\n\n\n载入程辑包：‘TH.data’\n\nThe following object is masked from ‘package:MASS’:\n\n
geyser\n\n& set.seed(1234)\n& library(coin)\n& oneway_test(response ~ trt, data = cholesterol, distribution = approximate(B=9999))\n\n\tApproximative K-Sample Fisher-Pitman Permutation Test\n\ndata:
response by\n\t trt (1time, 2times, 4times, drugD, drugE)\nchi-squared = 36.381, p-value & 2.2e-16\n\n& \u003C\u002Fcode\u003E此处需要注意两点，一个是此时的参考分布得自于数据9999次的置换，另外设置了随机数种子，确保每次结果的一致性。\u003Cp\u003E\u003Cb\u003E2.2 列联表的独立性\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003E通过chisq_test()或cmh_test()函数，我们可用置换检验判断两类别型变量的独立性。\u003Cbr\u003E\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& library(coin)\n载入需要的程辑包：survival\n& library(vcd)\n载入需要的程辑包：grid\n& Arthritis &- transform(Arthritis, Improved = as.factor(as.numeric(Improved)))\n& set.seed(134)\n& set.seed(1234)\n& chisq_test(Treatment ~ Improved, data = Arthritis, distribution = approximate(B=9999))\n\n\tApproximative Pearson Chi-Squared Test\n\ndata:
Treatment by Improved (1, 2, 3)\nchi-squared = 13.055, p-value = 0.0018\n\n& \n\u003C\u002Fcode\u003E\u003Cp\u003E经过9999次的置换，获得一个近似的卡方检验。在这个例子当中，将Improved从一个有序因子变成一个分类因子，这是因为如果你用有序因子，则coin()将生成一个线性也线性趋势检验，而不是卡方检验。Rstudio 提示如下：\u003Cbr\u003E\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& class(Arthritis$Improved)\n[1] \&factor\&\n& Arthritis &- transform(Arthritis, Improved = as.numeric(Improved))\n& class(Arthritis$Improved)\n[1] \&numeric\&\n& chisq_test(Treatment ~ Improved, data = Arthritis, distribution = approximate(B=9999))\nError in check(object) : \n
‘object’ does not represent a contingency problem\n& \n\u003C\u002Fcode\u003E\u003Cp\u003E\u003Cb\u003E2.3 数值变量间的独立性\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003Espearman_test()函数提供了两数值变量的独立性置换检验。\u003Cbr\u003E\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& states &- as.data.frame(state.x77)\n& set.seed(1234)\n& spearman_test(Illiteracy ~ Murder,data = state.x77, distribution = approximate(B=9999))\nError in model.frame.default(formula = ~Illiteracy, data = c(,
'data' must be a data.frame, not a matrix or an array\n& spearman_test(Illiteracy ~ Murder,data = state, distribution = approximate(B=9999))\nError in terms.formula(formula, data = data) : object 'state' not found\n& spearman_test(Illiteracy ~ Murder,data = states, distribution = approximate(B=9999))\n\n\tApproximative Spearman Correlation Test\n\ndata:
Illiteracy by Murder\nZ = 4.7065, p-value & 2.2e-16\nalternative hypothesis: true rho is not equal to 0\n\n& \n\u003C\u002Fcode\u003E\u003Cp\u003E基于9999次的重复的近似检验可知，独立性并不被满足。另外，细心的读者发现，我在输入的时候犯了两次错误，不难发现，在coin包中，可接受的数据类型是数据框。\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E2.4 两样本和K样本相关性检验\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& library(coin)\n& library(MASS)\n& wilcoxsign_test(U1~U2, data = UScrime, distribution = \&exact\&)\n\n\tExact Wilcoxon-Pratt Signed-Rank Test\n\ndata:
y by\n\t x (pos, neg) \n\t stratified by block\nZ = 5.9691, p-value = 1.421e-14\nalternative hypothesis: true mu is not equal to 0\n\n& \n\u003C\u002Fcode\u003E\u003Cp\u003E当处于不同组的观测已经被分配得当，或者使用了重复测量时，样本相关检验便可派上用场。对于两配对组的置换检验， 可使用wilcoxsign_test() 函数； 多于两组时， 使用friedman_test()函数。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E2.5 深入研究\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Cp\u003Ecoin包提供了一个置换检验的一般性框架，可以分析一组变量相对于其他任意变量，是否与第二组变量（可根据一个区组变量分层）相互独立。而且independence_test()函数可以让我们从置换角度来思考大部分传统检验，进而在面对无法用传统方法解决的问题时，使用户可以自己构建新的统计检验。当然，这种灵活性也是有门槛的：要正确使用该函数必须具备丰富的统计知识。更多函数细节请参阅包附带的文档（运行vignette(\&coin\&)即可）。vignette(\&包的名称\&)可以下载并用pdf阅读软件打开各个包的文档。\u003C\u002Fp\u003E\u003Cp\u003E下一节学习的lmPerm包，提供了线性模型（包括回归和方差分析）的置换方法。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E3 lmPerm包的置换检验\u003C\u002Fb\u003E\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E在lmPerm包中有lm()函数和aov()函数的修改版：lmp()和aovp()，进行置换检验，而非正态检验。参数也仅在原函数的基础上添加了perm=参数。perm参数一共有三个选项，意义分别如下：\u003C\u002Fp\u003E\u003Cul\u003E\u003Cli\u003EExact: 根据所有可能的排列组合生成精确检验；\u003Cbr\u003E\u003C\u002Fli\u003E\u003Cli\u003EProb: 从所有可能的排列中不断抽样，直至估计的标准差在估计的p值0.1之下，判停准则由可选的Ca参数控制；\u003Cbr\u003E\u003C\u002Fli\u003E\u003Cli\u003ESPR: 使用贯序概率比检验来判断何时停止抽样。\u003Cbr\u003E\u003C\u002Fli\u003E\u003C\u002Ful\u003E\u003Cbr\u003E\u003Cp\u003E为了在样本数过大时减小计算机运行时间，精确建议仅在样本数小于10才使用。因此，如果观测数大于10，perm=\&Exact\&将自动默认转为perm=\&Prob\&。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E3.1 简单回归和多项式回归\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& library(lmPerm)\n& set.seed(1234)\n& fitPerm &- lmp(weight ~ height, data = women, perm = \&Prob\&)\n[1] \&Settings:
unique SS : numeric variables centered\&\n& summary(fitPerm)\n\nCall:\nlmp(formula = weight ~ height, data = women, perm = \&Prob\&)\n\nResiduals:\n
Max \n-1.3 -0.7
3.1167 \n\nCoefficients:\n
Estimate Iter Pr(Prob)
&2e-16 ***\n---\nSignif. codes:
\n0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n\nResidual standard error: 1.525 on 13 degrees of freedom\nMultiple R-Squared: 0.991,\tAdjusted R-squared: 0.9903 \nF-statistic:
1433 on 1 and 13 DF,
p-value: 1.091e-14 \n\n& \n\u003C\u002Fcode\u003E\u003Ccode lang=\&text\&\u003E& fitPerm &- lmp(weight~height + I(height^2), data=women, perm=\&Prob\&)\n[1] \&Settings:
unique SS : numeric variables centered\&\n& summary(fitPerm)\n\nCall:\nlmp(formula = weight ~ height + I(height^2), data = women, perm = \&Prob\&)\n\nResiduals:\n
3Q \n-0....286151 \n
Max \n 0.597059 \n\nCoefficients:\n
Estimate Iter Pr(Prob)
&2e-16 ***\nI(height^2)
&2e-16 ***\n---\nSignif. codes:
\n0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n\nResidual standard error: 0.3841 on 12 degrees of freedom\nMultiple R-Squared: 0.9995,\tAdjusted R-squared: 0.9994 \nF-statistic: 1.139e+04 on 2 and 12 DF,
p-value: & 2.2e-16 \n\n& \n\u003C\u002Fcode\u003E\u003Cp\u003E以上分别为简单回归和多项式回归的代码，用置换检验来检验这些回归是非常容易的，换一个仅多一个Perm参数的函数即可。从summary()函数的输出可以看到增添的Iter列给出了达到判停标准所需要的迭代次数。\u003Cbr\u003E\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E3.2 多元回归\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& library(lmPerm)\n& set.seed(1234)\n& states &- as.data.frame(state.x77)\n& fitM &- lmp(Murder ~ Population + Illiteracy + Income + Frost, data = states, perm = \&Prob\&)\n[1] \&Settings:
unique SS : numeric variables centered\&\n& summary(fitM)\n\nCall:\nlmp(formula = Murder ~ Population + Illiteracy + Income + Frost, \n
data = states, perm = \&Prob\&)\n\nResiduals:\n
Max \n-4.746 -0.050
7.62104 \n\nCoefficients:\n
Estimate Iter Pr(Prob)
\nPopulation 2.237e-04
\nIlliteracy 4.143e+00 5000
0.0004 ***\nIncome
\n---\nSignif. codes:
\n0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n\nResidual standard error: 2.535 on 45 degrees of freedom\nMultiple R-Squared: 0.567,\tAdjusted R-squared: 0.5285 \nF-statistic: 14.73 on 4 and 45 DF,
p-value: 9.133e-08 \n\n& \n\u003C\u002Fcode\u003E\u003Cp\u003E与第八章的结果相比，Population不再显著。当两种方法的结果不一致时，需要更加谨慎的审视数据，可能是违反了相关前提假设或者存在离群点。\u003C\u002Fp\u003E\u003Cp\u003E\u003Cb\u003E3.3 单因素方差分析和协方差分析\u003C\u002Fb\u003E\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& library(lmPerm)\n& library(multcomp)\n载入需要的程辑包：mvtnorm\n载入需要的程辑包：survival\n载入需要的程辑包：TH.data\n载入需要的程辑包：MASS\n\n载入程辑包：‘MASS’\n\nThe following object is masked _by_ ‘.GlobalEnv’:\n\n
UScrime\n\n\n载入程辑包：‘TH.data’\n\nThe following object is masked from ‘package:MASS’:\n\n
geyser\n\n& set.seed(1234)\n& fit &- aovp(response ~ trt, data = cholesterol, perm = \&Prob\&)\n[1] \&Settings:
unique SS \&\n& anova(fit)\nAnalysis of Variance Table\n\nResponse: response\n
Df R Sum Sq R Mean Sq Iter
337.84 5000 & 2.2e-16 ***\nResiduals 45
\n---\nSignif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n& \u003C\u002Fcode\u003E\u003Cp\u003E以上为9.1节中单因素方差分析问题，即各种疗法对降低胆固醇的影响的置换检验，结果表明疗法的效果不全相同。\u003C\u002Fp\u003E\u003Ccode lang=\&text\&\u003E& fit &- aovp(weight ~ gesttime + dose, data = litter, perm = \&Prob\&)\n[1] \&Settings:
unique SS : numeric variables centered\&\n& anova(fit)\nAnalysis of Variance Table\n\nResponse: weight\n
Df R Sum Sq R Mean Sq Iter Pr(Prob)
\ngesttime
161.493 60 ***\ndose
45.708 31 *
\nResiduals 69
\n---\nSignif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\u003C\u002Fcode\u003E\u00}