R 语言实战（第二版）
## part 4 高级方法
-------------第13章 广义线性模型------------------
#前面分析了线性模型中的回归和方差分析，前提都是假设因变量服从正态分布#广义线性模型对非正态因变量的分析进行扩展：如类别型变量、计数型变量（非负有限值）#glm函数，对于类别型因变量用logistic回归，计数型因变量用泊松回归#模型参数估计的推导依据的是最大似然估计（最大可能性估计），而非最小二乘法#1.logistic回归library(AER)data("Affairs") #婚外情数据summary(Affairs)#将affairs转化为二值型因子(是否出轨)Affairs$ynaffair[Affairs$affairs>0] <- 1Affairs$ynaffair[Affairs$affairs==0] <- 0Affairs$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章 主成分分析和因子分析------------------------------------
#主成分分析（PCA）：一种降维方法，将大量相关变量转化为少量不相关变量（即主成分）#主成分是观测变量的线性组合#探索性因子分析（EFA）：一种发现多变量潜在结构的方法，通过寻找一组更小的、潜在的结构（即因子）来解释已观测到的、显式的变量间的关系#EFA是假设生成工具，常用于帮助理解众多变量间的关系，常用于社会科学的理论研究#R基础函数：princompfactanal#psych包：library(psych)principal #PCA分析fa #因子分析fa.parallel #含平行分析的碎石图factor.plot #绘制结果fa.diagram #载荷矩阵scree #碎石图#主成分和因子分析步骤（代码以PCA分析为例）：#a.数据预处理：确保无NAdata("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章 时间序列----------------------------------
#前面分析的数据大多是横向数据：一个给点时间点的测量值#纵向数据：随时间变化反复测量变量值，即时序数据#时序数据解决的问题：对数据的描述和预测#常用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章 聚类分析------------------------------------------
#一种数据归约，不是一个精确的定义，因此聚类方法众多#常见的两类聚类：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$sizefit.km$centersaggregate(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章 分类--------------------------
#有监督机器学习分类方法：逻辑回归、决策树、随机森林、支持向量机、神经网络等#将全部数据分为一个训练集和一个验证集，训练集用于建立预测模型，验证集用于测试模型的准确性。#数据分析的一般流程是通过训练集建立模型，基于验证集调节参数，在测试集上评价模型。如Rattle包（数据分析GUI）将3个数据集比例默认划分为70/15/15pkgs <- 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$cptableplotcp(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.svmsvm.pred <- predict(fit.svm,na.omit(df.validate)) #与随机森林不同，SVM不允许NAsvm.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章 处理缺失数据的高级方法-------------------------
#大部分情况下，在处理真实数据之前，面对缺失值，要么删除含有缺失值的实例，要么替换缺失值。#处理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,FVIM::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")#只利用了两个变量的可用观测（忽略其他变量），因此每次计算都只用了不同的数据子集，可能会有一些扭曲的结果，不建议使用#简单插补#均值、中位数、众数等具体某值来替换缺失值。#优点：不会减少分析过程可用的样本量。#缺点：对非随机数据会产生有偏结果。尤其是缺失值多时，会低估标准差、曲解变量相关性等，不建议使用。