|
R与ASReml-R统计分析教程（林元震）中国林业出版社 1-3章简单介绍了R的基本语法，然后第4章着重讲了各种统计方法，第5章讲R的绘图，最后一张讲ASReml-R这个包 1,install.packages(),library(),help(),example(),demo(),length(),attribute(),class(),mode(),dim(),names(),str(),head(), 2,rep,seq,paste,array,matrix,data.frame,list,c(),factor(), 3,缺失值处理（na.omit,na.rm=T），类型转换（as.numeric(),as.character(),as.factor(),as.logical()） 其中as.numeric()非常有用，在画图的时候经常需要加上，因为数据在处理的过程中经常被搞错成了字符串格式，而as.logical()可以进行分类，只有0，NA,NAN，NULL是FALSE 4,排序，合并，分割成子集，数据整合重构(reshape2和plyr包) 可以先了解一些R语言自带的数据包(见附录1)，然后试用一下aggregate函数，数据汇总，根据右边的因子来把左边的数据进行分割并处理一个函数 5,控制语句，自编函数 一、summary(), library(pastecs);options(digits=2);stat.desc(), library(psych);describe() 二、方差分析（analysis of variance,ANOVA）用来检验分组是否有显著差异 1，单因素+重复，数据框，df=data.frame(yield,treat) fit=aov(yield~treat,data=df) 可以用summary来查看这次分析结果summary(fit) 用TukeyHSD(fit)进行多重比较，或者duncan.test(fit,”treat”,alpha=0.05) 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) 3,双因素+重复，数据库首先要进行处理，把treat1和treat2合并成group来区分重复 Df$group=sapply(df,function(x) paste(df$treat1,df$treat2,sep=””)) 然后fit=aov(yield~treat1+treat2+group,data=df) 4，多元方差与此类似，不停的增加因子来区分变量及group 三、随机分组的检验 1，完全随机实验，等同于方差分析的单因素+重复（判断不同的处理是否有差异） 2，单因素随机区组实验，等同于方差分析的双因素无重复，其中（区组这个因素是人为控制的差异，不需要检验，主要检验我们的不同的处理是否有差异） 3，双因素随机区组实验，不等同于方差分析的双因素+重复，但是与之类似，其中重复这个变量与之前的group变量有点区别，这里是我们的区组，而不是treat1和treat2的简单组合，所以我们需要分析treat1和treat2处理间的差异，但同时不需要考虑区组的差异 fit=aov(yield~treat1*treat2+block,data=df)，如果treat1有2个水平，treat2有3个水平，那么之前的group应该是6个，但是我们的block是区组的个数，还是3个，数据是18个。 4,三因素随机区组实验 看下面的例子 本地图片，请重新上传 其中npk数据框里面有着N,P,K三个因素，每个因素都有两个水平，共8个group组合，分成了6个区组block，即为6个重复。但是每个group组合并没有包括所有的8个水平组合，只有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) 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列联表的比例是否变化 五、回归分析 1、简单线性回归，fit=lm(y~x)可以对此回归进行一系列分析，summary(fit)，round(fitted(fit),2) 预测值，round(residuals(fit),2)残差值，abline(fit)回归线 2、多项式回归fit=lm(y~x+I(x^2))以此类推 3、多元回归fit=lm(y~x1+x2+x3)以此类推，fit=lm(y~x1+x2+x1: x2)有交互项。 4、回归诊断，对fit对象plot可以输出四幅图par(mfrow=c(2,2)); plot(fit) 第一幅图是残差值与预测值的线性关系图，理论上残差值应该是随机分布在预测值的两端。 第二幅图是Q-Q图，判断残差值在标准正态分布下的概率，理论上应该是45度直线。 第三幅图是位置尺度图，判断同方差性，假设是方差不变，所以图中的点应该随机分布于水平线的两侧。 第四幅图是残差值的杠杆图，用来判断异常点，鉴别高杠杆点，离群点，强影响点，识别异常点。 5、广义线性模型 6、逻辑回归和泊松回归 六、概率分布 a) 分布+概率密度函数d+累计分布函数p+随机抽样r+分布检验ks.test(x,”pnorm”) b) 正态分布（norm），指数分布(exp)，二项分布，泊松分布，卡方分布(chisq)，伽马分布(gama)，贝塔分布(beta)，T分布，F分布，均匀分布（unif），韦伯分布(weibull),一般连续分布，一般离散分布。 c) 很复杂，见附录2 附录I，datasets（R自带数据包） 向量 因子 矩阵、数组 类矩阵 数据框 列表 类数据框 时间序列数据 本地图片，请重新上传 本地图片，请重新上传 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 http://zh./w/index.php?title=伽玛分布&variant=zh-cn 假设随机变量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)=s*a, Var(x)=s*a^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。
|
|
|
来自： 新用户3677sdB0 > 《基础》
