如何将负对数10分布应用于Python probplot函数以绘制p值的QQ图？

Question

我想把下面的代码从R语言转换成Python，使用的是scipy.stats.probplot这个库。

qqplot(-log10(ppoints(1000)), -log10(p_value))

这是一个Q-Q图，用来比较p值和均匀分布，采用的是负对数刻度。我想要的效果大概是这样的。（我知道还有其他库可以实现这个，但我现在只想用probplot来解决。）

probplot(-np.log10(p_values_data), dist="uniform", sparams=(0, 1), plot=plt)

这个方法不太对，因为x轴是均匀分布的。在这里，plt是通过import matplotlib.pyplot as plt引入的。我在网上找到了一些帖子，比如这个，但没有找到关于如何修改dist参数来适应-log10(uniform)的内容。

我该如何使用probplot来得到这个图呢？

这是对问题描述的进一步说明。

这里是数据生成的部分。

import numpy as np
from scipy.stats import chi2,probplot
from statsmodels.formula.api import ols
import matplotlib.pyplot as plt

def compute_p_with_chi2(x,y):
    model = ols('y ~ x', data=dict(y=y, x=x)).fit()
    t_stat = model.tvalues['x']
    p_value = 1-chi2.cdf(t_stat**2, 1)
    return p_value

def compute_pvalues(X_data,p_data):
  p_values = []
  for col in X_data.T:
      p_value = compute_p_with_chi2(col,p_data)
      p_values.append(p_value)
  return p_values

n = 100
p = 1000
X = np.random.binomial(2, 0.4, size=(n, p))
y = np.random.normal(size=n)

p_values = compute_pvalues(X,y)

当我对p值做直方图时，得到了预期的均匀分布。

plt.hist(p_values)

但是，使用probplot绘制Q-Q图时，我没有得到两条重叠的对角线。以下是我得到的结果。

probplot(-np.log10(p_values), dist="uniform", sparams=(0, 1), plot=plt)

我在这里附上了R语言中想要的输出，和上面的（第一）段代码。

我觉得这应该是个很简单的事情，但我似乎还是没搞明白。

Answer 1

听起来第一个参数需要是：

-np.log10(p_values_data)

这里的 p_values_data 是一些在 0 和 1 之间分布的值，而我们要将这些值与均匀分布进行比较。这和我之前的回答是一样的，只是换了一种可视化的方式。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import probplot

a, b = 1e-6, 1
p_values_data = np.logspace(-6, 0, 1000)
probplot(-np.log10(p_values_data), dist="uniform", sparams=(0, 6), plot=plt)
# plt.gca().set_xscale('log')  # if desired for some reason

与用户原始代码的主要区别在于，传递给均匀分布的参数应该反映出分布的范围——下限可以是 0，但在这种情况下，上限需要至少是 -log10(np.min(p_values_data))。如果你选择改变下限，请注意 SciPy 的均匀分布 是通过左端点和范围（端点之间的差值）来参数化的，而不是单独的左右端点。

如果其中一个 p_values_data 的值恰好是 0，那么 -log10(0) 是无穷大。在这种情况下，你需要明确说明你希望发生什么。

再次强调，如果这没有回答你的问题，请在提问时提供一个最小的可重现示例——在这种情况下，就是数据和你想要生成的图的示例。

Answer 2

这个dist参数接受一个对象，这个对象就像一个概率分布；更具体一点，它必须有一个ppf方法。如果你想把你的数据和对数均匀分布（也就是一个随机变量的对数是均匀分布的）进行比较，你可以这样做：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import probplot
from scipy import stats

a, b = 1e-6, 1
p_values_data = np.logspace(-6, 0, 100)
probplot(p_values_data, dist=stats.loguniform, sparams=(a, b), plot=plt)

如果你想让你的x轴是对数间隔的，你可以添加：

plt.gca().set_xscale('log')

如果这些内容没有回答你的问题，请在提问时提供一个最小可复现的例子——在这种情况下，就是你的数据和你想要生成的图的示例。

Answer 3

你可以手动处理这些数据，这样就能把它们和一个统一的分布进行比较：

transformed_data = -np.log10(p_values_data)

expected_quantiles = stats.uniform.ppf(np.linspace(0.001, 0.999, len(transformed_data)))

然后就是你之前提供的那个命令了