如何计算描述统计量

1. summary()函数。

最基本的就是summary函数。可以输出 最小值，最大值，四分位数，数值型变量的均值，以及因子向量和逻辑型向量的频数统计。
以下用mtcars数据集来举例

head(mtcars)
str(mtcars)
myvars <- c("mpg", "hp", "wt")
summary(mtcars[myvars])

2. sapply()函数。

使用sapply函数调用任意函数来实现。
使用方法：sapply(x, FUN, options)
定义一个mystats函数，后通过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来调用自己定义好的函数mystats
sapply(mtcars[myvars], mystats)

3. pastecs包中的stat.desc()函数。推荐使用这个。

stat.desc()这个函数是比较好用的计算描述统计量的函数，可以输出的结果比较多。
用法：stat.desc(x, basic=TRUE, desc=TRUE, norm=FALSE, p=0.95)
basic=TRUE，计算的基础指标包括：数据个数(nbr.val)、空值的个数(nbr.null)、缺失值的个数(nbr.na)、最小值(Min)、最大值(Max)、值域(range，即max-min)、所有非缺失值的总和(Sum)
desc=TRUE，计算：中位数(median)、平均数平均数(Mean)、平均数的标准误(SE.Mean)、平均数置信度为95%的置信区间(CI.Mean)、方差(Var)、标准差(std.dev)和变异系数(coef.var,其定义为标准差除以平均数)。
平均值置信度（95%）指的是在95%的置信度下计算出的平均值的允许误差,可以用平均值+或-这个数来计算置信区间的上限和下限.
norm=TRUE，计算偏度和峰度（以及他们的统计显著程度），s-w正态检验的结果。偏度系数skewness,其显著性判别依据(如果skew.2SE>1，则偏度显著不等于零)、峰度系数kurtosis、其显著性依据(与skew.2SE相同)、正态性的Shapiro-Wilk检验(Normaltest.W)及其概率(normtest.p)

install.packages("pastecs")
library(pastecs)
stat.desc(mtcars[myvars],basic=TRUE, desc=TRUE, norm=TRUE, p=0.95)

如何分组计算描述统计量

1. aggregate()函数。
   使用这个函数，可以计算平均数、标准差等返回值。但是每次只能计算一种。

data <-aggregate(df$decision.RESP, by=list(Subject=df$Subject,Intention=df$intention,Outcome=df$outcome), FUN=mean, na.rm=TRUE)
aggregate(mtcars[myvars], by=list(am=mtcars$am), mean)
aggregate(mtcars[myvars], by=list(am=mtcars$am), sd)

doBy包中的summaryBy()函数 推荐使用这个。和stat.desc搭配使用。

可以同时返回若干个统计量。
用法：summaryBy(formula, data=dataframe, FUN=function)
formula的格式：var1 + var2 + var3 + ... + varN ~ groupvar1 + groupvar2 + ... + groupvarN
～左侧的变量是需要分析的数值型变量，右侧的变量是类别型的分组变量。

function可以是任何内建或者用户自编的R函数。

install.packages("doBy")
library(doBy)
# 根据am的不同水平来 分别计算 am不同水平下 mpg hp wt的描述统计值。调用的函数是stat.desc。
# 这样可以输出stat.desc()函数的很多描述统计值。
summaryBy(mpg+hp+wt~am, data=mtcars, FUN=stat.desc)

直方图

par(mfrow=c(2,2)) # 画4幅图

par(mfrow=c(1,2))实现一页多图的功能
其中，通过设定函数par()的各个参数来调整我们的图形
mfrow=c(2,2) 就是画4幅图，
mfrow=c(3,5),是画15幅图，
例如 par(mfrow=c(2,3)) 一个图版显示2行，3列
作者：巴拉巴拉_9515
链接：https://www.jianshu.com/p/e5f86fd8d033
来源：简书
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

hist(mtcars$mpg)

hist(mtcars$mpg, 
     breaks=12,# 指定组的数量，也就是直方图中柱形的数量。
     # col="red", # 指定颜色
     xlab="Miles Per Gallon", #指定x轴的标题
     main="Colored histogram with 12 bins") #指定图形的标题

hist(mtcars$mpg,
     freq=FALSE, # 根据概率密度而不是频数绘制图形
     breaks=12,
     # col="red",
     xlab="Miles Per Gallon",
     main="Histogram, rug plot, density curve")
rug(jitter(mtcars$mpg,amount=0.01)) 
#轴须图(rug plot)：实际数据值的一种一维呈现方式。
# jitter 将向每个数据点添加一个小的随机值（一个±amount之间的均匀分布随机数），以避免重叠的点产生影响。
lines(density(mtcars$mpg), col="blue", lwd=2) # 添加一条蓝色、宽度2的 密度曲线。

x <- mtcars$mpg 
h<-hist(x, 
        breaks=12,
        # col="red",
        xlab="Miles Per Gallon",
        main="Histogram with normal curve and box") 
xfit<-seq(min(x), max(x), length=40) 
# dnorm是密度函数，也就是正态分布曲线函数 
# dnorm(x, mean = 0, sd = 1, log = FALSE)
yfit<-dnorm(xfit, mean=mean(x), sd=sd(x))
# h$mids[1:2] 含义是 h这一个list里面的mids变量的第一个和第二个数据，这里表示的是，横坐标每个柱形的中间点的坐标。
# diff(h$mids[1:2]) 第一个和第二个坐标的差值
yfit <- yfit*diff(h$mids[1:2])*length(x) 
lines(xfit, yfit, col="blue", lwd=2)
box() #将图形围绕起来的盒型

h[["mids"]]

箱线图

用法：boxplot(formula, data=dataframe)

par(mfrow=c(1,1))
boxplot(mtcars$mpg, main="Box plot", ylab="Miles per Gallon")
a<-boxplot.stats(mtcars$mpg)

a[["stats"]] 画图所需的数据：最小值、下四分位、中位数、上四分位、最大值。
a[["n"]] 数据的个数
a[["out"]] 离群点，是 范围±1.5IQR以外的值。IQR表示四分位距，即上四分位数与下四分位数的差值。
默认情况下，箱线图两条须的延伸极限不会超过 盒型各端±1.5IQR的范围。

boxplot(mpg~cyl, #y~a，将为类别型变量a的每个水平并列的生成数值型变量y的箱线图。
        # 如果有两个变量，则写为 y~a*b 表示为变量a和b所有水平的两两组合生成数值型变量y的箱线图。
        data=mtcars,
        main="Car Mileage Data",
        xlab="Number of Cylinders",
        ylab="Miles Per Gallon")

mtcars[c("cyl","am")]<-lapply(mtcars[c("cyl","am")],as.factor)
boxplot(mpg~cyl*am,
        data=mtcars,
        main="Car Mileage Data",
        xlab="Number of Cylinders",
        ylab="Miles Per Gallon")

boxplot(mpg~am*cyl,
        # 画图的顺序是，cyl的一个水平上，am的两个水平都画完。在画cyl的下一个水平
        data=mtcars,
        col=c("gold","darkgreen"),
        main="MPG Distribution by Auto Type",
        xlab="Auto Type",
        ylab="Miles Per Gallon")

小提琴图

是箱线图和核密度图的结合体。
用法：vioplot(x1, x2, ... , names=, col=)

install.packages("vioplot")
library(vioplot)
x1 <- mtcars$mpg[mtcars$cyl==4]
x2 <- mtcars$mpg[mtcars$cyl==6]
x3 <- mtcars$mpg[mtcars$cyl==8]
vioplot(x1, x2, x3,
        names=c("4 cyl", "6 cyl", "8 cyl"),
        col="gold")
# 图形中白点是中位数，黑色盒型的范围是上下四分位点，细黑线表示须。外部形状即为核密度估计。
title("Violin Plots of Miles Per Gallon", 
      ylab="Miles Per Gallon",
      xlab="Number of Cylinders")

点图

用法：dotchart(x, labels=)

dotchart(mtcars$mpg,
         labels=row.names(mtcars),# labels是由每个点的标签组成的向量
         cex=.7, # 指定标签的大小，包括labels，main和xlab
         main="Gas Mileage for Car Models",
         xlab="Miles Per Gallon")

x <- mtcars[order(mtcars$mpg),]
x$cyl <- factor(x$cyl)
x$color[x$cyl==4] <- "red"
x$color[x$cyl==6] <- "blue"
x$color[x$cyl==8] <- "darkgreen"
dotchart(x$mpg, 
         labels = row.names(x), 
         cex=.7, 
         groups = x$cyl,# 通过groups来选定一个因子，用以指定x中元素的分组方式。
         gcolor = "black",# 控制不同组 标签的颜色。这里是指 4 6 8三个数字的颜色。
         color = x$color,# 点和标签的颜色来自向量color
         pch=19, #在R语言中，点的样式由pch的取值决定。 
         main = "Gas Mileage for Car Models\ngrouped by cylinder", 
         xlab = "Miles Per Gallon")

频数表和列联表

install.packages("vcd")
library(vcd)
head(Arthritis)

table()

table()函数默认忽略缺失值
使用N个类别型变量（因子）创建一个N维列联表
用法：table(var1, var2, ..., varN)

一维列联表

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

prop.table()将频数转化为比例

prop.table(mytable) #prop.table()将频数转化为比例
prop.table(mytable)*100 #prop.table()将频数转化为比例，之后乘以100，变成百分比

二维列联表

mytable <- table(A, B)
A是行变量，B是列变量

mytable<-table(Arthritis$Improved,Arthritis$Treatment) 
x<-prop.table(mytable) 
# 不加数字的话，分母是 Improved和Treatment加起来的总人数。
# 得到的所有比例数字加起来等于1。也就是sum(x)等于1
x<-prop.table(mytable,1) 
# 1表示根据table()语句中的第1个变量(这里指Arthritis$Improved)来计算比例：
# 分别计算 Improved的三个水平的比例，分母是 Improved的三个水平的总数。
# 每一行的比例加起来等于1。apply(x,1,sum)等于1,1,1
x<-prop.table(mytable,2) 
# 2表示根据table()语句中的第2个变量(这里指Arthritis$Treatment)来计算比例：
# 分别计算 Treatment的两个水平的比例，分母是 Treatment 的两个水平的总数。
# 每一列的比例加起来等于1。apply(x,2,sum)等于1,1

mytable <- xtabs(~ A + B, data=mydata)

A是行变量，B是列变量

mytable <- xtabs(~ Treatment+Improved, data=Arthritis)
mytable

margin.table 生成频数

margin.table(mytable, 1)
# 1表示根据table()语句中的第1个变量(这里指Arthritis$Treatment)来计算频数
# 也就是计算 Treatment 的两个水平各自的频数。
margin.table(mytable, 2)
# 1表示根据 table()语句中的第2个变量(这里指Arthritis$Improved)来计算频数
# 也就是计算 Improved 的三个水平各自的频数。
margin.table(mytable)
# 不加数字的话，计算得到的是总人数。

addmargins()

addmargins(mytable) # 添加每行和每列的和
addmargins(mytable,1) # 添加每列的和
addmargins(mytable,2) # 添加每行的和
addmargins(prop.table(mytable,1),2)
addmargins(prop.table(mytable,2),1)

CrossTable() 生成二维列联表 推荐使用这个。可以生成多个参数。

install.packages("gmodels")
library(gmodels)
CrossTable(Arthritis$Treatment, Arthritis$Improved)
# 得到一个表格，包含数据总数、每种条件下的数据个数、以及比例。
# 这里比例包括 三种，分母分别是每一行的总数、每一列的总数、数据总数。

多维列联表

ftable()函数可以 以紧凑的方式输出多维列联表。推荐使用这个。

mytable <- xtabs(~ Treatment+Sex+Improved, data=Arthritis)
mytable
prop.table(mytable, c(1, 2))
# c(1,2)代表 Treatment 和 Sex，计算这两个变量2*2比例的列联表。
# 但是要根据 Improved 的不同水平来分别计算。
ftable(prop.table(mytable, c(1, 2)))
ftable(addmargins(prop.table(mytable, c(1, 2)),1))
ftable(addmargins(prop.table(mytable, c(1, 2)),2))

相关系数和显著性的计算

cor()函数用来计算相关系数;

用法：cor(x, use= , method= )
method有三种，pearson spearman kendall，默认是 pearson

cov()计算协方差。

states<- state.x77[,1:6]
cov(states) #协方差
cor(states) #相关系数
cor(states, method="spearman") # spearman等级相关
x <- states[,c("Population", "Income", "Illiteracy", "HS Grad")]
y <- states[,c("Life Exp", "Murder")]
cor(x,y) #计算得到4*2的相关矩阵

cor.test() 用来对相关系数的显著性进行检验

但是一次只能检验一种相关关系。
用法：cor.test(x, y, alternative = , method = )
alternative 取值有："two.side" "less"或"greater"。默认是 two.side，表示双尾检验（总体相关系数不等于!=0）。
less 表示，假设为 总体相关系数<0；greater 表示，假设为 总体相关系数>0
method有三种，pearson spearman kendall，默认是 pearson
cor.test(states[,3], states[,5])

如果要检验多种相关关系，需要使用 psych包中的corr.test()

用法：corr.test(x, y = NULL, use = "pairwise",method="pearson",adjust="holm",alpha=.05,ci=TRUE,minlength=5)

install.packages("psych")
library(psych)
corr.test(states, use="complete")
# 输出的结果里，对角线上方的数据 是经过多重比较矫正的结果。下方的是没有矫正的p值。
# use的取值 可以是 pairwise 和 complete，分别表示对缺失值执行成对删除 或 行删除。

pcor() 计算偏相关

用法：pcor(u, S)
u是数值向量，前两个数值为要计算相关系数的变量下标，其余的变量为要排除影响的变量下标。

install.packages("ggm")
library(ggm)
colnames(states)
pcor(c(1,5,2,3,6), cov(states)) #偏相关 等价于 pcor(c("Population","Murder","Income","Illiteracy","HS Grad"), cov(states))
# 在本例中，1,5对应的变量("Population","Murder")是计算相关的变量，2,3,6 ("Income","Illiteracy","HS Grad")是需要排除影响的变量。

pcor.test 偏相关系数的显著性检验

用法：pcor.test(r, q, n)
r是偏相关系数，q是要控制的变量数(以数值表示位置)，n为样本大小。
pcor.test(pcor(c(1,5,2,3,6), cov(states)),3,50)# 输出结果为 t值、df值和p值

t检验

独立样本t检验

使用方法：t.test(y ~ x, data) 或t.test(y1, y2)
默认参数是，假设方差不齐的t检验，双侧检验。
如果想要输出方差齐的t检验，需要加入 var.equal=TRUE

library(MASS) #使用MASS包中的数据集来做练习
w<-t.test(Prob ~ So, data=UScrime)
w<-t.test(Prob ~ So, data=UScrime,var.equal=TRUE)

独立样本的非参数检验

wilcoxon秩和检验，也被称作为 mann-whitney u检验

使用方法：wilcox.test(y ~ x, data) 或 wilcox.test(y1, y2)

with(UScrime, by(Prob, So, median))
wilcox.test(Prob ~ So, data=UScrime)

配对样本t检验

不需要方差齐性的假设
使用方法: t.test(y1, y2, paired=TRUE)

sapply(UScrime[c("U1","U2")], function(x)(c(mean=mean(x),sd=sd(x))))
with(UScrime, t.test(U1, U2, paired=TRUE))

非独立(配对)样本的非参数检验

wilcoxon符号秩和检验

使用方法：wilcox.test(y1, y2, paired=TRUE)

sapply(UScrime[c("U1","U2")], median)
with(UScrime, wilcox.test(U1, U2, paired=TRUE))

回归

对于OLS回归，通过使得 预测误差（残差）平方和最小 和 对影响变量的解释力度（R平方）最大，可获得模型参数。

简单线性回归

lm()函数

用法：myfit <- lm(formula, data)

formula的形式为：Y~ X1 + X2 + ... + Xk
~ 左边为响应变量(因变量)，~右边是各个预测变量(自变量)，预测变量之间用+符号隔开。
:符号表示变量的交互项。例如 y~x+z+x:z
*符号表示所有可能交互项的简洁方式。例如代码 y~xzw 等价于 y ~ x + z + w + x:z + x:w + z:w + x:z:w
^符号表示交互项达到某个次数。代码 y~(x+z+w)^2 就代表交互项就到2阶交互。等价于 y ~ x + z + w + x:z + x:w + z:w
.符号表示 除了因变量外的所有变量，这里不包含交互作用。如果一个数据框包含变量 x y z w。那么 代码 y~. 等价于 y ~ x + z + w
-符号表示从等式中移除某个变量。例如 y~(x+z+w)^2-x:w 等价于 y ~ x + z + w + x:z + z:w

-1 表示 删除截距项。y~x-1 拟合y在x上的回归，并强制直线通过原点。截距为0。

summary() 展示拟合模型的详细结果
coefficients() 列出拟合模型的模型参数（截距项和斜率β）
confint() 计算模型参数的置信区间（默认是95%）
fitted() 列出拟合模型的预测值
residuals() 列出拟合模型的残差：预测值-实际值
anova() 1.生成一个拟合模型的方差分析表 或2. 比较两个或多个拟合模型的方差分析表。

vcov() 列出拟合模型的协方差矩阵
AIC() 输出赤池信息统计量
plot() 生成评价拟合模型的诊断图
predict()用拟合模型对新的数据集预测响应变量值

fit <- lm(weight ~ height, data=women)
summary(fit)

结果里面的两个参数 Multiple R-squared 和 Adjusted R-squared 的值：
R-squared 表示模型可以解释预测变量的方差，也就是模型预测变量（预测值）解释响应变量（真实值）的程度。也是实际值和预测值之间相关系数的平方。Adjusted R-squared 与之类似，但是考虑了模型的参数数目。

women$weight
fitted(fit)
residuals(fit)
plot(women$height,women$weight,
     xlab="Height (in inches)",
     ylab="Weight (in pounds)")
abline(fit)

par(mfrow=c(2,2)) #两行两列的图展示在一起
plot(fit) #输出四张图 可以检验正态性 方差齐性 和异常值

线性回归的前提条件

1. 正态性

1.1 使用performance包 check_normality()

library(performance)
x<-check_normality(fit)
plot(x,type = "qq") #画qq图
plot(x, type = "pp") #画pp图
plot(x) #画残差的正态分布图

1.2 不使用包 直接画标准化残差 qq图

qqnorm(residuals(fit))
qqline(residuals(fit))

1.3 使用car包里的qqPlot() 画t分布下的学生化/标准化残差 qq图

library(car)
qqPlot(fit)

1.4 直接画直方图

hist(residuals(fit))

2. 检验自变量和因变量是否是线性关系

library(car)
crPlots(fit)

3. 检验方差齐性

3.1 直接画预测值和残差图 残差没有标准化

plot(fitted(fit),residuals(fit)) 
abline(h=0,lty=2)

3.2 使用performance包 check_heteroscedasticity()

x<-check_heteroscedasticity(fit)
plot(x)

3.3 使用car包中的 ncvTest()来做检验，同时用spreadLevelPlot()画图

library(car)
ncvTest(fit)
spreadLevelPlot(fit)

4. 误差的独立性

检验误差的序列相关性（看到一个资料说，仅仅时序类数据需要进行该诊断）。

使用car包中的durbinWatsonTest()函数来进行 Durbin-Watson 检验。

durbinWatsonTest(fit)

线性模型假设的综合检验

gvlma包中的gvlma()函数

gvlma()函数只适用于 lm()模型

install.packages("gvlma")
library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)

回归中的多重共线性问题

可用统计量 VIF(Variance Inflation Factor,方差膨胀因子) 来检测
vif>5 就认为存在严重的共线性。<5可以忽略

1. 可用car包中的vif()函数来计算得到

library(car)
vif(fit)
# sqrt(vif(fit)) > 2 # problem?

2. 可用performance包中的check_collinearity()函数来计算得到

check_collinearity(fit)

离群点检测

car包中的outlierTest()

outlierTest()函数可以求得最大标准化残差绝对值 Bonferroni 调整后的p值
该函数只根据单个最大（正或负）残差值的显著性来判断是否有离群点。若不显著则说明数据集中没有离群点；若显著，则你必须删除该离群点，然后在检测是否还有其他的离群点存在。

library(car)
outlierTest(fit)

强影响点检测

cutoff <- 4/(nrow(states)-length(fit$coefficients)-2) 
cutoff<-1
plot(fit, which=4, cook.levels=cutoff)
abline(h=cutoff, lty=2, col="red")

1. performance包中的check_outliers()

check_outliers(fit)

二项式回归

fit2 <- lm(weight ~ height + I(height^2), data=women)
summary(fit2)
fitted(fit2)
plot(women$height,women$weight,
     xlab="Height (in inches)",
     ylab="Weight (in lbs)")
lines(women$height,fitted(fit2))
par(mfrow=c(2,2)) 
plot(fit2)

多元线性回归

多元就是指 多个自变量。

states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
cor(states) #做多元回归之前，最好先计算一下变量之间的相关性。
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
summary(fit)
confint(fit)

检验误差的序列相关性（看到一个资料说，仅仅时序类数据需要进行该诊断）。

使用car包中的durbinWatsonTest()函数来进行 Durbin-Watson 检验。

library(car)
durbinWatsonTest(fit)
crPlots(fit)

有显著交互项的多元线性回归

fit <- lm(mpg ~ hp + wt + hp:wt, data=mtcars)
summary(fit)
install.packages("effects")
library(effects) 
# 用effects包的 effect()函数来画图展示 交互项的结果
# 用法：plot(effect(term,mod,,xlevels),multiline=TRUE)
# term是要画的项，mod是lm()拟合的模型，xlevels是一个list，用来指定变量要设定的常量值。
# multiline=TRUE 表示添加多条相应直线。multiline=FALSE 表示多个xlevels分开画多图。
plot(effect("hp:wt", fit,, list(wt=c(2.2,3.2,4.2))), multiline=TRUE)
plot(effect("hp:wt", fit,, list(wt=c(2.2,3.2,4.2))), multiline=FALSE)

模型比较

1. 用anova()函数进行模型比较

注意用anova()函数比较模型时，模型需要是嵌套模型。就比如下面的例子中，一个模型fit1包含了fit2。

states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
fit2 <- lm(Murder ~ Population + Illiteracy, data=states)
anova(fit2, fit1)
# 结果显示p>0.05，并不显著。说明两个模型fit1和fit2的拟合效果一样好。所以可以从模型中去掉"Income", "Frost"两个变量。

2. 用AIC(Akaike Information Criterion，赤池信息准则)来比较模型

AIC考虑了模型统计拟合度以及用来拟合的参数数目。AIC越小，代表模型拟合的越好，说明模型用较少的参数获得了足够的拟合度。

AIC方法不需要模型之间是嵌套关系。

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

怎么选择“最佳”模型

如何确定在模型中保留哪些变量 模型是最优的呢。anova()一般用于有2，3个模型比较的时候。如果模型特别多，该怎么办呢？

1. 逐步回归：向前逐步回归、向后逐步回归、向前向后逐步回归

使用MASS包中的stepAIC()函数

依据的是 精确AIC准则
下面是一个向后回归的例子

library(MASS)
states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
stepAIC(fit, direction="backward")

2. 全子集回归

全子集回归中，所有可能的模型都会被检验。

使用leaps包中的regsubsets()函数

可以通过R平方，调整R平方或Mallows Cp 统计量等准则来选择“最佳”模型

install.packages("leaps")
library(leaps) 
states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
leaps <-regsubsets(Murder ~ Population + Illiteracy + Income + Frost, data=states, nbest=4)
plot(leaps, scale="adjr2")
# 结果可以用car包中的subsets()函数绘制
library(car)
subsets(leaps, statistic="cp",
        main="Cp Plot for All Subsets Regression")
abline(1,1,lty=2,col="red")

交叉验证

通过交叉验证法，可以评价回归方程的泛化能力。

交叉验证：将一定比例的数据挑选出来作为训练样本（训练集），另外的样本作为保留样本（测试集）。先在训练样本上获取回归方程，然后在保留样本上做预测。
在k重交叉验证中，样本被分为k个子样本，轮流k-1个子样本组合作为训练集，另外1个样本作为测试集。
这样就可以得到k个预测方程，将k个预测集的预测表现结果求平均值。

1. R平方的k重交叉验证：bootstrap包中crossval()函数

方差分析

以下的概念中：大写字母表示组别因子，小写表示定量变量。
单因素方差分析ANOVA：y~A
含单个协变量的单因素协方差分析ANCOVA：y~x+A
双因素ANOVA：y~AB
含两个协变量的两因素协方差分析ANCOVA：y~x1+x2+AB
随机化区组：y~B+A（B是区组因子）
单因素组内ANOVA（单因素重复测量anova）：y~A+Error(Subject/A)
含单个组内因子（W）和单个组间因子（B）的重复测量anova：y~B*W+Error(Subject/W)
均衡设计（balanced design）是指 每个实验条件下的被试量相等。否则就是非均衡设计（unbalanced design）。
有多个（大于一个）因变量的设计被称为多元方差分析（MANOVA），若有协变量，则被称为 多元协方差分析（MANCOVA）。

表达式的顺序 在以下两种情况下，会对结果产生影响：

1. 因子不止一个，且是非均衡设计；
2. 存在协变量。

• 类型I方法（序贯型） 计算ANOVA效应：
  因为R中默认 使用类型I方法: y~A+B+A:B
  A对y的影响；A不做调整
  控制A时，B对y的影响；B根据A调整
  控制A和B的主效应时，A与B的交互效应。A:B交互项根据A和B调整。
• 类型II（分层型）：
  效应根据同水平或低水平的效应做调整。A根据B调整，B根据A调整，A:B交互项同时根据A和B调整。
• 类型III（边界型）：SPSS 和 SAS 默认使用的是类型III。所以如果要得到相同的结果，需要使用类型III。aov()函数默认使用的是类型I的方法，若要使用类型III，需要用car包中的Anova()函数，具体可参考 help(Anova,package="car")。
  每个效应根据模型其他各效应做相应调整。A根据B和A:B做调整，A:B交互项根据A和B调整。

样本大小越不平衡，效应项的顺序对结果的影响越大。一般来说，越基础性的效应越需要放在表达式前面。具体来讲：
首先是协变量，然后是主效应，接着是双因素的交互项，再接着是三因素的交互项，以此类推。
对于主效应，越基础性的效应越应放在表达式前面。

aov(formula, data=dataframe)

单因素方差分析 one-way anova

单因素方差分析 对应的 非参数检验 数据之间独立

Kruskal-Wallis检验

使用方法：kruskal.test(y ~ A, data)

states <- data.frame(state.region, state.x77)
kruskal.test(Illiteracy ~ state.region, data=states)
# 通过r语言实战第二版书中 作者写的wmc函数 做事后检验的多重比较
source("http://www.statmethods.net/RiA/wmc.txt")
wmc(Illiteracy ~ state.region, data=states, method="holm")
help(p.adjust) #查看其它多重比较矫正的方法
# p.adjust.methods 
# c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY",
#   "fdr", "none")

单因素重复测量方差分析 对应的 非参数检验 数据之间不独立

Friedman检验

使用方法：friedman.test(y ~ A | B, data)

RoundingTimes <-
  matrix(c(5.40, 5.50, 5.55,
           5.85, 5.70, 5.75,
           5.20, 5.60, 5.50,
           5.55, 5.50, 5.40,
           5.90, 5.85, 5.70,
           5.45, 5.55, 5.60,
           5.40, 5.40, 5.35,
           5.45, 5.50, 5.35,
           5.25, 5.15, 5.00,
           5.85, 5.80, 5.70,
           5.25, 5.20, 5.10,
           5.65, 5.55, 5.45,
           5.60, 5.35, 5.45,
           5.05, 5.00, 4.95,
           5.50, 5.50, 5.40,
           5.45, 5.55, 5.50,
           5.55, 5.55, 5.35,
           5.45, 5.50, 5.55,
           5.50, 5.45, 5.25,
           5.65, 5.60, 5.40,
           5.70, 5.65, 5.55,
           6.30, 6.30, 6.25),
         nrow = 22,
         byrow = TRUE,
         dimnames = list(1 : 22,
                         c("Round Out", "Narrow Angle", "Wide Angle")))
friedman.test(RoundingTimes)
## => strong evidence against the null that the methods are equivalent
##    with respect to speed
wb <- aggregate(warpbreaks$breaks,
                by = list(w = warpbreaks$wool,
                          t = warpbreaks$tension),
                FUN = mean)
wb
friedman.test(wb$x, wb$w, wb$t)
friedman.test(x ~ w | t, data = wb)