读书笔记(R语言)
作者:曾健明
公众号: 生信技能树
整理原因:在公众号这篇文章中看到如下一段话,自己最近又在练习用markdown格式写东西。
假如你接受了我的建议,把一本书看五遍,同一个领域的书籍至少看5本,那么希望你也记录一下读书笔记与我分享,欢迎来信交流,我的邮箱是 jmzeng1314@163.com
或者你把我的这个笔记重新编辑为markdown格式,也能发邮件给我获得认识我的资格哈!
R与ASReml-R统计分析教程(林元震)中国林业出版社
- 1-3章 R的基本语法
- 第4章 各种统计方法,
- 第5章 R的绘图,
- 第6章 ASReml-R这个包
语法重点:
install.packages(), library(), help(), example(), demo(), length(), attribute(), class(), mode(), dim(), names(), str(), head(), tail()
rep, seq, paste, array, matrix, data.frame, list, c(), factor()
(na.omit, na.rm = T)
缺失值处理,
as.numeric(), as.character(), as.factor(), as.logical()
##类型转换
as.numeric()非常有用,在画图的时候经常需要加上,因为数据在处理的过程中经常被搞错成了字符串格式;
as.logical()可以进行分类,只有0, NA, NAN, NULL是FALSE排序,合并,分割成子集,数据整合重构: reshape2, plyr包
可以先了解一些R语言自带的数据包(见附录1),然后试用一下aggregate函数,数据汇总,根据右边的因子来把左边的数据进行分割并处理一个函数控制语句,自编函数
统计分析
1. summary(), library(pastecs); options(digits = 2);stat.desc(), library(psych); describe()
2. 方差分析(analysis of variance,ANOVA)用来检验分组是否有显著差异
2.1 单因素+重复,数据框
df = data.frame(yield, treat)
fit = aov(yield ~ treat, data = df)
summary(fit) ##可以用summary来查看这次分析结果
TukeyHSD(fit) ##进行多重比较 方法1
duncan.test(fit, ”treat”, alpha = 0.05)##进行多重比较 方法2
2.2 双因素无重复,数据框
df=data.frame(yield,treat1,treat2)
fit=aov(yield~treat1+treat2,data=df)
## 这时候做多重比较就比较复杂了,
library(agricolae)
Duncan.test(fit,”treat1”,alpha=0.5)
Duncan.test(fit,”treat2”,alpha=0.5)
2.3 双因素+重复,数据库首先要进行处理,把treat1和treat2合并成group来区分重复
Df$group = sapply(df, function(x)paste(df$treat1, df$treat2, sep = ””)
fit = aov(yield ~ treat1 + treat2 + group, data = df)
2.4 多元方差与此类似,不停的增加因子来区分变量及group
3. 随机分组的检验
3.1 完全随机实验:
等同于方差分析的单因素+重复(判断不同的处理是否有差异)
3.2 单因素随机区组实验:
等同于方差分析的双因素无重复,其中(区组这个因素是人为控制的差异,不需要检验,主要检验我们的不同的处理是否有差异)
3.3 双因素随机区组实验:
不等同于方差分析的双因素+重复,但是与之类似,其中重复这个变量与之前的group变量有点区别,这里是我们的区组,而不是treat1和treat2的简单组合,所以我们需要分析treat1和treat2处理间的差异,但同时不需要考虑区组的差异
fit = aov(yield ~ treat1 * treat2 + block, data = df)
如果treat1有2个水平,treat2有3个水平,那么之前的group应该是6个,但是我们的block是区组的个数,还是3个,数据是18个。
3.4 三因素随机区组实验
看下面的例子:
图片
其中npk数据框里面有着N,P,K三个因素,每个因素都有两个水平,共8个group组合,分成了6个区组block,即为6个重复。但是每个group组合并没有包括所有的8个水平组合,只有4个而已,所以数据量也只有4个。这就是方差分析与随机区组分析最大的区别所在。
图片
4. 统计显著性检验(前提是符合各种概率分布模型)
4.1. t检验
i. 单样本,对一个多数据的向量x,看看是否是服从正态分布
qqnorm(x),qqline(x),正态QQ图,plot(density(x))核密度图 ,shapiro.test(x) 正态性检验,T.test(x,mu=8,alternative=”two.sided”)看看这个数据的均值与8的差异是否显著。
ii. 双样本检验是否显著差异,t.test(a,b)或者t.test(a~b)
4.2. 卡方检验
i. 是否符合一定的比例chisq.test(c(49,51),c(0.5,0.5)),看看扔100次硬币的正反面比例是否正常
ii. 图片
iii. P值为0.8415,所以显著的正常
a) 独立性检验
i. 2x2列联表或者2xc列联表独立性检验,主要是为了看某个处理是否改变了原来的标准比例,比如本来正常1:1的扔硬币,扔一百次是49:50的,但是现在换了一个小硬币,再扔一百次,结果是48:52,我们就想看看这个硬币是否改变了比例
图片
很明显可以看出比例未发生变化,同理可以扩展到RxC列联表的比例是否变化
5. 回归分析
5.1 简单线性回归,fit=lm(y~x)可以对此回归进行一系列分析,summary(fit),round(fitted(fit),2) 预测值,round(residuals(fit),2)残差值,abline(fit)回归线
5.2 多项式回归fit=lm(y~x+I(x^2))以此类推
5.3 多元回归fit=lm(yx1+x2+x3)以此类推,fit=lm(yx1+x2+x1: x2)有交互项。
5.4 回归诊断,对fit对象plot可以输出四幅图
par(mfrow=c(2,2))
plot(fit)
- 第一幅图是残差值与预测值的线性关系图,理论上残差值应该是随机分布在预测值的两端。
- 第二幅图是Q-Q图,判断残差值在标准正态分布下的概率,理论上应该是45度直线。
- 第三幅图是位置尺度图,判断同方差性,假设是方差不变,所以图中的点应该随机分布于水平线的两侧。
- 第四幅图是残差值的杠杆图,用来判断异常点,鉴别高杠杆点,离群点,强影响点,识别异常点。
5.5 广义线性模型
5.6 逻辑回归和泊松回归
6 概率分布
6.1 分布+概率密度函数d+累计分布函数p+随机抽样r+分布检验ks.test(x,”pnorm”)
6.2 正态分布(norm),指数分布(exp),二项分布,泊松分布,卡方分布(chisq),伽马分布(gama),贝塔分布(beta),T分布,F分布,均匀分布(unif),韦伯分布(weibull),一般连续分布,一般离散分布。
6.3 很复杂,见附录2
附录I : datasets(R自带数据包)
向量
euro #欧元汇率,长度为11,每个元素都有命名
landmasses #48个陆地的面积,每个都有命名
precip #长度为70的命名向量
rivers #北美141条河流长度
state.abb #美国50个州的双字母缩写
state.area #美国50个州的面积
state.name #美国50个州的全称
因子
state.division #美国50个州的分类,9个类别
state.region #美国50个州的地理分类
矩阵、数组
euro.cross #11种货币的汇率矩阵
freeny.x #每个季度影响收入四个因素的记录
state.x77 #美国50个州的八个指标
USPersonalExpenditure #5个年份在5个消费方向的数据
VADeaths #1940年弗吉尼亚州死亡率(每千人)
volcano #某火山区的地理信息(10米×10米的网格)
WorldPhones #8个区域在7个年份的电话总数
iris3 #3种鸢尾花形态数据
Titanic #泰坦尼克乘员统计
UCBAdmissions #伯克利分校1973年院系、录取和性别的频数
crimtab #3000个男性罪犯左手中指长度和身高关系
HairEyeColor #592人头发颜色、眼睛颜色和性别的频数
occupationalStatus #英国男性父子职业联系
类矩阵
eurodist #欧洲12个城市的距离矩阵,只有下三角部分
Harman23.cor #305个女孩八个形态指标的相关系数矩阵
Harman74.cor #145个儿童24个心理指标的相关系数矩阵
数据框
airquality #纽约1973年5-9月每日空气质量
anscombe #四组x-y数据,虽有相似的统计量,但实际数据差别较大
attenu #多个观测站对加利福尼亚23次地震的观测数据
attitude #30个部门在七个方面的调查结果,调查结果是同一部门35个职员赞成的百分比
beaver1 #一只海狸每10分钟的体温数据,共114条数据
beaver2 #另一只海狸每10分钟的体温数据,共100条数据
BOD #随水质的提高,生化反应对氧的需求(mg/l)随时间(天)的变化
cars #1920年代汽车速度对刹车距离的影响
chickwts #不同饮食种类对小鸡生长速度的影响
esoph #法国的一个食管癌病例对照研究
faithful #一个间歇泉的爆发时间和持续时间
Formaldehyde #两种方法测定甲醛浓度时分光光度计的读数
Freeny #每季度收入和其他四因素的记录
dating from #配对的病例对照数据,用于条件logistic回归
InsectSprays #使用不同杀虫剂时昆虫数目
iris #3种鸢尾花形态数据
LifeCycleSavings #50个国家的存款率
longley #强共线性的宏观经济数据
morley #光速测量试验数据
mtcars #32辆汽车在11个指标上的数据
OrchardSprays #使用拉丁方设计研究不同喷雾剂对蜜蜂的影响
PlantGrowth #三种处理方式对植物产量的影响
pressure #温度和气压
Puromycin #两种细胞中辅因子浓度对酶促反应的影响
quakes #1000次地震观测数据(震级>4)
randu #在VMS1.5中使用FORTRAN中的RANDU三个一组生成随机数字,共400组。
#该随机数字有问题。在VMS2.0以上版本已修复。
rock #48块石头的形态数据
sleep #两药物的催眠效果
stackloss #化工厂将氨转为硝酸的数据
swiss #瑞士生育率和社会经济指标
ToothGrowth #VC剂量和摄入方式对豚鼠牙齿的影响
trees #树木形态指标
USArrests #美国50个州的四个犯罪率指标
USJudgeRatings #43名律师的12个评价指标
warpbreaks #织布机异常数据
women #15名女性的身高和体重
列表
state.center #美国50个州中心的经度和纬度
类数据框
ChickWeight #饮食对鸡生长的影响
CO2 #耐寒植物CO2摄取的差异
DNase #若干次试验中,DNase浓度和光密度的关系
Indometh #某药物的药物动力学数据
Loblolly #火炬松的高度、年龄和种源
Orange #桔子树生长数据
Theoph #茶碱药动学数据
时间序列数据
airmiles #美国1937-1960年客运里程营收(实际售出机位乘以飞行哩数)
AirPassengers #Box & Jenkins航空公司1949-1960年每月国际航线乘客数
austres #澳大利亚1971-1994每季度人口数(以千为单位)
BJsales #有关销售的一个时间序列
BJsales.lead #前一指标的先行指标(leading indicator)
co2 #1959-1997年每月大气co2浓度(ppm)
discoveries #1860-1959年每年巨大发现或发明的个数
ldeaths #1974-1979年英国每月支气管炎、肺气肿和哮喘的死亡率
fdeaths #前述死亡率的女性部分
mdeaths #前述死亡率的男性部分
freeny.y #每季度收入
JohnsonJohnson #1960-1980年每季度Johnson & Johnson股票的红利
LakeHuron #1875-1972年某一湖泊水位的记录
lh #黄体生成素水平,10分钟测量一次
lynx #1821-1934年加拿大猞猁数据
nhtemp #1912-1971年每年平均温度
Nile #1871-1970尼罗河流量
nottem #1920-1939每月大气温度
presidents #1945-1974年每季度美国总统支持率
UKDriverDeaths #1969-1984年每月英国司机死亡或严重伤害的数目
sunspot.month #1749-1997每月太阳黑子数
sunspot.year #1700-1988每年太阳黑子数
sunspots #1749-1983每月太阳黑子数
treering #归一化的树木年轮数据
UKgas #1960-1986每月英国天然气消耗
USAccDeaths #1973-1978美国每月意外死亡人数
uspop #1790–1970美国每十年一次的人口总数(百万为单位)
WWWusage #每分钟网络连接数
Seatbelts #多变量时间序列。和UKDriverDeaths时间段相同,反映更多因素。
EuStockMarkets #多变量时间序列。欧洲股市四个主要指标的每个工作日记录,共1860条记录。
图片
图片
Warpbreaks这个数据集有3列变量,我们根据wool和tension这两个因子变量来分类对breaks这个数据变量求均值
图片
图片Airquality这个数据集有6个列变量,大气层,阳光,风,温度,月份,天数,虽然它们都是数据变量,但是我们可以把其中几个因子化来进行分类汇总计算,比如我们以month来作为因子,这样把数据分成了各个月份的,再对ozone和Temp进行分别求均值
图片
图片Chickwts这个数据有两列,不同的喂养环境下统计小鸡的重量,可以根据6中喂养环境来对各自的小鸡统计平均重量
图片
图片Esoph这个数据集有5个列变量,其中3个是因子,两个是数据,,同理做数据汇总
图片
图片
图片
这是一个时间序列数据,可以进行画图
还可以查看很多自己安装的包里面内置的数据
比如我安装一个ggplot2,它会自动下载几个相关的包一起安装
图片
data(package = "ggplot2")
可以查看这个包自带的数据集
图片
R还可以进行脚本运算,实习批量化处理数据
图片
附录二:各种统计分布函数
离散型
1. 二项分布Binomial distribution:binom
二项分布指的是N重伯努利实验,记为X ~ b(n,p),E(x)=np,Var(x)=np(1-p)
-
pbinom(q, size, prob)
, q是特定取值,比如pbinom(8, 20, 0.2)指第8次伯努利实验的累计概率。size指总的实验次数,prob指每次实验成功发生的概率 -
dbinom(x, size, prob)
, x同上面的q同含义。dfunction()对于离散分布来说结果是特定值的概率,对连续变量来说是密度(Density) -
rbinom(n, size, prob)
,产生n个b(size, prob)的二项分布随机数 -
qbinom(p, size, prob)
,quantile function 分位数函数。- 分位数:
若概率0<p<1,随机变量X或它的概率分布的分位数Za。是指满足条件p(X>Za)=α的实数。如t分布的分位数表,自由度f=20和α=0.05时的分位数为1.7247。 --这个定义指的是上侧α分位数 - α分位数:
实数α满足0 <α<1 时,α分位数是使P{X< xα}=F(xα)=α的数xα
双侧α分位数是使P{X<λ1}=F(λ1)=0.5α的数λ1、使 P{X>λ2}=1-F(λ2)=0.5α的数λ2。 - qbinom是上侧分位数,如qbinom(0.95,100,0.2)=27,指27之后P(x>=27)>=0.95。即对于b(100,0.2)为了达到0.95的概率至少需要27次重复实验。
- 分位数:
2. 负二项分布negative binomial distribution (帕斯卡分布)nbinom
掷骰子,掷到一即视为成功。则每次掷骰的成功率是1/6。要掷出三次一,所需的掷骰次数属于集合 { 3, 4, 5, 6, ... } 。掷到三次一的掷骰次数是负二项分布的随机变量。
- dnbinom(4, 3, 1/6) = 0.0334898,四次连续三次1的概率为这个数。
- 概率函数为f(k; r, p)=choose(k+r-1, r-1)p^r(1-p)^k, 当r=1时这个特例分布是几何分布
-
rnbinom(n, size, prob, mu)
其中n是需要产生的随机数个数,size是概率函数中的r,即连续成功的次数,prob是单词成功的概率,mu未知..(mu是希腊字母υ的读音)
3. 几何分布Geometric Distribution,geom
n次伯努利试验,前n-1次皆失败,第n次才成功的机率
-
dgeom(x, prob)
,注意这里的x取值是0:n,即dgeom(0,0.2)=0.2,以上的二项分布和负二项分布也是如此。 ngeom(n, prob)
4. 超几何分布Hypergeometric Distribution,hyper
它描述了由有限个(m+n)物件中抽出k个物件,成功抽出指定种类的物件的次数(不归还)。
- 概率:p(x) = choose(m, x) choose(n, k-x) / choose(m+n, k) for x = 0, ..., k.
- 当n=1时,这是一个0-1分布即伯努利分布,当n接近无穷大∞时,超几何分布可视为二项分布
-
rhyper(nn, m, n, k)
,nn是需要产生的随机数个数,m是白球数(计算目标是取到x个白球的概率),n是黑球数,k是抽取出的球个数 dhyper(x, m, n, k)
5. 泊松分布 Poisson Distribution,pois
p(x) = lambda^x exp(-lambda)/x!
for x = 0, 1, 2, .... The mean and variance are E(X) = Var(X) = λ. x ~ π(λ)
泊松分布的参数λ是单位时间(或单位面积)内随机事件的平均发生率.泊松分布适合于描述单位时间内随机事件发生的次数。如某一服务设施在一定时间内到达的人数,电话交换机接到呼叫的次数,汽车站台的候客人数,机器出现的故障数,自然灾害发生的次数等等.
rpois(n, lambda)
dpois(x,lambda)
连续型
6. 均匀分布 Uniform Distribution,unif
f(x) = 1/(max-min) for min <= x <= max
runif(n, min, max)
##生成16位数的随机数:
as.character(runif(1, 1000000000000000, 9999999999999999))
dunif(x, min, max) = 1 ## 恒定等于1/(max-min).
对于连续变量,dfunction的值是x去特定值代入概率密度函数得到的函数值。
7. 正态分布Normal Distribution,norm
f(x) = 1/(sqrt(2 pi) sigma) e^-((x - mu)^2/(2 sigma^2))
其中mu是均值,sigma是standard deviation标准差
理论上可以证明如果把许多小作用加起来看做一个变量,那么这个变量服从正态分布
rnorm(n, mean=0, sd=1) ##后两个参数如果不填则默认为0,1。
dnorm(x, mean, sd) ##sd是标准差。
画出正态分布概率密度函数的大致图形:
x <- seq(-3, 3, 0.1)
plot(x, dnorm(x)) ##plot中的x,y要有相关关系才会形成函数图。
qnorm(p, mean, sd) ##这个还是上侧分位数,如qnorm(0.05)=-1.644854,即x<=这个数的累计概率小于0.05
- 3sigma法则:对于正态分布的x,x取值在(mean-3sd,mean+3sd)几乎是在肯定的。
因为pnorm(3)-pnorm(-3)=0.9973002 - 用正太分布产生一个16位长的随机数字:
as.character(10^16*rnorm(1))
8. 伽玛分布Gamma Distribution,gamma
假设随机变量X为 等到第α件事发生所需之等候时间。
f(x)= 1/(s^a Gamma(a)) x^(a-1) e^-(x/s) for x >= 0, a > 0 and s > 0.
- Gamma分布中的参数α,称为形状参数(shape parameter),即上式中的s,β称为尺度参数(scale parameter)上式中的a
- E(x)=sa, Var(x)=sa^2. 当shape=1/2,scale=2时,这样的gamma分布是自由度为1的开方分
-
dgamma(x, shape, rate=1, scale = 1/rate)
, 请注意R在这里提供的rate是scale尺度参数的倒数,如果dgamma(0, 1, 2)
则表示dgamma(0, shape=1, rate=2)
,而非dgamma(0, shape=1, scale=2)
pgamma(q, shape, rate = 1, scale = 1/rate, lower.tail = TRUE, log.p = FALSE)
qgamma(p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE, log.p = FALSE)
rgamma(n, shape, rate = 1, scale = 1/rate)
9. 指数分布Exponential Distribution,exp
指数分布可以用来表示独立随机事件发生的时间间隔,比如旅客进机场的时间间隔、中文维基百科新条目出现的时间间隔等等。
记作X ~ Exponential(λ)
f(x) = lambda e^(- lambda x) for x >= 0
- 其中lambda λ > 0是分布的一个参数,常被称为率参数(rate parameter). E(x)=1/λ,Var(x)=1/λ^2
dexp(x, rate = 1, log = FALSE)
pexp(q, rate = 1, lower.tail = TRUE, log.p = FALSE)
qexp(p, rate = 1, lower.tail = TRUE, log.p = FALSE)
rexp(n, rate = 1)
假设在公交站台等公交车平均10分钟有一趟车,那么每小时候有6趟车,即每小时出现车的次数~ Exponential(1/6)
- 我们可以产生10个这些随机数看看rexp(10,1/6)
- 60/(rexp10,1/6)即为我们在站台等车的随机时间,如下:
[1] 6.443148 24.337131 6.477096 2.824638 15.184945 14.594903
[7] 7.133842 8.222400 42.609784 15.182827
- 可以看见竟然有一个42.6分钟的随机数出现,据说这种情况下你可以投诉上海的公交公司。
不过x符合指数分布,1/x还符合指数分布吗? - pexp(6,1/6)=0.6321206, 也就是说这种情况下只有37%的可能公交车会10分钟以内来。
按照以上分析一个小时出现的公交车次数应该不符合指数分布。
10. 卡方分布(non-central)Chi-Squared Distribution,chisq
它广泛的运用于检测数学模型是否适合所得的数据,以及数据间的相关性。数据并不需要呈正态分布
k个标准正态变量的平方和即为自由度为k的卡方分布。
E(x)=k,Var(x)=2k.
dchisq(x, df, ncp=0, log = FALSE)
pchisq(q, df, ncp=0, lower.tail = TRUE, log.p = FALSE)
qchisq(p, df, ncp=0, lower.tail = TRUE, log.p = FALSE)
rchisq(n, df, ncp=0)
- 其中df为degrees of freedom。
- ncp是non-centrality parameter (non-negative); ncp=0时是central卡方分布,ncp不为0时,表示这个卡方分布是由非标准正态分布组合而成,ncp=这些正态 分布的均值的平方和。
11. β分布Beta Distribution,beta
变量x仅能出现于0到1之间。
- 空气中含有的气体状态的水分。表示这种水分的一种办法就是相对湿度。即现在的含水量与空气的最大含水量(饱和含水量)的比值。我们听到的天气预告用语中就经常使用相对湿度这个名词。
- 相对湿度的值显然仅能出现于0到1之间(经常用百分比表示)。冬季塔里木盆地的日最大相对湿度和夏季日最小相对湿度。证实它们都符合贝塔分布
dbeta(x, shape1, shape2, ncp = 0, log = FALSE)
pbeta(q, shape1, shape2, ncp = 0, lower.tail = TRUE, log.p = FALSE)
qbeta(p, shape1, shape2, ncp = 0, lower.tail = TRUE, log.p = FALSE)
rbeta(n, shape1, shape2, ncp = 0)
- shape1,shape2是beta分布的两个参数。
- E(x)=s1/(s1+s2),var(x)=s1*s2/(s1+s2)^2 * (s1+s2+1)
12. t分布Student t Distribution,t
应用在当对呈正态分布的母群体的均值进行估计。当母群体的标准差是未知的但却又需要估计时,我们可以运用学生t 分布。
- 学生t 分布可简称为t 分布。其推导由威廉·戈塞于1908年首先发表,当时他还在都柏林的健力士酿酒厂工作。因为不能以他本人的名义发表,所以论文使用了学生 (Student)这一笔名。之后t 检验以及相关理论经由罗纳德·费雪的工作发扬光大,而正是他将此分布称为学生分布。
dt(x, df, ncp, log = FALSE)
pt(q, df, ncp, lower.tail = TRUE, log.p = FALSE)
qt(p, df, ncp, lower.tail = TRUE, log.p = FALSE)
rt(n, df, ncp)
- df是自由度
- ncp是non-centrality parameter delta,If omitted, use the central t distribution。ncp出现时表示分布由非标准的卡方分布构成。
13. F分布
一个F-分布的随机变量是两个卡方分布变量的比率。F-分布被广泛应用于似然比率检验,特别是方差分析中
df(x, df1, df2, ncp, log = FALSE)
pf(q, df1, df2, ncp, lower.tail = TRUE, log.p = FALSE)
qf(p, df1, df2, ncp, lower.tail = TRUE, log.p = FALSE)
rf(n, df1, df2, ncp)
- df1,df2是两个自由度,ncp同t分布中的ncp。