我的百分比数据应该适合什么样的曲线（或模型）？

机器算法验证 回归 造型 曲线拟合 百分比
2022-02-14 08:37:17

我正在尝试创建一个显示病毒拷贝和基因组覆盖率 (GCC) 之间关系的图。这是我的数据的样子：

病毒载量与 GCC

起初，我只是绘制了一个线性回归，但我的主管告诉我这是不正确的，并尝试使用 S 形曲线。所以我使用 geom_smooth 做到了这一点：

library(scales)
ggplot(scatter_plot_new, aes(x = Copies_per_uL, y = Genome_cov, colour = Virus)) +
    geom_point() +
    scale_x_continuous(trans = log10_trans(), breaks = trans_breaks("log10", function(x) 10^x), labels = trans_format("log10", math_format(10^.x))) +
        geom_smooth(method = "gam", formula = y ~ s(x), se = FALSE, size = 1) +
    theme_bw() +
    theme(legend.position = 'top', legend.text = element_text(size = 10), legend.title = element_text(size = 12), axis.text = element_text(size = 10), axis.title = element_text(size=12), axis.title.y = element_text(margin = margin (r = 10)), axis.title.x = element_text(margin = margin(t = 10))) +
    labs(x = "Virus copies/µL", y = "GCC (%)") +
    scale_y_continuous(breaks=c(25,50,75,100))

病毒载量与 GCC - geom_smooth

然而，我的主管说这也是不正确的，因为曲线看起来 GCC 可以超过 100%，但它不能。

我的问题是：显示病毒副本和 GCC 之间关系的最佳方式是什么？我想明确一点，A) 低病毒副本 = 低 GCC，并且 B) 在一定数量的病毒复制之后 GCC 平台期。

我研究了很多不同的方法——GAM、LOESS、逻辑、分段——但我不知道如何判断什么是我的数据的最佳方法。

编辑：这是数据：

>print(scatter_plot_new)  
Subsample   Virus   Genome_cov  Copies_per_uL
1   S1.1_RRAV   RRAV    100 92500
2   S1.2_RRAV   RRAV    100 95900
3   S1.3_RRAV   RRAV    100 92900
4   S2.1_RRAV   RRAV    100 4049.54
5   S2.2_RRAV   RRAV    96.9935 3809
6   S2.3_RRAV   RRAV    94.5054 3695.06
7   S3.1_RRAV   RRAV    3.7235  86.37
8   S3.2_RRAV   RRAV    11.8186 84.2
9   S3.3_RRAV   RRAV    11.0929 95.2
10  S4.1_RRAV   RRAV    0   2.12
11  S4.2_RRAV   RRAV    5.0799  2.71
12  S4.3_RRAV   RRAV    0   2.39
13  S5.1_RRAV   RRAV    4.9503  0.16
14  S5.2_RRAV   RRAV    0   0.08
15  S5.3_RRAV   RRAV    4.4147  0.08
16  S1.1_UMAV   UMAV    5.7666  1.38
17  S1.2_UMAV   UMAV    26.0379 1.72
18  S1.3_UMAV   UMAV    7.4128  2.52
19  S2.1_UMAV   UMAV    21.172  31.06
20  S2.2_UMAV   UMAV    16.1663 29.87
21  S2.3_UMAV   UMAV    9.121   32.82
22  S3.1_UMAV   UMAV    92.903  627.24
23  S3.2_UMAV   UMAV    83.0314 615.36
24  S3.3_UMAV   UMAV    90.3458 632.67
25  S4.1_UMAV   UMAV    98.6696 11180
26  S4.2_UMAV   UMAV    98.8405 12720
27  S4.3_UMAV   UMAV    98.7939 8680
28  S5.1_UMAV   UMAV    98.6489 318200
29  S5.2_UMAV   UMAV    99.1303 346100
30  S5.3_UMAV   UMAV    98.8767 345100
4个回答

（编辑考虑到下面的评论。感谢@BenBolker 和@WeiwenNg 的有用输入。）

将分数逻辑回归拟合到数据。它非常适合范围在 0 到 100% 之间的百分比数据，并且在生物学的许多领域理论上是合理的。

请注意，您可能必须将所有值除以 100 才能拟合它，因为程序经常期望数据介于 0 和 1 之间。正如 Ben Bolker 建议的那样，为了解决由二项式分布对方差的严格假设引起的可能问题，请使用取而代之的是拟二项分布。

我根据您的代码做了一些假设，例如您感兴趣的病毒有 2 种，它们可能显示不同的模式（即病毒类型和副本数之间可能存在交互作用）。

一、模型拟合：

dat <- read.csv('Book1.csv')
dat$logcopies <- log10(dat$Copies_per_uL)
dat$Genome_cov_norm <- dat$Genome_cov/100

fit <- glm(Genome_cov_norm ~ logcopies * Virus, data = dat, family = quasibinomial())
summary(fit)


Call:
glm(formula = Genome_cov_norm ~ logcopies * Virus, family = quasibinomial(), 
    data = dat)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.55073  -0.13362   0.07825   0.20362   0.70086  

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)  
(Intercept)          -5.9702     2.8857  -2.069   0.0486 *
logcopies             2.3262     1.0961   2.122   0.0435 *
VirusUMAV             2.6147     3.3049   0.791   0.4360  
logcopies:VirusUMAV  -0.6028     1.3173  -0.458   0.6510  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 0.6934319)

    Null deviance: 30.4473  on 29  degrees of freedom
Residual deviance:  2.7033  on 26  degrees of freedom

如果您相信 p 值，则输出并不表明这两种病毒存在显着差异。这与下面@NickCox 的结果形成对比，尽管我们使用了不同的方法。无论哪种方式，我对 30 个数据点都不是很有信心。

二、出图：

自己编写一种可视化输出的方法并不难，但似乎有一个 ggPredict 包可以为您完成大部分工作（不能保证，我自己没有尝试过）。代码将类似于：

library(ggiraphExtra)
ggPredict(fit) + theme_bw(base_size = 20) + geom_line(size = 2)

更新：我不再更普遍地推荐代码或 ggPredict 函数。在尝试之后，我发现绘制的点并不能完全反映输入数据，而是由于一些奇怪的原因而发生了变化（一些绘制的点高于 1 和低于 0）。所以我建议你自己编写代码，尽管那是更多的工作。

这与@mkt 没有不同的答案，但特别是图表不适合评论。我首先将Stata中的逻辑曲线（在记录预测变量之后）拟合到所有数据并得到这个图

在此处输入图像描述

一个方程是

100 invlogit(-4.192654 + 1.880951 log10( Copies))

现在，在病毒定义指标变量的最简单场景中，我分别为每种病毒拟合曲线。为了记录，这里是一个Stata脚本：

clear 
input id str9 Subsample   str4 Virus   Genome_cov  Copies_per_uL
1   S1.1_RRAV   RRAV    100 92500
2   S1.2_RRAV   RRAV    100 95900
3   S1.3_RRAV   RRAV    100 92900
4   S2.1_RRAV   RRAV    100 4049.54
5   S2.2_RRAV   RRAV    96.9935 3809
6   S2.3_RRAV   RRAV    94.5054 3695.06
7   S3.1_RRAV   RRAV    3.7235  86.37
8   S3.2_RRAV   RRAV    11.8186 84.2
9   S3.3_RRAV   RRAV    11.0929 95.2
10  S4.1_RRAV   RRAV    0   2.12
11  S4.2_RRAV   RRAV    5.0799  2.71
12  S4.3_RRAV   RRAV    0   2.39
13  S5.1_RRAV   RRAV    4.9503  0.16
14  S5.2_RRAV   RRAV    0   0.08
15  S5.3_RRAV   RRAV    4.4147  0.08
16  S1.1_UMAV   UMAV    5.7666  1.38
17  S1.2_UMAV   UMAV    26.0379 1.72
18  S1.3_UMAV   UMAV    7.4128  2.52
19  S2.1_UMAV   UMAV    21.172  31.06
20  S2.2_UMAV   UMAV    16.1663 29.87
21  S2.3_UMAV   UMAV    9.121   32.82
22  S3.1_UMAV   UMAV    92.903  627.24
23  S3.2_UMAV   UMAV    83.0314 615.36
24  S3.3_UMAV   UMAV    90.3458 632.67
25  S4.1_UMAV   UMAV    98.6696 11180
26  S4.2_UMAV   UMAV    98.8405 12720
27  S4.3_UMAV   UMAV    98.7939 8680
28  S5.1_UMAV   UMAV    98.6489 318200
29  S5.2_UMAV   UMAV    99.1303 346100
30  S5.3_UMAV   UMAV    98.8767 345100
end 

gen log10Copies = log10(Copies)
gen Genome_cov_pr = Genome_cov / 100
encode Virus, gen(virus)
set seed 2803 
fracreg logit Genome_cov_pr log10Copies i.virus, vce(bootstrap, reps(10000)) 

twoway function invlogit(-5.055519 + 1.961538 * x), lc(orange) ra(log10Copies)      ///
|| function invlogit(-5.055519 + 1.233273 + 1.961538 * x), ra(log10Copies) lc(blue) ///
|| scatter Genome_cov_pr log10Copies if Virus == "RRAV", mc(orange) ms(Oh)          ///
|| scatter Genome_cov_pr log10Copies if Virus == "UMAV", mc(blue) ms(+)             ///
legend(order(4 "UMAV" 3 "RRAV") pos(11) col(1) ring(0))                             ///
xla(-1 "0.1" 0 "1" 1 "10" 2 "100" 3 "10{sup:3}" 4 "10{sup:4}" 5 "10{sup:5}")        ///
yla(0 .25 "25" .5 "50" .75 "75" 1 "100", ang(h))                                    ///
ytitle(Genome coverage (%)) xtitle(Genome copies / {&mu}L) scheme(s1color)

这对一个很小的数据集进行了大力推动，但病毒的 P 值看起来支持同时拟合两条曲线。

Fractional logistic regression                  Number of obs     =         30
                                                Replications      =     10,000
                                                Wald chi2(2)      =      48.14
                                                Prob > chi2       =     0.0000
Log pseudolikelihood = -6.9603063               Pseudo R2         =     0.6646

-------------------------------------------------------------------------------
              |   Observed   Bootstrap                         Normal-based
Genome_cov_pr |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------+----------------------------------------------------------------
  log10Copies |   1.961538   .2893965     6.78   0.000     1.394331    2.528745
              |
        virus |
        UMAV  |   1.233273   .5557609     2.22   0.026     .1440018    2.322544
        _cons |  -5.055519   .8971009    -5.64   0.000    -6.813805   -3.297234
-------------------------------------------------------------------------------

在此处输入图像描述

另一种解决方法是使用贝叶斯公式，一开始可能有点繁重，但它往往更容易表达你的问题的细节，以及更好地了解“不确定性”在哪里是

Stan是一个 Monte Carlo 采样器，具有相对易于使用的编程接口，库可用于R 和其他，但我在这里使用 Python

我们和其他人一样使用 sigmoid：它具有生化动机，并且在数学上非常便于使用。这个任务的一个很好的参数化是：

import numpy as np

def sigfn(x, alpha, beta):
    return 1 / (1 + np.exp(-(x - alpha) * beta))

其中alpha定义 sigmoid 曲线的中点（即它与 50% 相交的位置）并beta定义斜率，接近零的值更平坦

为了展示它的样子，我们可以提取您的数据并绘制它：

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

df = pd.read_table('raw_data.txt', delim_whitespace=True)
df.columns = ['subsample', 'virus', 'coverage', 'copies']
df.coverage /= 100

x = np.logspace(-1, 6, 201)
plt.semilogx(x, sigfn(np.log(x), 5.5, 3), label='sigfn', color='C2')

sns.scatterplot(df.copies, df.coverage, hue=df.virus, edgecolor='none')

whereraw_data.txt包含您提供的数据，我将覆盖范围转换为更有用的内容。系数 5.5 和 3 看起来不错，并且给出的图与其他答案非常相似：

绘制数据和手动拟合

要使用 Stan“拟合”这个函数，我们需要使用它自己的语言定义我们的模型，该语言是 R 和 C++ 的混合体。一个简单的模型将类似于：

data {
    int<lower=1> N;  // number of rows
    vector[N] log_copies;
    vector<lower=0,upper=1>[N] coverage;
}
parameters {
    real alpha;
    real beta;
    real<lower=0> sigma;
}
model {
    vector[N] mu;
    mu = 1 ./ (1 + exp(-(log_copies - alpha) * beta));

    sigma ~ cauchy(0, 0.1);
    alpha ~ normal(0, 5);
    beta ~ normal(0, 5);

    coverage ~ normal(mu, sigma);
}

希望读起来没问题。当我们对模型进行采样时，我们有一个data块定义了我们期望的数据，parameters定义了被采样的东西，并model定义了似然函数。你告诉斯坦“编译”模型，这需要一段时间，然后你可以从它中提取一些数据。例如：

import pystan

model = pystan.StanModel(model_code=code)
model.sampling(data=dict(
    N=len(df),
    log_copies=np.log(df.copies),
    coverage=df.coverage,
), iter=10000, chains=4, thin=10)

import arviz
arviz.plot_trace(fit)

arviz使漂亮的诊断图变得容易，而打印拟合为您提供了一个漂亮的 R 风格参数摘要：

4 chains, each with iter=10000; warmup=5000; thin=10; 
post-warmup draws per chain=500, total post-warmup draws=2000.

        mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
alpha   5.51  6.0e-3   0.26   4.96   5.36   5.49   5.64   6.12   1849    1.0
beta    2.89    0.04   1.71   1.55   1.98   2.32   2.95   8.08   1698    1.0
sigma   0.08  2.7e-4   0.01   0.06   0.07   0.08   0.09    0.1   1790    1.0
lp__   57.12    0.04   1.76   52.9   56.1  57.58  58.51  59.19   1647    1.0

上的大标准偏差beta表示数据确实没有提供有关此参数的太多信息。还有一些在他们的模型拟合中给出 10+ 有效数字的答案有些夸大了

因为一些答案指出每种病毒可能需要自己的参数，所以我扩展了模型以允许alpha并beta因“病毒”而异。这一切都变得有点繁琐，但是这两种病毒几乎可以肯定具有不同的alpha值（即，对于相同的覆盖范围，您需要更多的 RRAV 副本/μL）并且显示这一点的图是：

数据和 MC 样本图

数据与以前相同，但我为 40 个后验样本绘制了一条曲线。 UMAV似乎相对确定，而RRAV可能遵循相同的斜率并需要更高的拷贝数，或者具有更陡峭的斜率和相似的拷贝数。大多数后部质量都需要更高的拷贝数，但这种不确定性可能解释了其他答案中的一些差异，发现不同的东西

我主要使用回答这个作为练习来提高我对 Stan 的了解，我在这里放了一个 Jupyter 笔记本，以防有人感兴趣/想要复制它。

试试sigmoid函数。这种形状有很多公式，包括逻辑曲线。双曲正切是另一种流行的选择。

鉴于这些图，我也不能排除一个简单的阶跃函数。恐怕您将无法区分阶跃函数和任意数量的 sigmoid 规范。您没有任何观察到您的百分比在 50% 范围内，因此简单的步骤公式可能是最简约的选择，其性能不比更复杂的模型差

其它你可能感兴趣的问题
上一篇为什么柯西分布如此有用？ 下一篇对于线性分类器，更大的系数是否意味着更重要的特征？