Welche Art von Kurve (oder Modell) sollte ich an meine prozentualen Daten anpassen?

15

Ich versuche, eine Abbildung zu erstellen, die die Beziehung zwischen Viruskopien und Genomabdeckung (GCC) zeigt. So sehen meine Daten aus:

Viruslast vs GCC

Zuerst habe ich nur eine lineare Regression gezeichnet, aber meine Vorgesetzten sagten mir, dass dies nicht korrekt sei, und versuchten es mit einer Sigmoidalkurve. Also habe ich das mit geom_smooth gemacht:

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

Viruslast vs GCC - geom_smooth

Meine Vorgesetzten sagen jedoch, dass dies auch falsch ist, da die Kurven den Eindruck erwecken, dass GCC über 100% hinausgehen kann, was nicht der Fall ist.

Meine Frage lautet: Wie lässt sich die Beziehung zwischen Viruskopien und GCC am besten veranschaulichen? Ich möchte klarstellen, dass A) niedrige Viruskopien = niedrige GCC, und dass B) nach einer bestimmten Menge von Viren die GCC-Plateaus kopiert.

Ich habe viele verschiedene Methoden recherchiert - GAM, LÖSS, logistisch, stückweise - aber ich weiß nicht, wie ich sagen soll, welche Methode für meine Daten die beste ist.

EDIT: das sind die Daten:

>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
Tee-Tee
quelle
6
Es scheint, als wäre eine logistische Regression am besten, da diese zwischen 0 und 100% liegt.
mkt - Setzen Sie Monica
1
Versuchen Sie (2) stückweise (lineare) Modelle.
User158565
3
Versuchen Sie, method.args=list(family=quasibinomial))die Argumente geom_smooth()in Ihren ursprünglichen ggplot-Code einzufügen.
Ben Bolker
4
PS Ich würde Sie ermutigen , Standardfehler nicht mit zu unterdrücken se=FALSE. Immer schön zu zeigen, wie groß die Unsicherheit tatsächlich ist ...
Ben Bolker
2
Sie haben nicht genügend Datenpunkte im Übergangsbereich, um mit einer Autorität zu behaupten, dass es eine glatte Kurve gibt. Ich könnte genauso gut eine Heaviside-Funktion an die Punkte anpassen, die Sie uns zeigen.
Carl Witthoft

Antworten:

6

Ein anderer Weg, dies zu erreichen, wäre die Verwendung einer Bayes'schen Formulierung. Es kann anfangs etwas schwer sein, aber es macht es in der Regel viel einfacher, Einzelheiten Ihres Problems auszudrücken und bessere Vorstellungen darüber zu bekommen, wo die "Unsicherheit" liegt. ist

Stan ist ein Monte-Carlo-Sampler mit einer relativ einfach zu bedienenden Programmierschnittstelle. Bibliotheken sind für R und andere verfügbar, aber ich verwende hier Python

Wir verwenden wie alle anderen auch ein Sigmoid: Es hat biochemische Motivationen und ist mathematisch sehr einfach zu handhaben. Eine schöne Parametrisierung für diese Aufgabe ist:

import numpy as np

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

Wo alphader Mittelpunkt der Sigmoidkurve definiert (dh wo er 50% kreuzt) und betadie Steigung definiert, sind Werte nahe Null flacher

Um zu zeigen, wie dies aussieht, können wir Ihre Daten abrufen und zeichnen mit:

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

Wo raw_data.txtenthält die Daten, die Sie angegeben haben, und ich habe die Berichterstattung in etwas Nützlicheres umgewandelt. Die Koeffizienten 5.5 und 3 sehen gut aus und geben eine Darstellung, die den anderen Antworten sehr ähnlich ist:

Plotdaten und manuelle Anpassung

Um diese Funktion mit Stan "anzupassen", müssen wir unser Modell mit einer eigenen Sprache definieren, die eine Mischung aus R und C ++ ist. Ein einfaches Modell wäre so etwas wie:

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

was hoffentlich OK lautet. Wir haben einen dataBlock, der die Daten definiert, die wir erwarten, wenn wir das Modell parametersabtasten, die Dinge, die abgetastet werden, und modeldie Wahrscheinlichkeitsfunktion definiert. Sie weisen Stan an, das Modell zu "kompilieren", was eine Weile dauert, und dann können Sie mit einigen Daten davon abtasten. beispielsweise:

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 Erleichtert das Ausdrucken von Diagnosediagnosen, und beim Drucken der Anpassung erhalten Sie eine schöne Zusammenfassung der Parameter im R-Stil:

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

Die große Standardabweichung von gibt an beta, dass die Daten nicht wirklich viele Informationen zu diesem Parameter liefern. Auch einige der Antworten, die 10+ signifikante Ziffern in ihren Modellanpassungen angeben, übertreiben die Dinge etwas

Da in einigen Antworten darauf hingewiesen wurde, dass jeder Virus möglicherweise seine eigenen Parameter benötigt, habe ich das Modell erweitert, um zuzulassen alphaund betaje nach "Virus" zu variieren. es wird alles ein bisschen fummelig, aber die beiden Viren haben mit ziemlicher Sicherheit unterschiedliche alphaWerte (dh Sie benötigen mehr Kopien / μl RRAV für die gleiche Abdeckung) und ein Diagramm, das Folgendes zeigt:

Plot von Daten und MC-Samples

Die Daten sind die gleichen wie zuvor, aber ich habe eine Kurve für 40 Proben des Seitenzahns gezeichnet. UMAVscheint relativ gut bestimmt zu sein, RRAVkönnte jedoch der gleichen Steigung folgen und eine höhere Kopienanzahl erfordern oder eine steilere Steigung und eine ähnliche Kopienanzahl aufweisen. Der Großteil der posterioren Masse benötigt eine höhere Kopienzahl, aber diese Unsicherheit könnte einige der Unterschiede bei anderen Antworten erklären, die andere Dinge finden

Ich meist Beantwortung dieses als Übung verwendete mein Wissen von Stan zu verbessern, und ich habe einen Jupyter Notebook diesen setzte hier , falls jemand interessiert ist / will dies replizieren.

Sam Mason
quelle
14

(Unter Berücksichtigung der folgenden Kommentare bearbeitet. Vielen Dank an @BenBolker & @WeiwenNg für hilfreiche Eingaben.)

Passen Sie eine partielle logistische Regression an die Daten an. Es eignet sich gut für prozentuale Daten, die zwischen 0 und 100% liegen und in vielen Bereichen der Biologie theoretisch gut begründet sind.

Beachten Sie, dass Sie möglicherweise alle Werte durch 100 teilen müssen, um sie anzupassen, da Programme häufig einen Datenbereich zwischen 0 und 1 erwarten. Um mögliche Probleme zu beheben, die durch die strengen Annahmen der Binomialverteilung in Bezug auf die Varianz verursacht werden, verwenden Sie a, wie von Ben Bolker empfohlen stattdessen Quasibinomialverteilung.

Ich habe einige Annahmen getroffen, die auf Ihrem Code basieren, z. B. dass es zwei Viren gibt, die Sie interessieren und die möglicherweise unterschiedliche Muster aufweisen (dh, es kann eine Wechselwirkung zwischen dem Virentyp und der Anzahl der Kopien geben).

Zunächst passte das Modell:

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

Wenn Sie den p-Werten vertrauen, deutet die Ausgabe nicht darauf hin, dass sich die beiden Viren signifikant unterscheiden. Dies steht im Gegensatz zu den Ergebnissen von @ NickCox, obwohl wir unterschiedliche Methoden angewendet haben. Mit 30 Datenpunkten wäre ich sowieso nicht sehr zuversichtlich.

Zweitens ist die Handlung:

Es ist nicht schwer, einen Weg zu finden, um die Ausgabe selbst zu visualisieren, aber es scheint ein ggPredict-Paket zu geben, das die meiste Arbeit für Sie erledigt (kann nicht dafür bürgen, ich habe es selbst nicht ausprobiert). Der Code sieht ungefähr so ​​aus:

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

Update: Ich empfehle den Code oder die ggPredict-Funktion nicht mehr allgemeiner. Nach dem Ausprobieren stellte ich fest, dass die eingezeichneten Punkte nicht genau die Eingabedaten widerspiegeln, sondern aus bizarren Gründen geändert wurden (einige der eingezeichneten Punkte lagen über 1 und unter 0). Daher empfehle ich, es selbst zu codieren, obwohl das mehr Arbeit ist.

mkt - Setzen Sie Monica wieder ein
quelle
7
Ich stimme dieser Antwort zu, möchte jedoch klarstellen, dass ich diese fraktionierte logistische Regression nennen würde. Ich denke, dieser Begriff würde allgemein anerkannt werden. Wenn die meisten Leute "logistische Regression" hören, denke ich an eine 0/1-abhängige Variable. Eine gute Antwort von Stackexchange zu dieser Nomenklatur ist hier: stats.stackexchange.com/questions/216122/…
Weiwen Ng
2
@teaelleceecee Sie müssen die Abdeckung offensichtlich zuerst durch 100 teilen.
Nick Cox
4
Verwenden Sie family=quasibinomial()diese Option, um die Warnung (und die zugrunde liegenden Probleme mit zu strengen Varianzannahmen) zu vermeiden. Nehmen Sie den Rat von @ mkt bezüglich des anderen Problems an.
Ben Bolker
2
Dies mag funktionieren, aber ich möchte die Leute warnen, dass Sie eine Voraussetzung haben sollten, bevor Sie eine Funktion anpassen, dass Ihre Daten tatsächlich dieser Funktion folgen sollten . Andernfalls schießen Sie ziemlich zufällig, wenn Sie eine passende Funktion auswählen, und Sie können sich von den Ergebnissen täuschen lassen.
Carl Witthoft
6
@CarlWitthoft Wir hören die Predigt, sind aber Sünder außerhalb des Gottesdienstes. Welche vorherige Prämisse hat Sie veranlasst, in anderen Kommentaren eine Heaviside-Funktion vorzuschlagen? Die Biologie ähnelt hier nicht dem Übergang an einer scharfen Schwelle. Nach meinem Verständnis ist die formale Theorie schwächer als die Daten. Ich stimme zu: Wenn die Leute denken, dass eine Sprungfunktion Sinn macht, sollten sie zu einer passen.
Nick Cox
11

Dies ist keine andere Antwort als @mkt, aber insbesondere Grafiken passen nicht in einen Kommentar. Ich habe zuerst eine logistische Kurve in Stata (nach dem Protokollieren des Prädiktors) an alle Daten angepasst und diese Grafik erhalten

Bildbeschreibung hier eingeben

Eine Gleichung lautet

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

Im einfachsten Szenario, in dem ein Virus eine Indikatorvariable definiert, passe ich die Kurven für jeden Virus separat an. Hier für den Datensatz ist ein Stata-Skript:

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) 

Dies drängt auf einen winzigen Datensatz, aber der P-Wert für den Virus scheint die Anpassung von zwei Kurven gemeinsam zu unterstützen.

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

Bildbeschreibung hier eingeben

Nick Cox
quelle
3

Probieren Sie die Sigmoid- Funktion aus. Es gibt viele Formulierungen dieser Form, einschließlich einer logistischen Kurve. Hyperbolischer Tangens ist eine weitere beliebte Wahl.

Angesichts der Handlungen kann ich auch eine einfache Sprungfunktion nicht ausschließen. Ich befürchte, Sie werden nicht in der Lage sein, zwischen einer Sprungfunktion und einer beliebigen Anzahl von Sigmoidspezifikationen zu unterscheiden. Sie haben keine Beobachtungen, bei denen Ihr Prozentsatz im Bereich von 50% liegt. Daher kann die einfache Schrittformulierung die sparsamste Wahl sein, die nicht schlechter abschneidet als komplexere Modelle

Aksakal
quelle
Es ist erwähnenswert, dass der hyperbolische Tangens eng mit der Sigmoidfunktion zusammenhängt, nämlich. . σ(x)=12(1+tanhx2)
JG
2
@JG "Sigmoid" ist für mich ein Oberbegriff für eine S-Kurve, aber Sie weisen zu Recht auf eine Verbindung zwischen zwei Spezifikationen eines Sigmoid hin
Aksakal
2

Hier sind die 4PL (4 parameter logistic) -Fits, sowohl eingeschränkt als auch nicht eingeschränkt, mit der Gleichung gemäß CA Holstein, M. Griffin, J. Hong, PD Sampson, "Statistische Methode zur Bestimmung und zum Vergleich der Nachweisgrenzen von Bioassays", Anal . Chem. 87 (2015) 9795 & ndash; 9801. Die 4PL-Gleichung ist in beiden Figuren gezeigt und die Parameterbedeutungen sind wie folgt: a = untere Asymptote, b = Steigungsfaktor, c = Wendepunkt und d = obere Asymptote.

Abbildung 1 beschränkt a auf 0% und d auf 100%:

Abb. 1 Eingeschränktes a & d

Abbildung 2 unterliegt keinen Einschränkungen für die 4 Parameter in der 4PL-Gleichung:

Abb. 2 Keine Einschränkungen

Das hat Spaß gemacht, ich mache keinen Vorwand, irgendetwas Biologisches zu wissen, und es wird interessant sein zu sehen, wie sich alles einpendelt!

Ed V
quelle
Danke, das ist wirklich hilfreich. Haben Sie sich nur gefragt, ob Sie dies in MATLAB mit der Fit-Funktion gemacht haben?
Teaelleceecee
1
Ich habe Igor Pro mit der in den Abbildungen gezeigten benutzerdefinierten Benutzerfunktion verwendet. Ich habe Igor Pro und seinen Vorgänger (Igor) seit 1988 verwendet, aber viele andere Programme können die Kurvenanpassung durchführen, z. B. Origin Pro und der sehr kostengünstige Kaleidagraph. Und es scheint, dass Sie R und (möglicherweise?) Zugang zu Matlab haben, von denen ich nichts weiß, außer dass sie äußerst fähig sind. Ich wünsche Ihnen viel Erfolg und hoffe, dass Sie beim nächsten Gespräch mit den Vorgesetzten gute Nachrichten erhalten! Vielen Dank auch für die Veröffentlichung der Daten!
Ed V
2

Ich habe die Daten aus Ihrem Streudiagramm extrahiert und meine Gleichungssuche ergab eine logistische Gleichung mit drei Parametern als guten Kandidaten: "y = a / (1,0 + b * exp (-1,0 * c * x))", wobei " x "ist die logarithmische Basis 10 für Ihr Diagramm. Die angepassten Parameter waren a = 9.0005947126706630E + 01, b = 1.2831794858584102E + 07 und c = 6.6483431489473155E + 00. Eine Anpassung der (log 10 x) Originaldaten sollte bei erneuter Anpassung zu ähnlichen Ergebnissen führen die ursprünglichen Daten unter Verwendung meiner Werte als anfängliche Parameterschätzungen. Meine Parameterwerte ergeben R-Quadrat = 0,983 und RMSE = 5,625 für die extrahierten Daten.

Handlung

BEARBEITEN: Nachdem die Frage so bearbeitet wurde, dass sie die tatsächlichen Daten enthält, ist hier eine grafische Darstellung unter Verwendung der obigen 3-Parameter-Gleichung und der anfänglichen Parameterschätzungen.

plot2

James Phillips
quelle
Es scheint ein Fehler in Ihrer Datenextraktion aufgetreten zu sein: Sie haben eine Reihe negativer Prozentwerte. Außerdem liegen Ihre Maximalwerte bei etwa 90% anstelle von 100% wie im ursprünglichen Diagramm. Möglicherweise haben Sie alles aus irgendeinem Grund um etwa 10% versetzt.
mkt - Setzen Sie Monica
Meh - das sind halb manuell extrahierte Daten, die Originaldaten werden benötigt. Dies ist normalerweise ausreichend für die Suche nach Gleichungen und natürlich nicht für die endgültigen Ergebnisse. Aus diesem Grund habe ich gesagt, dass ich meine Werte für den Parameter "Extract-O-Fit" als anfängliche Parameterschätzung für die Originaldaten verwenden soll.
James Phillips
Bitte beachten Sie, dass ich diese Antwort unter Verwendung der aktualisierten Daten aktualisiert habe, da die tatsächlichen Daten jetzt zum Beitrag hinzugefügt wurden.
James Phillips
Nur um es noch einmal zu wiederholen: Die Anwendung von z. B. einer Heaviside-Funktion kann ähnliche Fehlerwerte ergeben.
Carl Witthoft
1
@ JamesPhillips Ich werde versuchen, dies zu tun (Heaviside -> Fehlerbalken oder gleichwertig)
Carl Witthoft
2

Da ich meine große Klappe über Heaviside aufmachen musste, sind hier die Ergebnisse. Ich habe den Übergangspunkt auf log10 (viruscopies) = 2.5 gesetzt. Dann berechnete ich die Standardabweichungen der beiden Hälften des Datensatzes - das heißt, der Heaviside geht davon aus, dass die Daten auf beiden Seiten alle Ableitungen = 0 haben.

Rechte Seite std dev = 4,76
linke Seite std dev = 7,72

Da sich herausstellt, dass in jeder Charge 15 Proben enthalten sind, ist der Standard-Gesamtwert der Mittelwert oder 6,24.

Unter der Annahme, dass das in anderen Antworten angegebene "RMSE" insgesamt "RMS-Fehler" ist, scheint die Heaviside-Funktion mindestens so gut zu funktionieren wie die meisten "Z-Kurven" -Passungen (aus der fotografischen Antwortnomenklatur entlehnt), wenn nicht sogar besser als diese Hier.

bearbeiten

Nutzloses Diagramm, aber in Kommentaren angefordert:

Heaviside Kurvenanpassung

Carl Witthoft
quelle
Woukd Sie bitte ein Modell und Streudiagramm ähnlich wie in den anderen Antworten? Ich bin sehr gespannt auf diese Ergebnisse und Vergleiche. Bitte addieren Sie auch RMSE- und R-Quadrat-Werte zum Vergleich. Ich persönlich habe die Heaviside-Funktion noch nie benutzt und finde das sehr interessant.
James Phillips
@JamesPhillips Es gibt wirklich nichts zu zeigen - offensichtlich ist das Streudiagramm dasselbe; Alles, was ich getan habe, war, den Übergangspunkt manuell auszuwählen und den rohen Mittelwert für jeden Satz von Punkten (links und rechts) zu berechnen. Ich bin mir nicht sicher, ob eine große Bedeutung hat. R2
Carl Witthoft
Meine Absicht war es, ein Diagramm zu erstellen, das den in den anderen Antworten gemachten ähnlich ist, um einen direkten Vergleich mit diesen Antworten zu ermöglichen.
James Phillips
2
@ JamesPhillips Sie haben zwei Wünsche übrig. Wählen Sie mit Bedacht :-)
Carl Witthoft
Vielen Dank für die Handlung. Ich beobachte, dass in allen Plots in anderen Antworten die geplottete Gleichung der gekrümmten Form der Daten oben rechts folgt - dies ist nicht der Fall, so wie es die Heaviside-Funktion ist. Dies scheint Ihrer Behauptung visuell zu widersprechen, dass die Heaviside-Funktion genauso gut funktioniert wie die in den anderen Antworten angegebenen Gleichungen. Aus diesem Grund hatte ich zuvor die RMSE- und R-Quadrat-Werte angefordert und vermutet, dass die Heaviside-Funktion der Form nicht folgen würde der Daten in dieser Region und könnte schlechtere Werte für die Fit-Statistiken ergeben.
James Phillips