今日代码(200924)--缺失值处理

缺失值处理

对110个城市10年的数据进行缺失值处理。

knitr::opts_chunk$set(echo = T, message = FALSE, warning = FALSE)

导包

library(VIM)
library(mice)
library(readr)
library(psych)
library(fpc)
library(lattice)
library(MASS)

自定义函数


#统计行/列 缺失值函数
countNaN <- function(myline) {return(sum(is.na(myline)))
}#定位函数
#dataF只能传递矩阵或者数据框
myPosition <- function(dataF, flag = T) {posmtr <- matrix(0, ncol = 2)colnames(posmtr) <- c("rowNum", "colNum")numberCol <- dim(dataF)[2]numberRow <- dim(dataF)[1]for (item in c(1:numberCol)) {targetIndex <- which(dataF[, item])lenTar <- length(targetIndex)tempMatrix <- data.frame(rowNum = targetIndex, colNum = rep(item, lenTar))posmtr <- rbind(posmtr, tempMatrix)}posmtr <- posmtr[-1, ]return(posmtr)
}#回归函数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)
}

读取数据

mydata <- read.csv("../RawData0923.csv")dim(mydata) #1100 24head(mydata)
#str(mydata)

可以看到共有1100条数据,共有110个城市,十年的数据,且原本应该为数值型的变量，却显示为字符型，说明数据中可能存在字符型数据。

整理数据

将数据中有#DIV/0!的设置为缺失值

mydata2 <- mydatafor (item in c(2:length(mydata2))) {data_type <- class(mydata2[, item])if (data_type == "character") {temp_dataL <- mydata2[, item]temp_dataL[which(temp_dataL == "#DIV/0!")] = NAmydata2[, item] <- as.numeric(temp_dataL)}
}str(mydata)
str(mydata2)

查看整体的缺失值情况：

length(which(complete.cases(mydata2))) #524

可以看到总共有524个观测是完整的。

绘制缺失值图：

#png("../images/缺失值图1.png")
aggr(mydata2, prop = F, numbers = T)
#dev.off()

查看各个变量的缺失值情况:


cNa <- apply(mydata2, 2, countNaN)
cNa[which(cNa > 0)]

可以看到, 变量urbanRat缺失了398个数据, incomeRatUrbanRural缺失了154个数据, socSecEmployExpend缺失了89个数据, numCollegeStu缺失了78个数据.

考察变量间的相关性

#剔除年份和城市两个无关变量再计算相关系数
testdata <- mydata2[, -c(1, 2)]
#只保留完整数据行进行相关性分析
comdata01 <- testdata[which(complete.cases(testdata)== T), ]
#有524条完整记录
dim(comdata01) #524  22#计算相关系数
#cor(comdata01)#相关系数大于0.6的返回True
outcor <- cor(comdata01) > 0.6
write.csv(outcor, "../output/outcor.csv")
write.csv(cor(comdata01), "../output/realcor.csv")PosCor <- myPosition(outcor, T)new_cor_big_name <- data.frame(rowname = colnames(outcor)[PosCor[, 1]], colname = colnames(outcor)[PosCor[, 2]])

由结果可知，相关系数较高的变量为：
pro3industryGdp realGdpPerCapita;
nonAgriDevDegree realGdpPerCapita;
urbanRat realGdpPerCapita;
libCollect realGdpPerCapita;
advIndexIndusStruc pro3industryGdp
urbanRat pro3industryGdp
libCollect pro3industryGdp
numHospBed pro3industryGdp
numHospBed advIndexIndusStruc
urbanRat nonAgriDevDegree
libCollect urbanRat
libCollect amoforeCapUtil
socSecEmployExpend amoforeCapUtil

查看每年的缺失值数：


complete_data <- c()
for (item in c(2018:2009)) {dataYear <- subset(mydata2, year == item)#print(dim(dataYear))len_complete <- length(which(complete.cases(dataYear) == T))complete_data <- c(complete_data, len_complete)
}print(complete_data)

可以看到2018年只有41条数据是完整的，2016年只有42条是完整的，其他年份也分别有大量数据丢失。

再次查看每年具体有多少数据丢失：

missingdf <- as.data.frame(abs(is.na(mydata2)))
#dim(missingdf)
missingnum = c()
for (i in c(1:10)) {tempdata <- missingdf[(110*(i-1)+1):(110*i), ]num <- sum(tempdata)missingnum <- c(missingnum, num)
}
print(missingnum)

可以看到2018年有159个缺失值, 2016年有129个缺失值, 2015年有128个缺失值, 2014年有126个缺失值。

现在对城市数据的缺失值进行考察：

library(knitr)cityname <- names(table(mydata2$city))
missingnum <- c()
for ( i in c(1:length(cityname))) {tempdata <- subset(mydata2, city == cityname[i])num <- length(which(complete.cases(tempdata) == F))missingnum <- c(missingnum, num)
}
citydf <- data.frame(cityname = cityname, missing =missingnum)
#print(citydf)
kable(citydf[which(citydf$missing >= 1), ], caption = "有数据缺失的城市")

我们看到有些城市在10年中没有一条完整的数据,可能是因为这个城市的某个变量因为个各种原因没有被记录.

针对变量进行处理


missingdf <- as.data.frame(abs(is.na(mydata2)))
varmissing <- apply(missingdf, 2, sum)
varmissing

可以看到：
incomeRatUrbanRural缺失了154条数据
urbanRat缺失了398条数据
numCollegeStu缺失了78条数据
socSecEmployExpend缺失了89条数据

现在，我们对缺失值数大于等于50的变量进行删除：

delete_value_number <- which(varmissing >= 50)testdata2 <- mydata2[, -delete_value_number]
#str(testdata2)
length(which(complete.cases(testdata2))) #951

可以看到当我们删除这4列时,我们有951行完整数据，比之前的524行多了427行.

一元回归模型

varname <- colnames(mydata2)
mydataTemp <- mydata2
yearColNum <- which("year" == varname)#针对每个城市的每一个变量，我们用一元线性回归模型进行填补
#标准是：在10条数据中(以年份为自变量)，缺失数据必须小于等于4大于0，才能进行填补
#如果不满足该条件，那么我们就跳过，之后再处理for (i in c(1:length(cityname))) {for (y in c(3:length(varname))) {tempv <- mydataTemp[which(mydataTemp$city == cityname[i]), c(yearColNum, y)]numcom <- length(which(complete.cases(tempv)))if (numcom >= 6 & 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(yearColNum, y)] <- tempv} else {next} }
}

查看当前缺失值图：

aggr(mydataTemp, prop = F, numbers = T)

再次计算每个变量的缺失值数量：

missingdf2 <- as.data.frame(abs(is.na(mydataTemp)))
varmissing2 <- apply(missingdf2, 2, sum)
varmissing2

可以看到各个变量的缺失值数目均减少, 只有urbanRat变量含有100个以上的缺失值,而numCollegeStu有35个缺失值, socSecEmployExpend有25,nonAgriDevDegree有23个缺失值,ratIndexIndusStruc有23个缺失值,对于这些,我们需要进行进一步处理。

输出部分结果：

write.csv(mydataTemp, "../output/mydataTemp0923.csv")

主成分分析

创建完整的数据集：

tempdata <- mydataTemp[which(complete.cases(mydataTemp)) , ]
rownames(tempdata) <- paste(tempdata$city, tempdata$year, sep = "")#head(tempdata, 10)
#head(mydata2, 10)
compeledata1 <- tempdata[, -c(1, 2)]#summary(compeledata1)
#dim(compeledata1)#head(tempdata$city, 10)
head(paste(tempdata$city, tempdata$year, sep = ""), 10)

中心化标准化数据：

scaledata <- scale(compeledata1, center=T,scale=T)
fa.parallel(scaledata, fa = "pc")

通过碎石土，选取特征值大于1的主成分，即前前5个主成分。

主成分分析:

pc <- principal(scaledata, nfactors = 5, scores = T)
pc$loadings

主成分分析方法适用于变量之间存在较强相关性的数据,如果原始数据相关性较弱，运用主成分分析后不能起到很好的降维作用,从另一个角度看，主成分不仅可以进行降维，还可以对变量之间的相关性进行检验,我们观察因子载荷矩阵(在一个主成分中，对主成分影响较大的几个变量之间相关性较强,且有时候在理论上具有一些相似的特性)：
realGdpPerCapita与nonAgriDevDegree,ratIndexIndusStruc,amoforeCapUtil,
foreTradeCoef,libCollect和urbanRat有较强的相关性(RC1)
pro2industryGdp与pro3industryGdp,advIndexIndusStruc,socSecEmployExpend和numHospBed有较强的相关性(RC2)
gdpGrowthRate与proSocLab, incomeRatUrbanRural和dischaIndusWaste有较强的相关性(RC3)
proTechEduExpendGdp与PowConsumpPerGdpL有较强的相关性(RC4)
foreTradeCoef与numCollegeStu,PowConsumpPerGdpL有较强的相关性(RC5)

现在我们保存这个载荷矩阵:

write.csv(pc$loadings, "../output/因子载荷矩阵0923.csv")

三个城市群进行回归分析：


cityColNum <- which("city" == varname)
mydataTemp2 <- mydataTemp
numlib <- c(0, 41, 77, 110)
cityname2 <- mydataTemp[1:110, cityColNum]urbanRatNameIndex <- which(colnames(mydataTemp) == "urbanRat")
numCollegeStuNameIndex <- which(colnames(mydataTemp) == "numCollegeStu")
socSecEmployExpendNameIndex <- which(colnames(mydataTemp) == "socSecEmployExpend")
nonAgriDevDegreeNameIndex <- which(colnames(mydataTemp) == "nonAgriDevDegree")
ratIndexIndusStrucNameIndex <- which(colnames(mydataTemp) == "ratIndexIndusStruc")# realGdpPerCapita与nonAgriDevDegree,ratIndexIndusStruc,amoforeCapUtil,
# foreTradeCoef,libCollect和urbanRat有较强的相关性(RC1)
# pro2industryGdp与pro3industryGdp,advIndexIndusStruc,socSecEmployExpend和numHospBed有较强的相关性(RC2)
# gdpGrowthRate与proSocLab, incomeRatUrbanRural和dischaIndusWaste有较强的相关性(RC3)
# proTechEduExpendGdp与PowConsumpPerGdpL有较强的相关性(RC4)
# foreTradeCoef与numCollegeStu,PowConsumpPerGdpL有较强的相关性(RC5)# urbanRat变量含有100个以上的缺失值
# numCollegeStu有35个缺失值
# socSecEmployExpend有25个缺失值
# nonAgriDevDegree有23个缺失值
# ratIndexIndusStruc有23个缺失值for (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])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 ~ amoforeCapUtil + libCollect + foreTradeCoef + realGdpPerCapita,data = trainU)pre1 <- predict(lm1, testU)#print(pre1)#赋值tempdatalm[which(!complete.cases(tempdatalm$urbanRat)), urbanRatNameIndex] <- pre1    }#socSecEmployExpend与pro2industryGdp,advIndexIndusStruc和amoforeCapUtiltrainU <- tempdatalm[which(complete.cases(tempdatalm$numCollegeStu)), ]testU <- tempdatalm[which(!complete.cases(tempdatalm$numCollegeStu)), ]if (dim(testU)[1] != 0) {print(paste("我执行了2", i, sep = "-"))lm1 <- lm(numCollegeStu ~ PowConsumpPerGdpL + foreTradeCoef,data = trainU)pre1 <- predict(lm1, testU)#print(pre1)tempdatalm[which(!complete.cases(tempdatalm$numCollegeStu)), numCollegeStuNameIndex] <- pre1    }trainU <- tempdatalm[which(complete.cases(tempdatalm$socSecEmployExpend)), ]testU <- tempdatalm[which(!complete.cases(tempdatalm$socSecEmployExpend)), ]if (dim(testU)[1] != 0) {print(paste("我执行了3", i, sep = "-"))#pro2industryGdp与pro3industryGdp,advIndexIndusStruc,socSecEmployExpend和numHospBedlm1 <- lm(formula = socSecEmployExpend ~ numHospBed + advIndexIndusStruc + pro2industryGdp + pro3industryGdp,data = trainU)pre1 <- predict(lm1, testU)if (!all(complete.cases(pre1))) {print(pre1)}#print(pre1)tempdatalm[which(!complete.cases(tempdatalm$socSecEmployExpend)), socSecEmployExpendNameIndex] <- pre1    }trainU <- tempdatalm[which(complete.cases(tempdatalm$nonAgriDevDegree)), ]testU <- tempdatalm[which(!complete.cases(tempdatalm$nonAgriDevDegree)), ]if (dim(testU)[1] != 0) {print(paste("我执行了4", i, sep = "-"))lm1 <- lm(nonAgriDevDegree ~ amoforeCapUtil + libCollect + foreTradeCoef + realGdpPerCapita,data = trainU)pre1 <- predict(lm1, testU)#print(pre1)tempdatalm[which(!complete.cases(tempdatalm$nonAgriDevDegree)), nonAgriDevDegreeNameIndex] <- pre1    }trainU <- tempdatalm[which(complete.cases(tempdatalm$ratIndexIndusStruc)), ]testU <- tempdatalm[which(!complete.cases(tempdatalm$ratIndexIndusStruc)), ]if (dim(testU)[1] != 0) {print(paste("我执行了5", i, sep = "-"))lm1 <- lm(ratIndexIndusStruc ~ amoforeCapUtil + libCollect + foreTradeCoef + realGdpPerCapita,data = trainU)pre1 <- predict(lm1, testU)#print(pre1)#print(pre1)tempdatalm[which(!complete.cases(tempdatalm$ratIndexIndusStruc)), ratIndexIndusStrucNameIndex] <- pre1    }  mydataTemp2[which(mydataTemp2$city %in% cityname2[(numlib[i]+1):numlib[i+1]]), ] <- tempdatalm
}# [1] "我执行了1-1"
# [1] "我执行了4-1"
# [1] "我执行了5-1"
# [1] "我执行了1-2"
# [1] "我执行了2-2"
# [1] "我执行了1-3"
# [1] "我执行了2-3"
# [1] "我执行了3-3"
# [1] "我执行了4-3"
# [1] "我执行了5-3"

绘制缺失值图：

#第三次清理完毕
aggr(mydataTemp2, prop = F, numbers = T)

length(which(complete.cases(mydataTemp2))) #1077
dim(mydataTemp2)missingdf <- as.data.frame(abs(is.na(mydataTemp2)))
varmissing <- apply(missingdf, 2, sum)
varmissing

可以看到进行填补之后, 我们有1077条完整数据, 对于剩下的几个缺失值, 我们单独对其进行处理.现在, 我们直接对剩下的23条数据分别进行均值填补.

首先, 我们查看剩下的缺失值：

mydataTemp2[which(!complete.cases(mydataTemp2)), c(1:2, 5, 11:13, 16:23)]

均值填补:


cityname3 <- unique(mydataTemp2[which(!complete.cases(mydataTemp2)), "city"])
missingdf <- as.data.frame(abs(is.na(mydataTemp2)))
varmissing <- apply(missingdf, 2, sum)
lossingNum <- which(varmissing > 0)
mydataTemp3 <- mydataTemp2
for (itemCity in cityname3) {for (itemValue in lossingNum) {tempSerise <- mydataTemp3[which(mydataTemp3$city == itemCity), itemValue]if (sum(is.na(tempSerise)) > 0) {meanTemp <- mean(tempSerise, na.rm = T)tempSerise[is.na(tempSerise)] <- meanTempmydataTemp3[which(mydataTemp3$city == itemCity), itemValue] <- tempSerise}}
}

绘制缺失值图:

aggr(mydataTemp3, prop = F, numbers = T)

可以看到, 缺失值已经全部被填补.

输出填补结果：

write.csv(mydataTemp3, "../output/mydataTemp3.csv")