在日常科研工作中，缺失数据是一个很常见的问题。特别是在大型的数据集中，由于各种不可抗因素，数据缺失几乎是难以避免的。但这就带来一个问题：当我们面对缺失数据时，应该如何处理？直接删除含有缺失值的数据行似乎是一个简单且直接的方法，但这样会导致有效数据的损失。今天，我想为大家分享几种处理数据缺失的方法。请注意，这些方法各有利弊，最适合的方法应该基于具体的数据特点和研究目的来选择。
构建数据
首先，我们要读入一个30行、14列的生态数据集。这个数据集用于示范如何处理数据中的缺失值。通过随机抽样方法，我们在数据集的copy_SOC列中人为地产生了一些缺失值。
# 数据读入
test_data<- read.csv('H:/data/test_data.csv')
test_data$copy_SOC <- test_data$SOC
# 计算需要替换为NA的数据个数
num_na <- round(nrow(test_data) * 0.20)
# 随机选择 20%索引
random_indices <- sample(1:nrow(test_data), size=num_na)
# 替换选择的索引对应的数据为NA
test_data[random_indices,15] <- NA
colSums(is.na(test_data))
sites NPP ANPP Root.biomass SOC
0 0 0 0 0
TN pH Clay Silt Sand
0 0 0 0 0
Bulk.density total.PLFA Fungal.PLFAs Bacterial.PLFAs copy_SOC
0 0 0 0 6
看到其中copy_SOC列有6个缺失值
简单插补
均值/中位数插补
这是一个非常基础且常用的方法。适用于数据缺失是随机的情况。方法是直接用变量的均值或中位数替代缺失值。
# 使用列的均值、中位数或众数来填充缺失值。这是最简单的方法。
test_data$mean_copy_SOC <- test_data$copy_SOC
test_data$mean_copy_SOC[is.na(test_data$mean_copy_SOC)] <- mean(test_data$copy_SOC, na.rm = TRUE)
test_data$median_copy_SOC <- test_data$copy_SOC
test_data$median_copy_SOC[is.na(test_data$median_copy_SOC)] <- median(test_data$copy_SOC, na.rm = TRUE)
随机插补
直接从已有的观测值中随机选择一个值来替代缺失值。这种方法适用于数据缺失是完全随机的情况。
# 从已有的观测值中随机选择值来填充缺失值。
library(Hmisc)
test_data$Hmisc_copy_SOC <- test_data$copy_SOC
test_data$Hmisc_copy_SOC <- impute(test_data$Hmisc_copy_SOC, 'random')
# 当使用impute函数时，确保你的数据是数值型的，因为这个函数主要针对数值数据设计的。
# impute 也可以使用均值，中值进行插值
# impute(test_data$Hmisc_copy_SOC, 'mean')
# impute(test_data$Hmisc_copy_SOC, 'median')
基于模型的插补方法
线性回归插补
利用其他变量对有缺失值的变量进行线性回归预测，然后用预测值来替代缺失值。
# 使用已知的其他变量作为预测变量，进行线性回归，然后使用该回归模型来预测缺失值。
test_data$lm_copy_SOC <- test_data$copy_SOC
train_data <- test_data[!is.na(test_data$lm_copy_SOC),]
# 使用train_data建立线性模型
lm_fit <- lm(lm_copy_SOC ~ NPP+ANPP+Root.biomass+TN+pH+Clay+Silt+Sand+Bulk.density+total.PLFA+
Fungal.PLFAs+Bacterial.PLFAs,train_data )
# 对线性模型进行逐步回归，筛选变量
lm_fit2 <- step(lm_fit)
#模型总结
summary(lm_fit2)
Call:
lm(formula = lm_copy_SOC ~ NPP + ANPP + Root.biomass + TN + Clay +
Sand + Bulk.density + total.PLFA + Fungal.PLFAs + Bacterial.PLFAs,
data = train_data)
Residuals:
Min 1Q Median 3Q Max
-2.9593 -2.0936 0.2103 1.0633 4.2886
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -50.34984 29.42610 -1.711 0.1308
NPP -0.03340 0.01210 -2.760 0.0281 *
ANPP -0.34054 0.27252 -1.250 0.2516
Root.biomass 0.05054 0.04098 1.233 0.2573
TN 15.00918 1.48659 10.096 2.01e-05 ***
Clay 1.17952 1.16784 1.010 0.3461
Sand 0.65299 0.38771 1.684 0.1360
Bulk.density -9.35362 8.41716 -1.111 0.3032
total.PLFA -1.39401 0.90615 -1.538 0.1678
Fungal.PLFAs 2.88526 1.91431 1.507 0.1755
Bacterial.PLFAs 2.53241 1.72284 1.470 0.1850
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.586 on 7 degrees of freedom
Multiple R-squared: 0.9934, Adjusted R-squared: 0.984
F-statistic: 105.4 on 10 and 7 DF, p-value: 1.149e-06
# 删选除用于预测的数据集
predict_data <- test_data[is.na(test_data$lm_copy_SOC),names(coefficients(lm_fit2))[-1]]
# 使用该模型预测缺失值
predicted_values <- predict(lm_fit2, newdata = predict_data)
# 将预测的值插入到数据中的缺失位置
test_data$lm_copy_SOC[is.na(test_data$lm_copy_SOC)] <- predicted_values
这种方法首先使用其他已知变量建立线性模型，然后用该模型预测缺失值。
k-最近邻插补 (k-NN)
该方法通过查找整个数据集中与缺失值最接近的k个观测值来插补数据。
# 使用DMwR包的knnImputation函数，基于k-NN方法填充缺失值。
remotes::install_github("cran/DMwR")
library(DMwR)
test_data$DMwR_copy_SOC <- test_data$copy_SOC
knnImputation_data <- knnImputation(test_data)
test_data$DMwR_copy_SOC <- knnImputation_data$DMwR_copy_SOC
随机森林插补
随机森林是一种集成学习方法，可以用来处理缺失数据问题。
# 使用missForest包，该方法基于随机森林算法对缺失值进行插补。
library(missForest)
test_data$missForest_copy_SOC <- test_data$copy_SOC
result <- missForest(as.matrix(test_data))
result$OOBerror
test_data_missForest <- as.data.frame(result$ximp)
test_data$missForest_copy_SOC <- test_data_missForest$missForest_copy_SOC
多重插补 (Multiple Imputation)
多重插补是一个更为复杂的方法，但也是目前广泛被认为是处理缺失数据的最佳方法之一。
#有多种实现途径使用mice包进行多重插补。这是一种更复杂但被广泛接受的方法，它创建了多个数据集，并在每个数据集上进行分析。
library(mice)
test_data$mice_copy_SOC <- test_data$copy_SOC
# 进行插补
imputed_test_data <- mice(test_data[c(8:14,22)], m = 5, maxit = 50, method = 'pmm', seed = 10)
# m代表生成的数据集数量, 最大迭代50次， pmm 方法，也可以使用其他方法，具体有
# pmm any Predictive mean matching
# midastouch any Weighted predictive mean matching
# sample any Random sample from observed values
# cart any Classification and regression trees
# rf any Random forest imputations
# mean numeric Unconditional mean imputation
# norm numeric Bayesian linear regression
# norm.nob numeric Linear regression ignoring model error
# norm.boot numeric Linear regression using bootstrap
# norm.predict numeric Linear regression, predicted values
# lasso.norm numeric Lasso linear regression
# lasso.select.norm numeric Lasso select + linear regression
# quadratic numeric Imputation of quadratic terms
# ri numeric Random indicator for nonignorable data
# logreg binary Logistic regression
# logreg.boot binary Logistic regression with bootstrap
# lasso.logreg binary Lasso logistic regression
# lasso.select.logreg binary Lasso select + logistic regression
# polr ordered Proportional odds model
# polyreg unordered Polytomous logistic regression
# lda unordered Linear discriminant analysis
# 2l.norm numeric Level-1 normal heteroscedastic
# 2l.lmer numeric Level-1 normal homoscedastic, lmer
# 2l.pan numeric Level-1 normal homoscedastic, pan
# 2l.bin binary Level-1 logistic, glmer
# 2lonly.mean numeric Level-2 class mean
# 2lonly.norm numeric Level-2 class normal
# 2lonly.pmm any Level-2 class predictive mean matching
# 插补的数据
imputed_test_data$imp$mice_copy_SOC
# 选择第一个数据集
completed_test_data <- mice::complete(imputed_test_data)
test_data$mice_copy_SOC <- completed_test_data$mice_copy_SOC
数据插补效果展示
最后，我们可以使用散点图来直观地查看各种插补方法与原始数据之间的关系。
library(ggplot2)
#设置绘图主题
the <- theme_bw()+
theme(legend.position = "none",
axis.ticks = element_line(color = "black"),
axis.text = element_text(color = "black", size=13),
axis.title= element_text(color = "black", size=13),
axis.line = element_line(color = "black"),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank())
test_data %>%
dplyr::select(SOC, "mean_copy_SOC", "median_copy_SOC", "Hmisc_copy_SOC",
"lm_copy_SOC", "DMwR_copy_SOC", "missForest_copy_SOC",
"mice_copy_SOC") %>%
pivot_longer(cols = -1, ) %>%
ggplot(aes(x=value,y=SOC))+
geom_point() +
geom_smooth(method = 'lm',se=FALSE) +
stat_poly_eq(use_label(c( "R2", "P"), sep = "*\"; \"*"), formula = y ~ x)+
the+
labs(x= 'fited', y= 'real')+
facet_wrap(name~.,ncol=3)+
geom_abline(intercept = 0, slope = 1) # 1:1线
数据缺失是科研中常见的问题，但幸好我们有许多方法可以处理这个问题。本文介绍的方法只是其中的一部分，实际上还有许多其他的方法等待大家去探索和实践。希望这篇文章能对大家有所帮助！如果有任何问题或建议，欢迎留言交流。
4410
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
