R语言方法总结

1、计算描述性统计量1、summary():例： summary(mtcarsvars)summary() 函数提供了最小值、最大值、四分位数和数值型变量的均值，以及因子向量和逻 辑型向量的频数统计。range即最小2、apply() 函数或 sapply() 函数 计算所选择的任意描述性统计量。mean、 sd 、 var 、 min 、 max 、 median 、 length和 quantile 。函数 fivenum() 可返回图基五数总括( Tukey s five-number summary ，值、 下四分位数、中位数、上四分位数和最大值)。sapply()例： mystats -

2、 function(x, = FALSE) ifx - x!(x) m - mean(x)n - length(x) s - sd(x)skew - sum(x - m)A3/sA3)/n kurt - sum(x - m)A4/sA4)/n - 3 return(c(n = n, mean = m, stdev = s, skew = skew, kurtosis = kurt)sapply(mtcarsvars, mystats)3、 describe() ：Hmisc 包：返回变量和观测的数量、缺失值和唯一值的数目、平均值、 分位数，以及五个最大的值和五个最小的值。例： library(

3、Hmisc)describe(mtcarsvars)4、() ：pastecs 包若basic=TRUE(默认值)，则计算其中所有值、空值、缺失值的数量，以及最小值、最 值、值域，还有总和。若desc=TRUE(同样也是默认值)，则计算中位数、平均数、平均数的标准误、平均数置信 度为 95%的置信区间、方差、标准差以及变异系数。默认置若norm=TRUE(不是默认的)，则返回正态分布统计量，包括偏度和峰度(以及它们的统计 显着程度)和Shapiro - Wilk正态检验结果。这里使用了 p值来计算平均数的置信区间 信度为： 例： library(pastecs)(mtcarsvars)5、de

4、scribe() ： psych 包 计算非缺失值的数量、平均数、标准差、中位数、截尾均值、绝对中位差、最小值、最大值、值 域、偏度、峰度和平均值的标准误例： library(psych)describe(mtcarsvars)分组计算描述性统计量1、aggregate() ：例： aggregate(mtcarsvars, by = list(am = mtcars$am), mean)2、by() ：例： dstats - function(x)(c(mean=mean(x), sd=sd(x) by(mtcarsvars, mtcars$am, dstats) by(mtcars,var

5、s,mtcars$am,plyr:colwis(dstats)3、summaryBy() ：doBy 包例 library(doBy)summaryBy(mpg + hp + wt am, data = mtcars, FUN = mystats)4、() ：doBy 包( () 函数不允许指定任意函数，)例： library(psych)(mtcarsvars, mtcars$am)5、reshape 包分组 ：(重铸和融合)例： library(reshape)dstats - function(x) (c(n = length(x), mean = mean(x),sd = sd(x)d

6、fm - melt(mtcars, = c(mpg, hp,wt), = c(am, cyl)cast(dfm, am + cyl + variable ., dstats)频数表和列联表1、table() ：生成简单的频数统计表mytable - with(Arthritis, table(Improved)Mytable2、() ：频数转化为比例值(mytable)3、()*100 ：转化为百分比(mytable)*100二维列联表4、table(A,B)/xtabs(A+b,data=mydata)例： mytable - xtabs( Treatment+Improved, data=

7、Arthritis)5、() 和() ： 函数分别生成边际频数和比例(1: 行， 2：列 )行和与行比例(mytable, 1)(mytable, 1)列和与列比例(mytable, 2)(mytable, 2)(mytable)6、addmargins() ：函数为这些表格添加边际和 addmargins(mytable) admargins(mytable)addmargins(mytable, 1), 2) addmargins(mytable, 2, 1)() ： gmodels 包 例： library(gmodels) CrossTable(Arthritis$Treatment,

8、 Arthritis$Improved) 多维列联表1、table() 和 xtabs() ：都可以基于三个或更多的类别型变量生成多维列联表。2、ftable():例： mytable - xtabs( Treatment+Sex+Improved, data=Arthritis) mytable ftable(mytable)(mytable, 1)(mytable, 2)(mytable, 3)(mytable, c(1,3) ftable(mytable, c(1, 2) ftable(addmargins(mytable, c(1, 2), 3) gtable(addmargins(m

9、ytable, c(1, 2), 3) * 100 独立检验1、卡方独立性检验 ： ()例： library(vcd)mytable - xtabs(Treatment+Improved, data=Arthritis)(mytable)mytable - xtabs(Improved+Sex, data=Arthritis)(mytable)2、Fisher 精确检验： ()例： mytable - xtabs(Treatment+Improved, data=Arthritis) (mytable)3、Cochran-Mantel Haenszel 检验： ()例： mytable - x

10、tabs(Treatment+Improved+Sex, data=Arthritis) (mytable)相关性度量1、assocstats() ：例： library(vcd)mytable - xtabs(Treatment+Improved, data=Arthritis) assocstats(mytable)2、cor() ： 函数可以计算这三种相关系数，3、cov() ： 函数可用来计算协方差 例： states - , 1:6cov(states) cor(states) cor(states, method=spearman)x - states, c(Population,

11、 Income, Illiteracy, HS Grad)y - states, c(Life Exp, Murder)cor(x, y)4、pcor():偏相关 ggm包 例： library(ggm)pcor(c(1, 5, 2, 3, 6), cov(states)相关性的显着性检验1、()其中的 x 和 y 为要检验相关性的变量， alternative 则用来指定进行双侧检验或单侧检验(取 值为 、 less 或 greater ) ，而 method 用以指定要计算的相关类型( pearson 、kendall 或spearman)当研究的假设为总体的相关系数小于0时，请使用 al

12、ternative=less 。在研究的假设为总体的相关系数大于 0 时，应使用 alternative=greater 。在默认情 况下，假设为 alternative= (总体相关系数不等于 0)。例: (states, 3, states, 5)2、() : 可以为 Pearson 、 Spearman 或 Kendall 相关计算相关矩阵和显着性水平。 例: library(psych)(states, use = complete)3、() : psych 包t 检验1、(yx,data) (独立样本)例: library(MASS)(Prob So, data=UScrime)2、

13、(y1,y2,paired=TRUE)( 非独立 )例: library(MASS)sapply(UScrimec(U1, U2), function(x) (c(mean = mean(x),sd = sd(x) with(UScrime, (U1, U2, paired = TRUE)组间差异的非参数检验两组的比较:1、(yx,data) : 评估观测是否是从相同的概率分布中抽得例: with(UScrime, by(Prob, So, median)(Prob So, data=UScrime)2、(y1 ， y2,paried=TRUE) : 它适用于两组成对数据和无法保证正态性假设的

14、情境。例: sapply(UScrimec(U1, U2), median) with(UScrime, (U1, U2, paired = TRUE)多于两组的比较:1、(yA ， data): 各组独立(Illiteracy , data=states)2、(yA|B,data): 各组不独立 非参数多组比较 :1、npmc() : npmc包 例: class var - , c(Illiteracy) rm(class,var) library(npmc) summary(npmc(mydata), type = BF)aggregate(mydata, by = list(mydat

15、a$class), median)回归用一个或多个预测变量(也称自变量或解释变量) 来预测响应变量 或结果变量)的方法。1、 lm():拟合回归模型 lm(yx1+x2+x3,data)简单线性回归1、lm():( data 是数据框)例： fit - lm(weight height, data = women)summary(fit) women$weight fitted(fit) residuals(fit) plot(women$height, women$weight, main = Women Age 30-39, xlab = Height (in inches), ylab

16、= Weight (in pounds) 多项式回归例：fit2 - lm(weight height + I(heightA2), data = women) summary(fit2)plot(women$height, women$weight, main = Women Age 30-39, xlab = Height (in inches), ylab = Weight (in lbs) lines(women$height, fitted(fit2)2、 scatterplot()： 绘制二元关系图例： library(car)scatterplot(weight height,

17、data = women, spread = FALSE, = 2, pch = 19, main = Women Age 30-39, xlab = Height (inches), ylab = Weight (lbs.)多元线性回归1、scatterplotMatrix()： car 包scatterplotMatrix() 函数默认在非对角线区域绘制变量间的散点图， 和线性拟合曲线。对角线区域绘制每个变量的密度图和轴须图。例： fit - lm(Murder Population + Illiteracy + Income +Frost, data = states) 有交互项的多元线

18、性回归 例： fit - lm(mpg hp + wt + hp:wt, data = mtcars)summary(fit)1、 effect(): effects 包 ： 展示交互项的结果term 即模型要画的项， mod 为通过 lm() 拟合的模型， xlevels 是 定的常量值， multiline=TRUE 选项表示添加相应直线。例： library(effects)plot(effect(hp:wt, fit,xlevels=list(wt = c, , ), multiline = TRUE)也称因变量、效标变量并添加平滑 （ loess ）个列表，指定变量要设回归诊断 1、

19、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(heightA2), data = women-c(13, 15),) par(mfrow = c(2, 2

20、) plot(newfit) par(opar)gvlma 包提供了对所有线性模型假设进行检验的方法检验正态性：4、 qqPlot()： car 包： 学生化残差( studentized residual ，也称学生化删除残差或折叠化残 差)例： library(car)fit - lm(Murder Population + Illiteracy + Income + Frost, data = states) qqPlot(fit, labels = (states), = identify ,simulate = TRUE, main = Q-Q Plot)注： = identify选

21、项能够交互式绘图5、 fitted(): 提取模型的拟合值例： fitted(fit)“ Nevada”6、 residuals(): 二项式回归模型的残差例： residuals(fit)“ Nevada”7、 residplot()：生成学生化残差柱状图 (即直方图) ，并添加正态曲线、 核密度曲线和轴须图 它不需要加载 car 包例： residplot - function(fit, nbreaks=10) z - rstudent(fit)hist(z, breaks=nbreaks, freq=FALSE,xlab=Studentized Residual,main=Distrib

22、ution of Errors)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( Normal Curve, Kernel Density Curve),lty=1:2, col=c(blue,red), cex=.7)residplot(fit)误差的独立性8、 durbinWatsonTest()

23、：验证独立性例：durbinWatsonTest(fit)验证线性9、crPlots() : car包成分残差图也称偏残差图例：crPlots(fit)同方差性(car包的两个函数)10、ncvTest() :生成一个计分检验，零假设为误差方差不变，备择假设为误差方差随着拟合 值水平的变化而变化。若检验显着，则说明存在异方差性11、spreadLevelPlot() :添加了最佳拟合曲线的散点图，展示标准化残差绝对值与拟合值的关 系。例：library(car)n cvTest(fit) spreadLevelPlot(fit)线性模型假设的综合验证1、gvlma() : gvlma包：线性模

24、型假设进行综合验证，同时还能做偏斜度、峰度和异方差性的 评价例：library(gvlma)gvmodel 2就表明存在多重共线性问题例：vif(fit)sqrt(vif(fit) 2异常观测值1、outlierTest() : car包：求得最大标准化残差绝对值Bonferroni调整后的p值例：library(car)outlierTest(fit)高杠杆值点1、():观测点的帽子值大于帽子均值的2或3倍，即可以认定为高杠杆值点例：-fun ctio n(fit)p - len gth(coefficie nts(fit)n - length(fitted(fit)plot(hatvalu

25、es(fit), main = Index Plot of Hat Values)abline(h = c(2, 3) * p/n, col = red, lty = 2)identify(1:n, hatvalues(fit), names(hatvalues(fit)(fit)强影响点：Cook s D值大于4/(n-k -1)，则表明它是强影响点，其中n为样本量大小,预测变量数目。例：cutoff - 4/(nrow(states) - length(fit$coefficients) - 2)plot(fit, which = 4, = cutoff)abli ne(h = cutof

26、f, lty = 2, col = red)1、in flue ncePlot() : car包：离群点、杠杆值和强影响点的信息整合到一幅图形中例：in flue ncePlot(fit, = ide ntify, mai n = In flue nce Plot,sub = Circle size is proportial to Cooks Distance) 纵坐标超过 +2 或小于 ?2 的州可被认为是离群点， 水平轴超过或的州有高杠杆值(通常为预测值的组合)。 圆圈大小与影响成比例， 圆圈很大的点可能是对模型参数的估计造成的不成比例影响的强影响点 变量变换1、 powerTransf

27、orm():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,

28、 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、step

29、AIC() : 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, scal

30、e = adjr2)交叉验证1 、 crossval() 函 数： bootstrap 包 ： 实 现 k 重 交 叉 验 证 例： shrinkage - function(fit, k = 10) require(bootstrap)# define functions- function(x, y) lsfit(x, y)- function(fit, x) cbind(1, x) %*% fit$coef# matrix of predictorsx - fit$model, 2:ncol(fit$model)# vector of predicted valuesy - fit$mo

31、del, 1results - crossval(x, y, , , ngroup = k)r2 - cor(y, fit$A2r2cv - cor(y, results$A2cat(Original R-square =, r2, n)cat(k, Fold Cross-Validated R-square =, r2cv, n)cat(Change =, r2 - r2cv, n)2、shrinkage() ： 交叉验证 ;R 平方减少得越少，预测则越精确。例： fit - lm(Murder Population + Income + Illiteracy +Frost, data =

32、states)shrinkage(fit)相对重要性1、 scale() ：将数据标准化为均值为0、标准差为 1 的数据集，这样用 R 回归即可获得标准化的 回归系数。注意， scale() 函数返回的是一个矩阵，而 lm() 函数要求一个数据框zfit - lm(Murder Population + Income + Illiteracy +Frost, data = zstates)coef(zfit)2、 relweights()： 相对权重例： relweights - function(fit, .) R - cor(fit$model)nvar - ncol(R)rxx - R2

33、:nvar, 2:nvarrxy - R2:nvar, 1svd - eigen(rxx)evec - svd$vectorsev - svd$valuesdelta - diag(sqrt(ev)# correlations between original predictors and new orthogonal variableslambda - evec %*% delta %*% t(evec)lambdasq - lambdaA2# regression coefficients of Y on orthogonal variablesbeta - solve(lambda) %

34、*% rxyrsquare - colSums(betaA2)rawwgt - lambdasq %*% betaA2import - (rawwgt/rsquare) * 100lbls - names(fit$model2:nvar)rownames(import) - lblscolnames(import) - Weights# plot results barplot(t(import), = lbls, ylab = % of R-Square,xlab = Predictor Variables, main = Relative Importance of Predictor V

35、ariables,sub = paste(R-Square = , round(rsquare, digits = 3),.) return(import)# using relweights()fit - lm(Murder Population + Illiteracy + Income +Frost, data = states) relweights(fit, col = lightgrey) 方差分析1、 aov() =lm() 单因素方差分析2、 plotmeans(): 绘制带置信区间的图形 例： library(multcomp)attach(cholesterol) tabl

36、e(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 = Mean Plotnwith 95% 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 包： 多重均值比较 例： library(multcomp) par(mar = c(5, 4,