R语言方法总结

1、summary():

例： summary(mtcars[vars])

summary()函数提供了最小值、最大值、四分位数和数值型变量的均值，以及因子向量和逻 辑型向量的频数统计。

2、apply()函数或sapply()函数

计算所选择的任意描述性统计量。mean、 sd、 var、 min、 max、 median、 length、 range 和quantile。函数fivenum()可返回图基五数总括（Tukey’s five-number summary，即最小值、 下四分位数、中位数、上四分位数和最大值）。 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)/n kurt <- sum((x - m)^4/s^4)/n - 3 return(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值来计算平均数的置信区间（默认置信度为0.95：

例： library(pastecs)

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\ \"wt\"), id.vars = c(\"am\cast(dfm, am + cyl + variable ~ ., dstats)

频数表和列联表

1、table()：生成简单的频数统计表

mytable <- with(Arthritis, table(Improved)) Mytable

2、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) mytable

ftable(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\y <- states[, c(\"Life Exp\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包

t 检验

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\ 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\

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.region

var <- 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$weight fitted(fit) residuals(fit)

plot(women$height, women$weight, main = \"Women Age 30-39\ xlab = \"Height (in inches)\多项式回归

例：fit2 <- lm(weight ~ height + I(height^2), data = women)

summary(fit2)

plot(women$height, women$weight, main = \"Women Age 30-39\ xlab = \"Height (in inches)\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\

ylab = \"Weight (lbs.)\") 多元线性回归

1、scatterplotMatrix()：car包

scatterplotMatrix()函数默认在非对角线区域绘制变量间的散点图， 并添加平滑 （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()拟合的模型， xlevels是一个列表，指定变量要设定的常量值， multiline=TRUE选项表示添加相应直线。

例：library(effects)

plot(effect(\"hp:wt\ 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包：学生化残差（studentized residual，也称学生化删除残差或折叠化残差） 例：library(car)

fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)

qqPlot(fit, labels = row.names(states), id.method = \"identify\" ,simulate = TRUE, main = \"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 <- rstudent(fit)

hist(z, breaks=nbreaks, freq=FALSE, xlab=\"Studentized Residual\ main=\"Distribution of Errors\") rug(jitter(z), col=\"brown\")

curve(dnorm(x, mean=mean(z), sd=sd(z)), add=TRUE, col=\"blue\ lines(density(z)$x, density(z)$y, col=\"red\ legend(\"topright\

legend = c( \"Normal Curve\ lty=1:2, col=c(\"blue\

}

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 = \"Index Plot of Hat Values\") abline(h = c(2, 3) * p/n, col = \"red\

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, cook.levels = cutoff) abline(h = cutoff, lty = 2, col = \"red\")

1、influencePlot()：car包：离群点、杠杆值和强影响点的信息整合到一幅图形中 例：influencePlot(fit, id.method = \"identify\ sub = \"Circle size is proportial to Cook's Distance\") 纵坐标超过+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) # define functions

theta.fit <- function(x, y) { lsfit(x, y)

 }

theta.predict <- function(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(\"Original R-square =\

cat(k, \"Fold Cross-Validated R-square =\ cat(\"Change =\}

2、shrinkage()：交叉验证 ;R平方减少得越少，预测则越精确。 例：fit <- lm(Murder ~ Population + Income + Illiteracy + Frost, data = 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 <- R[2:nvar, 2:nvar] rxy <- R[2: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 <- lambda^2

# regression coefficients of Y on orthogonal variables beta <- solve(lambda) %*% rxy rsquare <- colSums(beta^2)

rawwgt <- lambdasq %*% beta^2 import <- (rawwgt/rsquare) * 100 lbls <- names(fit$model[2:nvar])

rownames(import) <- lbls

colnames(import) <- \"Weights\" # plot results

barplot(t(import), names.arg = lbls, ylab = \"% of R-Square\

xlab = \"Predictor Variables\ sub = paste(\"R-Square = \ ...)

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(trt), FUN = mean) aggregate(response, by = list(trt), FUN = sd) fit <- aov(response ~ trt) summary(fit) library(gplots)

plotmeans(response ~ trt, xlab = \"Treatment\ main = \"Mean Plot\\nwith 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, 6, 2))

tuk <- glht(fit, linfct = mcp(trt = \"Tukey\")) plot(cld(tuk, level = 0.05), col = \"lightgrey\") par(opar) 评估检验的假设条件

1、正态检验：library(car)

qqPlot(lm(response ~ trt, data = cholesterol), simulate = TRUE, main = \"QQ Plot\

2、方差齐性检验：bartlett.test(response ~ trt, data = cholesterol) 3、检测离群点：outlierTest() car包 library(car)

outlierTest(fit) 单因素协方差分析

例：data(litter, package = \"multcomp\") attach(litter) table(dose)

aggregate(weight, by = list(dose), FUN = mean) fit <- aov(weight ~ gesttime + dose) summary(fit)

1、effects() ：effects包 ：计算调整的均值 例： library(effects) effect(\"dose\

2、ancova() ：HH包 ：绘制因变量、协变量和因子之间的关系图 例：library(HH)

ancova(weight ~ gesttime + dose, data = litter)

3、interaction.plot() ：函数来展示双因素方差分析的交互效应 例：interaction.plot(dose, supp, len, type = \"b\ \"blue\"), pch = c(16, 18),

main = \"Interaction between Dose and Supplement Type\") 4、plotmeans()：gplots包 ：展示交互效应 例：library(gplots)

plotmeans(len ~ interaction(supp, dose, sep = \" \"), connect = list(c(1, 3, 5), c(2, 4, 6)), col = c(\"red\

main = \"Interaction Plot with 95% CIs\

xlab = \"Treatment and Dose Combination\") 5、interaction2wt()：HH包 ：可视化结果 例:library(HH)

interaction2wt(len ~ supp * dose) 6、colMeans():计算每列的平均值

7、nrow()/ncol :计算数组额行数和列数

8、mahalanobis()：用协方差来计算两点之间距离的方法 稳健多元方差分析

Wilks.test() :稳 健 单 因 素 MANOVA