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:
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))
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
la source
method.args=list(family=quasibinomial))
les arguments àgeom_smooth()
votre code ggplot d'origine.se=FALSE
. Toujours agréable de montrer aux gens à quel point l'incertitude est grande ...Réponses:
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:
où
alpha
dé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 platespour montrer à quoi cela ressemble, nous pouvons extraire vos données et les tracer avec:
où
raw_data.txt
contient 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: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:
qui, espérons-le, se lit bien. nous avons un
data
bloc qui définit les données que nous attendons lorsque nous échantillonnons le modèle,parameters
définissons les éléments qui sont échantillonnés etmodel
dé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: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:l'écart-type important
beta
indique 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 chosesparce 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
alpha
etbeta
varier en fonction du "virus". tout devient un peu compliqué, mais les deux virus ont presque certainement desalpha
valeurs 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:les données sont les mêmes qu'avant, mais j'ai tracé une courbe pour 40 échantillons de la partie postérieure.
UMAV
semble 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érentesJ'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.
la source
(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:
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: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.
la source
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.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
Une équation est
100
invlogit
(-4,192654 + 1,880951log10
(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:
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.
la source
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
la source
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%:
La figure 2 n'a aucune contrainte sur les 4 paramètres de l'équation 4PL:
C'était amusant, je ne prétends rien savoir de biologique et il sera intéressant de voir comment tout cela s'arrange!
la source
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.
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.
la source
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:
la source