我的百分比數據應該適合什麼樣的曲線（或模型）？
我正在嘗試創建一個顯示病毒拷貝和基因組覆蓋率（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 可以超過 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
另一種解決方法是使用貝葉斯公式，一開始可能有點繁重，但它往往更容易表達你的問題的細節，以及更好地了解“不確定性”在哪裡是
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')where
raw_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）並且顯示這一點的圖是：
數據與以前相同，但我為 40 個後驗樣本繪製了一條曲線。
UMAV似乎相對確定，而RRAV可能遵循相同的斜率並需要更高的拷貝數，或者俱有更陡峭的斜率和相似的拷貝數。大多數後部質量都需要更高的拷貝數，但這種不確定性可能解釋了其他答案中的一些差異，發現不同的東西我主要使用回答這個作為練習來提高我對 Stan 的了解，我在這裡放了一個 Jupyter 筆記本，以防有人感興趣/想要復制它。
