我是如何做生统作业的-2014试题
第一题
研究者希望分析一个地区的Acute Keshan disease。健康人群和病人的phosphorus content结果放在了“Keshan_disease.txt”里，第一行是病人的结果，第二行是健康人群的结果。我们需要知道这两者是否有显著性的差异。
a. 在试验之前，你需要描述一下这两组之间的差异。比如说，平均值，画个箱图
b. 假设两组数据都来自于正态总体，检验一下均值是否有差异。你还需要检验两者的方差是否相同？
c. 假设没有第a题的假设，使用非参数方法去检验一下。
d. 比较一下第b题和第c题的结果。
明确一下关键字： 描述性统计分析，方差检验，t检验，非参数检验。
请拿出《R语言实战》第二版，跟我一起翻书一起做。
第a题请翻到第7章：基本统计分析。
首先是用summary了解他们的均值和四分位数
Keshan_disease_data <- read.csv("Keshan_disease.csv")
summary(Keshan_disease_data)
# patient healthy
# Min. :1.600 Min. :1.280
# 1st Qu.:4.228 1st Qu.:2.045
# Median :4.730 Median :3.675
# Mean :4.652 Mean :3.355
# 3rd Qu.:5.385 3rd Qu.:4.228
# Max. :6.530 Max. :5.780我还需要知道一下他们的标准差
apply(Keshan_disease_data,2,sd)
# patient healthy
#1.058033 1.298421各种统计函数，数学函数请翻到R语言实战的高级数据管理。这里划重点，我觉得apply是必考的内容，具体说明点传送门，不懂好好去看，看了还是觉得不懂的，你只能请我喝咖啡了。
画个箱图超级简单,请翻到R语言实战的基本图形：
boxplot(Keshan_disease_data)第b题，首先要求我们检验方差是否相同, 用到是var.test。划重点，一定要先写你的假设，不然都不叫假设检验了。
H0: 健康的人样本的方差和病人样本的方差相同
Ha: 两组的方差不同
var.test(Keshan_disease_data$patient, Keshan_disease_data$healthy)
# F test to compare two variances
#data: Keshan_disease_data$patient and Keshan_disease_data$healthy
#F = 0.664, num df = 39, denom df = 39, p-value = 0.2055
#alternative hypothesis: true ratio of variances is not equal to 1
#95 percent confidence interval:
# 0.3511884 1.2554349
#sample estimates:
#ratio of variances
# 0.6639987由于P值> 0.05,我们有理由接受原假设(听了统计的文献课，我都不知道如何看待P值了)，因此方差相同。
然后我们就可以做t.test。不知道t.test怎么用怎么办，用help(test)，
Default S3 method:
t.test(x, y = NULL,
alternative = c(“two.sided”, “less”, “greater”),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95, …)
我们需要根据自己的要求修改参数，判断是否有差异，所以 alternative = c(“two.sided”), 方差相等，所以var.equal = TRUE。
还是记住，先把假设写了
H0: 两组的均值没有差异
Ha: 两组的均值存在差异
t.test(Keshan_disease_data$patient, Keshan_disease_data$healthy,
alternative = c("two.sided"), var.equal = TRUE)
# Two Sample t-test
# data: Keshan_disease_data$patient and Keshan_disease_data$healthy
# t = 4.8947, df = 78, p-value = 5.2e-06
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
# 0.7690202 1.8234798
# sample estimates:
# mean of x mean of y
# 4.65150 3.35525由于p=5.2e-06 小于0.01， 所以我们有把握拒绝原假设，认为两组的均值有差异（还是不知道如何看待P值）。
第c题的非参数检验，建议用wilcoxon rank sum或 binomial test， 那我就用wilcoxon rank sum吧。请翻到基本统计分析的组间差异的非参数检验，方法就是wilcox.test，不知道怎么用，请用help(wilcox.test).
注意：写下你的假设，我这里就不写了。
wilcox.test(Keshan_disease_data$patient, Keshan_disease_data$healthy,
alternative = c("two.sided"))
# Wilcoxon rank sum test with continuity correction
# data: Keshan_disease_data$patient and Keshan_disease_data$healthy
# W = 1252.5, p-value = 1.358e-05
# alternative hypothesis: true location shift is not equal to 0
# Warning message:
# In wilcox.test.default(Keshan_disease_data$patient, Keshan_disease_data$healthy, :
# 无法精確計算带连结的p值发现了一个warning message， 你是不是有点慌，不要担心，只是警告而已，又不是error。如果不想看到warning添加一个exact = FALSE。
第d题就可以说，非参数的符号秩方法算的P值比t检验大一点，但是依旧是否定原假设哦。
第二题
“anemia.txt”给出了30个病人的anemia数据，我们希望用线性回归建立reticulyte和lymphocyte之间的模型。
注：lymphocyte是响应变量（也就是y了），请给出线性模型和p值。
关键字：线性模型
请翻到第8章：回归，用lm()拟合回归模型,对了记得看162页的函数哦。
anemia_data <- read.csv("anemia.csv")
anemia_fit <- lm(lymphocyte~reticulyte, data=anemia_data)
summary(anemia_fit)
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -183.54 119.76 -1.533 0.137
# reticulyte 511.49 16.58 30.844 <2e-16 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# Residual standard error: 204.4 on 28 degrees of freedom
# Multiple R-squared: 0.9714, Adjusted R-squared: 0.9704
# F-statistic: 951.3 on 1 and 28 DF, p-value: < 2.2e-16所以p-value: < 2.2e-16,reticulyte的系数是511.49，那么模型就是
lymphocyte = 511.49 * reticulyte - -183.54
你要是觉得太单调了，可以对回归进行诊断哦，花一些图看看
par(mfrow = c(2,2))
plot(anemia_fit)
par(mfrow=c(1,1))
plot(anemia_data$reticulyte, anemia_data$lymphocyte)
abline(anemia_fit)你知道如何解释么？
第三题：
这次问题要用到datasets.txt， 数据来源于某个疾病的基因芯片。第一行是样品名（S1-S20）第一列是基因名（G1-G100)。数字都是表达量情况
a: G3基因在20个样本的密度图，（PDF哦），然后计算他的最小值，中位数，和方差。
b: 所有样本的基因分布的箱图。你需要检查一下离群点，然后删掉它重新来一遍。
这道题已经开始考验你的数据管理能力了，不太熟练的童鞋请到基本数据管理好好把代码敲一遍，当然我也写了专门的文章，点击传送门
A： 首先我们要用[]操作符提取数据
exp_data <- read.csv("datasets.csv", stringsAsFactors = F)
head(exp_data)
G3 <- exp_data[exp_data$X == "G3", ]
G3
# X S1 S2 S3 S4 S5 S6 S7 S8 S9
# 3 G3 9.881114 11.35708 10.27059 9.643856 10.37759 11.12393 10.55358 9.433481 9.793034
# S10 S11 S12 S13 S14 S15 S16 S17 S18
# 3 9.717762 10.13699 10.46827 9.600006 9.960798 10.02791 9.271463 10.46148 9.326767
# S19 S20
# 3 11.02746 11.05728结果G3是一个数据库， 我希望他是一个matrix或是一个vector，才可以让我方便作图，所以就要用到转置了。
G3.vector <- t(G3[2:length(G3)])那么问题来了，密度图，怎么画，急，在线求！！！请翻到120页核密度图
d <- density(G3.vector)
plot(d)min(G3.vector)
# [1] 9.271463
median(G3.vector)
# [1] 10.08245
var(G3.vector[,1])
# [1] 0.3880644PS: 核密度估计是用来估计随机变量概率密度函数的一种非参数方法。
B： 然后我们需要画箱图看看离群点（坑爹呢，输入法老提示利群店！！这广告也太明显吧！！）
boxplot(exp_data[,2:21])这个离群点还真不是一般的显著明显呀。先看这个数值大小
> exp_data$S7[exp_data$S7 > 80]
[1] 105
> exp_data$S8[exp_data$S8 > 80]
[1] 105估计是少打一个小数点把，没办法题目说删除，那我们也没有办法了，根据第72页做。
exp_data$S7[exp_data$S7 > 80] <- NA
exp_data$S8[exp_data$S8 > 80] <- NA
boxplot(exp_data[,2:21])第四题：
蛋白A被人是肺癌病人预兆的预测信号。因此他们从三个医院收集了病人和健康的人数据，包括蛋白A的表达量，还有有beta-actin（持家基因），治疗后的存活时间等。数据在lung_cancer.txt。可以看instruction.pdf看详细细节。
a. 单因素ANOVA检验是否数据来源对蛋白的A表达量有影响
b. 一些研究者认为蛋白A表达水平是否异常与status没有关系。基于Alpha医院数据，我们可以否定这个结论么
提示： logistic回归
关键字：ANOVA， logistic
第A题： 请记得写假设哦
H0 : 不同医院来源的数据的蛋白A表达水平没有差异
Ha： 不同医院来源的数据的蛋白A表达水平存在差异
这题我建议你可以先去看下我写的方差分析（良心之作），传送门点我。做方差分析脑子一定要画好一张表，如下
| 来源 | alpha | beta | gamma |
|---|---|---|---|
| x | 1 | 2 | 3 |
| … | … | … | … |
然后就可以做aov做方差分析了
PA <- read.table("lung_cancer.txt", header =T, sep = '\t', row.names = 1)
pa.fit <- aov(exp_A ~ as.factor(hospital), data=PA)有童鞋会问，为什么有一个as.factor呢？ 这是因为反差分析是研究连续型响应变量和离散型相应变量之间的关系。导入数据的时候，R默认不把数值型变量当做因子，所以需要自己把它转成因子了哦。
summary(pa.fit)
# Df Sum Sq Mean Sq F value Pr(>F)
# as.factor(hospital) 2 44670 22335 31.65 1.2e-12 ***
# Residuals 197 139010 706
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1p值小于0.01，说明不同医院来源的蛋白A表达量是有差异的。由于是单因素方差分析，所以不需要考虑到病人和健康的人之间的差异。
第B题：这题需要我们提出aplha医院的status数据做logistic回归分析。可以去看我的logistic回归的讲解哦（后台回复生统作业）。R语言实战翻到广义线性模型
记得写假设：
H0： 不同状态与蛋白是否异常表达无关
Ha: 不同状态与蛋白是否异常表达有关
PA.alpha <- subset(PA, hospital == 1)
PA.alpha.fit <- glm(abnormal ~ status, data =PA.alpha, family = binomial())
summary(PA.alpha.fit)
# Coefficients:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -0.4595 0.3687 -1.246 0.21268
# status 1.2186 0.4502 2.707 0.00679 **看看这个P值吧，明显小于0.01，所以拒绝Ho， 看来这只是一部分的不同意见而已呀。
第五题
不知不觉天就黑了，一遍做题目，一边写内容讲解，时间过得真快。同志们，这题厉害了，有难度哦。请听题：
NASH是一种常见的肝病， 发现于肥胖和糖尿病有关。于是建立老鼠模型（可怜的小白鼠呀），研究adipose tissue和liver的关系。野生型雄性C57B1/6老鼠喂食LFD或者HFD，21周。然后根据liver histology和food fed分成4组。相信描述在GDS4013.txt, GDS4013_sample.txt和mehtylation.txt
A: 读取数据，看看有没有NA, 有的话，你如何处理
B： 找到HF>LF的上调基因，列出前20个，用啥检验呢？
C: 你有来自另一个研究的正常状态肝组织的每个基因的甲基化水平数据。你想知道差异表达基因的甲基化水平和其他正常状态的基因，请随机找DEG相同数量的数据，用t-test检验
D： 提供第三题的power.(温馨小提示：这是两个独立样本，当你计算绝对效应值的使用混池方差哦)
关键字： power, t-test, apply, is.na
第A题：就是看看没有NA值了。 翻到R语言实战第18章，处理缺失数据的高级方法的识别缺失值。对了默认把字符当做因子，对于那么多数据会降低读取速度，而且基因名也不是因子，所以设置为FALSE。
GSD4013 <- read.csv("GDS4013.csv",stringsAsFactors = F)
table(is.na(GSD4013))
# FALSE
# 508459
all(complete.cases(GSD4013))
TRUEis.na结果全都是false， comple.case全部是true，说明没有缺失值。那如果有缺失值怎么办？
直接删除这一行数据，缺失较少
多重插补，用附近的值预测缺失值，缺失多
第B题：首先就比较LF和HF两组之间的差异，然后没有告诉我们是否服从正态分布，所以我们用非参数方法，由于两个组之间的数量不同。所以选择wilcoxon秩和检验。
这里要用apply了，请大家仔细看我的代码哦！
基本思路是写一个自定义秩和检验函数，能够返回p值，这个p值要有名字。
my.wilcox.test <- function(x){
x <- t(x)
LF <- as.numeric(x[2:10])
HF <- as.numeric(x[11:19])
fit <- wilcox.test(HF, LF, alternative = "greater", exact = FALSE )
pv <- fit$p.value
return(pv)
}
pvalues <- apply(GSD4013, 1, my.wilcox.test)然后排序后，取出前20个。
sort(pvalues)[1:20]
# [1] 0.0002061474 0.0002061474 0.0002061474 0.0002061474 0.0002061474 0.0002868168
# [7] 0.0003961341 0.0003961341 0.0003961341 0.0003961341 0.0003961341 0.0005398953
# [13] 0.0005431233 0.0005431233 0.0005431233 0.0007392324 0.0007392324 0.0007392324
# [19] 0.0007392324 0.0007392324这里我有一个疑问，要不要做多重试验校正呢。抱着试一试的态度，我发现结果没有一个显著性，估计是不做了。
第c题：这一题的目标是看一下正常组织中里那些差异表达的基因的甲基化程度和非差异表达基因的甲基化程序是否有差异。
如果正常组织两则差异的话， 说明了什么问题呢？
题目要用t检验，那么我就随机选择100个基因吧。既然是显著性，不是极显著，那么就0.05吧。随机选择用sample.
names(pvalues) <- GSD4013$X
set.seed(1993)
DEG_name <- sample(names(pvalues[pvalues < 0.05]), 100)
non_DEG_name <- sample(names(pvalues[pvalues > 0.05]), 100)
methy_DEG <- methy[methy$V1 %in% DEG_name, ]
methy_non_DEG <- methy[methy$V1 %in% non_DEG_name, ]然后看这两个方差是否相同
var.test(methy_DEG$V2, methy_non_DEG$V2)
# F test to compare two variances
# data: methy_DEG$V2 and methy_non_DEG$V2
# F = 1.2854, num df = 99, denom df = 99, p-value = 0.2133
# alternative hypothesis: true ratio of variances is not equal to 1
# 95 percent confidence interval:
# 0.8649028 1.9104758
# sample estimates:
# ratio of variances
# 1.2854480.2大于0.5, 认为两者相同吧。然后用t检验。
t.test(methy_DEG$V2, methy_non_DEG$V2, var.equal = TRUE)
# Two Sample t-test
# data: methy_DEG$V2 and methy_non_DEG$V2
# t = 3.7556, df = 198, p-value = 0.0002273
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
# 0.05336879 0.17138190
# sample estimates:
# mean of x mean of y
# 0.4963295 0.3839542p小于0.01,极显著哦，看来的确有差别。这是不是说“不是甲基化导致差异吧。”
最后一题，我们算power，用到是pwr.t.test() ， 请翻到226页
#install.packages("pwr") # 请安装这个包哦
library(pwr)
vars <- var(c(methy_DEG$V2, methy_non_DEG$V2))
ds = (mean(methy_DEG$V2) - mean(methy_non_DEG$V2)) / vars # 根据混池计算效应值
powers <- pwr.t.test(n=200, d=ds, sig.level = .05,
type="two.sample", alternative = "two.sided")
powers
# Two-sample t test power calculation
# n = 100
# d = 2.355184
# sig.level = 0.05
# power = 1
# alternative = two.sided
#
# NOTE: n is number in *each* group题目终于做完了，欢迎讨论，不是正确答案，提供一种思路。还有阅读原文，可以下载我百度云的代码。