Jenis kurva (atau model) apa yang harus saya masukkan ke data persentase saya?

15

Saya mencoba membuat gambar yang menunjukkan hubungan antara salinan virus dan cakupan genom (GCC). Seperti inilah data saya:

Viral load vs GCC

Pada awalnya, saya hanya merencanakan regresi linier tetapi pengawas saya mengatakan itu tidak benar, dan untuk mencoba kurva sigmoidal. Jadi saya melakukan ini menggunakan 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))

Viral load vs GCC - geom_smooth

Namun, penyelia saya mengatakan ini juga tidak benar karena kurva membuatnya terlihat seperti GCC bisa lebih dari 100%, yang tidak bisa.

Pertanyaan saya adalah: apa cara terbaik untuk menunjukkan hubungan antara salinan virus dan GCC? Saya ingin menjelaskan bahwa A) salinan virus rendah = GCC rendah, dan B) setelah sejumlah virus menyalin dataran tinggi GCC.

Saya telah meneliti banyak metode berbeda - GAM, LOESS, logistik, sebagian - tetapi saya tidak tahu bagaimana cara mengetahui metode apa yang terbaik untuk data saya.

EDIT: ini adalah datanya:

>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
penerima teh
sumber
6
Sepertinya regresi logistik akan menjadi yang terbaik, karena ini dibatasi antara 0 dan 100%.
mkt - Pasang kembali Monica
1
Coba (2) model sepotong-bijaksana (linier).
user158565
3
coba tambahkan method.args=list(family=quasibinomial))argumen ke geom_smooth()dalam kode ggplot asli Anda.
Ben Bolker
4
PS Saya akan mendorong Anda untuk tidak menekan kesalahan standar dengan se=FALSE. Selalu menyenangkan untuk menunjukkan kepada orang-orang seberapa besar sebenarnya ketidakpastian itu ...
Ben Bolker
2
Anda tidak memiliki cukup titik data di wilayah transisi untuk mengklaim dengan otoritas apa pun bahwa ada kurva yang mulus. Saya bisa dengan mudah memasukkan fungsi Heaviside ke poin yang Anda tunjukkan pada kami.
Carl Witthoft

Jawaban:

6

Cara lain untuk menggunakan ini adalah dengan menggunakan formulasi Bayesian, mungkin akan sedikit berat untuk memulainya tetapi cenderung membuatnya lebih mudah untuk mengungkapkan secara spesifik masalah Anda serta mendapatkan ide yang lebih baik tentang di mana "ketidakpastian" adalah

Stan adalah sampler Monte Carlo dengan antarmuka terprogram yang relatif mudah digunakan, perpustakaan tersedia untuk R dan lainnya, tetapi saya menggunakan Python di sini

kami menggunakan sigmoid seperti orang lain: ia memiliki motivasi biokimia serta secara matematis sangat nyaman untuk bekerja dengannya. parameterisasi yang bagus untuk tugas ini adalah:

import numpy as np

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

di mana alphamendefinisikan titik tengah dari kurva sigmoid (yaitu di mana ia melintasi 50%) dan betamendefinisikan kemiringan, nilai-nilai yang mendekati nol lebih rata.

untuk menunjukkan seperti apa ini, kami dapat menarik data Anda dan memplotnya dengan:

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')

di mana raw_data.txtberisi data yang Anda berikan dan saya mengubah cakupan menjadi sesuatu yang lebih berguna. koefisien 5.5 dan 3 terlihat bagus dan memberikan plot yang sangat mirip dengan jawaban lainnya:

data plot dan kecocokan manual

untuk "menyesuaikan" fungsi ini menggunakan Stan, kita perlu mendefinisikan model kita menggunakan bahasanya sendiri yang merupakan campuran antara R dan C ++. model sederhana akan menjadi seperti:

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);
}

yang mudah-mudahan berbunyi OK. kami memiliki datablok yang mendefinisikan data yang kami harapkan ketika kami mengambil sampel model, parametersmenentukan hal-hal yang disampel, dan modelmendefinisikan fungsi kemungkinan. Anda memberi tahu Stan untuk "mengkompilasi" model, yang membutuhkan waktu, dan kemudian Anda dapat mengambil sampel dengan beberapa data. sebagai contoh:

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 memudahkan plot diagnostik yang bagus, saat mencetak fit memberi Anda ringkasan parameter R-style yang bagus:

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

standar deviasi yang besar pada betamengatakan bahwa data benar-benar tidak memberikan banyak informasi tentang parameter ini. juga beberapa jawaban yang memberikan 10+ digit signifikan pada model mereka terlalu melebih-lebihkan

karena beberapa jawaban mencatat bahwa setiap virus mungkin memerlukan parameternya sendiri, saya memperluas model untuk memperbolehkan alphadan betabervariasi menurut "Virus". semuanya menjadi sedikit fiddly, tetapi kedua virus hampir pasti memiliki alphanilai yang berbeda (yaitu Anda memerlukan lebih banyak salinan / μL dari RRAV untuk cakupan yang sama) dan plot yang menunjukkan ini adalah:

plot data dan sampel MC

datanya sama seperti sebelumnya, tetapi saya telah menggambar kurva untuk 40 sampel posterior. UMAVtampaknya relatif ditentukan, sementara RRAVdapat mengikuti kemiringan yang sama dan membutuhkan jumlah salinan yang lebih tinggi, atau memiliki kemiringan yang lebih curam dan jumlah salinan yang serupa. sebagian besar massa posterior membutuhkan jumlah salinan yang lebih tinggi, tetapi ketidakpastian ini mungkin menjelaskan beberapa perbedaan dalam jawaban lain yang menemukan hal-hal yang berbeda.

Saya kebanyakan menggunakan menjawab ini sebagai latihan untuk meningkatkan pengetahuan saya tentang Stan, dan saya telah meletakkan buku catatan Jupyter ini di sini kalau-kalau ada yang tertarik / ingin meniru ini.

Sam Mason
sumber
14

(Diedit dengan mempertimbangkan komentar akun di bawah ini. Terima kasih kepada @BenBolker & @WeiwenNg untuk masukan yang bermanfaat.)

Sesuaikan regresi logistik fraksional dengan data. Sangat cocok untuk data persentase yang dibatasi antara 0 dan 100% dan dapat dibenarkan secara teoritis di banyak bidang biologi.

Perhatikan bahwa Anda mungkin harus membagi semua nilai dengan 100 agar sesuai, karena program sering memperkirakan data berkisar antara 0 dan 1. Dan seperti yang disarankan Ben Bolker, untuk mengatasi kemungkinan masalah yang disebabkan oleh asumsi ketat distribusi binomial mengenai varians, gunakan distribusi kuasibinomial sebagai gantinya.

Saya telah membuat beberapa asumsi berdasarkan kode Anda, seperti ada 2 virus yang Anda minati dan mereka mungkin menunjukkan pola yang berbeda (yaitu mungkin ada interaksi antara jenis virus dan jumlah salinan).

Pertama, modelnya pas:

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

Jika Anda mempercayai nilai-p, output tidak menyarankan bahwa kedua virus berbeda bermakna. Ini berbeda dengan hasil @ NickCox di bawah ini, meskipun kami menggunakan metode yang berbeda. Saya juga tidak akan terlalu percaya diri dengan 30 poin data.

Kedua, plot:

Tidak sulit untuk membuat kode cara memvisualisasikan output sendiri, tetapi tampaknya ada paket ggPredict yang akan melakukan sebagian besar pekerjaan untuk Anda (tidak dapat menjamin untuk itu, saya belum mencobanya sendiri). Kode akan terlihat seperti:

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

Pembaruan: Saya tidak lagi merekomendasikan kode atau fungsi ggPredict lebih umum. Setelah mencobanya saya menemukan bahwa titik-titik yang diplot tidak mencerminkan data input tetapi malah diubah karena alasan yang aneh (beberapa titik yang diplot berada di atas 1 dan di bawah 0). Jadi saya sarankan untuk mengkodekannya sendiri, meskipun itu lebih sulit.

mkt - Pasang kembali Monica
sumber
7
Saya mendukung jawaban ini, tetapi saya ingin membuat klarifikasi: Saya akan menyebut regresi logistik fraksional ini. Saya pikir istilah ini akan lebih dikenal luas. Ketika kebanyakan orang mendengar "regresi logistik", saya yakin mereka memikirkan variabel dependen 0/1. Satu jawaban Stackexchange yang bagus yang berhubungan dengan nomenklatur ini ada di sini: stats.stackexchange.com/questions/216122/…
Weiwen Ng
2
@teaelleceecee Anda jelas harus membagi cakupan dengan 100 terlebih dahulu.
Nick Cox
4
gunakan family=quasibinomial()untuk menghindari peringatan (dan masalah mendasar dengan asumsi varians yang terlalu ketat). Ikuti saran @ mkt untuk masalah lainnya.
Ben Bolker
2
Ini mungkin berhasil, tetapi saya ingin memperingatkan orang-orang bahwa Anda harus memiliki premis sebelum mencocokkan suatu fungsi yang data Anda sebenarnya harus mengikuti fungsi itu. Kalau tidak, Anda cukup banyak memotret secara acak ketika Anda memilih fungsi yang pas, dan Anda mungkin tertipu oleh hasilnya.
Carl Witthoft
6
@CarlWitthoft Kami mendengar khotbah tetapi berdosa di luar kebaktian. Apa premis sebelumnya yang mengarahkan Anda untuk menyarankan fungsi Heaviside di komentar lain? Biologi di sini tidak menyerupai transisi pada ambang yang tajam. Fakta penelitian di sini yang saya mengerti adalah bahwa teori formal lebih lemah daripada data. Saya setuju: jika orang berpikir bahwa fungsi langkah masuk akal, mereka harus cocok.
Nick Cox
11

Ini bukan jawaban yang berbeda dari @mkt tetapi grafik pada khususnya tidak akan cocok dengan komentar. Pertama-tama saya memasukkan kurva logistik di Stata (setelah mencatat prediktornya) ke semua data dan mendapatkan grafik ini

masukkan deskripsi gambar di sini

Persamaannya adalah

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

Sekarang saya memasukkan kurva secara terpisah untuk setiap virus dalam skenario paling sederhana dari virus yang mendefinisikan variabel indikator. Di sini untuk catatan adalah skrip 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) 

Ini mendorong keras pada dataset kecil tetapi nilai-P untuk virus terlihat mendukung pemasangan dua kurva secara bersama-sama.

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
-------------------------------------------------------------------------------

masukkan deskripsi gambar di sini

Nick Cox
sumber
3

Coba fungsi sigmoid . Ada banyak formulasi bentuk ini termasuk kurva logistik. Singgung hiperbolik adalah pilihan populer lainnya.

Mengingat plot, saya tidak bisa mengesampingkan fungsi langkah sederhana juga. Saya khawatir Anda tidak akan dapat membedakan antara fungsi langkah dan sejumlah spesifikasi sigmoid. Anda tidak memiliki pengamatan di mana persentase Anda berada dalam kisaran 50%, sehingga perumusan langkah sederhana bisa menjadi pilihan yang paling tidak tepat yang berkinerja tidak lebih buruk daripada model yang lebih kompleks

Aksakal
sumber
σ(x)=12(1+tanhx2)
2
@JG "sigmoid" adalah istilah umum untuk kurva-S, sejauh yang saya ketahui, tetapi Anda benar untuk menunjukkan tautan antara dua spesifikasi sigmoid
Aksakal
2

Berikut adalah 4PL (4 parameter logistik) yang cocok, baik terbatas maupun tidak terbatas, dengan persamaan seperti CA Holstein, M. Griffin, J. Hong, PD Sampson, "Metode Statistik untuk Menentukan dan Membandingkan Batas Deteksi Bioassays", Anal . Chem 87 (2015) 9795-9801. Persamaan 4PL ditunjukkan dalam kedua angka dan makna parameter adalah sebagai berikut: a = asimtot lebih rendah, b = faktor kemiringan, c = titik belok, dan d = asimtot atas.

Gambar 1 membatasi a hingga sama dengan 0% dan d sama dengan 100%:

Gbr. 1 Dibatasi a & d

Gambar 2 tidak memiliki kendala pada 4 parameter dalam persamaan 4PL:

2 Tidak ada kendala

Ini menyenangkan, saya tidak berpura-pura mengetahui sesuatu yang biologis dan akan menarik untuk melihat bagaimana semuanya beres!

Ed V
sumber
Terima kasih, ini sangat membantu. Hanya ingin tahu, apakah Anda melakukan ini di MATLAB dengan fungsi fit?
teaelleceecee
1
Saya menggunakan Igor Pro dengan fungsi pengguna yang ditentukan pengguna yang ditunjukkan pada gambar. Saya telah menggunakan Igor Pro dan pendahulunya (Igor) sejak tahun 1988, tetapi banyak program lain dapat melakukan pemasangan kurva, misalnya, Origin Pro dan Kaleidagraph yang sangat murah. Dan tampaknya Anda memiliki R dan (mungkin?) Akses ke Matlab, tidak ada yang saya ketahui kecuali bahwa mereka sangat mampu. Sukses terbaik dengan ini dan saya harap Anda mendapatkan kabar baik saat berikutnya Anda membahas hal-hal dengan penyelia! Juga, terima kasih telah mengirim data!
Ed V
2

Saya mengekstraksi data dari sebar Anda, dan pencarian persamaan saya menghasilkan persamaan tipe logistik 3-parameter sebagai kandidat yang baik: "y = a / (1.0 + b * exp (-1.0 * c * x))", di mana " x "adalah basis log 10 per plot Anda. Parameter yang dipasang adalah a = 9.0005947126706630E + 01, b = 1.2831794858584102E + 07, dan c = 6.6483431489473155E + 00 untuk data saya yang diekstraksi, kesesuaian (log 10 x) data asli akan menghasilkan hasil yang sama jika Anda menyesuaikan kembali data asli menggunakan nilai saya sebagai perkiraan parameter awal. Nilai parameter saya menghasilkan R-squared = 0,983 dan RMSE = 5,625 pada data yang diekstraksi.

merencanakan

EDIT: Sekarang pertanyaannya telah diedit untuk memasukkan data aktual, berikut adalah plot menggunakan persamaan 3-parameter di atas dan perkiraan parameter awal.

plot2

James Phillips
sumber
Tampaknya ada kesalahan dalam ekstraksi data Anda: Anda memiliki banyak nilai persentase negatif. Juga, nilai maksimum Anda adalah sekitar 90%, bukan 100% seperti pada plot aslinya. Anda mungkin memiliki segalanya diimbangi sekitar 10% untuk beberapa alasan.
mkt - Pasang kembali Monica
Meh - ini adalah data yang diekstraksi secara semi-manual, data asli diperlukan. Ini biasanya cukup untuk pencarian persamaan, dan tentu saja bukan untuk hasil akhir - itulah sebabnya saya mengatakan untuk menggunakan nilai parameter ekstrak-o-fit saya sebagai perkiraan parameter awal pada data asli.
James Phillips
Harap dicatat bahwa karena data aktual sekarang telah ditambahkan ke pos, saya telah memperbarui jawaban ini menggunakan data yang diperbarui.
James Phillips
Untuk mengulangi: penerapan, misalnya, fungsi Heaviside, dapat menghasilkan nilai kesalahan yang sama.
Carl Witthoft
1
@ JamesPhillips Saya akan berusaha melakukannya (Heaviside -> errorbars atau yang setara)
Carl Witthoft
2

Karena saya harus membuka mulut besar tentang Heaviside, inilah hasilnya. Saya mengatur titik transisi ke log10 (viruscopies) = 2.5. Lalu saya menghitung standar deviasi dari dua bagian dari kumpulan data - yaitu, Heaviside mengasumsikan data di kedua sisi memiliki semua turunan = 0.

Sisi RH std dev = 4.76
LH sisi std dev = 7.72

Karena ternyata ada 15 sampel dalam setiap batch, std dev keseluruhan adalah mean, atau 6,24.

Dengan asumsi "RMSE" yang dikutip dalam jawaban lain adalah "RMS error" secara keseluruhan, fungsi Heaviside akan muncul untuk melakukan setidaknya serta, jika tidak lebih baik daripada, sebagian besar "kurva-Z" (dipinjam dari nomenklatur respons fotografis) cocok sini.

sunting

Grafik tidak berguna, tetapi diminta dalam komentar:

Kurva Heaviside cocok

Carl Witthoft
sumber
Apakah Anda harap memposting model dan sebar serupa dengan apa yang dilakukan dalam jawaban lain? Saya paling penasaran melihat hasil ini dan membandingkan. Harap juga tambahkan nilai RMSE dan R-squared untuk perbandingan. Saya pribadi tidak pernah menggunakan fungsi Heaviside dan menemukan ini sangat menarik.
James Phillips
@JamesPhillips Sebenarnya tidak ada yang bisa ditampilkan - jelas scatterplotnya sama; semua yang saya lakukan adalah memilih titik transisi secara manual dan mengambil rata-rata mentah dari setiap set poin (kiri dan kanan). Saya tidak yakinR2memiliki banyak makna di sini.
Carl Witthoft
Makna saya adalah membuat plot yang mirip dengan yang dibuat di jawaban lain, untuk tujuan perbandingan langsung dengan jawaban itu.
James Phillips
2
@JamesPhillips Anda masih memiliki dua keinginan. Pilih dengan bijak :-)
Carl Witthoft
Terima kasih untuk plotnya. Saya mengamati bahwa di semua plot dalam jawaban lain, persamaan yang diplot mengikuti bentuk melengkung dari data di kanan atas - Anda tidak, seperti sifat fungsi Heaviside. Ini secara visual tampaknya bertentangan dengan pernyataan Anda bahwa fungsi Heaviside akan melakukan serta persamaan yang dipasang di jawaban lain - itulah sebabnya saya sebelumnya meminta nilai RMSE dan R-squared, saya menduga bahwa fungsi Heaviside tidak akan mengikuti bentuk dari data di wilayah ini dan mungkin menghasilkan nilai yang lebih buruk untuk statistik yang sesuai.
James Phillips