嘤嘤嘤还要学统计


基本操作

最基本

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
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
x1=c(171,175,159,155,152,158,154,164,168,166,159,164)
> x1
[1] 171 175 159 155 152 158 154 164 168 166 159 164
> x2=c(57,64,41,38,35,44,41,51,57,49,47,46);x2
[1] 57 64 41 38 35 44 41 51 57 49 47 46
> rbind(x1,x2)#按行合并
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
x1 171 175 159 155 152 158 154 164 168 166 159 164
x2 57 64 41 38 35 44 41 51 57 49 47 46
> cbind(x1,x2)#按列合并
x1 x2
[1,] 171 57
[2,] 175 64
[3,] 159 41
[4,] 155 38
[5,] 152 35
[6,] 158 44
[7,] 154 41
[8,] 164 51
[9,] 168 57
[10,] 166 49
[11,] 159 47
[12,] 164 46
#利用x1数据创建矩阵
> matrix(x1,nrow=3,ncol=4)
[,1] [,2] [,3] [,4]
[1,] 171 155 154 166
[2,] 175 152 164 159
[3,] 159 158 168 164

#创建行数列数发生变化的矩阵
> matrix(x1,nrow=4,ncol=3)
[,1] [,2] [,3]
[1,] 171 152 168
[2,] 175 158 166
[3,] 159 154 159
[4,] 155 164 164
#创建两个相同的矩阵
> A=B=matrix(1:12,nrow=3,ncol=4)
> A+B#矩阵加法
[,1] [,2] [,3] [,4]
[1,] 2 8 14 20
[2,] 4 10 16 22
[3,] 6 12 18 24
> A-B#矩阵减法
[,1] [,2] [,3] [,4]
[1,] 0 0 0 0
[2,] 0 0 0 0
[3,] 0 0 0 0

> A=matrix(1:9,nrow=3,ncol=3)
> B=matrix(1:9,nrow=3,ncol=3)
> A*B#矩阵对应元素的乘积
[,1] [,2] [,3]
[1,] 1 16 49
[2,] 4 25 64
[3,] 9 36 81
> A%*%B#矩阵的乘积
[,1] [,2] [,3]
[1,] 30 66 102
[2,] 36 81 126
[3,] 42 96 150


> A=matrix(1:16,nrow=4,ncol=4)
> diag(A)#获得矩阵对角线元素
[1] 1 6 11 16
> diag(diag(A))#利用对角线元素创建对角矩阵
[,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 0 6 0 0
[3,] 0 0 11 0
[4,] 0 0 0 16

> A=matrix(rnorm(16),4,4)
> solve(A)#求矩阵的逆
[,1] [,2] [,3] [,4]
[1,] -0.2144625 0.9291739 0.41134588 -0.4129970
[2,] -0.2990380 0.5107384 -0.03081756 0.6477244
[3,] -0.1700996 1.4089089 -0.49639658 0.4082339
[4,] -0.4334604 -0.3543499 -0.30915934 -0.5515849
> A=diag(4)+1
> A.e=eigen(A,symmetric=T)
> A.e#求特征根和特征向量
eigen() decomposition
$values
[1] 5 1 1 1

$vectors
[,1] [,2] [,3] [,4]
[1,] -0.5 0.8660254 0.0000000 0.0000000
[2,] -0.5 -0.2886751 -0.5773503 -0.5773503
[3,] -0.5 -0.2886751 -0.2113249 0.7886751
[4,] -0.5 -0.2886751 0.7886751 -0.2113249

> A=matrix(1:12,3,4)
> dim(A)#矩阵的维数
[1] 3 4
> nrow(A)#矩阵的行数
[1] 3
> ncol(A)#矩阵的列数
[1] 4
> rowSums(A)#矩阵按行求和
[1] 22 26 30
> rowMeans(A)#矩阵按行求均值
[1] 5.5 6.5 7.5
> colSums(A)#矩阵按列求和
[1] 6 15 24 33
> colMeans(A)#矩阵按列求均值
[1] 2 5 8 11

函数

apply()函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
> apply(A,1,sum)#矩阵按行求和
[1] 22 26 30
> apply(A,1,mean)#矩阵按行求均值
[1] 5.5 6.5 7.5
> A=matrix(rnorm(100),20,5)
> apply(A,2,var)#矩阵按列求方差
[1] 0.864552 1.153079 1.143116 1.264152 1.061721
> B=matrix(1:12,3,4)
#矩阵按列求函数结果
> apply(B,2,function(x,a)x*a,a=2)
[,1] [,2] [,3] [,4]
[1,] 2 8 14 20
[2,] 4 10 16 22
[3,] 6 12 18 24

数据框

#产生由X1和X2构建的数据框

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
> X=data.frame(x1,x2);X
x1 x2
1 171 57
2 175 64
3 159 41
4 155 38
5 152 35
6 158 44
7 154 41
8 164 51
9 168 57
10 166 49
11 159 47
12 164 46

#赋予数据框新的列标签
> X=data.frame('身高'=x1,'体重'=x2);X
身高 体重
1 171 57
2 175 64
3 159 41
4 155 38
5 152 35
6 158 44
7 154 41
8 164 51
9 168 57
10 166 49
11 159 47
12 164 46

读取

从剪切板读取

1
2
dat = read.table("clipboard",header=TRUE)
dat = read.table("clipboard")

从文本文件读取

1
2
dat = read.table("textdata.txt",header=TRUE)
dat = read.table("textdata.txt")

读取csv格式

1
X=read.csv("textdata.csv")

读取excel格式

  1. 下载包”readxl”
  2. 调用包library(readxl)
    3.读取文件
    1
    X=read_excel("data.xlsx")

定性变量分析

1
2
3
4
5
#将剪切板数据读入数据框中
> d=read.table("clipboard",header=T)
> head(d)#显示前6组数据
> attach(d)#绑定数据
> table(年龄)#以为列联表

单因素分析

1
2
3
4
#条形图
> barplot(table(年龄),col=1:7)
#饼图
> pie(table(结果))

两因素分析

1
2
3
4
#以性别分组的年龄条图
> barplot(table(年龄,性别),beside=T,col=1:7)
#以年龄分组的性别条图
> barplot(table(性别,年龄),beside=T,col=1:2)

三因素分析

1
2
3
4
#以年龄、性别排列的结果频数三维列联表
> ftable(年龄,性别,结果)
#以性别、年龄排列的结果频数三维列联表
> ftable(性别,年龄,结果)

当数据不使用时,记得解除绑定

1
> detach(d)

可视化

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
> X=read.table('clipboard',header=T);X
#按行做均值条形图
> barplot(apply(X,1,mean))
#修改横坐标位置
> barplot(apply(X,1,mean),las=3)

#按列做均值图条形
> barplot(apply(X,2,mean))
#按列做彩色均值图条形图
> barplot(apply(X,2,mean),col=1:8)

#按列做中位数条形图
> barplot(apply(X,2,median),col=1:8)
#按列做均值饼图
> pie(apply(X,2,mean))

# 垂直箱线图
> boxplot(X)
#水平箱线图
> boxplot(X,horizontal=T)

#简单星相图
> stars(X)
#带图例的星相图
> stars(X,key.loc=c(17,7))
#带图例度彩色星相图
> stars(X,key.loc=c(17,7),draw.segments=T)

#脸谱图
#加载aplpack包
> library(aplpack)
> faces(X)
#去掉变量1做脸谱图
> faaces(X[,-1])
#face(X,[,2:8])
#选择第1,5,6,9,18个观测的多元数据做脸谱图
> faces(X[C(1,5,6,9,18),])

#调和曲线图
> library(mvstats)#mvstats不是官方包
> plot.andrews(X)
> plot.andrews(X[C(1,5,6,9,18),])

多元相关及回归分析

1
2
3
x=c(171,175,159,155,152,158,154,164,168,166,159,164)#身高
y=c(57,64,41,38,35,44,41,51,57,49,47,46)#体重
plot(x,y)

画散点图
计算lxx=556.9 lyy=813 lxy=645.5
r=lxy/根号下(lxxlyy)=0.9593 高的正相关

相关系数计算函数cor()

cor(x,y=NULL,method=c(“pearson”,”kendall”,”spearman”))
x 数值向量、矩阵或数据框,y 空或数值向量、矩阵或数据框
method 计算方法,包括”pearson”,”kendall”或”spearman”三种,默认pearson

1
2
3
4
5
6
7
8
9
计算pearson相关系数:
> cor(x,y)

建立假设检验 H0:ρ=0,H1:ρ≠0,α=0.05
计算相关系数r的t值:tr=r-0/根号下(1-r^2/n-2)=10.74
> n=length(x)
> tr=r/sqrt((1-r^2)/(n-2));tr
计算t值和P值,作结论:
> cor.test(x,y)#相关系数假设检验
1
2
3
4
5
6
7
8
9
10
 Pearson's product-moment correlation

data: x and y
t = 10.743, df = 10, p-value = 8.21e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.8574875 0.9888163
sample estimates:
cor
0.9593031

由于p-value < 0.05,于是在α=0.05水准上拒绝H0,接受H1,可认为该人群身高与体重呈现正的线性关系。

建立直线回归方程

1
2
3
4
5
6
7
b=lxy(x,y)/lxy(x,x)#不过我这里报错没有lxy函数...lxy好像是离差平方和
a=mean(y)-b*mean(x)
c(a=a,b=b)

#散点图
plot(x,y)
lines(x,a+b*x)

1
2
3
4
#数据输入
d=read.table("clipboard",header=T);d
#拟合模型
m=lm(y~x,data=d);m

输出

1
2
3
4
5
6
Call:
lm(formula = y ~ x, data = d)

Coefficients:
(Intercept) x
-140.364 1.159

y估计=-140.364+1.159x
作回归直线

1
2
plot(y~x,data=d)#作散点图
abline(m)#添加回归线

模型的方差分析(ANOVA)

1
anova(fm)

1
summary(fm)#回归系数t检验

广义线性模型

Poisson分布族模型和拟Poisson分布族模型的使用方法为

1
2
glm(formula,family=poisson,data,...)
glm(formula,family=poisson(link=log),data,...)

(1)建立Poisson对数线性模型:

1
2
3
d=read.table("clipboard",header=T)
log=glm(y~x1+x2,family=poisson,data=d)#对数线性模型
summary(log)#检验结果

广义线性模型glm()的用法

1
glm(formula, family=gaussian, data, ...)

formula为公式,即为要拟合的模型,famliy为分布族,包括正态分布(gaussian)、二项分布(binomial)、泊松分布(poission)和伽玛分布(gamma),分布族还可以通过选项link=来指定使用的连接函数,data为可选择的数据框。

1
2
3
4
5
6
7
8
9
#(1)建立全变量logistic回归模型
d=read.table("clipboard",header=T)
logit<-glm(y~x1+x2+x3,family=binomial,data=d)#Logistic模型
summary(logit)#Logistic模型结果
#(2)逐步筛选变量logistic回归模型
#逐步筛选法变量选择
logit.step=step(logit)
#逐步筛选法变量选择结果
summar(logit.step)

1
2
d=read.table("clipboard",header=T)
anova(lm(Y~factor(A),data=d))#完全随机设计方差分析

判别分析

1
2
3
4
5
6
7
#线性判别分析
d=read_excel('data.xlsx','d')
boxplot(x1~G,d)
boxplot(x2~G,d)
t.test(x1~G,d)
t.test(x2~G,d)
summary(glm(G-1~x1+x2,family=binomial,d))#Logistic模型分析

线性判别分析函数lda的用法
lda(formula, data, …)
formula形如y~x1+x2+…的公式框架,data数据框

1
2
3
4
5
6
7
8
9
attach(d)
plot(x1,x2);text(x1,x2,G,adj=-0.5)#表示点所属类别G
library(MASS)
ld=lda(G~x1+x2);ld
lp=predict(ld)
G1=lp$class
data.frame(G,G1)
tab1=table(G,G1);tab1
sum(diag(prop.table(tab1)))

……没有实践写不下去