R语言入门与数据分析
本文转自Git-hub上一位大佬的。
本文仅作学习使用、不做任何商业用途
Creater: Fangyi Li
此笔记基于哔哩哔哩弹幕网视频：R语言入门与数据分析（up主id：基因学院）总结而成。
This note is created based on a video, Introduction to R programming language
and data analysis, on the website of Bilibili (id: 基因学院).
文章目录
- R语言入门与数据分析
- P1-6 Introduction of R
- P7-14 R的安装
- P15 R内置数据集
- P16-26 数据结构（R储存组织数据的方式）
- P28 获取数据
- P29-30 读入文件
- P31 写入文件
- P32 读写Excel文件
- P33 读写R格式文件
- P34-40 数据格式转换
- P 40 数据处理
- P41 R函数
- P42 选项参数
- P43 数学统计函数
- P44 描述性统计函数
- P45 频数统计函数
- P46 独立性检验函数
- P47 相关性分析函数
- P48 相关性检验函数
- P49 绘图函数
- P50 自定义函数
- P51 数据分析实战
- P52-53 线性回归（regression）
- P54 多元线性回归
- P55 回归诊断
- P56-57 方差分析（Analysis of Variance）
- P58 功效分析（Power Analysis）
- P59 广义线性模型
- P60 Logistic回归
- P61 主成分分析
- P62 因子分析
- P63 购物篮分析
- P64-65 后记
P1-6 Introduction of R
数据采集
数据挖掘（Data mining）
-
三个特点
-
数据整体性
-
不追求精确性
-
因果关系（x）
-
相关性（√）
-
通过算法发现新东西（规律）-> 相同数据不同结果 -> 结果可视化（数据去抽象化）
决策
R语言特点
-
编程开源（大量拓展包）
-
绘图功能
-
面向对象
-
不够规范、不支持多线程
-
函数后必须接括号
getwd()
# 看路径
df <- read.table("dataname.cvs", sep=",", header = T)
# df将数据文件命名
# read.table()读写Excel数据文件
# sep=","来表示空格用字符
table(df$supp, df$dose)
# table()用于制作新表格，$用于统计总量
# aggregate()测量均值、方差等
P7-14 R的安装
R的安装
# - 用 "<-" 赋值给一个局部变量
# - 用 "<<-" 赋值给一个全局变量
# - "<-"的快捷键：alt + -
mean()
# mean(1:5) -> 3 求1-5的平均数
# mean(1,2,3,4,5) -> 1 只输出第一个向量（1）的平均数，后面为无效输入
# mean(c(1,2,3,4,5)) -> 3 增加的c()帮助识别一串数字
ls()
# 列出所有变量
ls.str()
# 列出所有变量的详细信息
str()
# 括号里填写变量名，输出变量的信息
rm()
# 清除某个变量
rm(list = ls())
# 消除所有变量
history()
# 唤出历史
save.image()
# 保存工作空间（图片不会自动保存）
R拓展包安装
在线
# 通过访问镜像站点
# 拓展包的名字都需要被引号
install.package()
# eg. install.package("vcd")
install.packages()
# eg. install.packages(c("a", "b", "c"))
离线（源代码）
# 镜像站点下载源代码
# 需要注意包之间的依赖关系
require()
# 加载包（eg. require(vcd)）
detach() # eg. detach("package:vcd")
remove.packages() # eg. remove.packages("vcd")
# 去除包
installed.packages()
# 可用于拓展包移植，比如：
1.老设备
Rpack <- installed.packages() [,1] # Rpack赋值所有拓展包的名字
save(Rpack, file = "Rpack.Rdata") # 保存文件
2.新设备
for (i in Rpact) install.packages(i)
# ?函数名（一个问号查函数）eg. ?mean
# ??函数名（两个问好查拓展包）eg. ??heatmap
apropos(关键字)
# 查询括号内的关键字
RSiteSearch("关键字")
# 在R网址查询关键字 rseek.org
# R -> 大数据 -> 批量化、自动化
# Excel -> 鼠标操作方便
P15 R内置数据集
# - 存在datasets包中
# - data()输出所有包含的数据集 或 data(package = "包名")
# - 如果数据里面的某个变量被重新赋值，可以使用data("变量名")来赋值回数据集里的值
# 向量 -> 名称
# 因子 -> 分类
# 矩阵、数组 -> 不同东西的不同类型
# 类矩阵、数据框
data.frame()
# 制作数据框
data(数据集名称, package = "包名")
# 只加载拓展包里面的数据集，忽略包
P16-26 数据结构（R储存组织数据的方式）
1.向量（vector）
# - 用函数c()来创建向量（数组）
# - 逻辑型向量：TRUE(T) OR FALSE(F)
# - 输出等差数列
seq(from = 起始数, to = 结尾数, by = 间隔, length.out = 输出数量)
# - 重复输出
rep(数组, 次数)
# - 一个数组里面的向量必须是同一类型
> c(1, 2, "three")
"1" "2" "three"
# 导致前面的数字向量全部变成字符串
# - 向量化编程可以省略循环，比如：
> x <- c(1, 2, 3, 4, 5)
> y <- c(6, 7, 8, 9, 10)
> x * z + y
8 11 14 17 20
# - 使用c()来控制循环的次数
> x <- c(1, 2, 3, 4, 5)
> rep(x, c(2, 4, 1, 1, 3))
1 1 2 2 2 2 3 4 5 5 5
# - 向量索引
length(数组) # 数组数量
数组[位置] # 输出数组某个向量
数组[-位置] # 输出除了选中位置其余的所有向量
append() # 使用append不会mutate用来的向量
# - 向量运算
x ** y # x^y
x %% y # x/y的余数
x %/% y # x/y的整数部分（整除）
# 循环：
# 向量个数不同也可能可以相加，但长的向量个数必须是短的向量个数的整除数。
# 想要对比两个对象是否相等，应用两个等号==（一个等号相当于赋值）
abs() # 绝对值
sqrt() # 平方根
log(向量, base = log的底数) # 不加base默认底数是e
log10() # base是10的log
exp() # 以e为底的指数
ceiling() # 向上取整
floor() # 向下取整
trunc() # 向量内取整
# 四舍五入
round(向量, digits = 保留小数点位数)
signif(向量, digits = 总保留位数)
# 三角函数
sin()
cos()
tan()
# 统计函数
sum() # 求和
max() # 最大值
min() # 最小值
range() # 最大、小值
mean() # 均值
var() # 方差（variance）
sd() # 标准差（standard deviation）
prod() # 连乘的积
median() # 中位数
quantile(向量, 百分比combine) # 默认四分位数
# 索引值
which.max() # 向量中最大值的位置
which.min() # 向量中最小值的位置
which(向量==目标值) # 目标值在向量中的位置
2.矩阵与数组（matrix）
# （多维的向量）
matrix(向量, nrow = 行数, ncol = 列数, byrow = T/F) # 制作一个矩阵
# 向量的个数必须大于或等于matrix的格数（行数 x 列数 <= 向量数）
# Data length must be a sub-multiple or multiple of the number of rows
# byrow = T 时matrix先填行，再填列，反之亦然
dimnames(matrix) <- list(行名集, 列名集)
# 如何加行、列名
dim(向量) # 输出向量的维度，也可直接赋值来改变向量的维度
# > dim(x) # x = c(1, 2)
# NULL
# > dim(x) <- c(1, 2)
# > x
# [,1] [,2]
# [,1] 1 2
# 数组（更多维）
array(向量, 维数集, dimnames = 维度名的list)
# eg. array(1:6, c(1, 2, 3)指1行2列3层, dimnames = list(c("A1"),
# c("B1", "B2"),
# c("C1", "C2", "C3")))
# 矩阵索引
matrix[行数, 列数] # 列数和行数不一定只是一个数，可以是个集
matrix[行数] # 访问行数的第一列
matrix[行数,] # 访问行数的全列
matrix[-行数, 列数] # 访问除了制定行数的全部行数
matrix["行名", "列名"] # 如果matrix已有行、列名，也可直接用行、列名索引
# 矩阵运算
# 行列不一致的两个矩阵无法相加减
rowSums(matrix) # matrix每一行的和
colSums(matrix) # matrix每一列的和
rowMeans() # matrix每一行的平均数
colMeans() # matrix每一列的平均数
matrix1 * matrix2 # 内积
matrix1 %*% matrix2 # 外积
diag(matrix) # matrix的对角线值
t(matrix) # matrix的行列互换
3.列表（list）
list(第一对象名 = 第一个对象, 第二对象名 = 第二个对象, ...)
# 制作列表
# 访问单个对象
列表名[对象位置]
列表名["对象名"] # 这里访问需要加引号
# 访问多个对象
列表名[c(对象1位置, 对象2位置, ...)]
列表名[c("对象1名", "对象2名", ...)]
列表名$对象名 # 一般对象名会自动给
# ***mlist[1]和mlist[[1]]的区别***
mlist[1]
# 类型为list
# 要删减列表可以使用单括号负索引，如：mlist[-2]等于去除第二个list
mlist[[1]]
# 类型为一开始组成列表的元素，如vector、matrix等（也可能是列表）
# 所以要改或加mlist里面的元素，需要两个括号
# 删减列表可以把list里面的元素赋值NULL，如mlist[[1]] <- NULL
# 删减都会造成所有的数据在列表中的位置变化
4.数据框（数据集）
# - 行（row）表示观测（object）
# - 列（column）表示变量（dose）
# col1 col2 col3
# row1
# row2
# row3
# 矩阵必须为同一数据类型
# 数据框每一列必须同一类型，但每一行可以不同。
data.frame(变量1, 变量2, ...) # 制作数据框
# 访问数据框
数据集[列/"列名"] # 可以用c()一步查看多列，会输出一个带行的新数据框
数据集["行名",] # 一行 + 全列
数据集[, "列名"] # 全行 + 一列
数据集$列名 # 全行 + 一列
with(数据集, {列名}) # 直接输出全列
attach(数据集) # 把数据集的列名导入R，之后可以直接访问
detach(数据集) # 消除数据集的列名
5.因子（factor）
# 变量类型：
# - 名义型（nominal）：无顺序的/项目之间是独立的/比如名字
# - 有序型（ordinal）：项目之间有关联但不是连续的数量变化/比如poor, better, best
# - 连续型（continuous）：有顺序的/连续的数量变化/比如金额、增产率等
# 名义型（nominal）和有序型（ordinal）都属于因子
# 水平值（level）构成的向量为因子
# 因子的作用：计算概率和频数
# factors Lever 1 Level 2 Level 3
# * $ &
# * * * -> **** $$ &&&
# $ & &
# 因子型数据
factor(c(数据集), orderred = T/F, levels = c("level1", "level2", ...))
# 会输出数据集和levels的顺序
# 绘图（因子输出条形图????，向量输出散点图）
# 散点图
# ^ level
# |
# | ...
# | ...
# | ...
# |
# ----------------> index
# 条形图
# ^ 数量
# |
# | ---
# | --- | |
# | --- | | | |
# | | | | | | |
# ----------------> level
# cut(被分割的向量, (把0到100, 总level数量为10, 等分10份))
cut(1:100, c(seq(0, 100, 10)))
6.缺失数据
# - 测量值丢失
# - 没有值
# - 无效值
# NA表示缺失值（NA不一定等于零）
1 + NA
NA
NA == 0
NA
# 函数中关于缺失值的参数
na.rm = T/F # 是否移除缺失值（注：选了T后，函数计算相当于缺失值从未出现）
# 消除向量中的NA
na.omit(原向量) # 如果是数据集使用na.omit()，会把有NA的整行都删除掉
#其他缺失数据
NaN # 代表不可能的值（比如0/0），不存在的
is.nan() # 查询是否为NaN
Inf # 代表正无穷（1/0）
-Inf # 代表负无穷（-1/0）
is.infinite() # 查询是否为无穷
Na # 是存在的，但不知道是多少
7.字符串（string/str）
# 要加引号！
nchar() # 向量中每个元素（字符串）的长度
length() # 向量中元素的个数
paste()
# > paste(c("a", "b", "c"))
# "a" "b" "c" # 因为c()，三个向量不在一个维度，分别处理
# > paste("a", "b", "c")
# "a b c" # 默认是空格分割
# > paste("a", "b", "c", sep="-")
# "a-b-c" # sep为分割参数
substr(向量, start = 第一个字符, stop = 最后的字符) # 用于提取部分字符
toupper() # 全大写
tolower() # 全小写
# 用于分割字符串，返回的是一个列表不是向量
strsplit(字符串向量/集, "分隔符")
# 对向量1和向量2使用制定函数
outer(向量1, 向量2, FUN = 函数名)
8.日期和时间（date）
# 对时间序列的描述
# 利用前面的结果进行预测
"ts" # time series 时间数据的简称
# ts不是数据框
Sys.Date() # 查看系统时间
as.Date(向量, format = "%Y-%m-%d")
# 例子
a <- "2020-01-01"
class(a) # "character"
class(as.Date(a)) # "Date"
?strftime # 时间函数相关介绍
seq(as.Date("2017-01-01"), as.Date("2017-07-05"), by=5)
# by参数为相隔天数（从起始开始，每隔by天出一个数）
ts(data,
start = c(起始年, 起始年的第几个频率单位),
end = c(终止年, 终止年的第几个频率单位),
frequency = 频率数)
# Creating time-series objects
# frequency频率变化：
# 1 以年为单位
# 12 以月份为单位
# 4 以季度为单位
# 例子：
> sales <- round(runif(48, min=50, max=100))
> sales
[1] 56 75 75 85 51 73 63 79 73 88 57 81 72 85 58 60 87 77 55 86 65 71 50 91 64 63 59
[28] 60 64 82 75 81 60 75 95 50 75 87 79 85 57 50 62 86 52 72 98 97
> ts(sales, start = c(2010,5), end = c(2014,4), frequency = 1) #以年为单位
Time Series:
Start = 2014
End = 2017
Frequency = 1
[1] 56 75 75 85
# 因为frequency=1，所以起始为从2010开始算的第五年，即2014年。
# 结束为2014开始算的第四年，即2017年。
# 2014-2017，总4年，可以有四个数据指向时间。
# 因为可以fill的时间不够，多余的数据就被省略掉了。
ts(sales, start = c(2010,5), end = c(2014,4), frequency = 4) #以季度为单位
# Qtr1 Qtr2 Qtr3 Qtr4
# 2011 56 75 75 85
# 2012 51 73 63 79
# 2013 73 88 57 81
# 2014 72 85 58 60
# 2010开始的第五季度，即2010第一季度；2014开始的第四季度，即2014第四季度。
ts(sales, start = c(2010,5), end = c(2014,4), frequency = 12) #以月份为单位
ts
# Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
# 2010 56 75 75 85 51 73 63 79
# 2011 73 88 57 81 72 85 58 60 87 77 55 86
# 2012 65 71 50 91 64 63 59 60 64 82 75 81
# 2013 60 75 95 50 75 87 79 85 57 50 62 86
# 2014 52 72 98 97
P28 获取数据
键盘输入
# Creating a dataframe
patientID <- c(1, 2, 3, 4)
admdate <- c("10/15/2009", "11/01/2009", "10/21/2009", "10/28/20009")
age <- c(25, 34, 28, 52)
diabetes <- c("Type1", "Type2", "Type1", "Type1")
status <- c("Poor", "Improved", "Excellent", "Poor")
data <- data.frame(patientID, age, diabetes, status) # 数据框
data
# patientID age diabetes status
# 1 1 25 Type1 Poor
# 2 2 34 Type2 Improved
# 3 3 28 Type1 Excellent
# 4 4 52 Type1 Poor
data2 <- data.frame(patientID = character(0),
age = numeric(0),
diabetes = character(0),
status = character(0))
data2 <- edit(data2) # 可视编辑器
fix(data2)
读取外部文件
访问数据库
# ODBC = open database connectivity 开放数据库
install.packages("RODBC")
# 双向通信，读写和写入
P29-30 读入文件
?read.table()
x <- read.table("input.txt") # 默认空白分割
x <- read.table("input.csv", sep=",", header = T, skip = 5, nrows = 100,
na.strings = " ")
# 逗号分割
# 第一行是否为变量名称
# skip的行数
# nrows：总共导入多少行
# 自动转化缺失值
head(x, n = 6) # 提取最开始的六行（默认六行）
tail(x, n = 6) # 提取最末尾的六行
read.csv() # 读取csv文件
read.fwf() # fixed width file
x <- read.table("website", header = T)
install.packages("XML") # 网络读取数据的包
?readHTMLTable
install.packages("foreign") # 读取不同软件格式的数据的包
help(package = "foreign")
# Excel直接复制到剪贴板导入
x <- read.table("clipboard", header = T, sep = "\t") # Windows
x <- read.table(pipe("pbpaste")) # Mac
read.table(gzfile("input.txt.gz")) # gz file
?readLines
readLines("input.csv", n = 5) # n = 返回行数
?scan
scan("input.txt", what = list(numeric(0)))
scan("input.txt", what = list(X1 = character(3),
X2 = numeric(0),
X3 = numeric(0)))
# 数字表示读入的数量，0表示所有
P31 写入文件
?write
write.table(x, file = "x.txt")
write.table(x, file = "x.csv", sep = "\t", quote = F, append = F,
na = "NA", row.names = F)
# append = F清空源文件（T加在原文件后）
# row.names取消自动添加的行号
write.csv
P32 读写Excel文件
# Excel文件另存为csv格式
x <- read.csv()
# 复制剪贴板
x <- read.table(pipe("pbpaste"))
# Excel的格式：.xlsx
# XLConnect包，需要JAVA环境
install.packages("XLConnect") # 需要rJava包
x <- loadWorkbook("data.xlsx")
readWorksheet(x, 1) # 不同的工作表，数字为工作表位置
# 写入（四步法）
y <- loadWorkbook("file.xlsx", create = T)
createSheet(y, "Sheet 1")
writeWorksheet(y, data = x, sheet = "Sheet 1")
saveWorkbook(y)
# 写入（一步法）
writeWorksheetToFile("file.xlsx", data = x, sheet = "Sheet 1")
# xlsx包，需要JAVA环境
read.xlsx("data.xlsx", 1, startRow = 1, endRow = 100)
write.xlsx(x, file = "data.xlsx", sheetName = "Sheet 1", append = T/F)
# openxlsx包，不用JAVA环境
P33 读写R格式文件
# .RDS -> R data set
saveRDS(data, file = "data.RDS") # 储存为RDS文件至公共目录
x <- readRDS("data.RDS") # 读入RDS文件
# .R
# .RData
load("data.RData")
save(data1, data2, file = "data.RData")
save.image()
P34-40 数据格式转换
添加、删减、排序
数据框
is.data.frame(data)
new.data <- as.data.frame(data) # 其他数据类型转换为数据框
method(is)
method(as)
unlist(list)
数据取子集
who <- read.csc("WHO.csv", header = T) # 导入数据
who1 <- who[c(1:50), c(1:10)] # 取前50行，10列（想要的那部分）
who$Continent # 使用列名
who2 <- who[which(who$Continent == 7)] # 使用which函数来定位（逻辑判断）
who3 <- who[which(who$CountryID > 50 & who$CountryID <= 100)] # 多重逻辑
?subset
who4 <- subset(who, who$CountryID > 50 & who$CountryID <= 100) # data, logic
?sample # 随机提取部分数据
x <- 1:100
sample(x, n, replace = T/F) # n为输出数量，replace为是否可以重复
who[sample(who$CountryID, 30, replace = F),] # 随机抽取CountryID 30个子集
删除固定行
data[-1:-5] # 负索引
data$row <- NULL # 赋值NULL
数据框合并
合并列
data.frame(d1, d2) # 创造一个新的数据框
?cbind
cbind(d1, d2) # 在第一个数据框下改
合并行（d1和d2列名得相同）
?rbind
rbind(d1, d2)
# 矩阵也可以用cbind和rbind
去除重复项
取子集
duplicated(data)
data[duplicated(data),] # 输出重复部分
data[!duplicated(data),] # 输出非重复部分
unique(data) # 直接去除重复数据
翻转数据框
# Excel转置选项
# R
new.data <- t(data) # 行列翻转
rev() # 用于向量
data[rev(rownames(data)),] # 单翻转行，比如1->10变成10->1
列运算
new.data <- data.frame(c1 = datac1 * 2, c2 = datac2) # 组成一个新的数据框
transform(data, c1 = c1 * 2) # 直接修改原数据
transform(data, c3 = c1 * 2) # 在原数据上加一列新数据
数据框的排序
sort(data) # 数字从小到大，字符从a到z（只能用于向量，不能用于数据框）
data[sort(rownames(data)),] # 用sort来排序数据框的方法（输出新数据框）
order() # 返回索引值
data[order(data$target),] # 按照target列从小到大排列
data[order(-data$target),] # 按照target列从大到小排列
data[order(data$t1, data$t2),] # 两个排序条件
rank()
数据框数学运算
# for dataframe, matrix
apply(data, MARGIN = 1, FUN = sum()) # 行
apply(data, MARGIN = 2, FUN = mean()) # 列
lapply(list, FUN = length()) # for list
sapply() # for vector, matrix
tapply(data1, data2, FUN = length()) # for factor
apply系列函数->
分组计算
tapply() # 参数：vector / 返回值：vector
apply() # 参数：list, data.frame, array / 返回值：vector, matrix
多参数计算
mapply() # 参数：vector，不限个数 / 返回值：vector, matrix
循环迭代
lapply() # 参数：list, data.frame / 返回值: list
# 简化版
sapply() # 参数：list, data.frame / 返回值：vector, matrix
# 可设置返回值
vapply() # 参数：list, data.frame / 返回值：vector, matrix
# 递归版
rapply() # 参数：list / 返回值：list
环境空间遍历
eapply() # 参数：environment / 返回值：list
Centering and standardization 中心化与标准化
# 中心化（方差）：data set中的各项数据减去data set的均值
# 标准化（标准差）：中心化之后在除以data set的标准差，即data set中的各项数据减去data set的均值再除以data set的标准差
# 使不同类型数据大小接近
# 先中心化，数据差还是大的话再标准化
# 中心化：
x - mean(x)
# 标准化：
(x - mean(x)) / sd(x)
scale(x, center = T/F, scale = T/F)
# 参数center = T为中心化处理
# 参数scale = T为标准化处理
# 中心化和标准化后，绘制heat map区别会更加明显
reshape2包
x <- data.frame(k1 = c(NA,NA,3,4,5), k2 = c(1,NA,NA,4,5), data = 1:5)
y <- data.frame(k1 = c(NA,2,NA,4,5), k2 = c(NA,NA,3,4,5), data = 1:5)
merge(x, y, by = "k1")
merge(x, y, by = c("k1", "k2"))
install.packages("reshape2")
melt(data) # 可以把宽数据变成长数据，列名变为因子
melt(data, id.vars = ) # id.vars 是观测排序的参数
dcast(data, month+day ~variable) # 重铸数据，以month和day列为排序
dcast(data, month ~ variable, fun.aggregate = sum, na.rm = T)
# 单以month为排列参数会有很多行有相同优先值
# 所以fun.aggregate的函数用来把有着相同month的数据运行一个设定的函数
tidyr包（Tidy data）
install.packages("tidyr")
# 不能有两个相同行名和列名的unit
# 宽变长
gather(data, key = "Key", value = "Value", factor1, factor2, facotr3)
# factor其实就是列，把每一列变成一种因子
gather(data, key = "Key", value = "Value", 2:4) # 可以直接输入列的编号
# 长变宽
spread(data, key = "Key", value = "Value") # gather的反函数
# 会自动把key里面不同的因子变为不同的列名
df <- data.frame(x = c(NA, "a.b", "b.c", "a.d"))
separate(df, col = x, into = c("c1", "c2")) # 会自动识别分割符
df <- data.frame(x = c(NA, "a.b-c", "b-c", "a-d"), sep="-") # 有多个分割符使需要标明
unite(x, col = "united col", c1, c2, sep="-")
# separate的反函数，把两列连起来，需表明连置符
P 40 数据处理
dplyr包
install.packages("dplyr")
dplyr::filter(data, c2 > 10) # filter掉不符合后面逻辑的行
dplyr::distinct(rbind(data[1:10,], data[1:15,])) # distinct去掉重复的十行
dplyr::slice(data, 10:15) # 取出制定行
dplyr::sample_n(data, n) # 随机取出n行
dplyr::sample_frac(data, 0.1) # 按几率取出行
dplyr::arrange(data, c2) # 按照某一列的大小来排序（从小到大）
dplyr::arrange(data, desc(c2)) # 按照某一列相反的大小来排序（从大到小）
select()
统计函数
ummarise(data, avg = mean(c2))
summarise(data, sum = sum(c2))
# 链式操作符 %>%（管道）
# 用于实现将一个函数的输出传递给下一个函数，作为下一个函数的输入。
# 快捷键：ctrl + shift + M
# 相当于“然后”
# 例子：
head(data, 20) %>% tail(10)
# 先取出data的前二十行，再从这二十行里面取出最后的十行，所以最终输出为11-20行。
dplyr::group_by(data, c5) # 通过制定列的因子来进行分组，每一种因子为一组。
data %>% group_by(c5) # 另一种写法
data %>% group_by(c5) %>% summarise(avg = mean(c2)) %>% arrange((avg))
# 先分组，再算平均值，再排序
dplyr::mutate(data, new=c2+c3) # 加新的一列，并添加逻辑
双表格整合
a <- data.frame(x1 = c("A","B","C"), x2 = c(1,2,3))
b <- data.frame(x1 = c("A","B","D"), x2 = c(T,F,T))
dplyr::left_join(a, b, by="x1") # x1以第一个数据框为准
# x1 x2.x x2.y
# 1 A 1 TRUE
# 2 B 2 FALSE
# 3 C 3 NA
dplyr::right_join(a, b, by="x1") # x1以第二个数据框为准
# x1 x2.x x2.y
# 1 A 1 TRUE
# 2 B 2 FALSE
# 3 D NA TRUE # 因为a的x1没有D，所以x2.x是NA
dplyr::full_join(a, b, by="x1") # 把所有可能都列出来了
# x1 x2.x x2.y
# 1 A 1 TRUE
# 2 B 2 FALSE
# 3 C 3 NA
# 4 D NA TRUE
dplyr::semi_join(a, b, by="x1") # 列出a，b都有的可能
# x1 x2
# 1 A 1
# 2 B 2
dplyr::anti_join(a, b, by="x1") # 列出a，b的不同行
# x1 x2
# 1 C 3
intersect(data1, data2) # 取交集
dplyr::union_all(data1, data2) # 取并集
dplyr::union(data1, data2) # 取非冗余的并集
setdiff(data1, data2) # 取data1的补集
P41 R函数
# 例子：
state <- as.data.frame(state.x77, c("Murder", "Population", "Illiteracy", "Income", "Frost"))
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = state)
summary(fit)
# 函数的返回值
# 一些有返回，一些没有
# 绘图函数输出图形
# 向量
sum()
mean()
sd()
range()
median()
sort()
order()
# 矩阵或数据框
rbind()
cbind()
# 数字矩阵
heatmap()
# 可以使用help()来看不懂的函数
P42 选项参数
ls("package:base")
par()
# 规律
# 1. 输入控制部分
# 可以输出哪种类型的数据
# file：接一个文件
# data：一般指要输入一个数据框
# x：一个单独的对象，一般是向量，但也可能是矩阵或者列表
# x和y：函数需要两个输入变量
# x, y, z：函数需要三个输入变量
# formula：公式
# na.rm：是否删除缺失值
# ...表示参数可传递，或没有数量限制
# ~ 表示相关
# 2. 输出控制部分
# 3. 调节部分
# （1）根据名字判断选项的作用：
# color选项：控制颜色
# select 与选择有关
# font 与字体有关
# font.axis 坐标轴的字体
# lty = line type
# lwd = line width
# method 软件算法
# （2）选项接受哪些参数：
# main：字符串，不能是向量
# na.rm：TRUE or FALSE
# axis：side参数只能是1到4
# fig：包含四个元素的向量
# 不取值可能是NULL或者NA
# method参数与算法相关
P43 数学统计函数
# 概率函数
# R概率分布
# d 概率密度函数（density）
# p 分布函数（distribution）
# q 分布函数的反函数（quantum？）
# r 产生相同分布的随机数（random）
# 正态分布
dnorm()
pnorm()
qnorm()
rnorm()
# 离散分布
x <- rnorm(n=100, mean=15, sd=2)
qqnorm(x)
# 生产随机数
runif(n) # 生产n个0-1的随机数
runif(n, min = 1, max = 100) # min/max参数设置范围
round(runif(n, min = 1, max = 100)) # 使用round生成整数
dgamma()
# 随机数种子
set.seed(n) # 每次运行可以使之后的随机数相同，n为种子编码
P44 描述性统计函数
# summary()提供向量的各项数据，比如最大小值、平均值、四分位数、因子频数向量统计等
myvars <- mtcars[c("mpg", "hp", "wt", "am")]
summary(myvars)
fivenum(myvars$hp)
# 返回五个数：最小值、四分位数第一、中位数、四分位数第三、最大值
# Descriptive stats via describe (Hmisc)
install.packages("Hmisc")
library(Hmisc)
myvars <- c("mpg", "hp", "wt")
describe(mtcars{myvars})
# Descriptive stats via stat.desc (pastecs)
install.packages("pastecs")
library(pastecs)
stat.desc(mtcars{myvars})
stat.desc(mtcars, basic = T, desc = T, norm = T)
# basic = T 计算基本值
# desc = T 计算描述值
# norm = T 计算统计值
# Descriptive stats via describe (psych)
install.packages("psych")
library(psych)
describe(mtcars{myvars}, trim = 0.1)
# trim = 0.1 去除最低10%和最高10%的数据
Hmisc::describe(mtcars)
# 函数顺序：后入为主（后面载入的包函数会覆盖前面的包函数）
# 用::来表明函数来源
# Descriptive stats by group with aggregate
library(MASS)
aggregate() # 每次只能用一个统计函数
# Descriptive stats by group via summaryBy
install.packages("doBy")
library(doBy)
summaryBy(mpg+hp+wt ~ am, data=mtcars, FUN=mystats)
# Descriptive stats by group via describe.by (psych)
library(psych)
describeBy(mtcars[myvars], list(am=mtcars$am)) # 没办法使用自定义的函数
P45 频数统计函数
# frequency table
split(mtcars, as.factor(mtcars$cyl)) # 以cyl的因子分类
cut(mtcars$mpg, c(seq(10, 50, 10))) # 最后一个数为间隔数（10-20，20-30，30-40，40-50）
table(mtcars$cyl) # 频数的统计
prop.table(table(mtcars$cyl)) # 计算每种因子的比例
# 二维
library(vcd)
table(Arthritis$Treatment, Arthritis$Improved)
with(data=Arthritis, table(Treatment, Improved))
x <- xtabs(~ Treatment + Improved, data = Arthritis)
margin.table(x, 1)
prop.table(x, 1)
# 1 = 只添加行
# 2 = 只添加列
# 三维
y <- xtabs(~ Treatment + Improved + Sex, data = Arthritis)
ftable(y) # 平铺式的列联表
P46 独立性检验函数
# 根据频数信息，判断两类因子，彼此相关或相互独立的假设检验
# 独立性：变量之间是独立的，没有关系
# p-value (probability)：在原假设为真时，得到最大的或者超出所得到的检验统计量值的概率。
# 原假设：没有发生 / 备择假设：发生了
# 一般将p值定位到0.05，当p<0.05拒绝原假设（不靠谱），当p>0.05不拒绝原假设（靠谱）。
# p值设定的越小，越精确。
# 卡方检验（Chi-square test）
library(vcd)
mytable <- table(Arthritis$Treatment, Arthritis$Improved) # 变量的顺序不重要
chisq.test(mytable)
# Fisher's exact test
mytable <- xtabs(~Treatment+Improved, data=Arthritis)
fisher.test(mytable)
# Cochran-Mantel-Haenszel检验
mytable <- xtabs(~Treatment+Improved+Sex, data=Arthritis) # 变量的顺序很重要（因为要三个变量）
mantelhaen.test(mytable)
P47 相关性分析函数
相关性分析：对两个或多个具备相关性的变量元素进行分析，从而衡量两个变量因素的相关密切程度。
相关性的元素之间需要存在一定的联系或者概率才可以进行相关性分析（变量之间是否有关系）
正相关或负相关（相关系数）
Pearson相关系数
Spearman相关系数
Kendall相关系数
偏相关系数
多分格（polychoric）相关系数
多系列（polyserial）相关系数
# 使用cor()计算相关系数
cor(state,x77) # 默认method="pearson"
cor(state.x77, method="spearman")
# 正负号表示正相关或负相关
cov(state.x77) # 协方差
x <- state.x77[, c(1,2,3,6)]
y <- state.x77[, c(4,5)]
cor(x, y)
# partial correlations 偏相关系数
install.packages("ggm")
library(ggm)
# partial correlation of population(1) and murder rate(5), controlling
# for income(2), illiteracy rate(3), and HS graduation rate(6)
colnames(state.x77)
pcor(c(1,5,2,3,6), cov(state.x77))
P48 相关性检验函数
# Testing a correlation coefficient for significance
cor.test(state.x77[,3], state.x77[,5])
Confidence interval（置信区间）
样本统计量所构造的总体参数的估计区间（被测量参数的测量值的可信程度）
概率发生的范围
# Correlation matrix and tests of significance via corr.test
corr.test(states, use="complete")
library(psych)
corr.test(state.x77)
library(ggm)
x <- pcor(c(1,5,2,3,6), cov(state.x77))
pcor.test(x, 3, 50) # 3 -> 变量数，50 -> 样本数
# Two Sample t-test
library(MASS)
t.test(Prob ~ So, data=UScrime)
t.test(y ~ x, data = )
非参数检验（Nonparametric test）：在总体方差未知或知道甚少的情况下，利用样本对总体分布形态等进行推断的方法。
参数检验（Parametric test）:总体分布形式已知的情况下，对总体分布的参数（如均值、方差等）进行推断的方法。也就是数据分布已知，比如满足正态分布。
# dependent t test
sapply(UScrime[c("U1", "U2")], function(x)(c(mean=mean(x), sd=sd(x))))
with(UScrime, t.test(U1, U2, paired=T))
P49 绘图函数
R基础绘图系统
ls("package:graphics")
demo(graphics)
help(package=graphics)
（1）高级绘图：一步到位，直接绘制出图
（2）低级绘图：不能单独使用，得在高级绘图产生图形的基础上，对图形进行调整，比如加线、标题文字等
知道输入数据的格式很重要：
散点图：x和y两个坐标数据
直方图???：因子
热力图：数据矩阵
plot(women$height)
plot(womenheight, womenweight)
plot(as.factor(women$height))
plot(as.factor(mtcars$cyl)) # 直方图
plot(as.factor(mtcarscyl), mtcarscarb) # 箱线图
plot(mtcarscarb, as.factor(mtcarscyl)) # 散列图
plot(as.factor(mtcarscyl), as.factor(mtcarscarb)) # 脊柱图
plot(womenheight ~ womenweight) # 另一种写法
fit <-lm(height ~ weight, data=women)
plot(fit) # plot会自动检测输入的数据类型来进行不同的绘图
method(plot)
method(summary)
#S3系统：属性，泛型函数，方法
?par # parameter（用于微调）
par() # R绘图的默认设置
plot(as.factor(mtcars$cyl), col=c("red", "green", "blue")) # col为其中一项参数
lattice包
ggplot2包
grid包
P50 自定义函数
自己编写函数
输入不加括号的函数，输出源代码
- 函数名称
1.1 函数命令与功能相关
1.2 可以是字母与数字的结合，但必须是字母开头
-
函数声明
-
函数参数
-
函数体
myfun <- function(选项参数)
{
函数体
}
例子：计算偏度和峰度的函数
偏度（skewness）：统计数据分布偏斜方向和程度的度量，是统计数据分布非对称程度的数字特征。
峰度（peakedness；kurtosis）：峰态系数。表征概率密度分布曲线在平均值处峰值高低的特征数。
mystats <- function(x, na.omit=FALSE) {
if(na.omit)
x <- x[!is.na(x)] # x只取不包含NA的值
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, skewness=skew, kurtosis=kurt))
}
循环与向量化操作：函数内部通过循环实现向量化操作
if条件判断
score=70;if (score > 60) {print("Passed")} else {print("Failed")}
ifelse(score > 60, "Passed", "Failed") # 两种写法（ifelse自带输出功能）
for循环
for (i in 1:10) {print("Hello")}
while循环
i=1;while(i <= 10) {print("Hello");i=i+1;}
switch语句
循环三部分
条件判断，真或假
用于循环执行的结构
表达式
P51 数据分析实战
案例
小麦产量案例
房地产案例
性价比案例
药物试验案例
社会科学研究案例
量化投资案例（计量经济学）
反面案例
苹果的Think different？？？
P52-53 线性回归（regression）
回归：指那些用一个或多个预测变量（自变量、解释变量），来预测响应变量（因变量、效标变量、结果变量）的方法。
根据一大堆数据，找到它们的规律。
如何建立模型。抽象出数学公式，哪些因素与模型有关，需要利用多少样品，模型的准确率有多高，在实际运用中还是否有效。
线性回归 或 非线性回归
普通最小二乘回归法（OLS）
plot(womenheight, womenweight)
fit <- lm(formula = weight ~ height, data = women) # 线性回归分析（拟合）
summary(fit)
# 残差（residuals）越小，越接近线性关系
# R-squared 表示符合模型的数据占比
# 一般先看p-value是否<0.05，如果不小于一般说明模型不符合
# 再看R-squared，有多少数据符合模型
"~" #分割符号
左边为响应变量（y），右边是解释变量（x）
"+" # 分隔预测/解释变量
":" # 表示预测变量的交互项
比如，要通过x、z及x与z的交互项来预测y，代码为y ~ x + z + x:z
"*" # 表示所有可能交互项的简洁方式
代码 y ~ x * z * w 可展开为 y ~ x + z + w + x:z + x:w + z:w + x:z:w
"^" # 表示交互项达到某个次数
代码 y ~ (x + z + w)^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上的回归，并强制直线通过原点
I() # 从算术的角度来解释括号中的元素
比如，y~x+(z+w)^2 展开为 y~x+z+w+z:w
代码 y~x+ I((z+w)^2) 展开为 y~x+h （h是一个由z和w的平方和创建的新变量）
"function()" # 可以在表达式中用的数学函数
比如，log(y) ~ x + z + w 表示通过x、z和w来预测log(y)
线性拟合常用函数
summary() # 展示拟合模型的详细结果
coefficients() # 列出拟合模型的模型参数（y-intercept和slope）
confint(data, level = 0.5) # 提供模型参数的confidence interval（默认95%，用level参数调整）
fitted() # 列出拟合模型的预测值
residuals() # 列出拟合模型的残差值（残差（residuals）越小，越接近线性关系）
anova() # 生成一个拟合模型的方差分析表，或者比较两个或多个拟合模型的方差分析表
vcov() # 列出模型参数的协方差矩阵
AIC() # 输出赤池信息统计量
plot() # 生成评价拟合模型的诊断图
predict(model, new.data) # 用拟合模型对新的数据集预测响应变量值
abline() # 绘制拟合曲线
多项式
fit2 <- lm(weight ~ height + I(height^2), data = women)
summary(fit2)
lines(women$height, fitted(fit2), col="red")
lines(x, y) # 横坐标（x），纵坐标（y）（这里纵坐标用了预测值）
三次项
fit3 <- lm(weight ~ height + I(height^2) + I(height^3), data = women)
summary(fit3)
lines(women$height, fitted(fit3), col="blue")
不能拟合不足，也不能拟合过度
P54 多元线性回归
多个变量，权重也不同
states <- as.data.frame(state.x77, c("Murder", "Population", "Illiteracy", "Income", "Frost"))
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
summary(fit)
coef(fit)
qqplot(fit, labels = row.names(states), id.method = "identify", simulate = T, main = "Q-Q Plot")
fit2 <- lm(mpg ~ hp + wt + hp:wt, data = mtcars)
summary(fit2)
如何从众多模型中，选择最佳模型
比较模型,AIC值越小越好
m1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
m2 <- lm(Murder ~ Population + Illiteracy, data = states)
AIC(m1, m2)
逐步回归法（Backward stepwise selection）
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")
全子集回归法（All subsets regression）
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")
P55 回归诊断
fit <- lm(weight ~ height, data=women)
opar <- par(no.readonly = T)
par(mfrow=c(2,2)) # 表示横排竖排各显示两幅
plot(fit)
par(opar)
满足OLS模型统计假设
- 正态性：对于固定的自变量值，因变量值成正态分布
- 独立性：因变量之间互相独立
- 线性：因变量与自变量之间为线性关系
- 同方差性：因变量的方差不随自变量的水平不同而变化。也可称作不变方差。
fit2 <- lm(weight ~ height + I(height^2), data = women)
opar <- par(no.readonly = T)
par(mfrow=c(2,2))
plot(fit2)
par(opar)
抽样法验证
-
数据中有1000个样本，随机抽取500个数据进行回归分析
-
模型建好后，利用predict()，对剩下的500个样本进行预测，比较残差值
-
如果预测准确，说明模型可以，否则就需要调整模型
P56-57 方差分析（Analysis of Variance）
方差分析（Analysis of Variance / ANOVA），也称变异数分析
用于两个及两个以上样本均数差别的显著性检验。
从广义上来讲，方差分析属于回归分析的一种。但线性回归的因变量一般是连续型变量。而当自变量是因子时，研究关注的重点通常会从预测转向不同组之间差异的比较。
1.单因素方差分析ANOVA（组内，组间）
y ~ A
# 案例
library(multcomp)
attach(cholesterol)
table(trt)
aggregate(response, by=list(trt), FUN=mean)
fit <- aov(response ~ trt, data = cholesterol)
summary(fit)
plot(fit)
fit.lm <- lm(response ~ trt, data = cholesterol)
summary(fit.lm)
2.双因素方差分析ANOVA
y ~ A * B
# 案例
attach(ToothGrowth)
xtabs(~ supp+dose)
aggregate(len, by=list(supp, dose), FUN=mean)
aggregate(len, by=list(supp, dose), FUN=sd)
ToothGrowth$dose <- factor(ToothGrowth$dose)
fit <- aov(len ~ supp*dose, data=ToothGrowth)
summary(fit)
install.packages("HH")
interaction.plot(dose, supp, len, type = "b", col = c("red", "blue"),
pch = c(16,18),
main = "Interaction between Dose and Supplement Type")
3.协方差分析ANCOVA
y ~ x + A
y ~ x1 + x2 + A*B
# 单因素协方差案例
table(litter$dose)
attach(litter)
aggregate(weight, by=list(dose), FUN=mean)
fit <- aov(weight ~ gesttime+dose, data=litter) # gesttime为协变量
summary(fit)
4.多元方差分析MANOVA
# 案例
library(MASS)
attach(UScereal)
shelf <- factor(shelf)
aggregate(cbind(calories, fat, sugars), by=list(shelf), FUN=mean)
fit <- manova(cbind(calories, fat, sugars) ~ shelf, data=UScereal)
summary(fit)
summary.aov(fit)
5.多元协方差分析MANCOVA
?aov
# 顺序很重要
y ~ A + B + A:B
# 在R中，顺序为序贯型。效应根据表达式中先出现的效应做调整。
# A不做调整，B根据A做调整，A:B交互项根据A和B调整。
P58 功效分析（Power Analysis）
功效分析，可以帮助在给定置信度的情况下，判断检测到给定效应值时所需的样本量。
反过来，它也可以在给定置信度水平情况下，计算在某样本量内能检测到给定效应值的概率。
# 功效分析函数（pwr包中的函数）
install.packages("pwr")
pwr.2p.test() # 两比例（n相等）
pwr.2p2n.test() # 两比例（n不相等）
pwr.anova.test() # 平衡的单因素ANOVA
pwr.anova.test(k=2, f=25, sig.level=0.05, power=0.90)
pwr.chisq.test() # 卡方检验
pwr.f2.test() # 广义线性模型（Linear Models）
pwr.f2.test(u=3, f2=0.0769, sig.level=0.05, power=0.90)
pwr.p.test() # 比例（单样本）
pwr.r.test() # 相关系数
pwr.t.test() # t检验（单样本、两样本、配对）
pwr.t.test(d=.8, sig.level=.05, power=.9, type="two.sample", alternative="two.sided")
pwr.t.test(n=20, d=.5, sig.level=.01, type="two.sample", alternative="two.sided")
pwr.t2n.test() # t检验（n不相等的两样本）
功效分析理论基础
-
样本大小：实验设计中每种条件/组中观测的数目。
-
显著性水平（alpha）：由I型错误的概率来定义，也可以把它看作是发现效应不发生的概率。
-
功效：通过I减去II型错误的概率来定义，可以把它看作是真实效应发生的概率。
-
效应值：指的是在备择或研究假设下效应的量。效应值的表达式依赖于假设检验中使用的统计方法。
有其中三个，就可以求出第四个数
P59 广义线性模型
广义线性模型扩展了线性模型的框架，它包含了非正态因变量的分析。
?glm # Generalized Linear Models
data(breslow.dat, package="robust") # 加载数据集
names(breslow.dat)
summary(breslow.dat[c(6,7,8,10)]) # 统计数据集
泊松回归（Poisson regression）
用来为计数资料和列联表建模的一种回归分析。
假设因变量是泊松分布，并假设它平均值的对数可被未知参数的线性组合建模。
attach(breslow.dat)
# fit regression
fit <- glm(sumY ~ Base + Age + Trt, data = breslow.dat, family = poisson(link="log"))
summary(fit)
# intercept model parameter
coef(fit)
exp(coef(fit))
P60 Logistic回归
通过一系列连续型或类别型预测变量，来预测二值型结果变量时，Logistics回归很有用。
data(Affairs, package="AER")
summary(Affairs)
table(Affairs$affairs)
prop.table(table(Affairs$affairs))
prop.table(table(Affairs$gender))
# create binary outcome variable
Affairs$ynaffairs[Affairs$affairs > 0] <- 1
Affairs$ynaffairs[Affairs$affairs == 0] <- 0
Affairs$ynaffairs <- factor(Affairs$ynaffair, levels = c(0,1), labels = c("No", "Yes"))
table(Affairs$ynaffair)
# fit full model
attach(Affairs)
fit <- glm(ynaffairs ~ gender + age + yearsmarried + children +
religiousness + education + occupation + rating,
data = Affairs, family = binomial())
summary(fit)
# fit reduced model
fit2 <- glm(ynaffairs ~ age + yearsmarried + religiousness + rating,
data = Affairs, family = binomial())
summary(fit2)
# compare models
anova(fit, fit2, test = "Chisq")
coef(fit2)
exp(coef(fit2))
# calculate probability of extramarital affair by marital ratings
testdata <- data.frame(rating = c(1,2,3,4,5),
age = mean(Affair$age),
yearsmarried = mean(Affair$yearsmarried),
religiousness = mean(Affair$religiousness))
testdata$prob <- predict(fit2, newdata = testdata, type = "response")
testdata
# calculate probability of extramarital affair by age
testdata <- data.frame(rating = mean(Affairs$rating),
age = seq(17, 57, 10),
yearsmarried = mean(Affair$yearsmarried),
religiousness = mean(Affair$religiousness))
testdata$prob <- predict(fit2, newdata = testdata, type = "response")
testdata
P61 主成分分析
主成分分析（Principal Component Analysis / PCA）
一种数据降维技巧。
它能将大量相关变量转化为一组很少的不相关变量，这些无关的变量称为主成分。
主成分其实是对原始变量重新进行线性组合，将原来众多具有一定相关性的指标，重新组合为一组的新的相互独立的综合指标。
主成分分析公式
P
C
1
=
A
1
X
1
+
A
2
X
2
+
.
.
.
+
A
k
X
k
PC1 = A1X1 + A2X2 + ...+ AkXk
PC1=A1X1+A2X2+...+AkXk
主成分分析与因子分析步骤
-
数据预处理
-
选择分析模型
-
判断要选择的主成分/因子数目（可通过碎石图）
-
选择主成分/因子
-
旋转主成分/因子
-
解释结果（Optional）
-
计算主成分/因子得分（Optional）
# Principal components analysis of the Judge Ratings
library(psych)
USJudgeRatings
fa.parallel(USJudgeRatings, fa="pc", n.iter=100) # 碎石图（找出n.factors）
pc <- principal(USJudgeRatings, nfactors=1)
pc
# PCA Score
pc <- principal(USJudgeRatings, nfactors=1, scores=T) # score为是否需要计算主成分得分，默认F
pc$scores
# PCA Harman23.cor data
fa.parallel(Harman23.cor$cov, n.obs=302, fa="pc", n.iter=100, show.legend=F, main="Scree plot with parallel analysis")
# PCA of body measurement
pc <- principal(Harman23.cor$cov, nfactors=2, rotate="none")
pc
使成分不相关：正交旋转（去噪）
使成分相关：旋交旋转
P62 因子分析
探索性因子分析法（Exploratory Factor Analysis / EFA）
一系列用来发现一组变量的潜在结构的方法。
它通过寻找一组更小的、潜在的或隐藏的结构来解释已观测到的、显式的变量间的关系。
因子分析公式
X
i
=
A
1
F
1
+
A
2
F
2
+
.
.
.
+
A
k
F
k
+
U
i
Xi = A1F1 + A2F2 + ... + AkFk + Ui
Xi=A1F1+A2F2+...+AkFk+Ui
Xi：第i个可观测的变量
F：公共因子
Ui：Xi变量无法被公共因子解释的部分（误差）
主成分与因子分析比较
相同点
-
都对原始数据进行降维处理
-
都消除了原始指标的相关性对综合评价所造成的信息重复的影响
-
构造综合评价时所涉及的权数具有客观性
-
在信息损失不大的前提下，减少了评价工作量
-
用碎石图来判断主成分/因子数量
不同点
主成分分析：
-
用较少的变量表示原来的样本
-
目的是样本数据信息损失最小的原则下，对高维变量进行降维
-
参数估计，一般是求相关矩阵的特征值和相应的特征向量，取前几个计算主成分
-
应用：应用较少变量来解释各个样本的特征
因子分析：
-
用较少的因子表示原来的变量
-
目的是尽可能保存原变量的相互关系，寻找变量的公共因子
-
参数估计，指定几个因子，将其还原成相关系数矩阵，在和原样本相关矩阵最相似原则下，估计各个公因子的估计值
-
应用：找到具有本质意义的少量因子来归纳原来变量的特征
factanal()
options(digits = 2)
covariances <- ability.cov$cov
correlations <- cov2cor(covariances) # 转换成系数矩阵
fa.parallel(correlations, fa="both", n.obs=112, n.iter=100)
fa <- fa(correlations, nfactors=2, rotate="none", fm="pa")
# 正交旋转（PC1和PC2没有关系）
fa.varimax <- fa(correlations, nfactors=2, rotate="varimax", fm="pa")
# 旋交旋转（PC1和PC2有关系）
fa.promax <- fa(correlations, nfactors=2, rotate="promax", fm="pa")
factor.plot(fa.promax, labels = rownames(fa.promax$loadings))
fa.diagrams(fa.varimax, simple = F)
fa(correlations, nfactors=2, rotate="none", fm="pa", score = T)
fa$weights # 得分权重
P63 购物篮分析
# 数据加载
install.packages("arules")
data(Groceries)
inspect(Groceries)
fit <- apriori(Groceries, parameter = list(support=0.01, confidence=0.5))
summary(fit)
inspect(fit)
# lift指数越大，逻辑可能性越大
P64-65 后记
#白嫖不易，各自珍惜。
#希望世界和平。
287
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
