在R中，如何估算右删失缺失数据以遵循假设的分布？

Question

早些时候，已在 Stack Overflow here 中解释了左删失数据的随机插补以遵循假设的分布。使用 censlm 包可以轻松实现。

但是如果我想对右删失数据做类似的事情怎么办？

首先，左删失缺失数据的插补（省略图像注释代码）：

library(ggplot2)
library(dplyr)
set.seed(123)

# Simulate log-normally distributed biomarker data
original_data <- rlnorm(10000, 2.3, 0.4)

# Display minimum value
min(original_data)
#> [1] 2.142283

# Set the lower limit of quantification (= LLOQ)
lloq <- 8

# Erase values below lloq for plotting purposes
left_censored_data <- replace(original_data, original_data < lloq, NA) %>% na.omit()

# Display left-censored data (values below LLOQ erased)
ggplot() +
  geom_histogram(aes(left_censored_data), binwidth = 0.3) 

^{创建于 2023-08-26，使用 reprex v2.0.2}

可以使用 censlm 包轻松随机地估算缺失的左删失值，此处假设对数正态分布（省略图像注释代码）：

library(censlm)

# For imputation, replace values below LLOQ with the LLOQ value
obs <- replace(original_data, original_data < lloq, lloq)

# Compute (random) imputations with {censlm} to fill back the distribution below LLOQ
fit <- clm(log(obs) ~ 1, left = log(lloq))
imputed_data <- exp(imputed(fit))

# Combine the original and imputed data in long format for plotting purposes
overlayed <- data.frame(original_data, imputed_data)
overlayed <- stack(overlayed[, c(1,2)])
names(overlayed) <- c("values", "data_frame")

# Create an overlayed plot 
ggplot(overlayed) +
geom_density(aes(x=values, lty=data_frame, size=data_frame, color=data_frame), alpha=.5, bw = 1.0) +
scale_size_manual("type", values = c(0.5, 3), guide = "none")

^{创建于 2023-08-26，使用 reprex v2.0.2}

但是右删失数据怎么样？如果再次假设对数正态分布，如何计算高于量化上限 (= ULOQ) 的随机值？

# Right-censored data

# Set the UPPER limit of quantification (= ULOQ)
uloq <- 15

# Erase values above ULOQ for plotting purposes
right_censored_data <- replace(original_data, original_data > uloq, NA) %>% na.omit()

# Display right-censored data (values above ULOQ erased)
ggplot() +
  geom_vline(xintercept = uloq, linetype="dotted") +
  geom_histogram(aes(original_data), binwidth = 0.3, alpha = 0.0) +
  geom_histogram(aes(right_censored_data), binwidth = 0.3) +
  geom_text(aes(x = 16, y = 12,
    label = "Values above ULOQ (15.0) right-censored away"),
      hjust = 0, color="blue")

^{创建于 2023-08-26，使用 reprex v2.0.2}

一种非常艰苦的方法是水平翻转数据并再次使用censlm，但肯定存在更优雅的方法吗？

Answer 1

好吧，让我们更清楚地说明我的建议。您需要对数正态分布，因此您必须从不完整的数据中获取 μ、σ。我们需要两个估计值。

所以，刻薄是行不通的，因为我们没有正确的故事。 Stddev 不会起作用，因为我们没有正确的故事。但根据您所掌握的信息，您可以估计两个值：分布模式和半峰全宽 (FWHM)。 您可以从不完整的数据中获得这两个

我粗略地猜测众数为 8，半高宽为 14-5（右侧减去左侧）。 然后，我使用 Python 代码，将均值和 q25 的两个非线性方程替换为众数和 FWHM 方程。我得到了 μ、σ 以及分布形状和采样的一些估计。您只需找到良好的模式和 FWHM 估计器并将它们放入，它应该可以工作

代码，Python 3.10 Windows x64

import numpy as np
import scipy.optimize as opt

import matplotlib.pyplot as plt

LOG_2 = np.log(2.0)

def ff(variables):
    mode= 8. # taken from graph, to be replaced with proper estimator
    fwhm = 14.-5. # taken from graph, to be replaced with proper estimator

    μ, σ = variables

    # taken from wiki https://en.wikipedia.org/wiki/Log-normal_distribution
    mode_eq = np.exp(μ - σ*σ) - mode
    # taked from http://openafox.com/science/peak-function-derivations.html#lognormal
    fwhm_eq = np.exp((μ - σ*σ) + np.sqrt(2.0*σ*σ*LOG_2)) - np.exp((μ - σ*σ) - np.sqrt(2.0*σ*σ*LOG_2)) - fwhm

    return [mode_eq, fwhm_eq]

μ, σ = opt.fsolve(ff, (5,1) )
print(μ, σ)

rng = np.random.default_rng(135797537)
data = rng.lognormal(μ, σ, 200000)

fig, ax = plt.subplots()
ax.hist(data, bins=100)

print('mean: ' + str(np.mean(data)))
print('stdev: ' + str(np.std(data)))

代码产生输出

2.286994459923708 0.45557976057161764
mean: 10.930206083816197
stdev: 5.249043744738874

和图片