R 語言實戰（第二版）

## part 4 高級方法

-------------第13章 廣義線性模型------------------

#前面分析了線性模型中的回歸和方差分析，前提都是假設因變量服從正態分布
#廣義線性模型對非正態因變量的分析進行擴展：如類別型變量、計數型變量（非負有限值）
#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章 主成分分析和因子分析------------------------------------

#主成分分析（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章 時間序列----------------------------------

#前面分析的數據大多是橫向數據：一個給點時間點的測量值
#縱向數據：隨時間變化反復測量變量值，即時序數據
#時序數據解決的問題：對數據的描述和預測

#常用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$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章 分類--------------------------

#有監督機器學習分類方法：邏輯回歸、決策樹、隨機森林、支持向量機、神經網絡等
#將全部數據分為一個訓練集和一個驗證集，訓練集用於建立預測模型，驗證集用於測試模型的準確性。
#數據分析的一般流程是通過訓練集建立模型，基於驗證集調節參數，在測試集上評價模型。如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章 處理缺失數據的高級方法-------------------------

#大部分情況下，在處理真實數據之前，面對缺失值，要麽刪除含有缺失值的實例，要麽替換缺失值。
#處理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")
#只利用了兩個變量的可用觀測（忽略其他變量），因此每次計算都只用了不同的數據子集，可能會有一些扭曲的結果，不建議使用

#簡單插補
#均值、中位數、眾數等具體某值來替換缺失值。
#優點：不會減少分析過程可用的樣本量。
#缺點：對非隨機數據會產生有偏結果。尤其是缺失值多時，會低估標準差、曲解變量相關性等，不建議使用。

R 語言實戰-Part 4 筆記