我是靠谱客的博主 苗条睫毛膏,这篇文章主要介绍R 语言实战-Part 4 笔记R 语言实战(第二版),现在分享给大家,希望可以做个参考。


R 语言实战(第二版)

## part 4 高级方法

-------------第13章 广义线性模型------------------

复制代码
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
#前面分析了线性模型中的回归和方差分析,前提都是假设因变量服从正态分布 #广义线性模型对非正态因变量的分析进行扩展:如类别型变量、计数型变量(非负有限值) #glm函数,对于类别型因变量用logistic回归,计数型因变量用泊松回归 #模型参数估计的推导依据的是最大似然估计(最大可能性估计),而非最小二乘法 #1.logistic回归 library(AER) data("Affairs") #婚外情数据 summary(Affairs) #将affairs转化为二值型因子(是否出轨) Affairs$ynaffair[Affairs$affairs>0] <- 1 Affairs$ynaffair[Affairs$affairs==0] <- 0 Affairs$ynaffair <- factor(Affairs$ynaffair,levels = c(0,1),labels = c("No","Yes")) head(Affairs) table(Affairs$ynaffair) #将此ynaffair变量作为logistic回归的结果变量 fit.full <- glm(ynaffair~gender+age+yearsmarried+children+religiousness+education+occupation+rating, data = Affairs,family=binomial()) #family概率分布,binomial连接logit函数 summary(fit.full) #去除不显著的变量重新拟合模型,检验新模型是否拟合的好 fit.reduced <- glm(ynaffair~age+yearsmarried+religiousness+rating, data = Affairs,family = binomial()) summary(fit.reduced) #结果都显著 #两个模型嵌套(fit.reduced是fit.full的子集),可用anova进行比较,广义线性模型用卡方检验 anova(fit.reduced,fit.full,test = "Chisq") #卡方检验不显著,两个模型拟合一样,验证结果,可选择更简单的那个模型(fit.reduced)进行解释。 ##解释模型参数 #回归系数: coef(fit.reduced) exp(coef(fit.reduced)) #logistic回归用指数优势比(初始尺度)解释较好 #置信区间: exp(confint(fit.reduced)) #用predict函数评价预测变量对结果概率的影响 #创建一个数据探究某一预测变量(其他变量恒定)对结果概率的影响 #过度离势:观测到的因变量方差大于期望的二项分布的方差 #2.泊松回归 #一系列连续型或类别型预测变量来预测计数型结果变量 library(robust) data(breslow.dat) #癫痫数据 help(breslow.dat) summary(breslow.dat[c(6,7,8,10)]) #只关注这四个变量:两个协变量(base/age),重点关注药物治疗trt对癫痫发病数sumY的影响 opar <- par(no.readonly = T) par(mfrow=c(1,2)) attach(breslow.dat) hist(sumY,breaks = 20) boxplot(sumY~Trt) par(opar) #拟合泊松回归 fit <- glm(sumY~Base+Age+Trt,data=breslow.dat,family = poisson()) #family概率分布 summary(fit) #都显著 #解释模型参数 coef(fit) exp(coef(fit)) #因变量是以条件均值的对数形式来建模的,在初始尺度上解释回归系数比较容易 #注意:同logistic回归指数化参数相似,泊松模型中的指数化参数对因变量的影响是成倍增加的,而非线性相加。 #过度离势:因变量观测的方差大于泊松分布预测的方差,qcc包检验,略

-----------------------------第14章 主成分分析和因子分析------------------------------------

复制代码
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
#主成分分析(PCA):一种降维方法,将大量相关变量转化为少量不相关变量(即主成分) #主成分是观测变量的线性组合 #探索性因子分析(EFA):一种发现多变量潜在结构的方法,通过寻找一组更小的、潜在的结构(即因子)来解释已观测到的、显式的变量间的关系 #EFA是假设生成工具,常用于帮助理解众多变量间的关系,常用于社会科学的理论研究 #R基础函数: princomp factanal #psych包: library(psych) principal #PCA分析 fa #因子分析 fa.parallel #含平行分析的碎石图 factor.plot #绘制结果 fa.diagram #载荷矩阵 scree #碎石图 #主成分和因子分析步骤(代码以PCA分析为例): #a.数据预处理:确保无NA data("USJudgeRatings") #律师对美国高等法院法官的评价数据 help("USJudgeRatings") #b.选择因子模型:PCA/EFA(还需有估计因子模型方法,如最大似然估计) #简化11个变量信息,用PCA模型 #c.判断要选择的主成分/因子数目 #3种判别准则:特征值大于1准则(水平线),碎石检验(保留图形变化最大处之上的主成分),平行分析(虚线) fa.parallel(USJudgeRatings[,-1],fa="pc",n.iter = 100,show.legend = F) #三种准则表明选择一个主成分即可保留数据集大部分信息 #d.提取主成分/因子 principal(USJudgeRatings[,-1],nfactors = 1) #PC1包含了成分载荷(解释主成分含义),h2指成分公因子方差(方差解释度),u2指成分唯一性(即1-h2) #SS loadings指与主成分相关联的特征值,Proportion Var指每个主成分对整个数据集的解释程度 #e.旋转主成分/因子 #第二个数据例子:非原始数据,而是相关系数 data("Harman23.cor") #305个女孩的8个身体测量指标数据 head(Harman23.cor) fa.parallel(Harman23.cor$cov,n.obs = 302,fa="pc",n.iter = 100,show.legend = F) #建议两个主成分 principal(Harman23.cor$cov,nfactors = 2,rotate = "none") #当提取了多个成分时,对它们进行旋转可使结果(成分载荷阵)更具解释性 #正交旋转(成分保持不相关),如方差极大旋转 #斜交旋转(使成分变得相关) principal(Harman23.cor$cov,nfactors = 2,rotate = "varimax") #方差极大旋转 #PC变成RC,对载荷阵的列进行了去噪 #f.解释结果 #RC1第一主成分主要由前4个变量来解释(长度变量),RC2第二主成分主要由后4个变量来解释(容量变量) #g.计算主成分/因子得分 #获取每个观测在成分上的得分 principal(USJudgeRatings[,-1],nfactors = 1,scores = T) #非原始数据不能获得每个观测的主成分得分,只可计算主成分得分系数 rc <- principal(Harman23.cor$cov,nfactors = 2,rotate = "varimax") round(unclass(rc$weights),2) #unclass函数消除对象的类 #根据RC1和RC2的各个变量的系数计算PC1和PC2

-------------------------第15章 时间序列----------------------------------

复制代码
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
#前面分析的数据大多是横向数据:一个给点时间点的测量值 #纵向数据:随时间变化反复测量变量值,即时序数据 #时序数据解决的问题:对数据的描述和预测 #常用R包:stats,forecast,graphics,tseries等 #1.在R中生成时序对象 #时间序列对象:一种包括观测值、起始时间、终止时间和周期(如月、季度或年)的结构 scales <- c(18,33,23,56,78,90,51,24,34,56,78,23, 45,57,23,45,65,43,23,12,36,76,54,82) tscales <- stats::ts(scales,start=c(2003,1),frequency=12) #12表月度,1表年度,4表季度 tscales #获取对象信息 plot(tscales) start(tscales) end(tscales) frequency(tscales) #对象取子集 stats::window(tscales,start=c(2003,5),end=c(2004,6)) #2.时序的平滑化和季节性分解 #时序数据常有随机或误差成分,要撇开波动,画平滑曲线。采用简单移动平均方法,如居中移动(前后两个数取均值) library(forecast) par(mfrow=c(2,2)) ylim <- c(min(Nile),max(Nile)) plot(Nile) plot(forecast::ma(Nile,3)) #拟合一个简单的移动平均模型 plot(ma(Nile,7)) plot(ma(Nile,15)) #尝试多个不同的k值,找出最能反映数据规律的k,避免过平滑或欠平滑 #对于间隔大于1的时序数据(即存在季节性因子,如月度、季度数据),需要季节性分解为趋势因子(反映长期变化)、季节性因子(反映周期性变化)和随机因子(误差) #数据的分解通过两个模型:相加模型或相乘模型,具体分解步骤略 #3.指数预测模型 #用来预测时序未来值得最常用模型,短期预测能力较好。单指数、双指数、三指数模型选用的因子不同,略 #4.ARIMA预测模型 #预测值为最近的真实值和最近的预测误差组成的线性函数,比较复杂 #主要用于拟合具有平稳性的时间序列 #这些预测模型都是用的向外推断的思想,即假定未来条件和现在是相似的。 #实际上很多因素影响时序的趋势和模式,时间跨度越大,不确定越大。

-----------------------第16章 聚类分析------------------------------------------

复制代码
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
#一种数据归约,不是一个精确的定义,因此聚类方法众多 #常见的两类聚类:a.层次聚类:每个观测自成一类,每次两两合并,直至所有类并成一类。贪婪的分层算法,只能分给一个类 #b.划分聚类:首先指定类个数K,观测值随机分成K类,再重新聚合成类 #层次聚类算法:单联动(single linkage)、全联动(complete ~)、平均联动(average ~)、质心(centroid)、ward等 #划分聚类算法:K均值(K-means)、围绕中心点划分(PAM) #1.聚类一般步骤 #选择合适的变量 #缩放数据:如scale归一化 apply(data, 2, function(x){(x-mean(x)/sd(x))}) apply(data, 2, function(x){x/max(x)}) apply(data, 2, function(x){(x-mean(x))/mad(x)}) #mad平均绝对偏差 #寻找异常点 #计算距离:最常见欧氏(欧几里得)距离,常用于连续型数据。距离越大,异质性越大 #选择聚类算法:小样本量层次聚类合适,大样本量划分聚类 #获得聚类方法 #确定类数目 #获得最终聚类方案 #结果可视化 #解读类 #验证结果 data(nutrient,package = "flexclust") head(nutrient,4) d <- dist(nutrient) #计算距离 as.matrix(d)[1:4,1:4] #2.层次聚类 row.names(nutrient) <- tolower(row.names(nutrient)) nutrient.scaled <- scale(nutrient) d <- dist(nutrient.scaled) fit.average <- hclust(d,method = "average") #平均联动层次聚类 plot(fit.average,hang = -1,cex=0.8) #hang展示观测值标签在0下 #读图从下按层次往上读 #确定聚类数目 library(NbClust) devAskNewPage(ask = T) nc <- NbClust(nutrient.scaled,distance = "euclidean", min.nc = 2,max.nc = 15,method = "average") #平均联动聚类 table(nc$Best.n[1,]) table(nc$Best.nc[1,]) #聚类数2,3,5等赞同最多 barplot(table(nc$Best.nc[1,])) #一张张绘完图后,会给回聚类建议数目 devAskNewPage(ask = F) #关掉一张张展示 #获取最终聚类方案 clusters <- cutree(fit.average,k=5) table(clusters) aggregate(nutrient,by=list(cluster=clusters),median) aggregate(as.data.frame(nutrient.scaled),by=list(cluster=clusters),median) plot(fit.average,hang = -1,cex=0.8) rect.hclust(fit.average,k=5) #将5类框起来 #3.划分聚类 ##1)K means #绘图函数:根据图中弯曲选择适当类数 wssplot <- function(data,nc=15,seed=1234){ #nc考虑最大聚类数 wss <- (nrow(data)-1)*sum(apply(data,2,var)) #var方差函数 for (i in 2:nc) { set.seed(seed) wss[i] <- sum(kmeans(data,centers = i)$withinss) } plot(1:nc,wss,type = "b") } data(wine,package = "rattle") head(wine) df <- scale(wine[-1]) #确定聚类个数 wssplot(df) #三类之后下降速度减弱,因此三类可能拟合得好 library(NbClust) set.seed(1234) #每次随机选择k类,因此要用这个复现结果 devAskNewPage(ask = T) nc <- NbClust(df,min.nc = 2,max.nc = 15,method = "kmeans") #kmeans聚类 table(nc$Best.n[1,]) barplot(table(nc$Best.nc[1,])) #进行K均值聚类分析 set.seed(1234) fit.km <- kmeans(df,3,nstart = 25) #尝试25个初始配置 fit.km$size fit.km$centers aggregate(wine[-1],by=list(cluster=fit.km$cluster),mean) ##2)围绕中心点的划分 #K means聚类基于均值,对异常值敏感。围绕中心点划分(PAM)更为稳健 library(cluster) set.seed(1234) fit.pam <- pam(wine[-1],k=3,stand = T) #聚成3类,数据标准化 fit.pam$medoids #输出中心点 clusplot(fit.pam) #画出聚类方案

-----------------------第17章 分类--------------------------

复制代码
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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
#有监督机器学习分类方法:逻辑回归、决策树、随机森林、支持向量机、神经网络等 #将全部数据分为一个训练集和一个验证集,训练集用于建立预测模型,验证集用于测试模型的准确性。 #数据分析的一般流程是通过训练集建立模型,基于验证集调节参数,在测试集上评价模型。如Rattle包(数据分析GUI)将3个数据集比例默认划分为70/15/15 pkgs <- c("rpart","rpart.plot","party", #决策树模型及其可视化 "randomForest", #拟合随机森林 "e1071") #构造支持向量机 #install.packages(pkgs,depend=T) #1.数据准备 loc <- "http://archive.ics.uci.edu/ml/machine-learning-databases/" ds <- "breast-cancer-wisconsin/breast-cancer-wisconsin.data" url <- paste(loc,ds,sep = "") url <- "https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.data" breast <- read.table(url,sep = ",",header = F,na.strings = "?") names(breast) <- c("ID","clumpThickness","sizeUniformity","shapeUniformity","maginalAdhesion", "singleEpithelialCellSize","bareNuclei","blandChromatin", "normalNucleoli","mitosis","class") df <- breast[-1] df$class <- factor(df$class,levels = c(2,4),labels = c("benign","malignant")) #良性、恶性 set.seed(1234) train <- sample(nrow(df),0.7*nrow(df)) #训练集随机选70% df.train <- df[train,] df.validate <- df[-train,] #验证集随机选30% table(df.train$class) table(df.validate$class) #2.逻辑回归 #广义线性模型的一种,基础函数glm #拟合逻辑回归 fit.logit <- glm(class~.,data=df.train,family = binomial()) summary(fit.logit) #未过P值的,一般不会纳入最终模型 #对训练集外样本单元分类 prob <- predict(fit.logit,df.validate,type = "response") #预测恶性肿瘤的概率 logit.pred <- factor(prob>0.5,levels = c(FALSE,TRUE), labels = c("benign","malignant")) #评估预测准确性 logit.perf <- table(df.validate$class,logit.pred, dnn = c("Actual","Predicted")) #交叉表 logit.perf #准确率(76+118)/200 #逐步逻辑回归移除/增加多余变量(更小的AIC值) logit.fit.reduced <- step(fit.logit) #3.决策树 #对预测变量进行二元分离,构造一棵可用于预测新样本所属类别的树 ##1)经典决策树 #通常会得到一棵过大的树,出现过拟合,对于训练集外单元的分类性能较差,可用10折交叉验证法剪枝后的树进行预测 library(rpart) set.seed(1234) #生成树 dtree <- rpart(class~.,data = df.train,method = "class", parms = list(split="information")) dtree$cptable plotcp(dtree) #剪枝 dtree.pruned <- prune(dtree,cp=.0125) library(rpart.plot) #画出最终决策树 prp(dtree.pruned,type = 2,extra = 104,fallen.leaves = T) #对训练集外样本单元分类 dtree.pred <- predict(dtree.pruned,df.validate,type = "class") dtree.perf <- table(df.validate$class,dtree.pred, dnn = c("Actual","Predicted")) dtree.perf #准确率96% #对于水平数很多或缺失值很多的预测变量,决策树可能会有偏 ##2)条件推断树 #变量和分割的选取是基于显著性检验(置换检验)的,而非纯净度或同质性一类的度量 #可以不剪枝,生成过程更自动化 library(party) fit.ctree <- ctree(class~.,data = df.train) plot(fit.ctree) ctree.pred <- predict(fit.ctree,df.validate,type="response") ctree.pref <- table(df.validate$class,ctree.pred, dnn = c("Actual","Predicted")) ctree.pref #4.随机森林 #组成式的有监督学习方法:同时生成多个预测模型,将模型结果汇总以提升分类准确率 #所有大量决策树(randomForest包由传统决策树生成)预测类别中的众数类别即为随机森林所预测的这一样本单元的类别 #无法获得验证集时,是随机森林的优势。可计算变量的相对重要程度 #优点:分类准确率更高,可处理大规模问题(多样本单元、多变量),可处理含缺失值的训练集,可处理变量多于样本的数据 #缺点:分类方法较难理解和表达,需存储整个随机森林 library(randomForest) set.seed(1234) #生成森林(默认500棵树) fit.forest <- randomForest(class~.,data=df.train, na.action=na.roughfix, #NA替换(中位值/众数) importance=T) #变量重要性 fit.forest #给出变量重要性 importance(fit.forest,type = 2) #节点不纯度 #对训练集外样本点分类 forest.pred <- predict(fit.forest,df.validate) forest.perf <- table(df.validate$class,forest.pred, dnn = c("Actual","Predicted")) forest.perf #当预测变量间高度相关时,基于条件推断树的随机森林可能效果更好,party::cforest函数 #5.支持向量机SVM #可输出较准确的预测结果,模型基于较优雅的数学理论 #SVM旨在多维空间中找到一个能将全部样本单元分成两类的最优平面,这一平面使两类中距离最近的点的间距尽可能大 #在间距边界上的点称为支持向量 #N维空间(n个变量),最优超平面为N-1维 library(e1071) set.seed(1234) fit.svm <- svm(class~.,data=df.train) #默认将数据标准化(均值为0,标准差为1) fit.svm svm.pred <- predict(fit.svm,na.omit(df.validate)) #与随机森林不同,SVM不允许NA svm.perf <- table(na.omit(df.validate)$class,svm.pred, dnn = c("Actual","Predicted")) svm.perf #svm()默认通过径向基函数(RBF)将样本单元投射到高维空间 #svm拟合样本时,gamma和cost两个参数可能影响模型结果,两参数恒为正数 #选择调和参数 set.seed(1234) tuned <- tune.svm(class~.,data=df.train, gamma = 10^(-6:1), #尝试8个不同值(值越大训练样本到达范围越广) cost = 10^(-10:10)) #尝试21个成本参数(值越大惩罚越大) tuned #共拟合了168次,输出最优模型 #用这些参数拟合模型 fit.svm <- svm(class~.,data=df.train,gamma=0.01,cost=1) #评估交叉验证表现 svm.pred <- predict(fit.svm,na.omit(df.validate)) svm.perf <- table(na.omit(df.validate)$class,svm.pred, dnn = c("Actual","Predicted")) svm.perf #准确率有所提高 #同随机森林一样,SVM也可用于变量数远多于样本单元数的问题(如基因表达矩阵),但是分类准则较难理解和表述 #SVM本质上是一个黑盒子,在大量样本建模时不如随机森林 #6.选择预测效果最好的解 #分类器是否总能正确划分样本单元?最常用正确率来衡量 #预测一个分类器的准确性度量:敏感度、特异性、正例命中率、负例命中率、准确率 #评估二分类准确性 performance <- function(table,n=2){ if(!all(dim(table)==c(2,2))) stop("must be a 2X2 table!") #得到频数 tn=table[1,1] fp=table[1,2] fn=table[2,1] tp=table[2,2] #计算统计量 sensitivity=tp/(tp+fn) #敏感度 specificity=tn/(tn+fn) #特异性 ppp=tp/(tp+fp) #正例命中率 npp=tn/(tn+fn) #负例命中率 hitrate=(tp+tn)/(tp+tn+fp+fn) #正确率 #输出结果 result <- paste("sensitivity=",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.pref) #条件推断树 performance(forest.perf) #随机森林 performance(svm.perf) #支持向量机 #可从特异性和敏感度的权衡中提高分类性能。如可通过分类器的特异性来增加其敏感性 #ROC曲线可对一个区间内的阈值画出特异性和敏感性的关系,从而针对特定问题选择两者的最佳组合 #方法选择: #如果复杂、黑箱式的方法(如随机森林和SVM)相比简单方法(如逻辑回归和决策树)在预测效果上并没有显著提升,一般会选择较简单的方法

----------------第18章 处理缺失数据的高级方法-------------------------

复制代码
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
#大部分情况下,在处理真实数据之前,面对缺失值,要么删除含有缺失值的实例,要么替换缺失值。 #处理NA的步骤:识别缺失值,检查原因,删除或插补 #缺失值的分类: #①完全随机缺失(MCAR):该变量缺失值与其他变量/观测都不相关 #②随机缺失(MAR):该变量缺失值与其他变量相关,与它自己的未观测值不相关 #③非随机缺失(NMAR):非上述两类。 #大部分处理缺失值的方法都假定数据是随机缺失 #1.识别缺失值 #基础方法: data(sleep,package="VIM") sleep[complete.cases(sleep),] #列出没有缺失值的行 sleep[!complete.cases(sleep),] sum(is.na(sleep$Dream)) mean(is.na(sleep$Dream)) #该变量缺失值比例 mean(!complete.cases(sleep)) mean(is.na(sleep)) #为何与上不同结果? #complete.cases函数将NA和NaN视为缺失值,正负Inf视为有效值 #识别缺失值必须要is.na,myvar==NA无用 #2.探索缺失值模式 #图表查看缺失值: library(mice) md.pattern(sleep) library(VIM) aggr(sleep,pro=F,numbers=T) #比例和数值标签,默认分别T,F VIM::matrixplot(sleep) #热图的形式来展示缺失值(红色),其他颜色深浅代表值大小 #此函数还可图形交互,使变量重排 #缺失值散点图 VIM::marginplot(sleep[c("Gest","Dream")],pch=c(20), col = c("darkgray","red","blue")) #主体是变量完整数据散点图 #用相关性探索缺失值:变量间的缺失值以及缺失值和变量之间是否相关呢? x <- as.data.frame(abs(is.na(sleep))) head(sleep,5) head(x,5) y <- x[which(apply(x, 2, sum)>0)] #提取含缺失值的变量 head(y,5) cor(y) cor(sleep,y,use="pairwise.complete.obs") #表中相关系数并不是很大,表明数据是MCAR的可能性较小,更可能为MAR #三种流行的缺失值处理方法:推理法、行删除法、多重插补 #3.行删除(完整实例分析) mydata <- data[complete.cases(data),] mydata <- na.omit(data) options(digits = 1) cor(na.omit(sleep)) #为何有两位小数? cor(sleep,use = "complete.obs") #等同以上 fit <- lm(Dream~Span+Gest,data=na.omit(sleep)) summary(fit) fit <- lm(Dream~Span+Gest,data=sleep) #与变量相关的缺失值才删除 summary(fit) #4.多重插补MI #一种基于重复模拟的处理缺失值方法。通过生成3-10个数据集,每个模拟数据集中,缺失值用蒙特卡洛方法来填补 #基于mice包的处理过程: library(mice) imp <- mice(data,m) #m个插补数据集,默认为5 fit <- with(imp,analysis) #插补数据集的统计方法,如lm,glm,gam,nbrm等 pooled <- pool(fit) #m个统计分析平均结果 summary() data("sleep") imp <- mice(sleep,seed = 1234) fit <- with(imp,lm(Dream~Span+Gest)) pooled <- pool(fit) summary(pooled) imp #对象汇总 imp$imp$Dream #查看子成分,看到实际的插补值 complete(imp,action = 3) #观察m个插补数据集中的任意一个 #5.其他方法 #成对删除 cor(sleep,use = "pairwise.complete.obs") #只利用了两个变量的可用观测(忽略其他变量),因此每次计算都只用了不同的数据子集,可能会有一些扭曲的结果,不建议使用 #简单插补 #均值、中位数、众数等具体某值来替换缺失值。 #优点:不会减少分析过程可用的样本量。 #缺点:对非随机数据会产生有偏结果。尤其是缺失值多时,会低估标准差、曲解变量相关性等,不建议使用。

转载于:https://www.cnblogs.com/jessepeng/p/10657613.html

最后

以上就是苗条睫毛膏最近收集整理的关于R 语言实战-Part 4 笔记R 语言实战(第二版)的全部内容,更多相关R内容请搜索靠谱客的其他文章。

本图文内容来源于网友提供,作为学习参考使用,或来自网络收集整理,版权属于原作者所有。
点赞(81)

评论列表共有 0 条评论

立即
投稿
返回
顶部