Quel type de courbe (ou modèle) dois-je adapter à mes données de pourcentage?

15

J'essaie de créer une figure qui montre la relation entre les copies virales et la couverture du génome (GCC). Voici à quoi ressemblent mes données:

Charge virale vs GCC

Au début, je viens de tracer une régression linéaire mais mes superviseurs m'ont dit que c'était incorrect et d'essayer une courbe sigmoïdale. J'ai donc fait cela en utilisant 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))

Charge virale vs GCC - geom_smooth

Cependant, mes superviseurs disent que c'est incorrect aussi parce que les courbes donnent l'impression que GCC peut dépasser 100%, ce qu'il ne peut pas.

Ma question est: quelle est la meilleure façon de montrer la relation entre les copies de virus et GCC? Je tiens à préciser que A) faibles copies de virus = faible GCC, et que B) après une certaine quantité de virus copie les plateaux GCC.

J'ai recherché beaucoup de méthodes différentes - GAM, LOESS, logistique, par morceaux - mais je ne sais pas comment déterminer la meilleure méthode pour mes données.

EDIT: ce sont les données:

>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
la source
6
Il semble qu'une régression logistique serait préférable, car elle est limitée entre 0 et 100%.
mkt
1
Essayez (2) un modèle par morceaux (linéaire).
user158565
3
essayez d'ajouter method.args=list(family=quasibinomial))les arguments à geom_smooth()votre code ggplot d'origine.
Ben Bolker
4
PS Je vous encourage à ne pas supprimer les erreurs standard avec se=FALSE. Toujours agréable de montrer aux gens à quel point l'incertitude est grande ...
Ben Bolker
2
Vous n'avez pas suffisamment de points de données dans la région de transition pour affirmer avec autorité qu'il existe une courbe lisse. Je pourrais tout aussi bien adapter une fonction Heaviside aux points que vous nous montrez.
Carl Witthoft

Réponses:

6

une autre façon de procéder serait d'utiliser une formulation bayésienne, elle peut être un peu lourde au départ, mais elle a tendance à faciliter l'expression des spécificités de votre problème et à mieux comprendre où se situe «l'incertitude». est

Stan est un échantillonneur Monte Carlo avec une interface de programmation relativement facile à utiliser, des bibliothèques sont disponibles pour R et d'autres mais j'utilise Python ici

nous utilisons un sigmoïde comme tout le monde: il a des motivations biochimiques et il est mathématiquement très pratique de travailler avec. un bon paramétrage pour cette tâche est:

import numpy as np

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

alphadéfinit le milieu de la courbe sigmoïde (c'est-à-dire là où elle croise 50%) etbeta définit la pente, les valeurs plus proches de zéro sont plus plates

pour montrer à quoi cela ressemble, nous pouvons extraire vos données et les tracer avec:

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

raw_data.txtcontient les données que vous avez fournies et j'ai transformé la couverture en quelque chose de plus utile. les coefficients 5.5 et 3 sont agréables et donnent un tracé très similaire aux autres réponses:

données de tracé et ajustement manuel

pour "adapter" cette fonction en utilisant Stan, nous devons définir notre modèle en utilisant son propre langage qui est un mélange entre R et C ++. un modèle simple serait quelque chose comme:

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

qui, espérons-le, se lit bien. nous avons un databloc qui définit les données que nous attendons lorsque nous échantillonnons le modèle, parametersdéfinissons les éléments qui sont échantillonnés et modeldéfinit la fonction de vraisemblance. Vous dites à Stan de "compiler" le modèle, ce qui prend un certain temps, puis vous pouvez en échantillonner avec quelques données. par exemple:

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 facilite les jolis diagrammes de diagnostic, tandis que l'impression de l'ajustement vous donne un bon résumé des paramètres de style 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

l'écart-type important betaindique que les données ne fournissent vraiment pas beaucoup d'informations sur ce paramètre. certaines des réponses donnant plus de 10 chiffres significatifs dans leurs ajustements de modèle surestiment quelque peu les choses

parce que certaines réponses ont noté que chaque virus peut avoir besoin de ses propres paramètres, j'ai étendu le modèle pour autoriser alphaet betavarier en fonction du "virus". tout devient un peu compliqué, mais les deux virus ont presque certainement des alphavaleurs différentes (c'est-à-dire que vous avez besoin de plus de copies / μL de RRAV pour la même couverture) et un graphique montrant ceci est:

tracé des données et des échantillons MC

les données sont les mêmes qu'avant, mais j'ai tracé une courbe pour 40 échantillons de la partie postérieure. UMAVsemble relativement bien déterminé, alors queRRAV pouvant suivre la même pente et nécessiter un nombre de copies plus élevé, ou avoir une pente plus raide et un nombre de copies similaire. la majeure partie de la masse postérieure est sur le fait d'avoir besoin d'un nombre de copies plus élevé, mais cette incertitude pourrait expliquer certaines des différences dans d'autres réponses pour trouver des choses différentes

J'ai surtout utilisé la réponse à cela comme un exercice pour améliorer ma connaissance de Stan, et j'ai mis un cahier Jupyter ici au cas où quelqu'un serait intéressé / voudrait reproduire cela.

Sam Mason
la source
14

(Modifié en tenant compte des commentaires ci-dessous. Merci à @BenBolker et @WeiwenNg pour leur contribution utile.)

Ajustez une régression logistique fractionnée aux données. Il est bien adapté aux données en pourcentage limitées entre 0 et 100% et est bien justifié théoriquement dans de nombreux domaines de la biologie.

Notez que vous devrez peut-être diviser toutes les valeurs par 100 pour l'adapter, car les programmes s'attendent souvent à ce que les données varient entre 0 et 1. Et comme le recommande Ben Bolker, pour résoudre les problèmes possibles causés par les hypothèses strictes de la distribution binomiale concernant la variance, utilisez un distribution quasibinomiale à la place.

J'ai fait quelques hypothèses basées sur votre code, par exemple qu'il y a 2 virus qui vous intéressent et ils peuvent montrer des modèles différents (c'est-à-dire qu'il peut y avoir une interaction entre le type de virus et le nombre de copies).

Tout d'abord, le modèle convient:

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

Si vous faites confiance aux valeurs de p, la sortie ne suggère pas que les deux virus diffèrent de manière significative. Cela contraste avec les résultats de @ NickCox ci-dessous, bien que nous ayons utilisé différentes méthodes. Je ne serais pas très confiant de toute façon avec 30 points de données.

Deuxièmement, le tracé:

Il n'est pas difficile de coder un moyen de visualiser la sortie vous-même, mais il semble y avoir un package ggPredict qui fera la plupart du travail pour vous (je ne peux pas le garantir, je ne l'ai pas essayé moi-même). Le code ressemblera à quelque chose comme:

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

Mise à jour: je ne recommande plus le code ou la fonction ggPredict plus généralement. Après l'avoir essayé, j'ai constaté que les points tracés ne reflètent pas exactement les données d'entrée mais sont plutôt modifiés pour une raison bizarre (certains des points tracés étaient supérieurs à 1 et inférieurs à 0). Je recommande donc de le coder vous-même, bien que ce soit plus de travail.

mkt - Réintégrer Monica
la source
7
J'approuve cette réponse, mais j'aimerais apporter une précision: j'appellerais cette régression logistique fractionnée. Je pense que ce terme serait plus largement reconnu. Quand la plupart des gens entendent "régression logistique", je parie qu'ils pensent à une variable dépendante 0/1. Une bonne réponse de Stackexchange traitant de cette nomenclature est ici: stats.stackexchange.com/questions/216122/…
Weiwen Ng
2
@teaelleceecee Vous devez évidemment diviser la couverture par 100 en premier.
Nick Cox
4
utiliser family=quasibinomial()pour éviter l'avertissement (et les problèmes sous-jacents avec des hypothèses de variance trop strictes). Suivez les conseils de @ mkt sur l'autre problème.
Ben Bolker
2
Cela peut fonctionner, mais je voudrais avertir les gens que vous devez avoir une prémisse avant d'ajuster une fonction que vos données devraient en fait suivre cette fonction. Sinon, vous tirez à peu près au hasard lorsque vous choisissez une fonction appropriée, et vous pouvez être trompé par les résultats.
Carl Witthoft
6
@CarlWitthoft Nous entendons le sermon mais sommes des pécheurs en dehors du service. Quelle prémisse antérieure vous a amené à suggérer une fonction Heaviside dans d'autres commentaires? Ici, la biologie ne ressemble pas à une transition à un seuil précis. Si je comprends bien, la recherche ici est que la théorie formelle est plus faible que les données. Je suis d'accord: si les gens pensent qu'une fonction échelonnée a du sens, ils devraient en avoir une.
Nick Cox
11

Ce n'est pas une réponse différente de @mkt mais les graphiques en particulier ne rentreront pas dans un commentaire. J'adapte d'abord une courbe logistique dans Stata (après avoir enregistré le prédicteur) à toutes les données et j'obtiens ce graphique

entrez la description de l'image ici

Une équation est

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

Maintenant, j'adapte les courbes séparément pour chaque virus dans le scénario le plus simple de virus définissant une variable indicatrice. Voici pour l'enregistrement un script 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) 

Cela pousse dur sur un petit ensemble de données, mais la valeur P pour le virus semble favorable à l'ajustement conjoint de deux courbes.

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

entrez la description de l'image ici

Nick Cox
la source
3

Essayez la fonction sigmoïde . Il existe de nombreuses formulations de cette forme, y compris une courbe logistique. La tangente hyperbolique est un autre choix populaire.

Compte tenu des parcelles, je ne peux pas non plus exclure une fonction d'étape simple. Je crains que vous ne puissiez pas faire la différence entre une fonction pas à pas et un certain nombre de spécifications sigmoïdes. Vous n'avez aucune observation où votre pourcentage se situe dans la plage de 50%, donc la formulation d'étape simple peut être le choix le plus parcimonieux qui ne fait pas pire que les modèles plus complexes

Aksakal
la source
σ(X)=12(1+tanhX2)
2
@JG "sigmoid" est un terme générique pour une courbe en S, en ce qui me concerne, mais vous avez raison de pointer vers un lien entre deux spécifications d'un sigmoid
Aksakal
2

Voici les ajustements 4PL (logistique à 4 paramètres), à la fois contraints et non contraints, avec l'équation selon CA Holstein, M. Griffin, J. Hong, PD Sampson, «Méthode statistique pour déterminer et comparer les limites de détection des essais biologiques», Anal . Chem. 87 (2015) 9795-9801. L'équation 4PL est montrée dans les deux figures et les significations des paramètres sont les suivantes: a = asymptote inférieure, b = facteur de pente, c = point d'inflexion, et d = asymptote supérieure.

La figure 1 contraint a à 0% et d à 100%:

Fig. 1 Contrainte a & d

La figure 2 n'a aucune contrainte sur les 4 paramètres de l'équation 4PL:

Fig.2 Pas de contraintes

C'était amusant, je ne prétends rien savoir de biologique et il sera intéressant de voir comment tout cela s'arrange!

Ed V
la source
Merci, c'est vraiment utile. Je me demandais, avez-vous fait cela dans MATLAB avec la fonction fit?
teaelleceecee
1
J'ai utilisé Igor Pro avec la fonction utilisateur définie par l'utilisateur illustrée dans les figures. J'utilise Igor Pro et son prédécesseur (Igor) depuis 1988, mais de nombreux autres programmes peuvent faire l'ajustement de courbe, par exemple Origin Pro et le Kaleidagraph très bon marché. Et il semble que vous ayez un accès R et (éventuellement?) À Matlab, dont je ne connais rien, sauf qu'ils sont extrêmement capables. Meilleur succès avec cela et j'espère que vous obtiendrez de bonnes nouvelles la prochaine fois que vous discuterez des choses avec les superviseurs! Merci aussi d'avoir publié les données!
Ed V
2

J'ai extrait les données de votre nuage de points et ma recherche d'équation a révélé une équation de type logistique à 3 paramètres comme un bon candidat: "y = a / (1.0 + b * exp (-1.0 * c * x))", où " x "est la base de journal 10 par parcelle. Les paramètres ajustés étaient a = 9,0005947126706630E + 01, b = 1,2831794858584102E + 07, et c = 6,6483431489473155E + 00 pour mes données extraites, un ajustement des données originales (log 10 x) devrait donner des résultats similaires si vous ré-ajustez les données originales en utilisant mes valeurs comme estimations de paramètres initiaux. Mes valeurs de paramètres donnent R au carré = 0,983 et RMSE = 5,625 sur les données extraites.

terrain

MODIFIER: Maintenant que la question a été modifiée pour inclure les données réelles, voici un graphique utilisant l'équation à 3 paramètres ci-dessus et les estimations des paramètres initiaux.

plot2

James Phillips
la source
Il semble y avoir eu une erreur dans votre extraction de données: vous avez un tas de valeurs de pourcentage négatives. De plus, vos valeurs maximales sont d'environ 90% au lieu de 100% comme dans le tracé d'origine. Vous pouvez avoir tout compensé par environ 10% pour une raison quelconque.
mkt
Meh - il s'agit de données extraites semi-manuellement, les données d'origine sont requises. C'est généralement suffisant pour les recherches d'équations, et bien sûr pas pour les résultats finaux - c'est pourquoi j'ai dit d'utiliser mes valeurs de paramètres d'extraction-ajustement comme estimations de paramètres initiales sur les données originales.
James Phillips
Veuillez noter que comme les données réelles ont maintenant été ajoutées au message, j'ai mis à jour cette réponse en utilisant les données mises à jour.
James Phillips
Pour rappel: l'application, par exemple, d'une fonction Heaviside, peut produire des valeurs d'erreur similaires.
Carl Witthoft
1
@JamesPhillips J'essaierai de le faire (Heaviside -> barres d'erreur ou équivalent)
Carl Witthoft
2

Depuis que j'ai dû ouvrir ma grande bouche sur Heaviside, voici les résultats. J'ai défini le point de transition sur log10 (viruscopies) = 2,5. Ensuite, j'ai calculé les écarts-types des deux moitiés de l'ensemble de données - c'est-à-dire que le Heaviside suppose que les données de chaque côté ont toutes les dérivées = 0.

Côté droit std dev = 4,76
Côté gauche std dev = 7,72

Puisqu'il s'avère qu'il y a 15 échantillons dans chaque lot, le dev std global est la moyenne, ou 6,24.

En supposant que le «RMSE» cité dans d'autres réponses est globalement «l'erreur RMS», la fonction Heaviside semblerait faire au moins aussi bien, sinon mieux que la plupart des «courbes en Z» (empruntées à la nomenclature des réponses photographiques). ici.

Éditer

Graphique inutile, mais demandé dans les commentaires:

Ajustement courbe Heaviside

Carl Witthoft
la source
Souhaitez-vous s'il vous plaît poster un modèle et un diagramme de dispersion de la même manière que ce qui a été fait dans les autres réponses? Je suis très curieux de voir ces résultats et de les comparer. Veuillez également ajouter les valeurs RMSE et R au carré pour la comparaison. Personnellement, je n'ai jamais utilisé la fonction Heaviside et je trouve cela très intéressant.
James Phillips
R2
Mon sens était de faire une intrigue similaire à celles faites dans les autres réponses, à des fins de comparaison directe avec ces réponses.
James Phillips
2
@JamesPhillips, il vous reste deux souhaits. Choisissez judicieusement :-)
Carl Witthoft
Merci pour l'intrigue. J'observe que dans toutes les parcelles des autres réponses, l'équation tracée suit la forme courbe des données en haut à droite - ce n'est pas le cas, tout comme la nature de la fonction Heaviside. Cela semble visuellement contredire votre affirmation selon laquelle la fonction Heaviside ferait aussi bien que les équations affichées dans les autres réponses - c'est pourquoi j'avais précédemment demandé les valeurs RMSE et R au carré, je soupçonnais que la fonction Heaviside ne suivrait pas la forme des données dans cette région et pourrait donner des valeurs pires pour ces statistiques d'ajustement.
James Phillips