R语言方法总结

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

2、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)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 包：返回变量和观测的数量、缺失值和

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

4、p值来计算平均数的置信区间(默认置信度为 0.95:例: library(pastecs)stat.desc(mtcarsvars)5、describe() : psych 包计算非缺失值的数量、 平均数、 标准差、 中位数、 截尾均值、 绝对中位差、 最小值、 最大值、 值域、偏度、峰度和平均值的标准误例： library(psych)describe(mtcarsvars)分组计算描述性统计量1、aggregate() ：例： aggregate(mtcarsvars, by = list(am = mtcars$am), mean)2、by() ：例： dstats <- func

5、tion(x)(c(mean=mean(x), sd=sd(x) by(mtcarsvars, 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(mtcarsvars, mtc

6、ars$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)频数表和列联表*7-1

7、用十创建利愛建列联丸眈茁数临述table ( vatl,x上启ib s ( fa rmu 1丹*上闩ipropt出£生1工1曲HN/gd用1 rriargin . table trarffinslAddmargins (£ Ldh le ( t hl e)BE出N平类别壓变irt l因子）留建一亍贺整M监挖 懒携一个蠢式和一牛机絆或藪据忻创建一十Nit列嬴老 做阳!"小曲朮聲瞬血际艸盘将盅屮基札*肖莎赴阳衣 畑a抚-“足辽餌边年列衣i f R栏中羞目笳円创建一中聿崔俯-f盘列联表1、table():生成简单的频数统计表mytable <- with(Arth

8、ritis, table(lmproved)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+lmproved, data=Arthritis)5、margin.table()和prop.table():函数分别生成边际频数和比例 (1:行，2 :列)行和与行比例margin .table(mytable

9、, 1)prop.table(mytable, 1)列和与列比例margin .table(mytable, 2)prop.table(mytable, 2)prop.table(mytable)6、addmargins():函数为这些表格添加边际和 add marg ins(mytable)admarg in s(prop.table(mytable)addmarg in s(prop.table(mytable, 1), 2)addmarg in s(prop.table(mytable, 2, 1)7. crossTable() : gmodels 包 例: library(gmodel

10、s)CrossTable(Arthritis$Treatme nt, 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)ft

11、able(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=Arthr

12、itis)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)

13、mytable <- xtabs(Treatment+Improved, data=Arthritis) assocstats(mytable)2、cor() : 函数可以计算这三种相关系数，3、cov(): 函数可用来计算协方差例: states <- state.x77, 1:6cov(states)cor(states)cor(states, method="spearman")x <- states, c("Population", "Income", "Illiteracy", &quo

14、t;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、cor.test()其中的 x 和 y 为要检验相关性的变量， alternative 则用来指定进行双侧检验或单侧检验 (取 值为"two.side"、 "less"或"greater")，而 method用以指

15、定要计算的相关类型("pearson"、"kendall"或'spearman")当研究的假设为总体的相关系数小于0时，请使用alternative="less"。在研究的假设为总体的相关系数大于0时，应使用alternative="greater"。在默认情况下，假设为 alternative="two.side" (总体相关系数不等于 0)。例： cor.test(states, 3, states, 5)2、corr.test():可以为Pearson Spearman或

16、Kendall相关计算相关矩阵和显著性水平。例： library(psych)corr.test(states, use = "complete")3、pcor.test() ： psych 包t 检验1、 t.test(yx,data) (独立样本)例： library(MASS)t.test(Prob So, data=UScrime) 2、 t.test(y1,y2,paired=TRUE)( 非独立 )例： library(MASS)sapply(UScrimec("U1", "U2"), function(x) (c(mea

17、n = mean(x), sd = sd(x)with(UScrime, t.test(U1, U2, paired = TRUE)组间差异的非参数检验两组的比较：1、 wilcox.test(yx,data)：评估观测是否是从相同的概率分布中抽得例： with(UScrime, by(Prob, So, median)wilcox.test(Prob So, data=UScrime)2、wilcox.test(y1 ,y2,paried=TRUE):它适用于两组成对数据和无法保证正态性假设的情境。 例：sapply(UScrimec("U1", "U2&quo

18、t;), median)with(UScrime, wilcox.test(U1, U2, paired = TRUE)多于两组的比较：1、kruskal.test(yA ， data):各组独立例: states <- as.data.frame(cbind(state.region, state.x77)kruskal.test(Illiteracy state.regi on, data=states)2、 friedman.test(yA|B,data):各组不独立非参数多组比较：1、npmc():npmc 包例: class <- state.regionvar <

19、- state.x77, c("Illiteracy")mydata <- as.data.frame(cb in d(class, var)rm(class,var)library( npmc)summary (n pmc(mydata), type = "BF")aggregate(mydata, by = list(mydata$class), media n)回归用一个或多个预测变量(也称自变量或解释变量)来预测响应变量(也称因变量、效标变量或结果变量)的方法。1、lm():拟合回归模型 Im(yx1+x2+x3,data)表3讨拟合红性樓

20、型非帶有用的其抱噩熬£敌用世5umnAry()晨小掘A模!0的H却黠卑coe£fcier：ts )网岀扣合模室的根型倾（剛项和料咿confint()f刿出拟合模型的耐值Jualii ()anova()主X牛割音蟆塑的方養分析表或若比较曲牛韓史窖肱音模里时方差分折表vcovf)列岀槻3!擔u的肾方羞麵阵输岀赤池休且维讣址ploLQpredict()1011对新的数棍拿議测响憤童砂简单线性回归1、lm():( data 是数据框)例： fit <- lm(weight height, data = women)summary(fit) women$weight fitte

21、d(fit) residuals(fit)plot(women$height, women$weight, main = "Women Age 30-39", xlab = "Height (in inches)", ylab = "Weight (in pounds)")多项式回归例：fit2 <- lm(weight height + l(height2), data = women) summary(fit2)plot(women$height, women$weight, main = "Women Age 3

22、0-39", xlab = "Height (in inches)", ylab = "Weight (in lbs)")lines(women$height, fitted(fit2)2、 scatterplot()： 绘制二元关系图例： library(car) scatterplot(weight height, data = women, spread = FALSE, lty.smooth = 2, pch = 19, main = "Women Age 30-39", xlab = "Height (i

23、nches)", ylab = "Weight (lbs.)")多元线性回归1、scatterplotMatrix()： car 包scatterplotMatrix() 函数默认在非对角线区域绘制变量间的散点图， 并添加平滑 (loess) 和线性拟合曲线。对角线区域绘制每个变量的密度图和轴须图。例： fit <- lm(Murder Population + llliteracy + lncome +Frost, data = states)有交互项的多元线性回归例： fit <- lm(mpg hp + wt + hp:wt, data = mt

24、cars) summary(fit)1、 effect() : effects 包 ： 展示交互项的结果 term 即模型要画的项， mod 为通过 l m ()拟合的模型，xlevels 是一个列表，指定变量要设定的常量值， multiline=TRUE 选项表示添加相应直线。例： library(effects)plot(effect("hp:wt", fit,xlevels=list(wt = c(2.2, 3.2, 4.2), multiline = TRUE)回归诊断1、 confint() ： 求模型参数的置信区间例： fit <- lm(Murder P

25、opulation + llliteracy + lncome +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)plot (n ewfit)par(opar)carQJM供了大量函大大增舉

26、了拟合和评仰冋归模型的能力（卷见上、4 J-姦2包中的回归诊斷实用苗数£數目的帕Eg"分位SttUil翌tEQnTs; < i时保晋訂筈忤戲Du市in书弗松揃crPlons()诚分与我雄1划ncvTes-l:()对祁恒宜的逞签方普樹帘舒检验占pj：皂mdL尼丄白匕门QUtl i FirTpt-()Hnrr'CTrnnjttt/i 陀驗d.vFlot& ()脈甸的变最阳松i tlI u ence P J o t ()冋旧总昭用iCt±tt«rplot I)啊的RLMscat t-ezrpl o tKatr ix ()憎强的战卢.国审P

27、Tvif ti方基腥艸子gvlma包提供了对所有线性模型假设进行检验的方法检验正态性：4、qqPlot() : car包：学生化残差(studentized residual，也称学生化删除残差或折叠化残差) 例：library(car)fit <- lm(Murder Populati on + Illiteracy + In come + Frost, data = states)qqPlot(fit, labels = s(states), id.method = "identify" simulate = TRUE, main = "

28、Q-Q Plot")注：id.method = "identify"选项能够交互式绘图5、fitted():提取模型的拟合值例：fitted(fit) Nevada”6、residuals():二项式回归模型的残差例：residuals(fit) Nevada”7、residplot():生成学生化残差柱状图(即直方图)，并添加正态曲线、核密度曲线和轴须 图。它不需要加载car包例: residplot <- function(fit, nbreaks=10) z <- rstude nt(fit)hist(z, breaks=n breaks, fr

29、eq=FALSE,xlab="Studentized Residual",main="Distribution of Errors")rug(jitter(z), col="brow n")curve(d no rm(x, mean=mea n( z), sd=sd(z),add=TRUE, col="blue", lwd=2)lin es(de nsity(z)$x, den sity(z)$y,col="red", lwd=2, lty=2)lege nd("topright&qu

30、ot;,lege nd = c( "Normal Curve", "Kernel Den sity Curve"),lty=1:2, col=c("blue","red"), cex=.7)residplot(fit)误差的独立性8、durbinWatsonTest():验证独立性例：durbinWatsonTest(fit)验证线性9、crPlots() : car包成分残差图也称偏残差图例: crPlots(fit)同方差性(car包的两个函数)10、ncvTest():生成一个计分检验，零假设为误差方差不变，

31、备择假设为误差方差随着拟 合值水平的变化而变化。若检验显著，则说明存在异方差性11、spreadLevelPlot()：添加了最佳拟合曲线的散点图，展示标准化残差绝对值与拟合值的 关系。例：library(car)n cvTest(fit)spreadLevelPlot(fit)线性模型假设的综合验证1、gvlma() : gvlma包：线性模型假设进行综合验证，同时还能做偏斜度、峰度和异方差性的评价例: library(gvlma)gvmodel <- gvlma(fit)summary(gvmodel)多重共线性1、vif() : car包：函数提供VIF值， vif >2就表

32、明存在多重共线性问题例：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 <- le ngth(coefficie nts(fit)n <- length(fitted(fit)plot(hatvalues(fit), main = "In

33、 dex Plot of Hat Values")abl in e(h = c(2, 3) * p/n, col = "red", Ity = 2)identify(1:n, hatvalues(fit), names(hatvalues(fit)hat.plot(fit)强影响点：Cook' s D值大于4/(n-k -1),则表明它是强影响点， 其中n为样本量大小，k是 预测变量数目。例：cutoff <- 4/(nrow(states) - length(fit$coefficients) - 2)plot(fit, which = 4, co

34、ok .l evels = cutoff)abli ne(h = cutoff, lty = 2, col = "red")1、in flue ncePlot() ： car包：离群点、杠杆值和强影响点的信息整合到一幅图形中例：in flue ncePlot(fit, id.method = "ide ntify", main = "In flue nce Plot",sub = "Circle size is proportial to Cook's Dista nce")纵坐标超过+2或小于 2的州可被

35、认为是离群点，水平轴超过0.2或0.3的州有高杠杆值(通常为预测值的组合)。圆圈大小与影响成比例，圆圈很大的点可能是对模型参数的估计造成的不成比例影响的强影 响点变量变换*8-5竜见的娈按-2-1-0.500.512变換1叩1、 powerTransform():car包:函数通过 入的最大似然估计来正态化变量x。例：library(car)summary(powerTra nsform(states$Murder)2、boxTidwell():car包：通过获得预测变量幕数的最大似然估计来改善线性关系例：library(car)boxTidwell(Murder Populati on +

36、Illiteracy, data = states)模型比较1、 anova():基础包：比较两个嵌套模型的拟合优度例: fit1 <- lm(Murder Populati on + Illiteracy + In come +Frost, data = states)fit2 <- lm(Murder Populati on + Illiteracy, data = states)ano va(fit2, fit1)2、AIC() : AIC值越小的模型(可以不嵌套)要优先选择，它说明模型用较少的参数获得了足 够的拟合度。例: fit1 <- lm(Murder Popu

37、lati on + Illiteracy + In come +Frost, data = states)fit2 <- lm(Murder Populati on + 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

38、、 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) # de

39、fine functions theta.fit <- function(x, y) lsfit(x, y)theta.predict <- function(fit, x) cbind(1, x) %*% fit$coef# matrix of predictorsx <- fit$model, 2:ncol(fit$model)# vector of predicted valuesy <- fit$model, 1results <- crossval(x, y, theta.fit, theta.predict, ngroup = k)r2 <- c

40、or(y, fit$fitted.values)A2r2cv <- cor(y, results$cv.fit)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(Murd

41、er Population + Income + Illiteracy + Frost, data = states) shrinkage(fit)相对重要性1、scale():将数据标准化为均值为0、标准差为1的数据集，这样用 R回归即可获得标准化的回归系数。注意，scale()函数返回的是一个矩阵，而lm()函数要求一个数据框例: zstates <- as.data.frame(scale(states)zfit <- lm(Murder Population + Income + Illiteracy +Frost, data = zstates)coef(zfit)2、

42、relweights() ： 相对权重 例： relweights <- function(fit, .) R <- cor(fit$model) nvar <- ncol(R) rxx <- R2:nvar, 2:nvar rxy <- R2:nvar, 1 svd <- eigen(rxx) evec <- svd$vectors ev <- svd$values delta <- diag(sqrt(ev)# correlations between original predictors and new orthogonal var

43、iables lambda <- evec %*% delta %*% t(evec)lambdasq <- lambdaA2# regression coefficients of Y on orthogonal variables beta <- solve(lambda) %*% rxyrsquare <- colSums(betaA2) rawwgt <- lambdasq %*% betaA2 import <- (rawwgt/rsquare) * 100 lbls <- names(fit$model2:nvar) rownames(im

44、port) <- lbls colnames(import) <- "Weights"# plot resultsbarplot(t(import), names.arg = lbls, ylab = "% of R-Square",xlab = "Predictor Variables", main = "Relative Importance of Predictor Variables", sub = paste("R-Square = ", round(rsquare, di

45、gits = 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)table(trt) aggregate(response, by = list(

46、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)