Jaką krzywą (lub model) powinienem dopasować do danych procentowych?

15

Próbuję stworzyć postać, która pokazuje związek między kopiami wirusów a pokryciem genomu (GCC). Tak wyglądają moje dane:

Miano wirusa a GCC

Na początku po prostu nakreśliłem regresję liniową, ale moi przełożeni powiedzieli mi, że to nieprawda, i wypróbowałem krzywą sigmoidalną. Zrobiłem to za pomocą 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))

Obciążenie wirusowe vs GCC - geom_smooth

Jednak moi przełożeni twierdzą, że jest to również nieprawidłowe, ponieważ krzywe sprawiają, że wygląda na to, że GCC może przekroczyć 100%, czego nie może.

Moje pytanie brzmi: jaki jest najlepszy sposób na pokazanie związku między kopiami wirusów a GCC? Chcę wyjaśnić, że A) niska liczba kopii wirusa = niski GCC i że B) po określonej ilości wirusa kopiuje płaskowyże GCC.

Badałem wiele różnych metod - GAM, LOESS, logistyka, fragmentarycznie - ale nie wiem, jak powiedzieć, która metoda jest najlepsza dla moich danych.

EDYCJA: to są dane:

>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
teaelleceecee
źródło
6
Wygląda na to, że regresja logistyczna byłaby najlepsza, ponieważ jest ograniczona między 0 a 100%.
mkt - Przywróć Monikę
1
Wypróbuj (2) model częściowy (liniowy).
user158565
3
spróbuj dodać method.args=list(family=quasibinomial))argumenty do geom_smooth()oryginalnego kodu ggplot.
Ben Bolker
4
PS Zachęcam do nie tłumienia standardowych błędów se=FALSE. Zawsze miło jest pokazywać ludziom, jak duża jest niepewność ...
Ben Bolker
2
W regionie przejściowym nie ma wystarczającej liczby punktów danych, aby stwierdzić z jakąkolwiek władzą, że istnieje gładka krzywa. Równie dobrze mógłbym dopasować funkcję Heaviside do punktów, które nam pokazujesz.
Carl Witthoft

Odpowiedzi:

6

Innym sposobem na obejście tego byłoby użycie sformułowania bayesowskiego, może być nieco ciężkie na początku, ale zwykle znacznie ułatwia wyrażanie specyfiki twojego problemu, a także uzyskiwanie lepszych pomysłów na to, gdzie „niepewność” jest

Stan jest samplerem Monte Carlo ze stosunkowo łatwym w użyciu interfejsem programowym, biblioteki są dostępne dla R i innych, ale tutaj używam Pythona

używamy sigmoidu jak wszyscy: ma motywacje biochemiczne, a także jest matematycznie bardzo wygodny w pracy. niezłą parametryzacją tego zadania jest:

import numpy as np

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

gdzie alphaokreśla punkt środkowy krzywej sigmoidalnej (tj. gdzie przecina 50%) i betaokreśla nachylenie, wartości bliższe zeru są bardziej płaskie

aby pokazać, jak to wygląda, możemy pobrać twoje dane i wykreślić je za pomocą:

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

gdzie raw_data.txtzawiera dane, które podałeś, a ja przekształciłem zasięg w coś bardziej użytecznego. współczynniki 5.5 i 3 wyglądają ładnie i dają wykres bardzo podobny do innych odpowiedzi:

dane wykresu i dopasowanie ręczne

aby „dopasować” tę funkcję za pomocą Stana, musimy zdefiniować nasz model za pomocą własnego języka będącego mieszanką R i C ++. prosty model byłby mniej więcej taki:

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

co mam nadzieję, że brzmi OK. mamy datablok, który definiuje dane, których oczekujemy podczas próbkowania modelu, parametersdefiniuje rzeczy, które są próbkowane, i modeldefiniuje funkcję prawdopodobieństwa. Mówisz Stanowi, aby „skompilował” model, co zajmuje chwilę, a następnie możesz próbkować z niego przy użyciu niektórych danych. na przykład:

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 ułatwia tworzenie ładnych wykresów diagnostycznych, a drukowanie pasowania daje ładne podsumowanie parametrów w stylu 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

duże odchylenie standardowe betamówi, że dane naprawdę nie dostarczają zbyt wielu informacji na temat tego parametru. również niektóre z odpowiedzi, które wskazują ponad 10 cyfr znaczących w modelu, są nieco zawyżone

ponieważ w niektórych odpowiedziach zauważono, że każdy wirus może potrzebować własnych parametrów, rozszerzyłem model, aby zezwolić alphai betaróżnić się w zależności od „wirusa”. wszystko robi się trochę dziwnie, ale dwa wirusy prawie na pewno mają różne alphawartości (tj. potrzebujesz więcej kopii / μL RRAV dla tego samego zasięgu), a wykres pokazujący to:

wykres danych i próbek MC

dane są takie same jak poprzednio, ale narysowałem krzywą dla 40 próbek tylnej. UMAVwydaje się stosunkowo dobrze określony, chociaż RRAVmoże podążać tym samym nachyleniem i wymagać większej liczby kopii lub mieć bardziej strome nachylenie i podobną liczbę kopii. większa część tylnej masy wymaga większej liczby kopii, ale ta niepewność może wyjaśnić niektóre różnice w innych odpowiedziach dotyczących znalezienia różnych rzeczy

Najczęściej stosowane odpowiadając to jako ćwiczenie poprawić moją znajomość Stan, a ja umieścić Jupyter zeszyt to tu w przypadku gdy ktoś jest zainteresowany / chce powtórzyć to.

Sam Mason
źródło
14

(Edytowane z uwzględnieniem komentarzy poniżej. Podziękowania dla @BenBolker i @WeiwenNg za pomocne informacje).

Dopasuj ułamkową regresję logistyczną do danych. Jest dobrze dostosowany do danych procentowych, które są ograniczone od 0 do 100% i jest dobrze uzasadniony teoretycznie w wielu obszarach biologii.

Należy pamiętać, że może być konieczne podzielenie wszystkich wartości przez 100, aby je dopasować, ponieważ programy często oczekują, że dane będą zawierać się w przedziale od 0 do 1. I, jak zaleca Ben Bolker, aby rozwiązać ewentualne problemy spowodowane ścisłymi założeniami rozkładu dwumianowego dotyczącymi wariancji, użyj zamiast tego rozkład quasibinomial.

Podjąłem pewne założenia na podstawie twojego kodu, na przykład, że interesują Cię 2 wirusy, które mogą wykazywać różne wzorce (tj. Może istnieć interakcja między typem wirusa a liczbą kopii).

Po pierwsze, model pasuje:

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

Jeśli ufasz wartościom p, dane wyjściowe nie sugerują, że oba wirusy różnią się znacząco. Jest to sprzeczne z poniższymi wynikami @ NickCox, chociaż zastosowaliśmy różne metody. W każdym razie nie byłbym bardzo pewny siebie z 30 punktami danych.

Po drugie, spisek:

Nie jest trudno samodzielnie kodować sposób wizualizacji wyników, ale wydaje się, że istnieje pakiet ggPredict, który wykona większość pracy za ciebie (nie mogę za to ręczyć, sam tego nie próbowałem). Kod będzie wyglądał mniej więcej tak:

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

Aktualizacja: Nie polecam już bardziej ogólnie kodu ani funkcji ggPredict. Po wypróbowaniu okazało się, że wykreślone punkty nie odzwierciedlają dokładnie danych wejściowych, ale zostały zmienione z jakiegoś dziwnego powodu (niektóre z wykreślonych punktów były powyżej 1 i poniżej 0). Zalecam więc samodzielne kodowanie, ale to więcej pracy.

mkt - Przywróć Monikę
źródło
7
Popieram tę odpowiedź, ale chciałbym wyjaśnić: nazwałbym to ułamkową regresją logistyczną. Myślę, że ten termin byłby szerzej rozpoznawany. Kiedy większość ludzi słyszy „regresję logistyczną”, założę się, że myślą o zmiennej zależnej od 0/1. Jedna dobra odpowiedź Stackexchange dotycząca tej nomenklatury znajduje się tutaj: stats.stackexchange.com/questions/216122/...
Ng
2
@teaelleceecee Najwyraźniej najpierw musisz podzielić zasięg przez 100.
Nick Cox
4
użyj, family=quasibinomial()aby uniknąć ostrzeżenia (i podstawowych problemów związanych ze zbyt ścisłymi założeniami dotyczącymi wariancji). Skorzystaj z porady @ mkt na temat drugiego problemu.
Ben Bolker
2
To może zadziałać, ale chciałbym ostrzec ludzi, że powinieneś mieć przesłankę przed dopasowaniem funkcji, że twoje dane powinny być zgodne z tą funkcją. W przeciwnym razie będziesz strzelać losowo, gdy wybierzesz funkcję dopasowania, i możesz dać się zwieść wynikom.
Carl Witthoft
6
@CarlWitthoft Słyszymy kazanie, ale są grzesznikami spoza nabożeństwa. Jakie wcześniejsze przesłanki skłoniły cię do zasugerowania funkcji Heaviside w innych komentarzach? Biologia tutaj nie przypomina przejścia na ostrym progu. Faktem jest tutaj, jak rozumiem, fakt, że teoria formalna jest słabsza niż dane. Zgadzam się: jeśli ludzie myślą, że funkcja kroku ma sens, powinni ją dopasować.
Nick Cox,
11

To nie jest inna odpowiedź niż @mkt, ale w szczególności wykresy nie mieszczą się w komentarzu. Najpierw dopasowuję krzywą logistyczną w Stata (po zalogowaniu predyktora) do wszystkich danych i otrzymuję ten wykres

wprowadź opis zdjęcia tutaj

Równanie jest

100 invlogit(-4,192654 + 1,880951 log10( Copies))

Teraz dopasowuję krzywe osobno dla każdego wirusa w najprostszym scenariuszu definiowania zmiennej wskaźnikowej przez wirusa. Tutaj do nagrania jest skrypt 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) 

To mocno naciska na niewielki zestaw danych, ale wartość P dla wirusa wydaje się wspierać dopasowanie dwóch krzywych razem.

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

wprowadź opis zdjęcia tutaj

Nick Cox
źródło
3

Wypróbuj funkcję sigmoidalną . Istnieje wiele formuł tego kształtu, w tym krzywa logistyczna. Styczna hiperboliczna to kolejny popularny wybór.

Biorąc pod uwagę wykresy, nie mogę również wykluczyć prostej funkcji kroku. Obawiam się, że nie będziesz w stanie odróżnić funkcji kroku od dowolnej liczby specyfikacji sigmoid. Nie masz żadnych obserwacji, w których procent mieści się w zakresie 50%, więc prosty krok może być najbardziej oszczędnym wyborem, który nie ustępuje bardziej niż bardziej złożone modele

Aksakal
źródło
Warto zauważyć, że styczna hiperboliczna jest ściśle związana z funkcją sigmoidalną, a mianowicie. . σ(x)=12(1+tanhx2)
JG
2
@JG „sigmoid” jest, moim zdaniem, ogólnym terminem na krzywą S, ale masz rację, wskazując na związek między dwiema specyfikacjami
sigmoidu
2

Oto dopasowania 4PL (logistyka 4 parametrów), zarówno ograniczone, jak i nieograniczone, z równaniem według CA Holsteina, M. Griffina, J. Honga, PD Sampsona, „Statystyczna metoda określania i porównywania granic wykrywania testów biologicznych”, Anal . Chem. 87 (2015) 9795–9801. Równanie 4PL pokazano na obu rysunkach, a znaczenie parametru jest następujące: a = niższa asymptota, b = współczynnik nachylenia, c = punkt przegięcia, d = górna asymptota.

Rycina 1 ogranicza a do równości 0% id do równości 100%:

Ryc. 1 Ograniczony a d

Rycina 2 nie ma ograniczeń dotyczących 4 parametrów w równaniu 4PL:

Ryc. 2 Brak ograniczeń

To była świetna zabawa, nie udaję, że wiem coś biologicznego i ciekawe będzie, jak to wszystko się ułoży!

Ed V.
źródło
Dziękuję, to jest naprawdę pomocne. Zastanawiam się, czy zrobiłeś to w MATLAB z funkcją dopasowania?
teaelleceecee
1
Użyłem Igor Pro z funkcją użytkownika zdefiniowaną na rysunkach. Używam Igor Pro i jego poprzednika (Igor) od 1988 roku, ale wiele innych programów może dopasowywać krzywą, np. Origin Pro i bardzo niedrogi Kaleidagraph. I wygląda na to, że masz R i (prawdopodobnie?) Dostęp do Matlaba, o którym nic nie wiem poza tym, że są wyjątkowo zdolni. Powodzenia w tym i mam nadzieję, że przy następnej rozmowie z przełożonymi otrzymasz dobre wieści! Dziękujemy również za opublikowanie danych!
Ed V
2

Wyodrębniłem dane z twojego wykresu rozrzutu, a moje wyszukiwanie równań ujawniło 3-parametrowe równanie typu logistycznego jako dobrego kandydata: „y = a / (1.0 + b * exp (-1.0 * c * x))”, gdzie „ x ”to podstawa dziennika 10 na twoją działkę. Dopasowane parametry to = = 9.0005947126706630E + 01, b = 1,2831794858584102E + 07, oraz c = 6,6483431489473155E + 00 dla moich wyodrębnionych danych, dopasowanie oryginalnych danych (log 10 x) powinno dać podobne wyniki, jeśli ponownie dopasujesz oryginalne dane, używając moich wartości jako wstępnych oszacowań parametrów. Moje wartości parametrów dają R-kwadrat = 0,983 i RMSE = 5,625 na wyodrębnionych danych.

wątek

EDYCJA: Teraz, gdy pytanie zostało poddane edycji w celu uwzględnienia rzeczywistych danych, oto wykres wykorzystujący powyższe równanie 3-parametrowe i szacunkowe parametry początkowe.

działka 2

James Phillips
źródło
Wygląda na to, że wystąpił błąd podczas ekstrakcji danych: masz kilka ujemnych wartości procentowych. Ponadto maksymalne wartości wynoszą około 90% zamiast 100%, jak na oryginalnym wykresie. Z jakiegoś powodu możesz mieć wszystko przesunięte o około 10%.
mkt - Przywróć Monikę
Meh - są to częściowo ręcznie wyodrębnione dane, wymagane są oryginalne dane. Zwykle jest to wystarczające do wyszukiwania równań i oczywiście nie do ostatecznych wyników - dlatego powiedziałem, że używam moich wartości parametrów ekstrakcji o-fit jako wstępnych oszacowań parametrów na oryginalnych danych.
James Phillips
Pamiętaj, że ponieważ rzeczywiste dane zostały teraz dodane do postu, zaktualizowałem tę odpowiedź, używając zaktualizowanych danych.
James Phillips
Powtarzam: zastosowanie np. Funkcji Heaviside może dawać podobne wartości błędów.
Carl Witthoft
1
@JamesPhillips Spróbuję to zrobić (Heaviside -> paski błędów lub równoważne)
Carl Witthoft
2

Ponieważ musiałem otworzyć swoje wielkie usta na temat Heaviside, oto wyniki. Ustawiłem punkt przejścia na log10 (viruscopies) = 2,5. Następnie obliczyłem standardowe odchylenia dwóch połówek zbioru danych - to znaczy, że Heaviside zakłada, że ​​dane po obu stronach mają wszystkie pochodne = 0.

Od lewej strony std dev = 4,76 Po
lewej stronie std dev = 7,72

Ponieważ okazuje się, że w każdej partii jest 15 próbek, ogólny std dev jest średnią, czyli 6,24.

Zakładając, że „RMSE” ​​cytowany w innych odpowiedziach to ogólnie „błąd RMS”, wydaje się, że funkcja Heaviside działa co najmniej tak dobrze, a jeśli nie lepiej niż, większość „krzywej Z” (zapożyczonej z nomenklatury odpowiedzi fotograficznej) pasuje tutaj.

edytować

Bezużyteczny wykres, ale wymagany w komentarzach:

Dopasowana krzywa Heaviside

Carl Witthoft
źródło
Czy możesz zamieścić model i wykres punktowy podobnie jak w innych odpowiedziach? Najbardziej ciekawi mnie te wyniki i porównanie. Dodaj również RMSE i wartości R-kwadrat do porównania. Osobiście nigdy nie korzystałem z funkcji Heaviside i uważam to za bardzo interesujące.
James Phillips
R2
Miałem na myśli stworzenie fabuły podobnej do tej z innych odpowiedzi, w celu bezpośredniego porównania z tymi odpowiedziami.
James Phillips
2
@JamesPhillips, zostały ci dwa życzenia. Wybierz mądrze :-)
Carl Witthoft
Dziękuję za fabułę. Zauważyłem, że na wszystkich wykresach w innych odpowiedziach wykreślone równanie podąża za zakrzywionym kształtem danych w prawym górnym rogu - twoje nie, podobnie jak natura funkcji Heaviside. Wizualnie wydaje się to zaprzeczać twojemu twierdzeniu, że funkcja Heaviside zrobiłaby, a także równaniom zamieszczonym w innych odpowiedziach - dlatego wcześniej poprosiłem o wartości RMSE i R-kwadrat, podejrzewałem, że funkcja Heaviside nie będzie podążać za kształtem danych w tym regionie i może przynieść gorsze wartości dla tych pasujących statystyk.
James Phillips