R语言中级篇

SectionII基础统计与数学模型数据的产生dat1-rnorm(20)#####一组数的平均值dat1dat2-rnorm(10,mean=7,sd=5)dat21基础统计数据检验正态W检验方法shapiro.test()，提供了W的统计量和检验P值，我们当P0.05时候，我们可以认为这组数据是正态的，但是适用条件是所测的数据个数在3and5000之间格式test.nor-dat$biomassshapiro.test(test.nor)格式binom.test(x,n,p=0.5,alternative=c(“two.sided”,“less”,“greater”),conf.level=0.95)####数据是否是二项分布的binom.test(c(682,243),p=3/4)binom.test(682,682+243,p=3/4)#Thesame1基础统计求一组数据的95%置信区间xbar-mean(x)#####一组数的平均值sigma-sd(x)####标准差sem-sigma/sqrt(length(x))####平均标准误low-xbar+sem*qnorm(0.025)###置信区间下限high-xbar+sem*qnorm(0.975)####置信区间上限1基础统计随机抽样（randomsampling）sample(x,size,replace=FALSE)a-1:40asample(a,6,replace=FALSE)sample(c(H,T),10,replace=T)[1]THTHTHHHHT1基础统计概率统计（Probability）a-1:40asample(a,6)40个数字中，抽取想要的6个数字，由于是非放回抽样，每个抽取数字的概率分别为:1/40,1/39,1/38,1/37,1/36,1/351基础统计整体的概率：Pro-1/prod(40:35)如果是仅仅选取5个数字，那么就有不同组合，总体的概率为prod(5:1)/prod(40:35)或者是：1/choose(40,5)数据分布-密度图（Densities）x-seq(-4,4,0.1)plot(x,dnorm(x),type=l)1基础统计数据分布-密度图（Densities）x-0:50plot(x,dbinom(x,size=50,prob=.33),type=h)###prob解释##probabilityofsuccessoneachtrial1基础统计基础统计练习x-rnorm(20)#####产生一组随机树Exercise:##编辑一个能够返回所有基本统计量的函数var(x)mean(x)sd(x)1基础统计R语言中求导###############R语言求导y=D(expression(sin(x)+x^3),'x')x=2;eval(y);###就可以计算出f(2);值y=deriv(expression(sin(x)+x^3),'x',function.arg=T)y(2)1基础统计Apply函数的应用函数格式：apply(X,MARGIN,FUN,...)###MARGIN可以选择，1为按行，2为按列dat-matrix(1:6,2,3)#######创建一个2×3的矩阵apply(dat,1,mean)apply(dat,2,mean)apply(dat,1,sd)Ex:a-tapply(dat$biomass,dat$disturb,mean)1基础统计2、一般线性模型（GeneralLinearmodels）2.1一般线性模型（GeneralLinearmodels）一般线性回归模型generallinearmodels##Note1:适用的假设前提是1)变量的独立性，非独立的变量会有很多意想不到的结果出现；2）正态的残差分布；3）方差齐性；4）小的采样误差##Note2；线性模型lm:自变量是连续变量（即x是1,2,3,4,4.5这种）；ANOVA:变量是因子型的（比如加温、森林干扰类型等）；ANCOVA:既有因子变量又有连续变量2.1一般线性模型（GeneralLinearmodels）线性回归模型lm()lm(formula,data,subset,weights,na.action,method=qr,model=TRUE,x=FALSE,y=FALSE,qr=TRUE,singular.ok=TRUE,contrasts=NULL,offset,...)Ex1:dat-read.csv(data_eg1.csv,head=T)head(dat)lm.model-lm(biomass~num.sp,data=dat)summary(lm.model)###利用lm.model来看结果分析表2.1一般线性模型（GeneralLinearmodels）线性模型显示结果解读Call:lm(formula=biomass~num.sp,data=dat)Residuals:Min1QMedian3QMax-26.8425-8.8579-0.60877.897924.7994Coefficients:EstimateStd.ErrortvaluePr(|t|)(Intercept)1.834e+026.954e-01263.782e-16***num.sp6.356e-025.011e-0312.682e-16***---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1Residualstandarderror:10.19on598degreesoffreedomMultipleR-squared:0.212,AdjustedR-squared:0.2107F-statistic:160.9on1and598DF,p-value:2.2e-16模型公式模型残差分布分位数模型系数显著水平残差标准误模型的拟合解释能力模型整体评价2.1一般线性模型（GeneralLinearmodels）线性回归模型作图1）画出散点图plot(dat$num.sp,dat$biomass,main=Maintitle,xlab=speciesrichness,ylab=Biomass,col='blue',lwd=0.5)abline(lm.model,col='red',lty=2,lwd=2)添加由lm.model确定的回归线2）画出回归线对于回归线，我们也可以加指定的截距和斜率的线abline(a=a,b=slope)这个就可以调用summary里面的参数3）画出回归线放入这个回归模型的公式与一些参数r2=0.212text(x=400,170,y=1.834+6.356x)text(500,170,as.expression(substitute(r^2==rr,list(rr=round(r2,2)))))####添加公式和P值,round去小数点位数###把公式放上去的例子plot(1:10,1:10)text(4,9,expression(hat(beta)==(X^t*X)^{-1}*X^t*y))text(4,8.4,expression(hat(beta)==(X^t*X)^{-1}*X^t*y),cex=.75)2.1一般线性模型（GeneralLinearmodels）###这种结果的显示也可以是方差分析表的形式summary.aov(lm.model)##结果显示DfSumSqMeanSqFvaluePr(F)num.sp11669016689.5160.872.2e-16***Residuals59862039103.7---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.12.1一般线性模型（GeneralLinearmodels）对于要执行很多个模型，而且你还期望能够自动提取这些模型的参数，特别是在大量的随机模拟中，那么这些参数是以list的形式存储在summary(lm.model)中的，可以调用和存储的，如：res.summary-summary(lm.model)mode(res.summary)a-res.summary[[1]]aa-res.summary[[4]]####输出结果显示EstimateStd.ErrortvaluePr(|t|)(Intercept)183.424267370.695375275263.777380.000000e+00num.sp0.063560160.00501124312.683518.146262e-33r2-res.summary[[8]]###调取R2###或者是ls(res.summary)r2-res.summary$r.squaredr22.1一般线性模型（GeneralLinearmodels）一个模型的结果里面存储了不止这些显示的参数，还有很多的结果可以调用：其他的一些专门的看的函数，比如coef,effects,residuals,fitted模型预测predict()模型预测函数predict,主要传递的参数名字和模型一致，比如我们有模型lm.model和一列新的物种数new.sp，要预测基于这个物种数的生物量new.sp-c(90,45,26,79,357)new.spnew.bio-predict(lm.model,data.frame(num.sp=new.sp))new.bio2.1一般线性模型（GeneralLinearmodels）一个模型的构建完毕后，我们用什么方法来评价模型的好坏呢？有很多的评判标准，其中比如高的拟合R平方，相对较少的因子，低的AIC值（相对于同一的应变量Y值）fitvalue-fitted(lm.model)res-resid(lm.model)####获取残差。（残差即：模型的预测值与实际值之间的差值，即余下模型没有解释的部分）plot(fitvalue,res)###可看下残差的和预测值关系，来评判其模型拟合效果，较好的残差分布abline(h=0,col=red)####abline(h=0)，好的模型，残差的散点图是符合正态分布的###标准化的残差图fit2-rstandard(lm.model)plot(fit2,res)或者可以看自变量和x的残差图2.1一般线性模型（GeneralLinearmodels）###模型的残差正态检测shapiro.test(lm.model$resid)qqnorm(lm.model$resid)##如果一个模型不是很好，那么我们就可能会用glm模型来处理AIC(lm.model)###但是注意，AIC是只针对一组相同的y值才有比较意义。AIC：AkaikeInformationCriterion(1973)，后面会详细讲[1]4491.881####如果模型不够好，那么我们就准备对其进行转换，如log,sqrt,2.1一般线性模型（GeneralLinearmodels）较好的模型残差分布2.1一般线性模型（GeneralLinearmodels）残差标准残差2.2方差分析AnalysisofVariance（ANOVA）方差分析AnalysisofVariance（ANOVA）单因素方差分析dat$Increase_T-as.factor(dat$Increase_T)###注意，我们做的是水平处理的话，那么应该隐去相关的数值信息，做因子转换dis.aov-aov(biomass~Increase_T,data=dat)summary(dis.aov)##结果显示DfSumSqMeanSqFvaluePr(F)Increase_T2507.07253.5422.9032.728e-06***Residuals24265.6811.07---Signif.cod