代码记录

对经济数据集中的缺失值进行处理

• 前言

这个数据集中存在大量的缺失，主要原因是某几个年份的某些指标没找到，或者干脆就是某些指标很难找，导致该指标数据的大批量丢失，更有甚者,由于要查找的年份(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!")] = NA


mydata$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 <- mydata

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 <- lagrange(vecX, vecY, newX)
      print(newX)
      print(newY)
      print("+++++++")
      tempv[which(!complete.cases(tempv)), 2] <- newY
      mydataTemp[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] <- newY
      mydataTemp[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：100


write.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与amoforeCapUtil


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])
  #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和amoforeCapUtil
  trainU <- 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,ratIndexIndusStruc
  trainU <- tempdatalm[which(complete.cases(tempdatalm$regUrbanUnemployRat)), ]
  testU <- tempdatalm[which(!complete.cases(tempdatalm$regUrbanUnemployRat)), ]
  
  if (dim(testU)[1] != 0) {
    print(paste("我执行了3", i, sep = "-"))
    #numHospBed+dischaIndusWaste
    lm1 <- 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和proTechEduExpendGdp
  trainU <- 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和greenCoverRat
  trainU <- 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] <- meantemp

meantemp <- mean(mydataTemp3[which((mydataTemp3$city == "眉山市")), 16], na.rm = T)
mydataTemp3[which(!complete.cases(mydataTemp3$numHospBed) & (mydataTemp3$city == "眉山市")), 16] <- meantemp


mydataTemp3[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 <- 0
    for (i in 1:n) {
      Li <- 1
      for (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 <- pre
  return(newY)
}

• 处理感想

走一步看一步把。