¿Qué tipo de curva (o modelo) debo ajustar a mis datos porcentuales?

Question

Estoy tratando de crear una figura que muestre la relación entre las copias virales y la cobertura del genoma (GCC). Este es el aspecto de mis datos:

Al principio, me limité a trazar una regresión lineal, pero mis supervisores me dijeron que eso era incorrecto y que probara con una curva sigmoidal. Así que lo hice utilizando 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))

Sin embargo, mis supervisores dicen que esto también es incorrecto porque las curvas hacen parecer que GCC puede superar el 100%, lo cual no es posible.

Mi pregunta es: ¿cuál es la mejor manera de mostrar la relación entre las copias de virus y el CCG? Quiero dejar claro que A) bajas copias de virus = baja CCG, y que B) después de una cierta cantidad de copias de virus la CCG se estabiliza.

He investigado un montón de métodos diferentes - GAM, LOESS, logístico, a trozos - pero no sé cómo decir cuál es el mejor método para mis datos.

EDIT: estos son los datos:

>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

Answer 1

(Editado teniendo en cuenta los comentarios de abajo. Gracias a @BenBolker y @WeiwenNg por sus útiles aportaciones).

Ajuste una regresión logística fraccional a los datos. Se adapta bien a los datos porcentuales acotados entre 0 y 100% y está bien justificada teóricamente en muchas áreas de la biología.

Tenga en cuenta que es posible que tenga que dividir todos los valores por 100 para ajustarlos, ya que los programas suelen esperar que los datos oscilen entre 0 y 1. Y, como recomienda Ben Bolker, para abordar los posibles problemas causados por las estrictas suposiciones de la distribución binomial con respecto a la varianza, utilice en su lugar una distribución cuasibinomial.

He hecho algunas suposiciones basadas en tu código, como que hay 2 virus que te interesan y que pueden mostrar patrones diferentes (es decir, puede haber una interacción entre el tipo de virus y el número de copias).

En primer lugar, el ajuste del modelo:

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 se confía en los valores p, el resultado no sugiere que los dos virus difieran significativamente. Esto contrasta con los resultados de @NickCox más abajo, aunque utilizamos métodos diferentes. Yo no estaría muy seguro de ninguna de las dos cosas con 30 puntos de datos.

En segundo lugar, la trama:

No es difícil codificar una forma de visualizar el resultado por ti mismo, pero parece que hay un paquete ggPredict que hará la mayor parte del trabajo por ti (no puedo responder por él, no lo he probado yo mismo). El código será algo así:

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

Actualización: Ya no recomiendo el código ni la función ggPredict en general. Después de probarlo, descubrí que los puntos trazados no reflejan exactamente los datos de entrada, sino que se modifican por alguna extraña razón (algunos de los puntos trazados estaban por encima de 1 y por debajo de 0). Así que recomiendo que lo codifiques tú mismo, aunque eso es más trabajo.

Answer 2

Esta no es una respuesta diferente a la de @mkt pero los gráficos en particular no caben en un comentario. Primero ajusto una curva logística en Stata (después de registrar el predictor) a todos los datos y obtengo este gráfico

Una ecuación es

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

Ahora ajusto las curvas por separado para cada virus en el escenario más sencillo de que el virus defina una variable indicadora. Aquí, para que conste, hay un script de 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) 

Esto es un esfuerzo en un pequeño conjunto de datos, pero el valor P para el virus parece apoyar el ajuste de dos curvas conjuntamente.

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

Answer 3

Otra forma de hacerlo sería utilizar una formulación bayesiana, puede ser un poco pesado para empezar, pero tiende a hacer que sea mucho más fácil expresar los detalles de su problema, así como obtener mejores ideas de dónde está la "incertidumbre".

Stan es un muestreador de Monte Carlo con una interfaz programática relativamente fácil de usar, hay bibliotecas disponibles para R y otros pero aquí estoy usando Python

usamos una sigmoidea como todo el mundo: tiene motivaciones bioquímicas además de ser matemáticamente muy conveniente para trabajar. una buena parametrización para esta tarea es:

import numpy as np

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

donde alpha define el punto medio de la curva sigmoidea (es decir, donde cruza el 50%) y beta define la pendiente, los valores más cercanos a cero son más planos

para mostrar cómo se ve esto, podemos sacar sus datos y graficarlos con:

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

donde raw_data.txt contiene los datos que diste y transformé la cobertura a algo más útil. los coeficientes 5,5 y 3 se ven bien y dan un gráfico muy parecido a las otras respuestas:

para "ajustar" esta función usando Stan necesitamos definir nuestro modelo usando su propio lenguaje que es una mezcla entre R y C++. un modelo simple sería algo como

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

que espero que se lea bien. tenemos un data que define los datos que esperamos cuando muestreamos el modelo, parameters definir las cosas que se muestrean, y model define la función de probabilidad. Le dices a Stan que "compile" el modelo, lo que lleva un tiempo, y luego puedes muestrear a partir de él con algunos datos. por ejemplo:

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 hace que los gráficos de diagnóstico sean fáciles, mientras que la impresión del ajuste le da un buen resumen de parámetros al estilo de 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

la gran desviación estándar en beta dice que los datos realmente no proporcionan mucha información sobre este parámetro. también algunas de las respuestas que dan 10+ dígitos significativos en sus ajustes del modelo están exagerando un poco las cosas

Debido a que algunas respuestas señalaron que cada virus podría necesitar sus propios parámetros, amplié el modelo para permitir alpha et beta que varíe según el "Virus". todo se vuelve un poco complicado, pero es casi seguro que los dos virus tienen diferentes alpha (es decir, se necesitan más copias/μL de RRAV para la misma cobertura) y un gráfico que lo demuestre:

los datos son los mismos que antes, pero he dibujado una curva para 40 muestras de la posterior. UMAV parece relativamente bien determinado, mientras que RRAV podría seguir la misma pendiente y necesitar un mayor recuento de copias, o tener una pendiente más pronunciada y un recuento de copias similar. la mayor parte de la masa posterior está en necesitar un mayor recuento de copias, pero esta incertidumbre podría explicar algunas de las diferencias en otras respuestas encontrando cosas diferentes

Principalmente utilicé la respuesta a esto como un ejercicio para mejorar mi conocimiento de Stan, y he puesto un cuaderno Jupyter de esto aquí por si alguien está interesado/quiere replicar esto.

Answer 4

He extraído los datos de tu gráfico de dispersión, y mi búsqueda de ecuaciones ha dado como resultado una ecuación de tipo logístico de 3 parámetros como buen candidato: "y = a / (1,0 + b * exp(-1,0 * c * x))", donde "x" es el logaritmo de base 10 según tu gráfico. Los parámetros ajustados fueron a = 9,0005947126706630E+01, b = 1,2831794858584102E+07, y c = 6,6483431489473155E+00 para mis datos extraídos; un ajuste de los datos originales (log 10 x) debería dar resultados similares si se vuelven a ajustar los datos originales utilizando mis valores como estimaciones iniciales de los parámetros. Los valores de mis parámetros producen un R-cuadrado = 0,983 y un RMSE = 5,625 en los datos extraídos.

EDIT: Ahora que la pregunta ha sido editada para incluir los datos reales, aquí hay un gráfico utilizando la ecuación de 3 parámetros anterior y las estimaciones iniciales de los parámetros.

Answer 5

Prueba con sigmoide función. Hay muchas formulaciones de esta forma, incluyendo una curva logística. La tangente hiperbólica es otra opción popular.

Teniendo en cuenta los gráficos, tampoco puedo descartar una simple función escalonada. Me temo que no podrá diferenciar entre una función escalonada y cualquier número de especificaciones sigmoides. No tienes ninguna observación en la que tu porcentaje esté en el rango del 50%, por lo que la formulación escalonada simple puede ser la opción más parsimoniosa que no se comporta peor que los modelos más complejos