R语言方法总结

计算描述性统计量 :1、summary():例： summary(mtcars[vars])summary()函数提供了最小值、最大值、四分位数和数值型变量的均值，以及因子向量和逻辑型向量的频数统计。2、apply()函数或

sapply()函数计算所选择的任意描述性统计量。 mean、sd、和quantile。函数fivenum()可返回图基五数总括（下四分位数、中位数、上四分位数和最大值）。

var、min、max、median、length、rangeTukey’sfive-numbersummary，即最小值、sapply()例： mystats<-function(x,na.omit=FALSE){if(na.omit)x<-x[!is.na(x)]m<-mean(x)n<-length(x)s<-sd(x)skew<-sum((x-m)^3/s^3)/nkurt<-sum((x-m)^4/s^4)/n-3return(c(n=n,mean=m,stdev=s,skew=skew,kurtosis=kurt))}sapply(mtcars[vars],mystats)3、describe()：Hmisc包：返回变量和观测的数量、缺失值和唯一值的数目、平均值、分位数，以及五个最大的值和五个最小的值。例： library(Hmisc)describe(mtcars[vars])4、stat.desc()：pastecs包若basic=TRUE（默认值），则计算其中所有值、空值、缺失值的数量，以及最小值、最值域，还有总和。若desc=TRUE（同样也是默认值），则计算中位数、平均数、平均数的标准误、平均数置信度为95%的置信区间、方差、标准差以及变异系数。若norm=TRUE（不是默认的），则返回正态分布统计量，包括偏度和峰度（以及它们的统计显着程度）和 Shapiro–Wilk正态检验结果。这里使用了 p值来计算平均数的置信区间（默认置信度为例： library(pastecs)

大值、0.95：stat.desc(mtcars[vars])5、describe()：psych包计算非缺失值的数量、平均数、标准差、中位数、截尾均值、绝对中位差、最小值、最大值、值域、偏度、峰度和平均值的标准误例： library(psych)describe(mtcars[vars])分组计算描述性统计量1、aggregate()：例：aggregate(mtcars[vars],by=list(am=mtcars$am),mean)2、by()：例： dstats<-function(x)(c(mean=mean(x),sd=sd(x)))by(mtcars[vars],mtcars$am,dstats)by(mtcars[,vars],mtcars$am,plyr::colwis(dstats))3、summaryBy()：doBy包例library(doBy)summaryBy(mpg+hp+wt~am,data=mtcars,FUN=mystats)4、describe.by()：doBy包（describe.by()函数不允许指定任意函数，）例：library(psych)describe.by(mtcars[vars],mtcars$am)5、reshape包分组：（重铸和融合）例：library(reshape)dstats<-function(x)(c(n=length(x),mean=mean(x),sd=sd(x)))dfm<-melt(mtcars,measure.vars=c("mpg","hp","wt"),id.vars=c("am","cyl"))cast(dfm,am+cyl+variable~.,dstats)频数表和列联表1、table()：生成简单的频数统计表mytable<-with(Arthritis,table(Improved))Mytable2、prop.table()：频数转化为比例值prop.table(mytable)3、prop.table()*100：转化为百分比prop.table(mytable)*100二维列联表4、table(A,B)/xtabs(~A+b,data=mydata)例：mytable<-xtabs(~Treatment+Improved,data=Arthritis)5、margin.table()和prop.table()：函数分别生成边际频数和比例(1:行，2：列)行和与行比例margin.table(mytable,1)prop.table(mytable,1)列和与列比例margin.table(mytable,2)prop.table(mytable,2)prop.table(mytable)6、addmargins()：函数为这些表格添加边际和addmargins(mytable)admargins(prop.table(mytable))addmargins(prop.table(mytable,1),2)addmargins(prop.table(mytable,2,1)7.crossTable()：gmodels包例：library(gmodels)CrossTable(Arthritis$Treatment,Arthritis$Improved)多维列联表1、table()和xtabs()：都可以基于三个或更多的类别型变量生成多维列联表。2、ftable():例：mytable<-xtabs(~Treatment+Sex+Improved,data=Arthritis)mytableftable(mytable)margin.table(mytable,1)margin.table(mytable,2)margin.table(mytable,3)margin.table(mytable,c(1,3))ftable(prop.table(mytable,c(1,2)))ftable(addmargins(prop.table(mytable,c(1,2)),3))gtable(addmargins(prop.table(mytable,c(1,2)),3))*100独立检验1、卡方独立性检验 ：chisq.test()例：library(vcd)mytable<-xtabs(~Treatment+Improved,data=Arthritis)chisq.test(mytable)mytable<-xtabs(~Improved+Sex,data=Arthritis)chisq.test(mytable)2、Fisher精确检验：fisher.test()例：mytable<-xtabs(~Treatment+Improved,data=Arthritis)fisher.test(mytable)3、Cochran-Mantel—Haenszel检验：mantelhaen.test()例：mytable<-xtabs(~Treatment+Improved+Sex,data=Arthritis)mantelhaen.test(mytable)相关性度量1、assocstats()：例：library(vcd)mytable<-xtabs(~Treatment+Improved,data=Arthritis)assocstats(mytable)2、cor()：函数可以计算这三种相关系数，3、cov()：函数可用来计算协方差例：states<-state.x77[,1:6]cov(states)cor(states)cor(states,method="spearman")x<-states[,c("Population","Income","Illiteracy","HSGrad")]y<-states[,c("LifeExp","Murder")]cor(x,y)4、pcor()：偏相关 ggm包例：library(ggm)pcor(c(1,5,2,3,6),cov(states))相关性的显着性检验1、cor.test()其中的 x 和y为要检验相关性的变量， alternative则用来指定进行双侧检验或单侧检验（取值为"two.side"、"less"或"greater"），而method用以指定要计算的相关类型（ "pearson"、"kendall"或"spearman"）当研究的假设为总体的相关系数小于 0时，请使用 alternative="less"。在研究的假设为总体的相关系数大于 0时，应使用 alternative="greater"。在默认情况下，假设为alternative="two.side"（总体相关系数不等于 0）。例：cor.test(states[,3],states[,5])2、corr.test()：可以为Pearson、Spearman或Kendall相关计算相关矩阵和显着性水平。例：library(psych)corr.test(states,use="complete")3、pcor.test()：psych包检验1、t.test(y~x,data)（独立样本）例：library(MASS)t.test(Prob~So,data=UScrime)2、t.test(y1,y2,paired=TRUE)(非独立)例：library(MASS)sapply(UScrime[c("U1","U2")],function(x)(c(mean=mean(x),sd=sd(x))))with(UScrime,t.test(U1,U2,paired=TRUE))组间差异的非参数检验两组的比较：1、wilcox.test(y~x,data)：评估观测是否是从相同的概率分布中抽得例：with(UScrime,by(Prob,So,median))wilcox.test(Prob~So,data=UScrime)2、wilcox.test(y1，y2,paried=TRUE)：它适用于两组成对数据和无法保证正态性假设的情境。例：sapply(UScrime[c("U1","U2")],median)with(UScrime,wilcox.test(U1,U2,paired=TRUE))多于两组的比较：1、kruskal.test(y~A，data):各组独立kruskal.test(Illiteracy~state.region,data=states)2、friedman.test(y~A|B,data):各组不独立非参数多组比较 :1、npmc() :npmc包例：class<-state.regionvar<-state.x77[,c("Illiteracy")]rm(class,var)library(npmc)summary(npmc(mydata),type="BF")aggregate(mydata,by=list(mydata$class),median)回归用一个或多个预测变量（也称自变量或解释变量）来预测响应变量（也称因变量、效标变量或结果变量）的方法。1、lm():拟合回归模型 lm(y~x1+x2+x3,data)简单线性回归1、lm(): （data是数据框）例：fit<-lm(weight~height,data=women)summary(fit)women$weightfitted(fit)residuals(fit)plot(women$height,women$weight,main="WomenAge30-39",xlab="Height(ininches)",ylab="Weight(inpounds)")多项式回归例：fit2<-lm(weight~height+I(height^2),data=women)summary(fit2)plot(women$height,women$weight,main="WomenAge30-39",xlab="Height(ininches)",ylab="Weight(inlbs)")lines(women$height,fitted(fit2))2、scatterplot()：绘制二元关系图例：library(car)scatterplot(weight~height,data=women,spread=FALSE,lty.smooth=2,pch=19,main="WomenAge30-39",xlab="Height(inches)",ylab="Weight(lbs.)")多元线性回归1、scatterplotMatrix()scatterplotMatrix()

：car包函数默认在非对角线区域绘制变量间的散点图，

并添加平滑

（loess）和线性拟合曲线。对角线区域绘制每个变量的密度图和轴须图。例：

fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)有交互项的多元线性回归例：fit<-lm(mpg~hp+wt+hp:wt,data=mtcars)summary(fit)1、effect():effects包 ：展示交互项的结果term即模型要画的项， mod为通过lm()拟合的模型，值， multiline=TRUE 选项表示添加相应直线。

xlevels

是一个列表，指定变量要设定的常量例：library(effects)plot(effect("hp:wt",fit,xlevels=list(wt=c(2.2,3.2,4.2))),multiline=TRUE)回归诊断1、confint()：求模型参数的置信区间例：fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)confint(fit)2、plot()：生成评价模型拟合情况的图形例：

fit<-lm(weight~height,data=women)par(mfrow=c(2,2))plot(fit)3、lm():删除观测点例：newfit<-lm(weight~height+I(height^2),data=women[-c(13,15),])par(mfrow=c(2,2))plot(newfit)par(opar)gvlma

包提供了对所有线性模型假设进行检验的方法检验正态性：4、qqPlot()：car

包：学生化残差（

studentizedresidual，也称学生化删除残差或折叠化残差）例：library(car)fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)qqPlot(fit,labels=s(states),id.method="identify",simulate=TRUE,main=注：id.method="identify"选项能够交互式绘图5、fitted():提取模型的拟合值例：fitted(fit)[“Nevada”]6、residuals():二项式回归模型的残差例：residuals(fit)[“Nevada”]7、residplot()：生成学生化残差柱状图（即直方图），并添加正态曲线、核密度曲线和轴须图。

"Q-QPlot")它不需要加载car包例：residplot<-function(fit,nbreaks=10){z<-rstudent(fit)hist(z,breaks=nbreaks,freq=FALSE,xlab="StudentizedResidual",main="DistributionofErrors")rug(jitter(z),col="brown")curve(dnorm(x,mean=mean(z),sd=sd(z)),add=TRUE,col="blue",lwd=2)lines(density(z)$x,density(z)$y,col="red",lwd=2,lty=2)legend("topright",legend=c("NormalCurve","KernelDensityCurve"),lty=1:2,col=c("blue","red"),cex=.7)}residplot(fit)误差的独立性8、durbinWatsonTest()

：验证独立性例：durbinWatsonTest(fit)验证线性9、crPlots()：car包成分残差图也称偏残差图例：crPlots(fit)同方差性 （car包的两个函数）10、ncvTest() ：生成一个计分检验，零假设为误差方差不变，备择假设为误差方差随着拟合值水平的变化而变化。若检验显着，则说明存在异方差性11、spreadLevelPlot()：添加了最佳拟合曲线的散点图，展示标准化残差绝对值与拟合值的关系。例：library(car)ncvTest(fit)spreadLevelPlot(fit)线性模型假设的综合验证1、gvlma()：gvlma包：线性模型假设进行综合验证，同时还能做偏斜度、峰度和异方差性的评价例：library(gvlma)gvmodel<-gvlma(fit)summary(gvmodel)多重共线性1、vif() ：car包：函数提供 VIF值， vif>2就表明存在多重共线性问题例：vif(fit)sqrt(vif(fit))>2异常观测值1、outlierTest()

：car

包：求得最大标准化残差绝对值

Bonferroni

调整后的

p值例：library(car)outlierTest(fit)高杠杆值点1、hat.plot():观测点的帽子值大于帽子均值的 2或3倍，即可以认定为高杠杆值点例：hat.plot<-function(fit){p<-length(coefficients(fit))n<-length(fitted(fit))plot(hatvalues(fit),main="IndexPlotofHatValues")abline(h=c(2,3)*p/n,col="red",lty=2)identify(1:n,hatvalues(fit),names(hatvalues(fit)))}hat.plot(fit)强影响点 ：Cook’sD值大于预测变量数目。

4/(n-k-1)，则表明它是强影响点，其中

n为样本量大小，

k是例：cutoff<-4/(nrow(states)-length(fit$coefficients)-2)plot(fit,which=4,cook.levels=cutoff)abline(h=cutoff,lty=2,col="red")1、influencePlot()：car包：离群点、杠杆值和强影响点的信息整合到一幅图形中例：influencePlot(fit,id.method="identify",main="InfluencePlot",sub="CirclesizeisproportialtoCook'sDistance")纵坐标超过+2或小于?2的州可被认为是离群点，水平轴超过 0.2或0.3的州有高杠杆值（通常为预测值的组合）。圆圈大小与影响成比例，圆圈很大的点可能是对模型参数的估计造成的不成比例影响的强影响点变量变换1、powerTransform():car包:函数通过 λ的最大似然估计来正态化变量 x。例：library(car)summary(powerTransform(states$Murder))2、boxTidwell():car包：通过获得预测变量幂数的最大似然估计来改善线性关系例：library(car)boxTidwell(Murder~Population+Illiteracy,data=states)模型比较1、anova()：基础包：比较两个嵌套模型的拟合优度例：fit1<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)fit2<-lm(Murder~Population+Illiteracy,data=states)anova(fit2,fit1)2、AIC()：AIC值越小的模型 (可以不嵌套)要优先选择，它说明模型用较少的参数获得了足够的拟合度。例：fit1<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)fit2<-lm(Murder~Population+Illiteracy,data=states)AIC(fit1,fit2)变量选择1、stepAIC()：MASS包：逐步回归模型例：library(MASS)fit1<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)stepAIC(fit,direction="backward")2、regsubsets()：leaps包：全子集回归例：library(leaps)leaps<-regsubsets(Murder~Population+Illiteracy+Income+Frost,data=states,nbest=4)plot(leaps,scale="adjr2")交叉验证1、crossval()函数：bootstrap包：实现k重交叉验证例：shrinkage<-function(fit,k=10){require(bootstrap)#definefunctionstheta.fit<-function(x,y){lsfit(x,y)}theta.predict<-function(fit,x){cbind(1,x)%*%fit$coef}#matrixofpredictorsx<-fit$model[,2:ncol(fit$model)]vectorofpredictedvaluesy<-fit$model[,1]results<-crossval(x,y,theta.fit,theta.predict,ngroup=k)r2<-cor(y,fit$fitted.values)^2r2cv<-cor(y,results$cv.fit)^2cat("OriginalR-square=",r2,"\n")cat(k,"FoldCross-ValidatedR-square=",r2cv,"\n")cat("Change=",r2-r2cv,"\n")}2、shrinkage()：交叉验证 ;R平方减少得越少，预测则越精确。例：fit<-lm(Murder~Population+Income+Illiteracy+Frost,data=states)shrinkage(fit)相对重要性1、scale()：将数据标准化为均值为 0、标准差为 1的数据集，这样用注意， scale()函数返回的是一个矩阵，而 lm()函数要求一个数据框

R回归即可获得标准化的回归系数。zfit<-lm(Murder~Population+Income+Illiteracy+Frost,data=zstates)coef(zfit)2、relweights()：相对权重例：relweights<-function(fit,...){R<-cor(fit$model)nvar<-ncol(R)rxx<-R[2:nvar,2:nvar]rxy<-R[2:nvar,1]svd<-eigen(rxx)evec<-svd$vectorsev<-svd$valuesdelta<-diag(sqrt(ev))correlationsbetweenoriginalpredictorsandneworthogonalvariableslambda<-evec%*%delta%*%t(evec)lambdasq<-lambda^2regressioncoefficientsofYonorthogonalvariablesbeta<-solve(lambda)%*%rxyrsquare<-colSums(beta^2)rawwgt<-lambdasq%*%beta^2import<-(rawwgt/rsquare)*100lbls<-names(fit$model[2:nvar])rownames(import)<-lblscolnames(import)<-"Weights"#plotresultsbarplot(t(import),names.arg=lbls,ylab="%ofR-Square",xlab="PredictorVariables",main="RelativeImportanceofPredictorVariables",sub=paste("R-Square=",round(rsquare,digits=3)),...)return(import)}#usingrelweights()fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)relweights(fit,col="lightgrey")方差分析1、aov()=lm()单因素方差分析2、plotmeans():绘制带置信区间的图形例：library(multcomp)attach(cholesterol)table(trt)aggregate(response,by=list(trt),FUN=mean)aggregate(response,by=list(trt),FUN=sd)fit<-aov(response~trt)summary(fit)library(gplots)plotmeans(response~trt,xlab="Treatment",ylab="Response",main="MeanPlot\nwith95%CI")detach(cholesterol)多重比较1、TukeyHSD()：对各组均值差异的成对检验例：TukeyHSD(fit)par(las=2)par(mar=c(5,8,4,2))plot(TukeyHSD(fit))par(opar)2、glht()：multcomp包：多重均值比较例