R语言和社会统计分析—个人总结（后续更新和优化结构）
基本语法
包的安装和下载
# 获取包含R包的库位置
> .libPaths()
# 查看已经安装的包
> library()
# 安装包
> install.packages("packagename")
# 加载包
> library(packagesname)
# 查看已经加载的包
> (.packages())
# 卸载加载的包（注意不是删除包）
> detach("package:packagename")
# 删除包
> remove.packages("packagename")
向量
# 创建向量
> a <- c(1, 4)
> a <- rep(x, y)
> a <- seq(1, 10, 2)
# 命名向量
> names(a) <- c('a', 'b')
# 向量相加（**保证两个向量的长度是相一致的**）
矩阵
# 创建矩阵
> a <- maxtrix(c(1: 12), ncol=, nrow=)
# 命名矩阵
> colnames(a) <- c('a', 'b', 'c')
> rownames(a) < c('a', 'b', 'c')
# 矩阵运算（行列相等）
> %*%
# 矩阵的行列运算
> mean(dat[1,]) # 对列计算平均值
> apply(dat,1,mean) # 对dat的第一行计算平均值
# 逻辑判断
> dat[dat[1,]>2]
> which(dat[1]>2)
因子
# 使用factor()函数
> f <- factor(x=charactor(), levels, labels=levels, exclude = NA, ordered = is.ordered(x), namax = NA)
'''
x 为创建因子的数据，是一个向量；
levels：因子数据的水平，默认是x中不重复的值;
labels：标识某水平的名称，与水平一一对应，以方便识别，默认取levels的值；
exclude：从x中剔除的水平值，默认为NA值；
ordered：逻辑值，因子水平是否有顺序（编码次序），若有取TRUE，否则取FALSE；
nmax：水平个数的限制。
'''
数组
# 创建数组
> array(1:12,dim=c(2,3,2)) # 使用较少，不做叙述
列表
# 创建列表
> > list_data <- list("Red", "Green", c(21,32,11), TRUE, 51.23, 119.1)
# 列表索引命名
> names(list_data) <- c("1st Quarter", "A_Matrix", "A Inner list")
# 列表索引
> print(list_data[1])
> print(list_data$A_Matrix)
# 添加元素
> list_data[4] <- NULL # 直接赋值
# 合并列表
> merged.list <- c(list1,list2)
# lapply循环
> lapply(list,mean) # 列表有用
> apply(x,mean) # 对向量有用
# 列表可以被转换为一个向量，以便能用于进一步操纵向量的元素, 所有关于向量的算术运算可以在列表被转换为矢量之后被应用。要做到这一点转换，使用`unlist()` 函数。它以列表作为输入，并产生一个向量
数据框
# 创建数据框
> data.frame(col1, col2, ..., row.name=NULL, check.rows = FALSE, check.names=TRUE, stringsAsFactors = default.stringsAsFactors())
其中，row.name用于指定各行（样本）的名称，默认没有名称，使用从1开始自增的序列来标识每一行；
check.rows用于用来检查行的名称和数量是否一致，默认为FALSE；
check.names来检查变量（列）的名称是否唯一且符合语法，默认为TRUE；
用来描述是否将字符型向量自动转换为因子，默认转换，若不改变的话使用stringsAsFactors = FALSE来指定即可。
# 数据框的引用
（1）df1["score"] #仍为一个数据框, 也是一个列表
（2）df1[,"score"] #返回的是向量
（3）df1[3] #同（1）
（4）df1[,3] #同（2）
（5）df1[c(1,3)] #返回第1列和第3列的数据
（6）df1[c(1,3),] #返回第1行和第3行的数据
（7）df1[c(1,3),c(2,3)] #返回第1行和第3行与第2列和第3列交叉处的数据
（8）df1$name #以因子的形式返回name列
（9）df1[["name"]] #以因子的形式返回name列
（10）df1[[1]][1] #返回第1分量的第一个元素值，王宏
（11）df1[['name']][1] #返回name分量第一个元素值：王宏
（12）df1$name[1] #返回name分量第一个元素值：王宏
# 增加样本量或者数据
# 可以使用rbind()函数和cbind()函数将新行或新列添加到数据框变量中
> rbind() # 增加行
> cbind() # 增加列
# 修改某列的第一个值
> df1$name[1] <- "王宏伟" #将王宏的值修改为王宏伟
> df1[1,2] <- "女" #将第一行第2列的值修改为“女”
> df1[[1]][2]<-"马兰兰" #将第一列第二个值改为“马兰兰”
# 删除行或列
> df1<-df1[-2,] #删除第2行数据
> df1<-df1[,-3] #删除第3列的数据
> df1<-df1[-c(1,3),] #删除第1行和第3行的数据
# 数据框连接
> merge(data1, data2, by='ID')
读取和保存数据
'''读取数据'''
# 读取csv
> data <- read.csv('.\\统计学\\example\\ch1\\table1_1.csv')
> head(data,6) # 读取前 6行的数据
# 读取 Excel数据
> library(xlsx) #需要安装 xlsx 包
> data <- read.xlsx("file",n) # n 为要导入工作表的序号
# 读取 spss数据
> library(foregin) # 已经默认安装
> data <- read.spss("file",use.value.labels=TRUE,as.data.frame=TRUE)
# 读取 R格式数据
> data <- load('.\\统计学\\example\\ch1\\example1_1.RData')
'''保存数据'''
# 保存 R格式数据
> save(data,file = '.\\...\\name.Rdata')
# 保存 csv格式数据
> write.csv(data,file = '.\\...\\name.csv')
# 保存 xlsx格式
> library(xlsx)
> write.xlsx(data, "data.xlsx",sheet.name="sheet1")
条件语句
if语句
x <- 30L # R语言中，在正整数后加 L来表示整型数据（正整数）
if(is.integer(x)) {
print("X is an Integer")
}
1234
[1] "X is an Integer"
1
if…else语句
y <- list('a', 'v', 'd')
if('a' %in% y){ # %in% 运算符 检查元素是否在向量中
print('a is in list')
}else{ # 注意这里的 else语句并不在if的花括号当中
print('a is not in list')
}
1234567
[1] "a is in list"
1
x <- c("what","is","truth")
if("Truth" %in% x) {
print("Truth is found the first time")
} else if ("truth" %in% x) {
print("truth is found the second time")
} else {
print("No truth found")
}
123456789
[1] "truth is found the second time"
1
switch语句
# 创建一个函数，输入的值和选择的函数类型来输出结果。
> centre <- function(x, type) {
switch(type,
mean = mean(x),
median = median(x),
trimmed = mean(x, trim = .1))
}
centre(c(1,2,4,5),'mean')
12345678910111213
3
1
循环语句
while循环
ant <- 2
while(ant<5){
print('hello')
ant = ant + 1
}
12345
[1] "hello"
[1] "hello"
[1] "hello"
123
for循环
> v <- LETTERS[1:4] # LETTERS为26个大写字母向量。
for(i in v){
print(i)
}
12345
[1] "A"
[1] "B"
[1] "C"
[1] "D"
1234
repeat循环
> i <- 1
> sum <- 0
repeat
{
sum = sum + i
if( i >= 100) #如果已循环加到了100，则使用break跳出repeat循环
break
i <- i + 1
}
print(sum)
12345678910
[1] 5050
1
next语句
R语言存在next语句，当我们想跳过循环的当前迭代而不终止它时便可使用next。 遇到next时，R解析器跳过本次迭代，并开始循环的下一次迭代。
> k <- LETTERS[1:6]
for ( i in k) {
if (i == "D") {
next
}
print(i)
}
1234567
[1] "A"
[1] "B"
[1] "C"
[1] "E"
[1] "F"
12345
数据探索
集中程度
# 均值：mean()
# 中位数：median()
# 分位数：quantile()
分散程度
# 极差，方差，标准差
# range(),var(),sd()
# 标准误和偏度系数，峰度系数
数据摘要表
# summary()
找出重数（向量）
> a <- c(1,2,3,4,6)
> table(a)
缺失值处理
> complete.cases() # 返回一个逻辑向量,输出的是每行的缺失情况
> no.omit() # 直接返回去除含有缺失值所在行后的数据
# which(is.na(data), arr.ind=T)
> subset(data, data[, 3] == '女')
排序
> sort(x,decreasing=False) # 返回升序向量
> order(x,decreasing=False) # 返回升序后的下标, decreasing是针对类别型向量的
> order(-x) # 降序排列
三角函数
> sin() # sin(a) a表示弧度
> cos()
> tan()
> asin()
> acos()
> atan()
判断数据类型
> class(a)
> character(a)
数据格式转换
> 简易，不做叙述
绘图基础
- 主体
| 特征 | 参数 |
|---|---|
| 符号种类 | pch |
| 符号大小 | cex |
| 符号填充色 | bg |
| 线条线型 | lty |
| 线条宽度 | lwd |
| 颜色 | col |
- 坐标轴
| 特征 | 参数 |
|---|---|
| 刻度位置 | at |
| 刻度长度和方向 | tcl |
| 刻度横坐标范围 | xlim |
| 刻度纵坐标范围 | ylim |
| 刻度文字内容 | label |
| 刻度文字颜色，大小，字体 | col.axis，cex.axis，font.axis |
- 坐标标题
| 特征 | 参数 |
|---|---|
| 标题横坐标内容 | xlab |
| 标题纵坐标内容 | ylab |
| 标题文字颜色 | col.lab |
| 标题文字大小 | cex.lab |
| 标题文字字体 | font.lab |
| 特征 | 参数 |
|---|---|
| 主标题内容 | main |
| 副标题内容 | sub |
| 主标题文字，大小，颜色 | col.main，cex.main，font.main |
| 副标题文字，大小，颜色 | col.sub，cex.sub，font.sub |
- 图形尺寸，边界，布局
> par(mforw=c(行数，列数), mar=c（n1, n2, n3, n4))
> layout(matrix(c(1, 2, 3, 4), 2, 2, byrow=TRUE))
| 特征 | 参数 |
|---|---|
| 图形的宽和高 | pin |
| 以数值向量为边界大小，下，左，上，右，英寸 | mai |
| 以数值变量为边界代销，下，左，上，右，英分，默认c（5, 4, 4, 2）+0.1 | mar |
绘图工具
- 条形图
'''假设x为一个向量'''
> barplot(x, horiz=TRUE) # 生成平行条形图
> barplot(x) # 生成垂直条形图
> barplot(x, main='', xlab='', ylab='')
'''条形图的微调:par'''
> cex.names=? # 设置字体大小
> names.arg=c('') # 自己定义标签文本
- 堆砌条形图
'''假设x为一个矩阵'''
> barplot(x, main='', xlab='', ylab='',col=c(''), legend=rownames(x))
- 分组条形图
'''假设x为一个矩阵'''
> barplot(x, main='', xlab='', ylab='',col=c(''), legend=rownames(x), beside=TRUE)
# 条形图不一定要基于计数数据或频率数据，可以使用数据整合函数并将结果传递给barplot函数
- 其他图形：箱线图，核密度图，饼图，小提琴图，直方图，点图
基本统计分析
###本章中所使用的方法过多，因此一般情况下只介绍一种常用的###
本章描述：对于各个向量的一些基本统计分析的指标进行计算
描述性统计分析—连续型变量统计
描述：分析连续型变量的中心趋势，变化性和分布形状。
主要方法（也是后期用的比较多的方法）
> myvars <- c('mpg', 'wt', 'hp')
> summary(myvars)
'''返回最小值min，最大值max，四分位数，数值型变量的均值'''
分组计算描述性统计量
> myvars <- c('mpg', 'wt', 'hp')
> aggregate(myvars, by=list(am=mtcars$am), mean)
频数表和列联表—类别型变量统计
- 一维列联表（一个向量内独自比较）
> mytable <- with(iris, table(iris$Species))
> prop.table(mytable) # 把频数转换成比例值
'''table()函数默认忽略缺失值，要在频数表中视NA为一个有效的类别，请设定参数useNA = "ifany"'''
- 二维列联表（两个向量之间相互比较）
> mytable <- table(A, B)
> mytable <- xtabs(~A+B, data=mydata)
> prop.table(mytable, 1) # 把频数转换成比例值,其中的1表示的第几个变量
> library(gmodels) # 第二种方法
> CrossTable(A, B)
- 三维列联表
> mytable <- xtabs(A + B + C, data = mydata)
> ftable(mytable) # 生成多层索引
> margin.table(mytable, 1) # 生成边际频数
> margin.table(mytable, c(1,3)) # 生成第一个变量乘以第三个变量的边际频数
> ftable(prop.table(mytable, c(1,3))) # 生成第一个变量乘以第三个变量的改善比例
- 独立性检验：检验类别型变量的独立性的方法
- 卡方独立性检验
> chisq.test(mytable) # 传入的数据是一个二维列联表
'''如果两种要素之间不存在某种关系，则p > 0.05;如果两种要素之间存在某种关系, 则p < 0.01'''
- Fisher精确检验
> fisher.test(mytable) # mytable的行列必须大于2
t检验—找到两个向量是否相关（两个向量不全为类别型变量）
- 独立样本下的t检验
> t.test(A~B, data=mydata)
'''如果两种要素之间不存在某种关系，则p < 0.001'''
- 非独立性下的t检验
> t.test(A, B, paired=TRUE)
回归
本章描述：用一个或多个预测变量来预测响应变量,变量都是连续型数值
OLS回归
- 用lm（）拟合回归模型
> lm(y~x1+x2+x3......xn) # 右边为响应变量，左边为预测变量
'''查看两个变量之间是否有关系看p值和0.001的大小'''
| 符号 | 用途 |
|---|---|
| ~ | 分隔符号，左边为预测变量，右边为解释变量 |
| + | 分隔预测变量 |
| ： | 表示预测变量的交互项 |
| * | 表示所有可能交互项的简洁方式 |
| ^ | 表示交互项达到某个次数 |
| . | 表示包含除因变量之外的所有变量 |
| - | 减号，表示从等式中移除某个变量 |
| -1 | 删除截距项，如y~x-1拟合y在x上的回归，并强制直线通过原点 |
| I（） | 从算术的角度来解释括号中的元素，如y~x + I（z + w）等于y~x + h |
| function | 可以在表达式中用的数学函数 |
| 函数 | 用途 |
|---|---|
| summa() | 展示拟合模型的详细结果 |
| coefficients() | 列出拟合模型的模型参数（截距想和斜率） |
| confint() | 提供模型参数的置信区间（默认95%） |
| residuals() | 列出拟合模型的残差值 |
| anova() | 生成一个拟合模型的方差分析表，或者比较两个或更多拟合模型的方差分析表 |
| vcvo() | 列出模型参数的协方差矩阵 |
| AIC() | 输出赤池信息统计量 |
| plot() | 生成评价拟合模型的诊断图 |
| predict() | 用拟合模型对新的数据集预测响应变量值 |
| fitted() | 列出拟合模型的预测值 |
-
简单线性回归
-
多项式线性回归
-
多元线性回归
# 监测二变量之间的关系
> cor(A, B, data=mydata)
回归诊断—用来告诉模型是否合适
- 一般方法
> fit <- lm(A~B， data=mydata)
> par(mfrow=c(2,2))
> plot(fit)
- 改进方法
- 正态性
> library(car)
> fit <- lm(A~B， data=mydata)
> qqplot(fit, id.method='identify', simulate=TRUE, mian='Q-Q plot')
# 如果不符合，使用变量变换
> summary(powerTransform(A)) # 使用最大似然值来正态化变量x²
- 误差独立性
> durbinWatsonTest(fit) # 检验p值
- 线性
# 通过成分残差图或者偏残差图，查看因变量和自变量是否具有线性关系
> library(car)
> crPlots(fit)
# 如果不符合，使用变量变换
> boxTidwell(A, B, data=mydata) # 获得预测变量幂数的最大似然值估计来改善线性关系
- 同方差性
> library(car)
> ncvTest(fit) # 检验p值大小
> spreadLevelPlot(fit) # 绘制图形
- 总和验证方法
> library(gvlma)
> gvmodel <- gvlma(fit)
> summary(gvmodel)
- 多重共线性—p值成立，但是回归系数不显著, 是由于某些定量原因改变导致的，比如四舍五入等等
> library(car)
> vif(fit)
> sqrt(vif(fit)) > 2
异常观测值
- 离群点
> library(car)
> outlierTest(fit)
-
高杠杆值点
-
强影响点
> library(car)
> avPlots(fit, ask=FALSE, id.method="identify")
- 总和统计三点
> library(car)
> influencePlot(fit, id.method="identify", main='', sub='')
选择最佳的回归模型
> anvoa(fit_1, fit_2) # 模型比较，AIC值较小的优先选择
# 变量选择
> library(MASS)
> fit <- lm(A, B, data=mydata)
> stepAIC(fit, direction='backward') # 向后回归
交叉验证—验证预测目标在真实世界的效果
# 最优目标：预测误差（残差）平方和最小和对响应变量的解释度（R平方）最大，可获得最优模型
相对重要性
> zfit <- lm(A, B, data=mydata)
> coef(zfit) # 计算权重值
方差分析
本章描述：在对于名义（类别）型因子和有序（有序类别）型因子作为预测变量建模，当包含的因子是解释变量的时候，重点通常会从预测转移到组别差异的分析
ANVOA模型拟合
- aov函数
| 设计 | 表达式 |
|---|---|
| 单因素 | y ~ A |
| 含单个协变量的单因素 | y ~ x + A |
| 双因素 | y ~ B + A |
| 含两个协变量的双因素 | y ~ x1 + x2 + A*B |
| 随机化区组 | y ~ B + A（B是区组因子） |
| 单因素组内 | y ~ A + Error |
| 含单个组内因子（w）和单个组内因子（B）的重复测量 | y ~ B * W + Error |
单因素方差分析
> fit <- aov(A ~ B)
> summary(fit)
# 绘制图像
> library(gplots)
> plotmeans(A~B)
多重比较
'''各种类别型因子x对于变量y有差异的时候，但是具体不知道哪种不同，因此使用多重比较来查看'''
> TuKeyHSD(fit) # 查看p值
# 绘制图像
> par(las=2)
> par(mar=c(5, 8, 4, 2))
> plot(TuKeyHSD(fit))
- 评估条件的假设条件(方差诊断)
# 正态性检验
> library(car)
> qqplot(lm(A~B, data=mydata), simulate=TRUE, mian='Q-Q plot', labels=FALSE)
# 方差齐性
> bartlett.test(y~x, data=mydata)
# 监测离群点
> library(car)
> outlierTest(fit)
单因素协方差分析
协变量：包含一个或者多个变量的协变量（不可控因素）
> fit <- aov(y~ x+a)
# 查看未调整的均值
> aggregate(y, by=list(x), FUN=mean)
# 去除协变量后的均值
> library(effects)
> effect("a", fit)
# 多重比较
-
评估检验的假设条件
检验回归斜率的同质性（在交互作用下）
> fit <- aov(y~x*a, data=mydata)
> summary(fit)
- 结果可视化
> library(HH)
> ancova(y~x + a, data=mydata)
双因素方差分析
# 结果可视化
interaction.plot(y,x)
重复测量方差分析
多元方差分析
> y <- rbind(A, B ,C, data=mydata)
# 检验多元正态性
> center <- colMeans(y)
> n <- nrow(y)
> p <- ncol(y)
> cov <- cov(y)
> d <- mahalanobis(y, center, cov)
> coord <- qqplot(qchisq(ppoints(n), df=p),
d, main='', ylab='')
# 检验多元离群点
> library(mvoutlier)
> outliers <- aq.plot(y)
> outliers
# 稳健多元方差分析
> library(rrcow)
> wilks.testly(y, x, method="mcd")
功效分析
本章描述：功效分析可以在给定置信度的情况下，判断监测到给定效应值时所需样本量
假设检验速览
定义：在统计假设检验中，首先要对总体分布参数设定一个假设（零假设，H0），然后从总体分布中抽样，通过样本计算所得的统计量来对总体参数进行判断，如果零假设为真，则可以拒绝原假设，接受对立面（H1），一般以p值的0.05为边界。
pwr包做功效分析
可更改参数：样本大小，显著性水平，功效和效应值
| 函数 | 功效计算的对象 |
|---|---|
| pwr.2p.test() | 两比例（n相等） |
| pwr.2p2n.test() | 两比例（n不相等） |
| pwr.anvoa.test() | 平衡的单因素 ANVOA |
| pwr.chisq.test() | 卡方检验 |
| pwr.f2.test() | 广义线性模型 |
| pwr.p.test() | 比例（单样本） |
| pwr.r.test() | 相关系数 |
| pwr.t.test() | t检验（单样本，两样本，配对） |
| pwr.t2n.test() | t检验（n个不相等的两样本） |
- t检验
- 方差分析
- 线性模型
- 相关性
- 比例检验
- 卡方检验
- 在新情况选择合适的效应值
图形可视化
中级绘图
散点图
# 一般方法
> plot(A, B, main='', xlab='', ylab='') # 散点图
> abline(lm(A~B), col='', lwd='', lty='') # 虚线图
> lines(lowess(A, B), col='', lwd='', lty='') # 直线图
# 改进方法
> library(car)
> scatterplot(A ~ B|C, data=mydata, lwd='', lty='',
main='', xlab='', ylab='',
legend.plot = TRUE,
id.method = "identify",
labels = row.names(matcars)
boxplots='xy'
)
- 散点图矩阵
> library(car)
> scatterplotMattix(~A+B, data=mydata, spread=FALSE, smoother.args=list(lty=''), main='')
- 高密度散点图
> library(hexbin)
> with(mydata, {
bin <- hexbin(x, y, xbins='')
plot(bin, main='')
})
- 三维散点图
> library(scatterplot3d)
> attach(matcars)
> s3d <- scatterplot3d(A, B, C,
pch=16, highlight.3d=TRUE,
type='h',
main='')
> fit <- lm(A~B+C)
> s3d$plane3d(fit)
- 旋转三维散点图
> library(car)
> with(mtcars,
scatter3d(A, B, C))
- 气泡图
> attach(mtcars)
> r <- sqrt(disp/pi)
> symbols(A, B, circle=r, inches=0.30, fg='', bg='', main='', xlab='', ylab='')
> text(A, B, rownames(mtcars), cex=0.6)
> detach(mtcars)
折线图
type格式
| 类型 | 图形外观 |
|---|---|
| p | 只有点 |
| l | 只有线 |
| o | 实心点和线 |
| b,c | 线连接点 |
| s,S | 阶梯线 |
| h | 直方图式的垂直线 |
| n | 不生成任何点和线 |
###下面两种图有些花里胡哨，就不做详解了蛤
相关图
马赛克图
重抽样和自助法
本章描述：是对于数据抽样存在某些缺点，比如混合分布，样本量过小等等因素，此时因用到重抽样和自助法
用coin包做置换检验
imperm包的置换检验
广义线性模型
本章描述：在线性模型中如果因变量为正态分布（甚至连续型变量）不合理,此时应该要用到广义线性模型
广义线性模型和glm()参数
'''用来拟合广义线性模型，形式与lm类似'''
> glm(formula, family=family(link=function), data=)
Logistic回归
环境：通过一系列连续型和/或类别型变量来预测二值型结果变量时，logistics是一个非常有用的工具
> fit <- glm(A~B)
> summary(fit)
# 获得模型系数
> coef(fit)
> exp(coef(fit)) # 指数化系数
过度离势
# 残差偏差/残差自由度
> deviance(fit.reduced)/df.residual(fit.reduced)
# 在模型中直接检验过度离势
> family = "binomial" # 第一次
> family = "quasibinomial" # 第二次
泊松回归
环境：当通过一系列连续型和/或类别型预测变量来预测计数型结果变量时，泊松回归是一个非常有用的工具
> family = poisson()
# 获得模型系数
> coef(fit)
> exp(coef(fit)) # 指数化系数
# 在模型中直接检验过度离势
> family = "binomial" # 第一次
> family = "quasibinomial" # 第二次
主成分分析和因子分析
本章描述：对于信息度过于复杂的多变量数据进行降维压缩
paych包使用
| 函数 | 描述 |
|---|---|
| principle（） | 含多种可选的方差旋转方法的主成分分析 |
| fa（） | 可用主轴，最小残差，加权最小平方或最大似然值估计的因子分析 |
| fa.parallel（） | 含平行分布的碎石图 |
| factor.plot（） | 绘制因子分析或主成分分析的结果 |
| fa.diagram（） | 挥着因子分析或主成分分析的载荷矩阵 |
| scree（） | 因子分析和主成分分析的碎石图 |
主成分分析（PCA）
- 判断主成分的个数
最常见的方法是基于特征值的方法，每个主成分都与相关系数矩阵的特征值相关联，第一主成分和最大的特征值相关联，第二主成分和第二大的特征值相关联。Cattell碎石检验则绘制了特征值与主成分数之间的图像。在图形变化最大处之上的主成分都可以保留—最后看的是横坐标。
第二种是依据初始矩阵相同大小的随机数据矩阵来判断要提取的特征值，若基于真实数据的某个特征值大于一组随机数据矩阵相应的平均特征值，那么该主成分可以保留。
> library(psych)
> fa.parallel(USJudgeRatings[,-1], fa="pc", n.iter=100, show.legend=FALSE, main='') # 去除了CONT这一列
- 提取主成分
> principal(r, nfactors=, rotate=, scores=)
+ r是相关矩阵或原始数据矩阵
+ nfactors设定主成分数（默认为1）
+ rotate指定旋转的方法（默认最大方差旋转（varimax））
+ scores设定是否需要计算主成分得分（默认不需要）
# 输出结果
> PC1栏包含了成分载荷，指的是预测变量与主成分的相关系数，如果不止一个主成分，则会有PC2, PC3
- 主成分旋转
目的：尽可能对成分去噪
类型：正交旋转（变得不相关），斜交旋转（变得相关）
- 获取主成分得分
探索性因子分析（EFA）
- 判断需要提取的公共因子数
> library(psych)
> cov <- ability.cov$cov # ability.cov是相关矩阵
> cor <- cov2cor(cov)
> fa.parallel(cor, n.obs=112, fa='both', n.iter=100, main='')
- 提取公共因子
> fa(r, nfactors=, n.obs=, rotate=, score=, fm=)
+ n.obs是观测数（输入相关数矩阵时需要填写）
+ fm设定因子化方法（默认是极小残差法）
- 可视化结果
> fa.daigram(fa.promax, simple=FALSE)
聚类分析
本章描述：数据归约
层次聚类分析
- 定义：每一个观测值自成一类，这些类两两合并，直到所有的类被聚成一类为止
# 聚类分析的一般步骤
'''
1. 选择合适的变量
2. 缩放数据
3. 寻找异常点
4. 计算距离
5. 选择聚类算法
6. 获得一种或多种聚类方法
7. 确定类的数目
8. 获得最终的聚类解决方案
9. 结果可视化
10. 解读类
11. 验证结果
'''
- 层次聚类方法
| 聚类方法 | 两类之间的距离定义 |
|---|---|
| 单联动（single） | 一个类的点和另一个类中的点的最小距离 |
| 全联动（complete） | 一个类的点和另一个类中的点的最大距离 |
| 平均联动（average） | 一个类中的点和另一个类的点的平均距离 |
| 质心（centroid） | 两类中质心之间的距离，对单个的观测值来说，质心就是变量的值 |
| ward法（ward） | 两个类之间所有变量的方差分析的平方和 |
# 以平均联动聚类为例
> data("nutrient", package="flexclust") # 此处以flexclust包内的nutrient数据集来做
> row.names(nutrient) <- tolower(row.names(nutrient)) # 把索引都变成小写
> nutrient.scaled <- scale(nutrient) # 标准化数据到均值为0，方差为1
> d <- dist(nutrient.scaled ) # 产生距离矩阵
> fit.average <- hclust(d, method='average') # 运用平均联动方法
> plot(fit.average, hang=-1, cex=.8, main='') # 可视化
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-IQ3n8UuP-1602567128520)(E:\R笔记&作业\聚类分析图.png)]
- NbClust包
# 选择聚类个数
# 使用的是NbClust包中的26个评判标准得到的推荐聚类个数
> library(NbClust)
> devAskNewPage(ask = TRUE)
> nc <- NbClust(nutrient.scaled, distance="enclidean", min.nc=2, max.nc=15, method='average')
> table(nc$Beat.n[1,])
> barplot(table(nc$Beat.n[1,]))
# 获取最终的聚类方案
> clusters <- cutree(fit.average, k=5) # 把树状图分成五类
> table(clusters) # 分配情况
> aggregate(nutrient, by=list(cluster=clusters), median) # 描述聚类，获取中位数
> plot(fit.average, hang=-1, cex=.8, main='')
> rect.hclust(fit.average, k=5) # 叠加五类的解决方案
划分聚类分析
定义：首先定义类的K值，然后观测值被随机分成K类，再重新形成聚合的类
- K均值聚类
原理：把每个数据点分配到离他最近的中心点，重新计算每类中的点到该类中心点距离的平均值，分配每个数据到他最近的中心点（对于异常值很敏感）
# 确定聚类分布个数
> data(wine, package='rattle')
> df <- scale(wine[-1])
> wssplot(df)
> library(NbClust)
> set.seed(1234)
> devAskNewPage(ask = TRUE)
> nc <- NbClust(df, min.nc=2, max.nc=15, method='kmeans')
> table(nc$Beat.n[1,])
> barplot(table(nc$Beat.n[1,]), xlab='', ylab='', main='')
- 围绕中心点的分布
原理：随机选择K个观测值，计算观测值到各个中心的距离/相异性，把每个观测值分配到最近的中心点，计算每个中心点到每个观测值距离的总和，选择一个个该类中不是中心的点，并和中心点互换，重新把每个点分配到距他最近的中心点，再次计算总成本，如果总成本比前一次少，把新的点作为中心点
# 对葡萄酒数据使用基于质心的划分方法
> library(cluster)
> set.seed(1234)
> fit.pam <- pam(wine[-1], k=3, stand=TRUE)
> fit.pam$medoids
> clusplot(fit.pam, mian='')
- 避免不存在的类
# 发现不存在的类
> library(fMultivar)
> set.seed(1234)
> df <- rnorm(1000, rho=.5) # 从相关系数为0.5的二元正态分布抽取1000个观测值
> df <- as.data.frame(df)
> plot(df, main='')
# 确定聚类个数
> wssplot(df)
> library(NbClust)
> nc <- NbClust(df, min.nc=2, max.nc=15, method='kmeans')
> dev.new()
> barplot(table(nc$Beat.n[1,]), xlab='', ylab='', main='')
分类
数据准备（划分数据集）
# 对原始数据进行处理
# 假设第一列是类别型向量
> df <- breast[-1]
> df$class <- factor(df$class, levels=c(2,4), labels=c('A', 'B'))
> ste.seed(1234)
> train <- sample(nrow(df), 0.7*nrow(df))
> df.train <- df[train,]
>df.validate <- df[-train,]
原理： 是一种广义线性模型，可以根据一组数值变量预测二元输出
> fit.logit <- glm(class~.,data=mydata, family=binomial()) # 拟合逻辑回归
> summary(fit.logit) #检查模型
> prob <- predict(fit.logit, data.validata, type='response')
> logit.pred <- factor(prob > .5, levels=c(FALSE,TRUE), labels=c('a', 'b')) # 对训练集集外样本单元分类
> logit.perf <- table(df.validata$class, logit.pred, dnn=c('Actual', 'Predicted')) # 评估预测准确性
> logit.perd
决策树
- 经典决策树
> library(rpart)
> set.seed(1234)
> dtree <- rpart(class~., data=df.train, method='class', parms=list(split="information") ) # 生成树
> plotcp(dtree) # 画出交叉验证误差与复杂度参数的关系图
> dtree.pruned <- prune(detree, cp=.0125) # 剪枝
> library(rpart.plot)
> prp(dtree.pruned, type=2, extra=104, fallen.leaves=TRUE, main='') # 可视化
> dtree.pred <- predict(dtree data.validata, type='class')
> dtree.perf <- table(df.validata$class, dtree.pred, dnn=c('Actual', 'Predicted')) # 对训练集外样本单元分类
> dtree.perf
- 条件推断树
# 是在显著性检验后才可用
> library(party)
> fit.ctree <- ctree(class~., data=df.train)
> plot(fit.ctree, main='')
> ctree.pred <- predict(fit.ctree data.validata, type='response')
> ctree.perf <- table(df.validata$class, logit.pred, dnn=c('Actual', 'Predicted')) # 评估预测准确性
> ctree.perf
随机森林
原理：在随机森林中，会生成多个预测模型，并将模型的结果汇总以提升分类准确率
> library(randomForest)
> set,seed(1234)
> fit.forest <- randomForest(class~., data=df.train, na.action=na.roughfix, importance=TRUE)
> fit.forest # 生成森林
> importance(fit.forest, type=2) # 给出变量重要性
> forest.pred <- predict(fit.forest, df.validate)
> forest.perf <- table(df.validate$class, forest.pred, dnn=c('A', 'B'))
> forest.perf # 对训练集外样本点分类
支持向量机—用量少
原理：SVM旨在在多维空间中找到一个能将全部样本单元分成两类的最优平面，这一平面应使两类中距离最近的间距（margin）尽可能大，在间距边界上的点被称为支持向量（support vector他们决定间距），分割的超平面位于间距的中间。
> library(e1071)
> set.seed(1234)
> fit.svm <- svm(class~., data=df.train)
> fit.svm
> svm.pred <- predict(fit.svm, na.omit(df.validate))
> svm.perf <- table(na.omit(df.validate)$class, svm.pred, dnn=c('A', 'B'))
> svm.perf
'''由于方差较大的预测变量通常是对SVM的生成影响较大，svm函数默认在生成模型前对每个变量标准化，使其均值为0，标准差为1，与随机森林算法不同的是，SVM在预测新样本的时不允许有缺失值出现'''
选择调和参数—带RBF核的SVM模型
> set.seed(1234)
> tuned <- tune.svm(class~., data=df.train, gamma=10^(-6, 1), cost=10^(-10, 10)) # 变换参数
> tuned # 输出最优模型
> fit.svm <- svm(class~., data=df.train, gamma =.01, cost=1) # 用这些参数拟合模型
> fit.svm
> svm.pred <- predict(fit.svm, na.omit(df.validate)) # 评估交叉验证表现
> svm.perf <- table(na.omit(df.validate)$class, svm.pred, dnn=c('A', 'B'))
> svm.perf
选择预测结果最好的解
预测准确性度量
| 统计量 | 解释 |
|---|---|
| 敏感度 | 正类的样本单元被成功预测的概率，也叫召回率或正例覆盖率 |
| 特异性 | 负类的样本单元被成功预测的概率，也叫负例覆盖率 |
| 正例命中率 | 被预测为正类的样本单元中，预测正确的样本单元占比，也叫精确度（precision） |
| 负例命中率 | 被预测为负类的样本单元中，预测正确的样本单元占比 |
| 准确率 | 被分类正确的样本单元所占比重，也叫ACC |
- 评估二分类准确性
> performance <- function(table, n=2){
if(!all(dim(table) == c(2, 2)))
stop('must be a 2 x 2 table')
tn = table[1, 1]
fp = table[1, 2]
fn = table[2, 1]
tp = table[2, 2]
sensitivity = tp/(tp+fn)
specificity = tn/(tn+fp)
ppp = tp/(tp+fp)
npp = tn/(tn+fn)
hitrate = (tp+tn)/(tp+tn+fp+fn)
result <- paste(
"\nsensitivity = ", round(sensitivity, n),
"\nspecificity = ", round(specificity, n),
"\npositive predictive value = ", round(ppp, n),
"\nnegative predictive value = ", round(npp, n),
"\naccuracy = ", round(hitrate, n), "\n",sep=""
)
cat(result)
+ }
- 验证
> performance(logit.perf) # 逻辑回归
> performance(dtree.perf) # 经典决策树
> performance(ctree.perf) # 条件决策树
> performance(forest.perf) # 随机森林
> performance(svm.perf) # SVM支持向量机
###在此拓展介绍一种高级绘图法
ggplot2高级绘图
| 函数 | 添加 | 选项 |
|---|---|---|
| geom_bar() | 条形图 | color，fill，alpha |
| geom_boxplot() | 箱线图 | color，fill，alpha |
| geom_density() | 密度图 | color，fill，alpha |
| gemo_histogram() | 直方图 | color，fill，alpha |
| gemo_hline() | 水平图 | color，alpha |
| gemo_jitter() | 抖动图 | color，size，alpha，shape |
| gemo_line() | 线图 | colorvalpha，linetype，size |
| gemo_point() | 散点图 | color，alpha |
| gemo_rug() | 地毯图 | color，side |
| gemo_smooth() | 拟合曲线 | method，formula，color，fill，linetype，size |
| geom_text() | 文字注解 | 参见帮助 |
| geom_violin() | 小提琴图 | color，fill，alpha |
| gemo_vline() | 垂线 | color，alpha |
后续会对一些地方更新的。
大家加油！！！！
扫一扫
5万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
