Adattare una curva sigmoidale a punti con ggplot

Aug 25 2020

Ho un semplice dataframe per le misurazioni della risposta da un trattamento farmacologico a varie dosi:

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)

Posso tracciare questi punti utilizzando il codice seguente e ottenere le basi di base della visualizzazione che vorrei:

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

La prossima cosa che vorrei fare con questi dati è aggiungere una curva sigmoidale al grafico, che si adatti ai punti tracciati per ogni farmaco. Successivamente, vorrei calcolare l'EC50 per questa curva. Mi rendo conto che potrei non avere l'intera gamma della curva sigmoidale nei miei dati, ma spero di ottenere la migliore stima possibile con quello che ho. Inoltre, il punto finale per il farmaco_1 non segue l'andamento previsto di una curva sigmoidale, ma questo in realtà non è inaspettato poiché le soluzioni in cui si trova il farmaco possono inibire le risposte ad alte concentrazioni (ogni farmaco è in una soluzione diversa). Vorrei escludere questo punto dai dati.

Mi sto bloccando nella fase di adattamento di una curva sigmoidale ai miei dati. Ho esaminato alcune altre soluzioni per adattare le curve sigmoidali ai dati, ma nessuna sembra funzionare.

Un post che è molto vicino al mio problema è questo: (sigmoide) curve fitting glm in r

Sulla base di esso, ho provato:

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

Questo dà il seguente errore e sembra che l'impostazione predefinita sia la stampa di linee rette:

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

Ho anche provato la soluzione da questo collegamento: Adattare una curva sigmoidale a questi dati ossi-Hb

In questo caso, ottengo il seguente errore:

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

e nessuna riga viene aggiunta al grafico.

Ho provato a cercare entrambi questi errori ma non riesco a trovare una ragione che abbia senso con i miei dati.

Qualsiasi aiuto sarebbe molto apprezzato!

Risposte

2 BenBolker Aug 25 2020 at 06:27

Come ho detto in un commento, lo userei solo geom_smooth()per un problema molto semplice; appena vado nei guai lo uso nlsinvece.

La mia risposta è molto simile a quella di @ Duck, con le seguenti differenze:

Mostro sia gli accoppiamenti non ponderati che quelli ponderati (varianza inversa).
Per far funzionare gli adattamenti ponderati, ho dovuto utilizzare il nls2pacchetto, che fornisce un algoritmo leggermente più robusto
Uso SSlogis()per ottenere la selezione iniziale dei parametri automatica (autoaccensione)
Faccio tutte le previsioni al di fuori di ggplot2, poi le inseriscogeom_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)
)

Credo che l'EC50 sia equivalente al xmidparametro ... notare le grandi differenze tra stime ponderate e non ponderate ...

1 Duck Aug 25 2020 at 05:27

Suggerirei il prossimo approccio che è vicino a ciò che desideri. Ho anche provato con un'impostazione per i tuoi dati usando la binomialfamiglia, ma ci sono alcuni problemi sui valori tra 0 e 1. In tal caso avresti bisogno di una variabile aggiuntiva per determinare le rispettive proporzioni. Il codice nelle righe seguenti utilizza un'approssimazione non lineare per tracciare l'output.

Inizialmente, i dati:

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

In un minimo quadrato non lineare, è necessario definire i valori iniziali per la ricerca dei parametri ideali. Usiamo il codice successivo con la funzione di base nls()per ottenere quei valori iniziali: