R語言方法總結(jié)

1、計(jì)算描述性統(tǒng)計(jì)量:1、summary():例： summary(mtcarsvars)summary()函數(shù)提供了最小值、最大值、四分位數(shù)和數(shù)值型變量的均值，以及因子向量和邏輯型向量的頻數(shù)統(tǒng)計(jì)。2、apply()函數(shù)或sapply()函數(shù)計(jì)算所選擇的任意描述性統(tǒng)計(jì)量。mean、 sd、 var、 min、 max、 median、 length、 range和quantile。函數(shù)fivenum()可返回圖基五數(shù)總括（Tukeys five-number summary，即最小值、下四分位數(shù)、中位數(shù)、上四分位數(shù)和最大值）。sapply() 例： mystats <- function(x

2、, 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/s3)/n kurt <- sum(x - m)4/s4)/n - 3 return(c(n = n, mean = m, stdev = s, skew = skew, kurtosis = kurt) sapply(mtcarsvars, mystats)3、describe()： Hmisc包：返回變量和觀測的數(shù)量、缺失值和唯一值的數(shù)目、平均

3、值、分位數(shù)，以及五個(gè)最大的值和五個(gè)最小的值。 例： library(Hmisc) describe(mtcarsvars)4、stat.desc()：pastecs包若basic=TRUE（默認(rèn)值），則計(jì)算其中所有值、空值、缺失值的數(shù)量，以及最小值、最 大值、值域，還有總和。若desc=TRUE（同樣也是默認(rèn)值），則計(jì)算中位數(shù)、平均數(shù)、平均數(shù)的標(biāo)準(zhǔn)誤、平均數(shù)置信度為95%的置信區(qū)間、方差、標(biāo)準(zhǔn)差以及變異系數(shù)。若norm=TRUE（不是默認(rèn)的），則返回正態(tài)分布統(tǒng)計(jì)量，包括偏度和峰度（以及它們的統(tǒng)計(jì)顯著程度）和ShapiroWilk正態(tài)檢驗(yàn)結(jié)果。這里使用了p值來計(jì)算平均數(shù)的置信區(qū)間（默認(rèn)置信度為0

4、.95：例： library(pastecs) stat.desc(mtcarsvars)5、describe()：psych包計(jì)算非缺失值的數(shù)量、平均數(shù)、標(biāo)準(zhǔn)差、中位數(shù)、截尾均值、絕對中位差、最小值、最大值、值域、偏度、峰度和平均值的標(biāo)準(zhǔn)誤 例： library(psych) describe(mtcarsvars)分組計(jì)算描述性統(tǒng)計(jì)量1、aggregate()：例：aggregate(mtcarsvars, by = list(am = mtcars$am), mean)2、by()：例： dstats <- function(x)(c(mean=mean(x), sd=sd(x)

5、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()函數(shù)不允許指定任意函數(shù)，）例：library(psych)describe.by(mtcarsvars, mtcars$am)5、reshape包分組：（重鑄和融合）例：library(res

6、hape)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)頻數(shù)表和列聯(lián)表1、table()：生成簡單的頻數(shù)統(tǒng)計(jì)表mytable <- with(Arthri

7、tis, table(Improved)Mytable2、prop.table()：頻數(shù)轉(zhuǎn)化為比例值prop.table(mytable)3、prop.table()*100：轉(zhuǎn)化為百分比prop.table(mytable)*100二維列聯(lián)表4、table(A,B)/xtabs(A+b,data=mydata)例：mytable <- xtabs( Treatment+Improved, data=Arthritis)5、margin.table()和prop.table()：函數(shù)分別生成邊際頻數(shù)和比例 (1:行，2：列)行和與行比例margin.table(mytable, 1)pr

8、op.table(mytable, 1)列和與列比例margin.table(mytable, 2)prop.table(mytable, 2)prop.table(mytable)6、addmargins()：函數(shù)為這些表格添加邊際和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

9、$Treatment, Arthritis$Improved)多維列聯(lián)表1、table()和xtabs()：都可以基于三個(gè)或更多的類別型變量生成多維列聯(lián)表。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,

10、c(1, 2)ftable(addmargins(prop.table(mytable, c(1, 2), 3)gtable(addmargins(prop.table(mytable, c(1, 2), 3) * 100獨(dú)立檢驗(yàn)1、卡方獨(dú)立性檢驗(yàn) ：chisq.test()例：library(vcd)mytable <- xtabs(Treatment+Improved, data=Arthritis)chisq.test(mytable)mytable <- xtabs(Improved+Sex, data=Arthritis)chisq.test(mytable)2、Fish

11、er精確檢驗(yàn)：fisher.test() 例：mytable <- xtabs(Treatment+Improved, data=Arthritis) fisher.test(mytable)3、Cochran-MantelHaenszel檢驗(yàn)：mantelhaen.test() 例：mytable <- xtabs(Treatment+Improved+Sex, data=Arthritis) mantelhaen.test(mytable)相關(guān)性度量1、assocstats()： 例：library(vcd)mytable <- xtabs(Treatment+Impro

12、ved, data=Arthritis)assocstats(mytable)2、cor()：函數(shù)可以計(jì)算這三種相關(guān)系數(shù)，3、cov()：函數(shù)可用來計(jì)算協(xié)方差例：states <- state.x77, 1:6cov(states)cor(states)cor(states, method="spearman")x <- states, c("Population", "Income", "Illiteracy", "HS Grad")y <- states, c("

13、Life Exp", "Murder")cor(x, y)4、pcor()：偏相關(guān) ggm包例：library(ggm)pcor(c(1, 5, 2, 3, 6), cov(states)相關(guān)性的顯著性檢驗(yàn)1、cor.test()其中的x和y為要檢驗(yàn)相關(guān)性的變量， alternative則用來指定進(jìn)行雙側(cè)檢驗(yàn)或單側(cè)檢驗(yàn)（取值為"two.side"、 "less"或"greater"） ，而method用以指定要計(jì)算的相關(guān)類型（"pearson"、"kendall"或

14、"spearman"）當(dāng)研究的假設(shè)為總體的相關(guān)系數(shù)小于0時(shí)，請使用alternative="less"。在研究的假設(shè)為總體的相關(guān)系數(shù)大于0時(shí)，應(yīng)使用alternative="greater"。在默認(rèn)情況下，假設(shè)為alternative="two.side"（總體相關(guān)系數(shù)不等于0）。 例：cor.test(states, 3, states, 5)2、corr.test()：可以為Pearson、 Spearman或Kendall相關(guān)計(jì)算相關(guān)矩陣和顯著性水平。例：library(psych)corr.test(sta

15、tes, use = "complete")3、 pcor.test()：psych包t 檢驗(yàn)1、t.test(yx,data)（獨(dú)立樣本）例：library(MASS)t.test(Prob So, data=UScrime)2、 t.test(y1,y2,paired=TRUE)(非獨(dú)立) 例：library(MASS)sapply(UScrimec("U1", "U2"), function(x) (c(mean = mean(x), sd = sd(x)with(UScrime, t.test(U1, U2, paired =

16、 TRUE)組間差異的非參數(shù)檢驗(yàn)兩組的比較：1、wilcox.test(yx,data) ：評估觀測是否是從相同的概率分布中抽得例：with(UScrime, by(Prob, So, median)wilcox.test(Prob So, data=UScrime)2、 wilcox.test(y1，y2,paried=TRUE)：它適用于兩組成對數(shù)據(jù)和無法保證正態(tài)性假設(shè)的情境。例：sapply(UScrimec("U1", "U2"), median)with(UScrime, wilcox.test(U1, U2, paired = TRUE)多于

17、兩組的比較：1、 kruskal.test(yA，data):各組獨(dú)立例：states <- as.data.frame(cbind(state.region, state.x77)kruskal.test(Illiteracy state.region, data=states)2、 friedman.test(yA|B,data):各組不獨(dú)立非參數(shù)多組比較:1、 npmc() :npmc包例：class <- state.regionvar <- state.x77, c("Illiteracy")mydata <- as.data.frame(

18、cbind(class, var)rm(class,var)library(npmc)summary(npmc(mydata), type = "BF")aggregate(mydata, by = list(mydata$class), median)回歸用一個(gè)或多個(gè)預(yù)測變量（也稱自變量或解釋變量）來預(yù)測響應(yīng)變量（也稱因變量、效標(biāo)變量或結(jié)果變量）的方法。1、 lm(): 擬合回歸模型 lm(yx1+x2+x3,data) 簡單線性回歸1、 lm(): （data是數(shù)據(jù)框） 例：fit <- lm(weight height, data = women)summary

19、(fit)women$weightfitted(fit)residuals(fit)plot(women$height, women$weight, main = "Women Age 30-39", xlab = "Height (in inches)", ylab = "Weight (in pounds)")多項(xiàng)式回歸例：fit2 <- lm(weight height + I(height2), data = women)summary(fit2)plot(women$height, women$weight, main

20、 = "Women Age 30-39", xlab = "Height (in inches)", ylab = "Weight (in lbs)")lines(women$height, fitted(fit2)2、 scatterplot() ：繪制二元關(guān)系圖例：library(car)scatterplot(weight height, data = women, spread = FALSE, lty.smooth = 2, pch = 19, main = "Women Age 30-39", xlab

21、 = "Height (inches)", ylab = "Weight (lbs.)")多元線性回歸1、 scatterplotMatrix()：car包scatterplotMatrix()函數(shù)默認(rèn)在非對角線區(qū)域繪制變量間的散點(diǎn)圖， 并添加平滑 （loess）和線性擬合曲線。對角線區(qū)域繪制每個(gè)變量的密度圖和軸須圖。例：fit <- lm(Murder Population + Illiteracy + Income + Frost, data = states)有交互項(xiàng)的多元線性回歸例：fit <- lm(mpg hp + wt + hp

22、:wt, data = mtcars)summary(fit)1、 effect() : effects包 ：展示交互項(xiàng)的結(jié)果 term即模型要畫的項(xiàng)， mod為通過lm()擬合的模型， xlevels是一個(gè)列表，指定變量要設(shè)定的常量值， multiline=TRUE選項(xiàng)表示添加相應(yīng)直線。 例：library(effects)plot(effect("hp:wt", fit,xlevels=list(wt = c(2.2, 3.2, 4.2), multiline = TRUE)回歸診斷1、 confint()：求模型參數(shù)的置信區(qū)間 例：fit <- lm(Murde

23、r Population + Illiteracy + Income + Frost, data=states)confint(fit)2、 plot()：生成評價(jià)模型擬合情況的圖形例：fit <- lm(weight height, data = women)par(mfrow = c(2, 2)plot(fit)3、 lm() : 刪除觀測點(diǎn) 例：newfit <- lm(weight height + I(height2), data = women-c(13, 15),)par(mfrow = c(2, 2)plot(newfit)par(opar)gvlma包提供了對所有

24、線性模型假設(shè)進(jìn)行檢驗(yàn)的方法 檢驗(yàn)正態(tài)性：4、qqPlot()：car包：學(xué)生化殘差（studentized residual，也稱學(xué)生化刪除殘差或折疊化殘差）例：library(car)fit <- lm(Murder Population + Illiteracy + Income + Frost, data = states)qqPlot(fit, labels = s(states), id.method = "identify" ,simulate = TRUE, main = "Q-Q Plot")注：id.method

25、= "identify"選項(xiàng)能夠交互式繪圖5、fitted():提取模型的擬合值 例：fitted(fit)“Nevada”6、residuals():二項(xiàng)式回歸模型的殘差 例：residuals(fit)“Nevada”7、 residplot()：生成學(xué)生化殘差柱狀圖（即直方圖），并添加正態(tài)曲線、核密度曲線和軸須圖。它不需要加載car包例：residplot <- function(fit, nbreaks=10) z <- rstudent(fit) hist(z, breaks=nbreaks, freq=FALSE, xlab="Studen

26、tized Residual", main="Distribution 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( "Norma

27、l Curve", "Kernel Density Curve"), lty=1:2, col=c("blue","red"), cex=.7)residplot(fit)誤差的獨(dú)立性8、 durbinWatsonTest() ：驗(yàn)證獨(dú)立性例：durbinWatsonTest(fit)驗(yàn)證線性9、crPlots()：car包成分殘差圖也稱偏殘差圖 例：crPlots(fit)同方差性 （car包的兩個(gè)函數(shù)）10、ncvTest() ：生成一個(gè)計(jì)分檢驗(yàn)，零假設(shè)為誤差方差不變，備擇假設(shè)為誤差方差隨著擬合值水平的變化而變化。若檢驗(yàn)

28、顯著，則說明存在異方差性11、spreadLevelPlot()：添加了最佳擬合曲線的散點(diǎn)圖，展示標(biāo)準(zhǔn)化殘差絕對值與擬合值的關(guān)系。 例：library(car) ncvTest(fit) spreadLevelPlot(fit)線性模型假設(shè)的綜合驗(yàn)證1、 gvlma() ：gvlma包：線性模型假設(shè)進(jìn)行綜合驗(yàn)證，同時(shí)還能做偏斜度、峰度和異方差性的評價(jià) 例：library(gvlma) gvmodel <- gvlma(fit) summary(gvmodel)多重共線性1、 vif() ：car包 ：函數(shù)提供VIF值， >2就表明存在多重共線性問題 例：vif(fit)sqrt(v

29、if(fit) > 2異常觀測值1、 outlierTest() ：car包 ：求得最大標(biāo)準(zhǔn)化殘差絕對值Bonferroni調(diào)整后的p值例：library(car)outlierTest(fit)高杠桿值點(diǎn)1、 hat.plot() :觀測點(diǎn)的帽子值大于帽子均值的2或3倍，即可以認(rèn)定為高杠桿值點(diǎn) 例：hat.plot <- function(fit) p <- length(coefficients(fit) n <- length(fitted(fit) plot(hatvalues(fit), main = "Index Plot of Hat Value

30、s") abline(h = c(2, 3) * p/n, col = "red", lty = 2) identify(1:n, hatvalues(fit), names(hatvalues(fit)hat.plot(fit)強(qiáng)影響點(diǎn) ：Cooks D值大于4/(n-k -1)，則表明它是強(qiáng)影響點(diǎn)，其中n 為樣本量大小， k 是預(yù)測變量數(shù)目。 例：cutoff <- 4/(nrow(states) - length(fit$coefficients) - 2)plot(fit, which = 4, cook.levels = cutoff)abline

31、(h = cutoff, lty = 2, col = "red")1、 influencePlot()：car包：離群點(diǎn)、杠桿值和強(qiáng)影響點(diǎn)的信息整合到一幅圖形中 例：influencePlot(fit, id.method = "identify", main = "Influence Plot", sub = "Circle size is proportial to Cook's Distance")縱坐標(biāo)超過+2或小于-2的州可被認(rèn)為是離群點(diǎn)，水平軸超過0.2或0.3的州有高杠桿值（通常為預(yù)測值的組

32、合）。圓圈大小與影響成比例，圓圈很大的點(diǎn)可能是對模型參數(shù)的估計(jì)造成的不成比例影響的強(qiáng)影響點(diǎn)變量變換1、powerTransform():car包:函數(shù)通過 的最大似然估計(jì)來正態(tài)化變量。例：library(car)summary(powerTransform(states$Murder)2、 boxTidwell():car包：通過獲得預(yù)測變量冪數(shù)的最大似然估計(jì)來改善線性關(guān)系 例：library(car)boxTidwell(Murder Population + Illiteracy, data = states)模型比較1、 anova()：基礎(chǔ)包：比較兩個(gè)嵌套模型的擬合優(yōu)度 例：fit1

33、<- lm(Murder Population + Illiteracy + Income + Frost, data = states)fit2 <- lm(Murder Population + Illiteracy, data = states)anova(fit2, fit1)2、 AIC()：AIC值越小的模型(可以不嵌套)要優(yōu)先選擇，它說明模型用較少的參數(shù)獲得了足夠的擬合度。 例：fit1 <- lm(Murder Population + Illiteracy + Income + Frost, data = states)fit2 <- lm(Murde

34、r 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 Pop

35、ulation + Illiteracy + Income + Frost, data = states, nbest = 4)plot(leaps, scale = "adjr2")交叉驗(yàn)證1、 crossval() 函 數(shù)：bootstrap 包 ：實(shí) 現(xiàn) k 重 交 叉 驗(yàn) 證 例：shrinkage <- function(fit, k = 10) require(bootstrap) # define functions theta.fit <- function(x, y) lsfit(x, y) theta.predict <- functi

36、on(fit, x) cbind(1, x) %*% fit$coef # matrix of predictors x <- fit$model, 2:ncol(fit$model) # vector of predicted values y <- fit$model, 1 results <- crossval(x, y, theta.fit, theta.predict, ngroup = k) r2 <- cor(y, fit$fitted.values)2 r2cv <- cor(y, results$cv.fit)2 cat("Origin

37、al R-square =", r2, "n") cat(k, "Fold Cross-Validated R-square =", r2cv, "n") cat("Change =", r2 - r2cv, "n")2、 shrinkage()：交叉驗(yàn)證 ;R平方減少得越少，預(yù)測則越精確。 例：fit <- lm(Murder Population + Income + Illiteracy + Frost, data = states)shrinkage(fit)相對重要性1

38、、 scale()：將數(shù)據(jù)標(biāo)準(zhǔn)化為均值為0、標(biāo)準(zhǔn)差為1的數(shù)據(jù)集，這樣用R回歸即可獲得標(biāo)準(zhǔn)化的回歸系數(shù)。注意， scale()函數(shù)返回的是一個(gè)矩陣，而lm()函數(shù)要求一個(gè)數(shù)據(jù)框 例：zstates <- as.data.frame(scale(states)zfit <- lm(Murder Population + Income + Illiteracy + Frost, data = zstates)coef(zfit)2、 relweights() ：相對權(quán)重 例：relweights <- function(fit, .) R <- cor(fit$model)

39、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 variables lambda <- evec %*% delta %*% t(evec) lambdasq <- lambda2 # regress

40、ion coefficients of Y on orthogonal variables beta <- solve(lambda) %*% rxy rsquare <- colSums(beta2) rawwgt <- lambdasq %*% beta2 import <- (rawwgt/rsquare) * 100 lbls <- names(fit$model2:nvar) rownames(import) <- lbls colnames(import) <- "Weights" # plot results barp

41、lot(t(import), names.arg = lbls, ylab = "% of R-Square", xlab = "Predictor Variables", main = "Relative Importance of Predictor Variables", sub = paste("R-Square = ", round(rsquare, digits = 3), .) return(import)# using relweights()fit <- lm(Murder Populati

42、on + Illiteracy + Income + Frost, data = states)relweights(fit, col = "lightgrey")方差分析1、 aov() =lm() 單因素方差分析2、plotmeans():繪制帶置信區(qū)間的圖形例：library(multcomp)attach(cholesterol)table(trt)aggregate(response, by = list(trt), FUN = mean)aggregate(response, by = list(trt), FUN = sd)fit <- aov(resp

43、onse trt)summary(fit)library(gplots)plotmeans(response trt, xlab = "Treatment", ylab = "Response", main = "Mean Plotnwith 95% CI")detach(cholesterol)多重比較1、 TukeyHSD()：對各組均值差異的成對檢驗(yàn) 例：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,