代碼記錄

---

對經濟數據集中的缺失值進行處理

• 前言

這個數據集中存在大量的缺失，主要原因是某幾個年份的某些指標沒找到，或者干脆就是某些指標很難找，導致該指標數據的大批量丟失，更有甚者,由于要查找的年份(2017-2018)較近，所以缺失值巨多。

• 代碼

library(VIM) library(mice) library(readr) library(psych) library(fpc)workl <- "C:\\Users\\goatbishop\\Desktop\\data" setwd(workl) getwd()mydata <- read.csv("data01.csv", stringsAsFactors = F) head(mydata) str(mydata)####整理數據#### #將數據中有#DIV/0!的設置為缺失值 mydata$incomeRatUrbanRural[which(mydata$incomeRatUrbanRural == "#DIV/0!")] = NA mydata$advIndexIndusStruc[which(mydata$advIndexIndusStruc == "#DIV/0!")] = NAmydata$incomeRatUrbanRural <- as.numeric(mydata$incomeRatUrbanRural) mydata$advIndexIndusStruc <- as.numeric(mydata$advIndexIndusStruc)####繪制缺失值圖####png("缺失值圖1.png") aggr(mydata) dev.off()####考察相關性#### #剔除年份和城市兩個無關變量再計算相關系數 testdata <- mydata[, -c(1, 2)] #只保留完整數據行進行相關性分析 comdata01 <- testdata[which(complete.cases(testdata)== T), ]#有301條完整記錄 dim(comdata01) #總共710條記錄 dim(testdata)#計算相關系數 cor(comdata01) #相關系數大于0.6的返回True outcor <- cor(comdata01) > 0.6 write.csv(outcor, "outcor.csv") write.csv(cor(comdata01), "realcor.csv")#根據相關性分析，相關系數在0.6以上的為：#realGdpPerCapita與pro3industryGdp, nonAgriDevDegree, advIndexIndusStruc,urbanRat #pro3industryGdp和urbanRat與amoforeCapUtil #pro3industryGdp與urbanRat, amoforeCapUtil #urbanRat與numCollegeStu #socSecEmployExpend與amoforeCapUtil#但僅僅是作為參考，之后在回歸中可能會用到####現在對2017年和2018年的數據進行缺失值處理,驗證這兩年數據大量缺失的說法#### data2017 <- subset(mydata, year == 2017) data2018 <- subset(mydata, year == 2018) dim(data2017) dim(data2018)##缺失值為71個，全部缺失 length(which(complete.cases(data2017) == T))#僅有8個完整數據行 length(which(complete.cases(data2018) == T))#也就是說2*71-8 = 134個數據丟失，占了全部丟失數據的32.76%#再次驗證 missingdf <- as.data.frame(abs(is.na(mydata))) dim(missingdf) missingnum = c() for (i in c(1:10)) {tempdata <- missingdf[(71*(i-1)+1):(71*i), ]num <- sum(tempdata)missingnum <- c(missingnum, num) } missingnum #138 268 74 46 46 56 27 36 38 56 #可以看到2017年-2018年有大量的缺失值 #尤其是2017年####現在對城市數據進行考察####cityname = names(table(mydata$city)) missingnum = c() for ( i in c(1:length(cityname))) {tempdata <- subset(mydata, city == cityname[i])num <- length(which(complete.cases(tempdata) == F))missingnum <- c(missingnum, num) }citydf <- data.frame(cityname = cityname, missing =missingnum)#我們看到有些城市10年中沒有完整數據 #可能是因為這個城市的某個變量因為個各種原因沒有記錄，沒有展示出來####針對變量進行處理#### missingdf <- as.data.frame(abs(is.na(mydata))) varmissing <- apply(missingdf, 2, sum) varmissing #可以看到InvestOutputRate缺失值為105個， #urbanRat缺失值為122個 #regUrbanUnemployRat缺失值為115個 #socSecEmployExpend缺失值為100個 #這些變量的缺失值大于等于100 #剔除掉這些變量，再查看完整數據行的數量 testdata2 <- mydata[, -c(1, 2, 6, 13, 14, 17)] str(testdata2) length(which(complete.cases(testdata2))) #可以看到我們有 491行完整數據，比之前的301行####拉格朗日插值法(放棄)####varname <- colnames(mydata) mydataTemp <- mydatafor (i in c(1:length(cityname))) {for (y in c(3:length(varname))) {tempv <- mydataTemp[which(mydataTemp$city == cityname[i]), c(1, y)]numcom <- length(which(complete.cases(tempv)))if (numcom >= 7 & numcom <= 9) {newX <- tempv[which(!complete.cases(tempv)), 1]vecX <- tempv[which(complete.cases(tempv)), 1]vecY <- tempv[which(complete.cases(tempv)), 2]newY <- lagrange(vecX, vecY, newX)print(newX)print(newY)print("+++++++")tempv[which(!complete.cases(tempv)), 2] <- newYmydataTemp[which(mydataTemp$city == cityname[i]), c(1, y)] <- tempv} else {next} } }#由于多重插補法會涉及k階多項式，這會導致過擬合的線性 #并且在樣本邊界處會出現急速上升和急速下降的趨勢 #又由于我們的很多缺失值都處于2017至2018年所以不適合用多項式擬合 #這里，我們就建立簡單的一元回歸方程####一元回歸模型(OK)####varname <- colnames(mydata) mydataTemp <- mydata#針對每個城市的每一個變量，我們用一元線性回歸模型進行填補 #標準是：在10條數據中(以年份為自變量)，缺失數據必須小于3大于0，才能進行填補 #如果不滿足該條件，那么我們就跳過，之后再處理for (i in c(1:length(cityname))) {for (y in c(3:length(varname))) {tempv <- mydataTemp[which(mydataTemp$city == cityname[i]), c(1, y)]numcom <- length(which(complete.cases(tempv)))if (numcom >= 7 & numcom <= 9) {newX <- tempv[which(!complete.cases(tempv)), 1]vecX <- tempv[which(complete.cases(tempv)), 1]vecY <- tempv[which(complete.cases(tempv)), 2]newY <- myregression(vecX, vecY, newX)if (min(vecY) > 0 & min(newY) < 0 ) {print("預測異常...")next}print(newX)print(newY)print("+++++++")tempv[which(!complete.cases(tempv)), 2] <- newYmydataTemp[which(mydataTemp$city == cityname[i]), c(1, y)] <- tempv} else {next} } }#查看缺失值圖 #保存當前圖形參數設置(沒必要) opar <- par(no.readonly = T) aggr(mydata)png("一元回歸處理后的情況.png") aggr(mydataTemp) dev.off()#查看缺失值情況 #計算每個變量的缺失值數量 missingdf2 <- as.data.frame(abs(is.na(mydataTemp))) varmissing2 <- apply(missingdf2, 2, sum) varmissing2#可以看到下面3個變量的缺失值數量不小于100 #urbanRat：102 #regUrbanUnemployRat：101 #socSecEmployExpend：100write.csv(mydataTemp, "mydataTemp.csv")####主成分分析#####創建完整的數據集 tempdata <- mydataTemp[which(complete.cases(mydataTemp)) , ] rownames(tempdata) <- paste(tempdata$city, tempdata$year, sep = "")compeledata1 <- tempdata[, -c(1, 2)]summary(compeledata1) dim(compeledata1) #head(compeledata1)#中心化標準化數據 scaledata <- scale(compeledata1, center=T,scale=T)fa.parallel(scaledata, fa = "pc") #通過碎石土，選取特征值大于1的主成分，即前前5個主成分#主成分分析 pc <- principal(scaledata, nfactors = 5, scores = T) pc$loadings #好的，現在我們保存這個載荷矩陣 write.csv(pc$loadings, "因子載荷矩陣.csv")#主成分分析方法適用于變量之間存在較強相關性的數據 #如果原始數據相關性較弱，運用主成分分析后不能起到很好的降維作用 #從另一個角度看，主成分不僅可以進行降維，還可以對變量之間的相關性進行檢驗 #我們觀察因子載荷矩陣(在一個主成分中，對主成分影響較大的幾個變量之間相關性較強,且有時候在理論上具有一些相似的特性)： #urbanRat與ratIndexIndusStruc和greenCoverRat有較強的相關性(RC1) #socSecEmployExpend與pro2industryGdp,advIndexIndusStruc和amoforeCapUtil有較強的相關性(RC2) #regUrbanUnemployRat與numHospBed和dischaIndusWaste有較強的相關性(RC5) #PowConsumpPerGdpL與InvestOutputRate和proTechEduExpendGdp有較強的相關性(RC3)#pc$scores write.csv(pc$scores, "主成分得分.csv")####三個城市群進行回歸分析#####我們將城市分為3個城市群，分別缺失值與其相關的變量進行多元回歸分析#citys1 <- read.table("citys1.txt") #citys2 <- read.table("citys2.txt") #citys3 <- read.table("citys3.txt")#dim(citys1) #dim(citys2) #dim(citys3)mydataTemp2 <- mydataTemp numlib <- c(0, 27, 55, 71) cityname2 <- mydataTemp[1:71, 2]#urbanRat與ratIndexIndusStruc和greenCoverRat有較強的相關性(RC1) #socSecEmployExpend與pro2industryGdp,advIndexIndusStruc和amoforeCapUtil有較強的相關性(RC2) #regUrbanUnemployRat與numHospBed和dischaIndusWaste有較強的相關性(RC5) #PowConsumpPerGdpL與InvestOutputRate和proTechEduExpendGdp有較強的相關性(RC3) #foreTradeCoef與ratIndexIndusStruc和greenCoverRat有較強的相關性(RC1)#realGdpPerCapita與pro3industryGdp, nonAgriDevDegree, advIndexIndusStruc,urbanRat #pro3industryGdp與urbanRat, amoforeCapUtil #urbanRat與numCollegeStu #socSecEmployExpend與amoforeCapUtilfor (i in c(1:3)) {tempdatalm <- mydataTemp2[which(mydataTemp2$city %in% cityname2[(numlib[i]+1):numlib[i+1]]), ]#print((numlib[i]+1):numlib[i+1])#urbanRat與ratIndexIndusStruc和greenCoverRat有較強的相關性#且在前面的相關系數計算中我們得知：#pro3industryGdp和realGdpPerCapita與其相關系數較高trainU <- tempdatalm[which(complete.cases(tempdatalm$urbanRat)), ]testU <- tempdatalm[which(!complete.cases(tempdatalm$urbanRat)), ]if (dim(testU)[1] != 0) {print(paste("我執行了1", i, sep = "-"))lm1 <- lm(urbanRat ~ ratIndexIndusStruc+greenCoverRat+pro3industryGdp+realGdpPerCapita,data = trainU)pre1 <- predict(lm1, testU)#print(pre1)tempdatalm[which(!complete.cases(tempdatalm$urbanRat)), 13] <- pre1 }#socSecEmployExpend與pro2industryGdp,advIndexIndusStruc和amoforeCapUtiltrainU <- tempdatalm[which(complete.cases(tempdatalm$socSecEmployExpend)), ]testU <- tempdatalm[which(!complete.cases(tempdatalm$socSecEmployExpend)), ]if (dim(testU)[1] != 0) {print(paste("我執行了2", i, sep = "-"))lm1 <- lm(socSecEmployExpend ~ pro2industryGdp+advIndexIndusStruc+amoforeCapUtil,data = trainU)pre1 <- predict(lm1, testU)#print(pre1)tempdatalm[which(!complete.cases(tempdatalm$socSecEmployExpend)), 17] <- pre1 }#regUrbanUnemployRat與numHospBed和dischaIndusWaste#因為numHospBed和dischaIndusWaste同樣存在缺失，也是我們需要填補的對象#這里，我們試著用全部數據進行回歸，并采用逐步回歸法選取出參與回歸的自變量#proSocLab, incomeRatUrbanRural,realGdpPerCapita,InvestOutputRate,ratIndexIndusStructrainU <- tempdatalm[which(complete.cases(tempdatalm$regUrbanUnemployRat)), ]testU <- tempdatalm[which(!complete.cases(tempdatalm$regUrbanUnemployRat)), ]if (dim(testU)[1] != 0) {print(paste("我執行了3", i, sep = "-"))#numHospBed+dischaIndusWastelm1 <- lm(formula = regUrbanUnemployRat ~ realGdpPerCapita + proSocLab + InvestOutputRate + ratIndexIndusStruc + incomeRatUrbanRural,data = trainU[, -c(1, 2)])pre1 <- predict(lm1, testU)if (!all(complete.cases(pre1))) {print(pre1)}#print(pre1)tempdatalm[which(!complete.cases(tempdatalm$regUrbanUnemployRat)), 14] <- pre1 }#PowConsumpPerGdpL與InvestOutputRate和proTechEduExpendGdptrainU <- tempdatalm[which(complete.cases(tempdatalm$PowConsumpPerGdpL)), ]testU <- tempdatalm[which(!complete.cases(tempdatalm$PowConsumpPerGdpL)), ]if (dim(testU)[1] != 0) {print(paste("我執行了4", i, sep = "-"))lm1 <- lm(PowConsumpPerGdpL ~ InvestOutputRate+proTechEduExpendGdp,data = trainU)pre1 <- predict(lm1, testU)#print(pre1)#print(pre1)tempdatalm[which(!complete.cases(tempdatalm$PowConsumpPerGdpL)), 19] <- pre1 }#foreTradeCoef與ratIndexIndusStruc和greenCoverRattrainU <- tempdatalm[which(complete.cases(tempdatalm$foreTradeCoef)), ]testU <- tempdatalm[which(!complete.cases(tempdatalm$foreTradeCoef)), ]if (dim(testU)[1] != 0) {print(paste("我執行了5", i, sep = "-"))lm1 <- lm(foreTradeCoef ~ greenCoverRat+ratIndexIndusStruc,data = trainU)pre1 <- predict(lm1, testU)#print(pre1)#print(pre1)tempdatalm[which(!complete.cases(tempdatalm$foreTradeCoef)), 23] <- pre1 } mydataTemp2[which(mydataTemp2$city %in% cityname2[(numlib[i]+1):numlib[i+1]]), ] <- tempdatalm }#[1] "我執行了1-1" #[1] "我執行了2-1" #[1] "我執行了3-1" #[1] "我執行了1-2" #[1] "我執行了3-2" #[1] "我執行了5-2" #[1] "我執行了3-3" #[1] "我執行了4-3" #[1] "我執行了5-3"#aggr(mydataTemp)#第三次清理完畢 png("缺失值圖3.png") aggr(mydataTemp2) dev.off()length(which(complete.cases(mydataTemp2))) dim(mydataTemp2) #還有25行數據的缺失值未處理mydataTemp2[which(!complete.cases(mydataTemp2)), c(1:2, 15:16, 20, 22)] #可以看到在這些缺失值中： #眉山市的libCollect和numHospBed幾乎全部丟失(且和其他變量相關性很小) #需要注意的是numHospBed與其他變量相關性較強，但是之前我們利用這個變量作為多元回歸方程的自變量 #對其他缺失變量進行填補，所以為了防止混亂，就沒有在上面的多元回歸中進行填補 #資陽市的numCollegeStu有5年丟失(均值填補) #不同城市2017-2018年的dischaIndusWaste缺失(均值填補)write.csv(mydataTemp2, "第2次回歸處理后.csv")####均值處理####mydataTemp3 <- mydataTemp2#用資陽市其他年份的numCollegeStu進行填補 meantemp <- mean(mydataTemp3[which((mydataTemp3$city == "資陽市")), 22], na.rm = T) mydataTemp3[which(!complete.cases(mydataTemp3$numCollegeStu) & (mydataTemp3$city == "資陽市")), 22] <- meantemp#對dischaIndusWaste缺失的值進行均值填補namecitymissing <- mydataTemp3[which(!complete.cases(mydataTemp3$dischaIndusWaste)), 2] namecitymissingU <- unique(namecitymissing)for(ik in c(1:length(namecitymissingU))) {meantemp <- mean(mydataTemp3[which((mydataTemp3$city == namecitymissingU[ik])), 20], na.rm = T)mydataTemp3[which(!complete.cases(mydataTemp3$dischaIndusWaste) & (mydataTemp3$city == namecitymissingU[ik])), 20] <- meantemp }#利用2018年數據對眉山市的libCollect和numHospBed進行填補 #因為只有1年，所以均值就是2018年的值 meantemp <- mean(mydataTemp3[which((mydataTemp3$city == "眉山市")), 15], na.rm = T) mydataTemp3[which(!complete.cases(mydataTemp3$libCollect) & (mydataTemp3$city == "眉山市")), 15] <- meantempmeantemp <- mean(mydataTemp3[which((mydataTemp3$city == "眉山市")), 16], na.rm = T) mydataTemp3[which(!complete.cases(mydataTemp3$numHospBed) & (mydataTemp3$city == "眉山市")), 16] <- meantempmydataTemp3[which(!complete.cases(mydataTemp3)), ]#繪制缺失值圖png("最后一次缺失值圖.png") aggr(mydataTemp3) dev.off()write.csv(mydataTemp3, "mydataTemp3.csv")####密度聚類(放棄, 大量城市的數據擁擠在一起，有的則極其分散)#### dbscandata <- as.data.frame(pc$scores) dim(dbscandata) head(dbscandata)db <- dbscan(dbscandata, eps = 0.5, MinPts = 6) plotcluster(dbscandata, db$cluster) #plot(db, dbscandata)#db$cluster table(db$cluster)classdf <- data.frame(city = tempdata$city, class = db$cluster) usefuldata <- mydataTemp length(which(complete.cases(usefuldata) == T))dim(usefuldata)for (i in c(1:10)) {temp <- classdf$city[which(classdf$class == i)]temprtable <- table(temp)okdata <- temprtable[which(temprtable >= 4)]if (length(okdata) > 2) {print(okdata)tempdf <- usefuldata[which(usefuldata$city %in% names(okdata)),]#tempdf2 <- tempdf[, which(apply(is.na(tempdf), 2, sum) != 0)]assign(paste("tempdf", i, sep = ""), tempdf)print(i)} }#mice的method參數詳解 #"pmm"表示用預測的均值匹配 #"logreg"表示用邏輯回歸擬合 #“polyreg"表示多項式擬合 #“polr“表示采用比例優勢模型擬合等####函數定義####lagrange <- function(vectorX, vectorY, newX) {print(vectorX)print(vectorY)n <-length(vectorX)newY <- c()for (xnum in c(1:length(newX))) {x <- newX[xnum]lagr <- 0for (i in 1:n) {Li <- 1for (j in 1:n) {if (i!=j)Li <- Li*(x - vectorX[j])/(vectorX[i] - vectorX[j])}lagr <- Li*vectorY[i] + lagr}newY <- c(newY, lagr)}return(newY) }myregression <- function(vectorX, vectorY, newX) {print(vectorX)print(vectorY)lmtemp <- lm(vectorY ~ vectorX)print(lmtemp$coef)pre <- predict(lmtemp, data.frame(vectorX = newX))newY <- prereturn(newY) }

• 處理感想

走一步看一步把。

總結

以上是生活随笔為你收集整理的今日代码(200708)--缺失值处理的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。