Ajustar una curva sigmoidea a puntos con ggplot

Question

Tengo un marco de datos simple para las mediciones de respuesta de un tratamiento farmacológico en varias dosis:

drug <- c("drug_1", "drug_1", "drug_1", "drug_1", "drug_1", 
  "drug_1", "drug_1", "drug_1", "drug_2", "drug_2", "drug_2", 
        "drug_2", "drug_2", "drug_2", "drug_2", "drug_2")

conc <- c(100.00, 33.33, 11.11, 3.70, 1.23, 0.41, 0.14, 
        0.05, 100.00, 33.33, 11.11, 3.70, 1.23, 0.41, 0.14, 0.05)

mean_response <- c(1156, 1833, 1744, 1256, 1244, 1088, 678, 489, 
        2322, 1867, 1333, 944, 567, 356, 200, 177)

std_dev <- c(117, 317, 440, 200, 134, 38, 183, 153, 719,
      218, 185, 117, 166, 167, 88, 50)

df <- data.frame(drug, conc, mean_response, std_dev)

Puedo trazar estos puntos usando el siguiente código y obtener la base básica de la visualización que me gustaría:

p <- ggplot(data=df, aes(y=mean_response, x= conc, color = drug)) +
  geom_pointrange(aes(ymax = (mean_response + std_dev), ymin = (mean_response - std_dev))) +
  scale_x_log10()

p

Lo siguiente que me gustaría hacer con estos datos es agregar una curva sigmoidea al gráfico, que se ajuste a los puntos graficados para cada medicamento. Después de eso, me gustaría calcular la EC50 para esta curva. Me doy cuenta de que es posible que no tenga todo el rango de la curva sigmoidea en mis datos, pero espero obtener la mejor estimación posible con lo que tengo. Además, el punto final para el fármaco_1 no sigue la tendencia esperada de una curva sigmoidea, pero esto en realidad no es inesperado, ya que las soluciones en las que se encuentra el fármaco pueden inhibir las respuestas a altas concentraciones (cada fármaco está en una solución diferente). Me gustaría excluir este punto de los datos.

Me quedo atascado en el paso de ajustar una curva sigmoidea a mis datos. He examinado otras soluciones para ajustar curvas sigmoidales a los datos, pero ninguna parece funcionar.

Una publicación que está muy cerca de mi problema es la siguiente: ajuste de curva (sigmoidea) glm en r

Basado en eso, intenté:

p + geom_smooth(method = "glm", family = binomial, se = FALSE)

Esto da el siguiente error y parece predeterminado para trazar líneas rectas:

`geom_smooth()` using formula 'y ~ x'
Warning message:
Ignoring unknown parameters: family 

También probé la solución de este enlace: Ajustar una curva sigmoidea a estos datos de oxi-Hb

En este caso, aparece el siguiente error:

Computation failed in `stat_smooth()`:
Convergence failure: singular convergence (7) 

y no se agregan líneas al gráfico.

He intentado buscar ambos errores, pero parece que no puedo encontrar una razón que tenga sentido con mis datos.

¡Cualquier ayuda será muy apreciada!

Answer 1

Como dije en un comentario, solo lo usaría geom_smooth()para un problema muy fácil; tan pronto como tengo problemas, utilizo nlsen su lugar.

Mi respuesta es muy similar a la de @ Duck, con las siguientes diferencias:

• Muestro ajustes ponderados tanto no ponderados como (varianza inversa).
• Para que los ajustes ponderados funcionen, tuve que usar el nls2paquete, que proporciona un algoritmo un poco más robusto
• Utilizo SSlogis()para obtener la selección de parámetros inicial automática (autoencendido)
• Hago toda la predicción fuera de ggplot2, luego la introduzco engeom_line()

p1 <- nls(mean_response~SSlogis(conc,Asym,xmid,scal),data=df,
          subset=(drug=="drug_1" & conc<100)
        ## , weights=1/std_dev^2  ## error in qr.default: NA/NaN/Inf ...
          )

library(nls2)
p1B <- nls2(mean_response~SSlogis(conc,Asym,xmid,scal),data=df,
            subset=(drug=="drug_1" & conc<100),
            weights=1/std_dev^2)

p2 <- update(p1,subset=(drug=="drug_2"))
p2B <- update(p1B,subset=(drug=="drug_2"))

pframe0 <- data.frame(conc=10^seq(log10(min(df$conc)),log10(max(df$conc)), length.out=100))
pp <- rbind(
    data.frame(pframe0,mean_response=predict(p1,pframe0),
               drug="drug_1",wts=FALSE),
    data.frame(pframe0,mean_response=predict(p2,pframe0),
               drug="drug_2",wts=FALSE),
    data.frame(pframe0,mean_response=predict(p1B,pframe0),
               drug="drug_1",wts=TRUE),
    data.frame(pframe0,mean_response=predict(p2B,pframe0),
               drug="drug_2",wts=TRUE)
)

library(ggplot2); theme_set(theme_bw())
(ggplot(df,aes(conc,mean_response,colour=drug)) +
 geom_pointrange(aes(ymin=mean_response-std_dev,
                     ymax=mean_response+std_dev)) +
 scale_x_log10() +
 geom_line(data=pp,aes(linetype=wts),size=2)
)

Creo que el EC50 es equivalente al xmidparámetro ... observe las grandes diferencias entre las estimaciones ponderadas y no ponderadas ...

Answer 2

Sugeriría el siguiente enfoque que se acerque a lo que desea. También probé con una configuración para sus datos usando binomialfamilia, pero hay algunos problemas con los valores entre 0 y 1. En ese caso, necesitaría una variable adicional para determinar las proporciones respectivas. El código de las siguientes líneas utiliza una aproximación no lineal para esbozar su salida.

Inicialmente, los datos:

library(ggplot2)
#Data
df <- structure(list(drug = c("drug_1", "drug_1", "drug_1", "drug_1", 
"drug_1", "drug_1", "drug_1", "drug_1", "drug_2", "drug_2", "drug_2", 
"drug_2", "drug_2", "drug_2", "drug_2", "drug_2"), conc = c(100, 
33.33, 11.11, 3.7, 1.23, 0.41, 0.14, 0.05, 100, 33.33, 11.11, 
3.7, 1.23, 0.41, 0.14, 0.05), mean_response = c(1156, 1833, 1744, 
1256, 1244, 1088, 678, 489, 2322, 1867, 1333, 944, 567, 356, 
200, 177), std_dev = c(117, 317, 440, 200, 134, 38, 183, 153, 
719, 218, 185, 117, 166, 167, 88, 50)), class = "data.frame", row.names = c(NA, 
-16L))

En unos mínimos cuadrados no lineales, es necesario definir valores iniciales para la búsqueda de parámetros ideales. Usamos el siguiente código con la función base nls()para obtener esos valores iniciales:

#Drug 1
fm1 <- nls(log(mean_response) ~ log(a/(1+exp(-b*(conc-c)))), df[df$drug=='drug_1',], start = c(a = 1, b = 1, c = 1)) #Drug 2 fm2 <- nls(log(mean_response) ~ log(a/(1+exp(-b*(conc-c)))), df[df$drug=='drug_2',], start = c(a = 1, b = 1, c = 1))

Con este enfoque inicial de parámetros, esbozamos la gráfica utilizando geom_smooth(). Usamos nuevamente nls()para encontrar los parámetros correctos:

#Plot
ggplot(data=df, aes(y=mean_response, x= conc, color = drug)) +
  geom_pointrange(aes(ymax = (mean_response + std_dev), ymin = (mean_response - std_dev))) +
  geom_smooth(data = df[df$drug=='drug_1',],method = "nls", se = FALSE, formula = y ~ a/(1+exp(-b*(x-c))), method.args = list(start = coef(fm1), algorithm='port'), color = "tomato")+ geom_smooth(data = df[df$drug=='drug_2',],method = "nls", se = FALSE,
              formula = y ~ a/(1+exp(-b*(x-c))),
              method.args = list(start = coef(fm0),
                                 algorithm='port'),
              color = "cyan3")

La salida: