导语：
尽可能地挖掘出数据中有用的信息，是对数据提供者最大的尊重。因为很多时候数据提供者（有时是自己）是付出了巨大的牺牲才获得这些数据。
完成一篇高质量的文章（如CNS，PNAS或者单个领域的顶刊），需要一个（1）有新意的idea， （2）与该idea相对应的完美的实验或数据收集，（3）深入的统计分析【挖掘出数据中的规律】，通常3者缺一不可。此文主要关注第3步：深入的统计分析。
什么是“深入分析”呢？
这里，我们可以简单理解为“找到规律”，最好是一类现象通用的规律，规律可以是：
- 研究对象的时空变化规律，即研究对象哪里高（where, 空间分布规律），什么时候高（when, 时间分布规律），哪些类型的对象高（which 或者 who, 需要某些统计检验验证）；
- 为什么会在这里高，这个时候高，这些类型的高 (why)，即找到哪些因素（温度？盐度？气压？光照？竞争？）造成了这种状况及各因素的影响程度（重要性分析）；
- 最终我们还需要知道这些影响因素以什么样的形式定量地影响研究对象(how)？直线关系，抛物线关系，S型曲线关系，线性整流函数式曲线关系（满足一定阈值后开始影响）还是其他类型的曲线关系，该直线/曲线上的关键点在哪？
天行有常，不为尧存，不为桀亡。虽然大自然的运行有其内在的规律，但这个“规律”很多时候并不那么明显，而需要我们通过已有的现象或数据去慢慢揭示，科研工作可以看成一点一点，一步一步揭示大自然/宇宙中的规律的过程，即上文中的[2]和[3]。
PNAS案例分析
下面将以北京大学医学部的任爱国教授于2011年发表在PNAS上的一篇文章为例（下图），介绍如何应用开源的（免费）R语言实现该文章中的统计分析：
（1）Mann-Whitney U 检验
（2）Pearson卡方检验
（3）无条件逻辑回归[Unconditional Logistic Regression]
介绍此3个分析方法的过程中，作者也将简要介绍如何通过R语言做出漂亮的分析图。
Mann–Whitney U test (也叫 Wilcoxon rank-sum检验， Mann–Whitney–Wilcoxon (MWW)检验 或者 Wilcoxon–Mann–Whitney检验) ，作者应用此方法检验处理组与对照组的中位数是否有显著性差异（非参检验）。代码中为模块2第28行，(),就这一行，是不是非常简单。
模块2中1-26行为随机生成2个元组（array）, 做核密度估计（KDE）分布图（用于视觉感官上看两者的差异），最后代入()函数检验两者的中位数是否有显著性差异。
注意事项：Mann-Whitney U Test使用的各个数的秩（Rank，或者顺序），应用Excel实施该检验的步骤可参考如下链接： http://www.real-statistics.com/non-parametric-tests/mann-whitney-test/
第一步：导入数据分析与作图所需要库设置工作目录
library(ks) #plot kde2
library(Amelia) #可视化数据缺失情况
library(pscl) #pR2 for goodness of fitting
library(ROCR) #ROC曲线
cat("\014") # clear the console
rm(list=ls()); # clear past variables
#***设置工作工作目录(working directory)***#
#*********** 方便后续读写文件 ************#
setwd("XX/Desktop/GlobalStat");
第二步： Mann-Whitney U Test
#*********** 产生两个随机数序列 ***********#
#***** 用()保证结果的可重复性 ****#
(1);
x <- runif(15, min=0, max=10);
print(x);
(2);
y <- runif(15, min=2, max=12);
print(y);
#********两个样本的分布 Kernel Density Estimation *******#
fhat <- kde(x=x);
fhat2 <- kde(x=y);
#***作图：将生成的图片保存为, 300 dpi***#
jpeg("", width = 6, height = 6, units = "in", res = 300, type = "cairo");
plot(fhat, col="blue", xlim = c(-10,20), xaxs = "i", ylim = c(0,), yaxs = "i");
par(new=TRUE); #画图于此画布
plot(fhat2, col="red",xlim = c(-10,20), xaxs = "i", ylim = c(0,), yaxs = "i", main="");
title(main = "核密度估计", family="SimHei") #中文字体
grid();
legend("topleft", legend=c("VarX", "VarY"), col=c("blue", "Red"), lty=1, cex=1); #添加图例
() #作图完成后一定不要忘记该代码
#***** Wilcoxon Rank Sum Test or Mann-Whitney U test *****#
(x,y);
第三步：Pearson卡方检验
卡方检验：作者应用该方法检验了处理组与对照组的不同群体间（如孕妇年龄段，怀孕早期是否有被动吸烟史，是否感染感冒等）的比例是否有显著性差异（Table 1）。R语言中卡方检验的代码也非常简单，()[chi square test的缩写]模块3 第5行和第11行
注意事项：
卡方检验只适用于分类数据，如性别 （男、女） 或颜色 （红、绿、蓝） 等， 不适用于数值数据，例如身高、体重、价格等（但转化为范围后，如>180cm 与≤180cm, 适用该方法）。
计算步骤与该检验的详细介绍与应用可参考“数学乐”网站的如下链接：https://www.shuxuele.com/data/chi-square-test.html
下面代码中使用的数据选自文中的Table 1.
#******************** 卡方检验(Maternal Age) *********************#
MFever <- as.table(rbind(c(25, 50), c(4, 46)));
dimnames(MFever) <- list(Treat = c("Cases", "Controls"),
MaternalAge = c("Yes","No"));
(XsqFever <- (MFever)); # Print test summary
#************************ 卡方检验(Fever) ************************#
M <- as.table(rbind(c(39, 29, 32), c(32, 23, 45)));
dimnames(M) <- list(Treat = c("Cases", "Controls"),
MaternalAge = c("<25","25-29", ">=30"));
(XsqAge <- (M)); # Print test summary
XsqAge$observed # observed counts (same as M)
XsqAge$expected # expected counts under the null
XsqAge$residuals # Pearson residuals
XsqAge$stdres # standardized residuals
第5行与第11行代码将输出
Pearson's Chi-squared test with Yates' continuity correction
data: MFever
X-squared = 9.4308, df = 1, p-value = 0.002134
Pearson's Chi-squared test
data: M
X-squared = 3.5773, df = 2, p-value = 0.1672
第四步： 逻辑斯蒂回归
（无条件）逻辑斯蒂回归：作者应用此方法探究了不同因素与污染物含量和神经管缺陷（NTDs）的关联，即用于辅助确定某个（些）因素是否与NTDs相关。
模型拟合部分仅为模块4第23-24行：
logmodel <- glm(Survived ~.,family=binomial(link='logit'),
data=train);
其余代码主要为数据的预处理和拟合结果优劣评估和可视化，是不是也非常简单。
注意事项：
无条件逻辑斯蒂回归中的“无条件”是相对于“有条件”而言。在某些二分类算法中，如垃圾邮件检测和癌症检测，通常垃圾邮件或癌症案例的比例较低（5%左右或稍高一些），在建模过程中我们需要对这一类别加以比较高的权重，以使两个分类相对均衡。该文中，不同分类的比例差异不大，所以作者使用的无条件逻辑回归。
下面的案例使用了Kaggle的泰坦尼克数据， 里面包含了乘客的相关信息与是否生还的结果，是用于构建二分类模型的典型案例。以下代码参考自如下链接：How to perform a Logistic Regression in R
#******* (unconditional) logistic regression *******#
#****** 使用Kaggle 的泰坦尼克数据 *********#
# https://www.kaggle.com/c/titanic/data #
traindata <- read.csv('',header=T,(""));
#* 查看数据完整性/是否有缺失
sapply(traindata,function(x) sum((x)));
sapply(traindata, function(x) length(unique(x)));
missmap(traindata, main = "Missing values vs observed");
data <- subset(traindata,select=c(2,3,5,6,7,8,10,12));
data$Age[(data$Age)] <- mean(data$Age,);
#***** 查看Sex与Embarked是否为分类变量(factor) *****#
(data$Sex);
(data$Embarked);
contrasts(data$Sex); #查看dummyfied之后的结果
#********** 去除Embarked记录缺失行（2行）**********#
data <- data[!(data$Embarked),];
rownames(data) <- NULL;
#*************** 下面为模型拟合部分 **************#
train <- data[1:800,];
test <- data[801:889,];
logmodel <- glm(Survived ~.,family=binomial(link='logit'),
data=train);
summary(logmodelk);
#***** 查看模型偏差 (Analysis of Deviance) *****#
anova(logmodel, test="Chisq");
#***** 计算模型McFadden R2来评估拟合效果 ******#
pR2(logmodel);
#检验模型预测能力
Pred <- predict(logmodel,newdata=subset(test,select=c(2,3,4,5,6,7,8)),
type='response');
Pred <- ifelse(Pred > );
misClasificError <- mean(Pred != test$Survived);
print(paste('Accuracy',1-misClasificError));
#*** ROC曲线 ***#
p <- predict(logmodel, newdata=subset(test,select=c(2,3,4,5,6,7,8)),
type="response");
pr <- prediction(p, test$Survived);
prf <- performance(pr, measure = "tpr", x.measure = "fpr");
#********* 计算AUC(Area under curve) *********#
auc <- performance(pr, measure = "auc");
auc <- auc@[[1]];
print(auc);
AUCText <- paste0("AUC = ", round(auc,3));
#**作图：将生成的图片保存为, 300 dpi, 6in × 6in **#
jpeg("", width = 6, height = 6, units = "in",
res = 300, type = "cairo");
# mgp参数: The margin line (in mex units) for the
# [1]坐标轴标题, [2]坐标轴刻度标签 and [3]坐标轴刻度.
par(mar = c(4,4,3,3)+0.1, mgp = c(2.5, 1, 0));
plot(prf, col="blue",xlim = c(0,1),ylim = c(0,1),
xaxt="n", yaxt="n")
grid()
plot(prf, col="blue", main = "",
xlim = c(0,1),ylim = c(0,1),
xaxt="n",yaxt="n", add=TRUE)
title(main = "ROC Curve")
xtk <- seq(0, 1, by = )
ytk <- seq(0, 1, by = )
axis(side=1, at = xtk)
axis(side=2, at = ytk)
text(0.8, 0.8, AUCText, pos = 3, cex = 1.0)
()
如果您喜欢这篇文章，希望您能花一秒时间留下您的小手印 Thanks(･ω･)ﾉ
转播
