Jimmy那些事儿

R语言_统计分析

发表 | 更新 | 分类 ， ， | 阅读次数

运用R语言完成数据的统计分析。

包括统计函数、交叉表、相关性等。

统计函数

基本统计量

名称 含义 示例
mean(x) 均值
median(x) 中位数
sd(x) 标准差
var(x) 方差
mad(x) 绝对中位差
quantile(x,) 分位数，probs是0-1之间的数值向量 quantile(x,probs) quantile(x,c(.3,.7))
range(x) 求值域 rang(c(1,2,3,4))，返回(1,4)
sum(x)
diff(x,) 滞后分差； diff(x,lag=n) 1)diff(c(1,5,23,29),lag=2)
min(x)
max(x)
scale(x,) 为数据对象x按列进行中心化或标准化 scale(x,center=TRUE,scale=TRUE)

scale

默认情况下，scale对矩阵或数据框的指定列进行均值=0，标准差=1的标准化；

中心化 center=T

标准化 center=T,scale=T

若对任意均值，SD为标准差，M为均值

newdata <- scale(mydata)*SD+M

mad( ) - 绝对中位差

1.求原样本S1中位数，记为m1
2.原样本S1各元素减去m1，得到的值去绝对值，得到新的样本S2
3.求新样本S2的中位数m2
4.返回的值=1.4826*m2

概率函数

名称 作用 全称 示例 说明
d 概率函数 density
p 分布函数 distribution function pnorm(x,mean=a,sd=b) 生成mean=a,sd=b的x【左侧的面积值】；pnorm(1.96)=0.975
q 分位数函数 quantile function qnorm(0.x,mean=a,sd=b) 生成0.x分位数，mean=a，sd=b的情况下；默认为标准正态
r 生成随机数 random rnorm(x,mean=a,sd=b) 生成x个mean=a,sd=b的随机数；默认情况下a=0,b=1，且mean/sd可省略
概率函数 含义
beta Beta分布
binom 二项分布
nbinom 负二项分布
multinom 多项分布
norm 正态分布
lnorm 对数正态分布
unif 均匀分布
t t分布
f F分布
logis Logistic分布
exp 指数分布
chisq （非中心）卡方分布
pois 泊松分布
geom 几何分布
hyper 超几何分布
cauchy 柯西分布
signrank Wilcoxon符号秩分布
wilcox Wilcoxon秩和分布
weibull Weibull分布

通用函数

1
2
3
## 截尾
trim=0.n  # 在统计函数中加入trim表示截去最高和最低各n%，计算中间(1-n)%的值
c<-matrix(rnorm(12),1,mean,trim=0.2)


交叉表

  • what ：一种分类汇总的表格 Cross Tabulations
  • why ：用于查看变量的频数、占比

多结果列联表

1
2
3
4
5
6
7
8
9
10
11
12
13
14
library(gmodels)
CrossTable(x,y, digtis=3, max.width=5, expected=F, prop.r=T, prop.c=T, prop.t=T, prop.chisq=TRUE, chisq=F, fisher=F, mcnemar=F, missing.include=F, format=c("SAS","SPSS"), dnn=NULL,...) 
# x，一个矩阵或向量。如果出现y，则x必须是一个向量
# max.width= ，如果想要一个单柱立式表，可设为1
# expected=F/T, 如果为T，chisq将被设置为TRUE和Chi-Square单元计数将被纳入
# prop.r - 行比例 ; prop.c - 列比例;  prop.t - 表比例; 
# prop.chisq - 卡方贡献的每个单元格
# format="SAS" / "SPSS" ,默认输出格式为"SAS"
# 独立性检验  chisq=F / fisher=F / mcnemar= F
CrossTable(data$col, data$col); / with(data, CrossTable(col1, col2))
CrossTable(HR_data$sales, HR_data$salary, prop.r=T,prop.c=F, prop.t=F, prop.chisq=F)  # 只有prop.r=T ,表示对每一行求比例

单结果列联表

  • 表格紧凑化 - ftable(table)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
计算频数
table(object1, object2, object3, [useNA=ifany])
# 若出现useNA=ifany，则将NA视为一个有效类别
# x <- with(data, table(col))
xtabs(formula, data=matrix/df )  # 从矩阵/数据框中提取
xtabs(~col1+col2, data=frame_name)
-----------------------------------------------------------------------------------
计算占比
prop.table(table, margins) 
# 一般都需指定所求变量的比例（默认为表比例）
## prop.table(table) --对所有变量求比例
## prop.table(table, 1) --对第一个变量求比例
## porp.table(table, 2） --对第二个变量求比例
## prop.table(table, c(1,2)) -对第一、第二个变量求比例；当维度大于2时
## prop.table(table)*100 -- 得到百分比形式
-----------------------------------------------------------------------------------
计算边际和
margin.table(table, margins)
添加新的结果
addmargins(table, margins)
# addmargins(table) -对所有变量创建边际和
# addmargins(prop.table(table, 1), 2) -对table表中第一个变量求比例，并对第二个变量创建边际和


相关性

  • what ： 描述定量变量之间的关系；衡量两个样本之间的相关性；若为0，则表示两者不相关，即独立
    • 相关系数
      • 线性相关程度-Pearson积差相关系数
      • 分级定序变量的相关程度-Spearman等级相关系数
      • 非参数等级相关程度-Kendall’s Tau相关系数
    • 协方差：一个变量的偏离程度，会对另一个变量的偏离程度产生多大的影响；
      • 如果其中一个大于自身的期望值，另外一个也大于自身的期望值，那么两个变量之间的协方差就是正值。
      • 如果X 与Y 是统计独立的，那么二者之间的协方差就是0
      • 但是反过来并不成立，即如果X 与Y 的协方差为0，二者并不一定是统计独立的。取决于协方差的相关性η
    • 方差：当两个变量是相关的特殊情况
    • 偏相关：控制一个变量或多个定量变量时，另外两个变量之间的关系
  • why ：确定变量之间的影响程度的强弱，以得出变量的重要性[为后续分析提供变量选择]

相关系数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
cor(matrix/df, use="everything", method="pearson")
# 先选择定量变量的数据集 data<- select(data, )
# use= ，表示对缺失值NULL的处理方式
## use="everything"(默认) , 结果为missing;
## use="all.obs",结果报错;  
## use="complete.obs"--行删除; 
## use="pairwise.complete.obs"--成对删除
# method=， 表示选择的相关系数种类
## method="pearson"（默认）
## method-"spearman"
## method="kendall"
cor(matr/df) - 得到一个两两计算相关性的方阵
cor(HR_comma[ , 1:8]) # 只可对定量变量做相关性
# 求一组变量与另一组（非两两）的相关矩阵
x <- states[ , c("col", "col")]
y <- states[ , c("col", "col", "col"]
cor(x, y)
     
                 
-----------------------------------------------------------------------------------              
cov(matrix/df, use="everything", method="pearson")

–

1
2
3
4
5
6
7
8
偏相关
library(ggm)
pcor(vector, S)
# vector，一个数值型向量, 数值代表变量的序数；前两个数值表示要计算的变量下标，其余为要排除的变量
# S，变量的协方差阵
colnames(data)
pcor(c(1,2,4,5), cov(data))

相关图 - 可视化

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
library(corrgram)
corrgram(dataX, order=TRUE, lower.panel=, upper.panel=, text.panel=, diag.panel=, main="")
# dataX，是一行一个观测值的数据框，即长格式  或者是一个 相关性矩阵
## [reshape2包] - melt(data, [id="" / c("", "")] 
# order=TRUE, 相关矩阵使用【主成分分析法】对变量进行重排序
# lower/upper.panel=,控制非对角线的设置
## NULL - 不显示
## panel.shade - 用阴影深度表示相关性大小
## panel.pie - 用饼图填充比例表示相关性大小（顺时针，正相关）
## panel.pts - 画一个散点图
## panel.ellipse - 画一个置信椭圆+平滑曲线
## panel.conf - 画出相关性+置信区间（用数值表示，而非图片）
# text.panel=panel.txt - 输出变量名
# cex.labels=num，变量名的大小
# diag.panel=，控制主对角线的设置
## panel.minmax - 输出变量的最大最小值和变量名
## panel.ednsty - 输出核密度曲面和变量名
corrgram(mtcars /cor(data), order=TRUE, lower.panel=panel.pei, upper.panel=panel.conf, text.panel=panel.txt, main="")
-----------------------------------------------------------------------------------  
library(corrplot)
corrplot(corr, method="number", order="FPC" tl.col="black", title="",tl.srt=45，mar=c(1,1,3,1)...)
# corr，是一个相关性矩阵，且必须是方形；
# method="", 可视化的图形，circle（默认），square, ellipse, number, shade, color, pie
# type=""，图形矩阵的位置，full（默认，全部完整矩阵），lower/upper,显示下三角/上三角
# order=""，相关矩阵的排序方法，original(默认,原来的顺序), FPC(第一主成分), hclust(层次聚类次序), AOE(特征向量角的顺序)
# add=FALSE(默认)，不添加到现有的图形中（即创建一个新的图形），若为TRUE，则添加到现有的图形
# mar=c(1,1,3,1)  # 用来调节边界大小，以放下标题
# tl. 文本标签
## tl.cex= # 文本标签的大小
## tl.col= 文本标签的颜色
corrplot(corr,method="circle")


类别变量


总体均值对比


回归

  • what ：建立因变数Y与自变数 X之间关系的模型
  • why ：
    • 了解变量之间是否相关？相关的方向与强度？
    • 并建立数学模型以观察特定变量来预测因变量；（可以给出自变量来估计因变量的条件期望）
  • how ：分析步骤
  1. 确定因变量与自变量，选择回归方程
1
2
3
4
因变量 - 二值型数值/类别型变量
# Logistic回归
因变量
  1. 整理数据格式
1
2
3
4
5
6
7
8
Logistic回归
赋值 - 二分法 
# data$col[data$col cond] <- 1
# data$col[data$col cond] <- 0
因子化
# data$col <- as.factor(data$col)
  1. 拟合模型
1
2
Logistic回归 
# fit <- glm(formula, data=, family=binomail)
  1. 查看统计量，并重新选择自变量
1
summary(fit)
  1. 重新拟合模型，并查看统计量
1
2
3
Logistic回归 - 
# fit2 <- glm(formula, data=, family=binomal)
# summary(fit2)
  1. 模型选择 - 嵌套对比
1
2
Logistic回归 - 卡方检验(若p值<0.05，说明多余的变量对方程有显著影响，不能舍去)
# anova(fit, fit2, test="Chisq")
  1. 解释模型
1
2
Logisctic回归
# 指数化 - exp(coef(fit))
  1. 深层次分析
  • 交叉验证：- 评价回归方程的泛化能力；即，回归方程对新观测样本的预测能力
  • 模型评价：评价预测变量对结果概率的影响【在选择最佳模型之后】
  • 过度离势 [ Logistic回归 ]
    • [what] - 过度离势：观测到的响应变量的方差 > 期望的二项分布的方差(σ²=npq)；
    • [why] - 过度离势会导致标准误检验和不精确的显著性检验；若出现过度离势，可用类二项分布来代替，即glm(,,family=quasibinomial)
    • [how] - 比较模型的残差偏差与残差自由度的比值；若比值非常接近1，则没有过度离势；若比1大很多，则存在过度离势；deviance(fit) / df.residual(fit)
  1. 预测模型
1
2
Logistic回归
# predict(fit, newdata=, type="response")


确定统计量




参考


表达式

  • 常见表达式
表达式 作用 示例 解释
～ 分隔符号 y～x+z+w 左边-因变量/响应变量 右边-自变量/解释变量
+ 分隔预测变量
： 预测变量的交互项 y～x+z+x:z
* 所有可能交互项的表达方式 y～xzw—y～x+z+w+x:z+x:w+z:w
^ 交互项的某个次数 y～(x+z+w)^2—y～x+z+w+x:z+x:w+z:w 交互项最高次为2次
. 包含除因变量之外的所有变量 y～.—y～x+z+w 当一个数据框包含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() [大写i ] 从算术（而非表示式）的角度来解释括号中的元素 y～x+I((z+w)^2) 表示的是x+(z+w)²，而非x+z+w+z:w
function 可以在表达式中运用的数学函数 log(y)～x+z+w
mpg ~ wt \ cyl 表示按条件（cyl）绘图； 即按cyl的水平为【分组条件】分别绘制mpg和wt的关系图


  • 研究设计表达式
表达式 作用 示例
y～A 单因素ANOVA
y～x+A 含单个协变量的单因素ANCOVA
y～A * B 双因素ANOVA 展开为 ~A+B+A:B
y～x1+x2+A*B 含两个协变量的双因素ANCOVA
y～B+A（B是区组因子） 随机化区组
y～A + Error(Subject/A) 单因素组内ANOVA
y～B*W+Error(Subject/W) 含单个组内因子(W)和单个组间因子(B)的重复测量ANOVA 展开为 ~B+W+B:W

1.小写字母，定量变量

2.大写字母，组别因子（若不转换为factor，则默认为定量协变量）

3.Subject，被试者独有的标志变量

4.Error(Subject/A)，表示组内因子