R语言方法总结

1、R语言方法总结WTD standardization office WTD 5AB- WTDK 08- WTD 2C计 算 描 述 性 统 计 量1、summary():例 ：summary(mtcarsvars) summary。函数提供了最小值、最大值、四分位数和数值型变量的均值，以及因子向量和逻辑型向量的频数统计。2、apply。函数或sapplyO函数计算所选择的任意描述性统计量。mean、sd、varx min、max. median, length.range和 quantileo 函数 fivenum()可返回图基五数总括(Tukey's five-number summ

2、ary,即最小值、下四分位数、中位数、上四分位数和最大值)。sapplyO例： mystats <- function(x, = FALSE) ifx<- x!(x)m <- mean(x)n <- length(x)s <- sd(x)skew <- sum(x - m)A3/sA3)/nkurt <- sum(x - m)A4/sA4)/n - 3return(c(n = n, mean = m, stdev = s, skew = skew, kurtosis = kurt)sapply(mtcarsvars, mystats)3、describ

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

4、态检验结果。这里使用了 p值来计算平均数 的置信区间(默认置信度为：伍ij ： library(pastecs)(mtcarsvars)5、describe() ： psych 包计算非缺失值的数量、平均数、标准差、中位数、截尾均值、绝对中位差、最小值、最大值、值域、偏度、峰度和平均值的标准误例：libraiy(psych)describe(mtcars vars)分组计算描述性统计量1、aggregate()：例：aggregate(mtcarsvars, by = list(am = mtcarsSam), mean)2、by()：例：dstats <- function(x)(c(

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

6、stats <- function(x) (c(n = length(x), mean = mean(x),sd = sd(x)dfm <- nielt(mtcars, = c("mpg", "hp",nwt"), = c("am", Hcyl")cast(dfm, am + cyl + variable dstats)频数表和列联表1、table()：生成简单的频数统计表mytable <- with(Arthritis, table(Improved)Mytable2、()：频数转化为比例值(

7、mytable)3、()*100 ：转化为百分比(mytable)* 100二维列联表4、table(A,B)/xtabs(-A+b,data=mydata)例：mytable <- xtabs( Treatment+Improved. data=Arthritis)5、。和0 :函数分别生成边际频数和比例(1：行，2 ：列)行和与行比例(mytable, 1)(mytable, 1)列和与列比例(mytable, 2)(mytable, 2)(mytable)6. addmarginsO ：函数为这些表格添加边际和addmargins(mytable)admargins(mytable

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

9、mytable, 3)(mytable, c( 1,3)ftable(mytable, c(l, 2)ftable(addmargins(mytable, c(l, 2), 3)gtable(addmargins(mytable, c(l, 2), 3) * 100独立检验1、卡方独立性检验：()例：library(vcd)mytable <- xtabs(-Treatment+Improved, data=Arthritis)(mytable)mytable <- xtabs(-Improved+Sex, data=Arthritis)(mytable)2、Fisher精确检验：

10、0例：mytable <- xtabs(-Treatment+Improved. data=Ailhritis)(mytable)3、Cochran-Mantel一Haenszel 检验：()例：mytable <- xtabs(-Treatment+Improved+Sex, data=Arthritis)(mytable)相关性度量1、assocstats():例：library(vcd)mytable <- xtabs(-Treatment+Improved, data=Arthritis) assocstats(mytable)2、cor():函数可以计算这三种相关系

11、数，3、cov():函数可用来计算协方差例:states <- , 1:6cov(states)cor(states)cor(states, method=HspearmanH)x <- states, cC'PopuIation"，"Income", "Illiteracy”，"HS Grad”)y <- states, c("Life Exp' "Murder")cor(x, y)4、pcor():偏相关ggm包例：library(ggm)pcor(c( 1, 5, 2, 3,

12、 6), cov(states)相关性的显着性检验Is 0其中的X和y为要检验相关性的变量，alternative则用来指定进行双侧检验或单侧检验(取值为“less”或“greater”)，而method用以指定要计算的相关类型("pearson”、“kendall"或“spearman”）当研究的假设为总体的相关系数小于。时，请使用alternative=less”。在研究的假设为总体的相关系数大于。时，应使用altemative="greater”。在默认情况下，假设为如ernative=”“(总体相关系数不等于0)。例:(states, 3, states,

13、 5)2、():可以为Pearson. Speannan或Kendall相关计算相关矩阵和显着性水平。例：library(psych)(states, use = "complete")3、()： psych 包t检验1、(yx,data)(独立样本)例：library(MASS)(Prob So, data=UScrime)2、(yl,y2,paired=TRUE)俳独立)例：library(MASS)sapply(UScrimec("Ur nU2"), function(x) (c(mean = mean(x),sd = sd(x)with(UScr

14、ime, (UI, U2, paired = TRUE)组间差异的非参数检验两组的比较：1、(y-x,data):评估观测是否是从相同的概率分布中抽得例：with(UScrime, by(Prob, So, median)(Prob So, data=UScrime)2、(yl, y2,paried=TRUE):它适用于两组成对数据和无法保证正态性假设的情境。例：sapply(UScrimec("Ur "U2"), median)with(UScrime, (UI, U2, paired = TRUE)多于两组的比较：1、(yA, data):各组独立(Illit

15、eracy ，data=states)2、(yAIB,data):各组不独立非参数多组比较：1、npmc() :npmc 包例：class <-vai- <- , c("Illiteracy")rm(class,var)libraiy(npmc)sunimary(npmc(mydata), type = "BF”)aggregate(mydata, by = list(mydata$class), median)回归用一个或多个预测变量(也称自变量或解释变量)来预测响应变量(也称因变 量、效标变量或结果变量)的方法。1、lm():拟合回归模型 Im(y

16、-xl+x2+x3,data)简单线性回归lx lm(): (data 是数据框)例:fit <- lm(weight height, data = women)summary(fit)women$weightfitted(fit)residuals(fit)plot(women$height, women$weight, main = HWomen Age 30-39”,xlab = HHeight (in inches)11, ylab = MWeight (in pounds)")多项式回归 例：fit2 <- lm(weight height +1(height八

17、2), data = women)summary(fit2)plot(women$height, womenSweight, main = HWomen Age 30-39”， xlab = "Height (in inches)11, ylab = HWeight (in lbs)”)lines(women$height, fitted(fit2)2、scatterplot():绘制二元关系图l ： library (car)scatterplot(weight height, data = women, spread = FALSE,=2, pch = 19, main = HW

18、omen Age 30-39", xlab = "Height (inches)”, ylab = "Weight (lbs.)")多元线性回归1、scatterplotMatrix() : car 包scatterplotMatrix。函数默认在非对角线区域绘制变量间的散点图，并添加平滑(loess)和线性拟合曲线。对角线区域绘制每个变量的密度图和轴须图。根lj ： fit <- lm(Murder Population + Illiteracy + Income +Frost, data = states)有交互项的多元线性回归例：fit &l

19、t;- Im(mpg hp + wt + hp:wt, data = mtcars)summary(fit)1、effect(): effects包：展示交互项的结果term即模型要画的项，mod为通过lm()拟合的模型，xlevels是一个列表，指定变量要设定的常量值，multiline=TRUE选项表示添加相应直线。例：library(effects)plot(effect(nhp:wtH, fit,xlevels=list(wt = c,),multiline = TRUE)回归诊断1、confintO :求模型参数的置信区间伍ij ： fit <- lm(Murder Popul

20、ation + 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 +1(height八2), data = women-c(13, 15)J)par(mfrow = c(2, 2)plot(newfit)par(opar)gvlma包提供了对所有线性模型假设进行检验的

21、方法检验正态性：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”选项能够交互式绘图5、fitted。:提取模型的拟合值

22、例：fitted(fit)“Nevada”6、residuals。:二项式回归模型的残差例：residuals(fit)“Nevada”7. residplotO ：生成学生化残差柱状图(即直方图)，并添加正态曲线、核密度曲线 和轴须图。它不需要加载car包 例：residplot <- function(fit, nbreaks=10) z <- rstudent(fit)hist(z, breaks=nbreaks, freq=FALSE, x!ab=HStudentized Residual", main=,'Distribution of ErrorsH)

23、 rug(jitter(z), col=nbrownM) curve(dnorm(x, mean=mean(z), sd=sd(z),add二TRUE, col="blue", lwd=2) lines(density(z)$x, density(z)$y, col=nred* lwd=2, lty=2)legend(ntoprightH,legend = c( "Normal Curve", "Kernel Density Curve11), lty=l:2, coI=c(HblueH/,redH), cex=.7)residplot(fit

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

25、验证1、gvlmaO :gvlma包：线性模型假设进行综合验证，同时还能做偏斜度、峰度和异方差性的评价 例：library(gvlma)gvniodel <- gvlma(fit)summary(gvmodel)多重共线性1、vif() ：car包：函数提供VIF值， 阿 >2就表明存在多重共线性问题例：vif(fit)sqrt(vif(fit) > 2异常观测值1、outlleiTest。： car包：求得最大标准化残差绝对值Bonfewoni调整后的p值例：library(car)outlierTest(fit)高杠杆值点1、0：观测点的帽子值大于帽子均值的2或3倍，即可

26、以认定为高杠杆值点例:<-function(fit)p <- length(coefficients(fit)n <- length(fitted(fit)plot(hatvalues(fit), main = HIndex Plot of Hat Values11)abline(h = c(2, 3) * p/n, col = nred, Ity = 2)identify( 1 :n, hatvalues(fit)5 names(hatvakies(fit)(fit)强影响点：Cook*D值大于4/(n-k-l),则表明它是强影响点，其中n为样本量大小,k是预测变量数目。例

27、:cutoff <- 4/(nrow(states) - length(fit$coefficients) - 2)plot(fit, which = 4, = cutoff)abline(h = cutoff, Ity = 2, col = Hredn)1、iiifluencePlot。： car包：离群点、杠杆值和强影响点的信息整合到一幅图形中例：influencePlot(fit, = "identify”，main = "Influence Plot",sub = HCircle size is proportial to Cook's Di

28、stance*1)纵坐标超过+2或小于?2的州可被认为是离群点，水平轴超过或的州有高杠杆值(通常为预测值的组合)。圆圈大小与影响成比例，圆圈很大的点可能是对模型参数的估计造成的不成比例影响的强影响点变量变换1、powerTransform():car包:函数通过人的最大似然估计来正态化变量例：library(car)summary(powerTransform(states$Murder)2、boxTidwen():car :通过获得预测变量靠数的最大似然估计来改善线性关系例：library(car)boxTidwell(Murder Population + Illiteracy, data

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

30、ome +Frost, data = states)fit2 <- lm(Murder Population + Illiteracy, data = states)AIC(fitl,fit2)变量选择1、stepAICO : MASS包:逐步回归模型例：library(MASS)fitl <- lm(Murder Population + Illiteracy + Income +Frost, data = states)stepAIC(fit, direction = "backward”)2、regsubsets() : leaps 包:全子集回归例：library

31、Qe叩s)leaps <- regsubsets(Murder Population + Illiteracy +Income + Frost, data = states, nbest = 4)plot(Ieaps, scale = Hadjr2H)交叉验证1、crossval()函数：bootstrap包：实现k重交叉验证例:shrinkage <- function(fit, k = 10) require( bootstrap)# define functions<-function(x, y) lsfit(x, y)<-function(fit, x) cbi

32、nd(l, x) %*% fit$coef)# matrix of predictorsx <- fit$model, 2:ncol(fitSniodel)# vector of predicted valuesy <- fitSniodelf, 1results <- crossval(x, y, ngroup = k)r2 <- cor(y5 fit$A2r2cv <- cor(y, results$A2cat("Original R-square =", r2, "n")cat(k, "Fold Cross-

33、Validated R-square ="，r2cv, ”rT)cat("Change =", r2 - r2cv, nii")2、shrinkage。：交叉验证;R平方减少得越少，预测则越精确°例：fit <- lm(Murder Population + Income + Illiteracy +Frost, data = states) shrinkage(fit)相对重要性1、scale():将数据标准化为均值为0、标准差为1的数据集，这样用R回归即可获得 标准化的回归系数。注意，scale。函数返回的是一个矩阵，而lm()函数要

34、求一个数据 框zfit <- hn(Murder Population + Income + Illiteracy + Frost, data = zstates)coef(zfit)2、relweights():相对权重例：relweights <- function(fit,.) R <- cor(fitSmodel) nvar v- ncol(R) rxx <- R2:nvar, 2:nvar rxy <- Rnvar, 1 svd <- eigen(rxx) evec <- svdSvectors ev <- svd$vakies del

35、ta v- diag(sqrt(ev) # correlations between original predictors and new orthogonal variables lambda <- evec %*% delta %*% t(evec) lambdasq V- lambda八2# regression coefficients of Y on orthogonal variables beta <- solve(lambda) %*% rxy rsquare <- colSums(betaA2) rawwgt <- lambdasq %*% beta

36、A2 import <- (rawwgt/rsquare) * 100 Ibis <- names(fitSmodel2:nvar) rownames(import) <" Ibis colnames(import) v- "Weights”# plot resultsbarplot(t(import), = Ibis, ylab = "% of R-Squarefxlab = "Predictor Variables", main = "Relative Importance of Predictor Variab

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

38、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 = "Response11, main = "Mean Plotiiwith 95% CI")detach(cholesterol)多重比较1、TukeyHSDO :对各组均值差异的成对检验例：TukeyHSD(fit)par(las = 2)par (mar = c(5, 8, 4, 2) plot(TukeyHSD(fit) par(opar)2、glht() : multcomp :多重均值比较例：library(mu