Anpassen einer Sigmoidkurve an Punkte mit ggplot

Question

Ich habe einen einfachen Datenrahmen für die Ansprechmessungen einer medikamentösen Behandlung in verschiedenen Dosen:

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)

Ich kann diesen Punkt mit dem folgenden Code zeichnen und die grundlegende Grundlage für die Visualisierung erhalten, die ich möchte:

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

Das nächste, was ich mit diesen Daten machen möchte, ist, dem Diagramm eine Sigmoidkurve hinzuzufügen, die zu den aufgezeichneten Punkten für jedes Medikament passt. Anschließend möchte ich die EC50 für diese Kurve berechnen. Mir ist klar, dass ich möglicherweise nicht den gesamten Bereich der Sigmoidkurve in meinen Daten habe, aber ich hoffe, mit dem, was ich habe, die bestmögliche Schätzung zu erhalten. Auch der Endpunkt für Arzneimittel_1 folgt nicht dem erwarteten Trend einer Sigmoidkurve, aber dies ist tatsächlich nicht unerwartet, da die Lösungen, in denen sich das Arzneimittel befindet, Reaktionen bei hohen Konzentrationen hemmen können (jedes Arzneimittel befindet sich in einer anderen Lösung). Ich möchte diesen Punkt aus den Daten ausschließen.

Ich stecke beim Schritt des Anpassens einer Sigmoidkurve an meine Daten fest. Ich habe mir einige andere Lösungen zum Anpassen von Sigmoidkurven an Daten angesehen, aber keine scheint zu funktionieren.

Ein Beitrag, der meinem Problem sehr nahe kommt, ist folgender: (Sigmoid-) Kurvenanpassung glm in r

Basierend darauf habe ich versucht:

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

Dies führt zu folgendem Fehler und scheint standardmäßig gerade Linien zu zeichnen:

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

Ich habe auch die Lösung über diesen Link ausprobiert: Anpassen einer Sigmoidkurve an diese Oxy-Hb-Daten

In diesem Fall wird folgende Fehlermeldung angezeigt:

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

und dem Plot werden keine Linien hinzugefügt.

Ich habe versucht, diese beiden Fehler nachzuschlagen, kann jedoch keinen Grund finden, der für meine Daten sinnvoll ist.

Jede Hilfe wäre sehr dankbar!

Answer 1

Wie ich in einem Kommentar sagte, würde ich nur geom_smooth()für ein sehr einfaches Problem verwenden; Sobald ich in Schwierigkeiten gerate, benutze ich nlsstattdessen.

Meine Antwort ist @ Duck's sehr ähnlich, mit den folgenden Unterschieden:

• Ich zeige sowohl ungewichtete als auch (inverse Varianz) gewichtete Anpassungen.
• Um die gewichteten Anpassungen zum Laufen zu bringen, musste ich das nls2Paket verwenden, das einen etwas robusteren Algorithmus bietet
• Ich verwende SSlogis(), um eine automatische (selbststartende) anfängliche Parameterauswahl zu erhalten
• Ich mache die ganze Vorhersage außerhalb von ggplot2und speise sie dann eingeom_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)
)

Ich glaube, der EC50 entspricht dem xmidParameter ... beachten Sie die großen Unterschiede zwischen gewichteten und ungewichteten Schätzungen ...

Answer 2

Ich würde den nächsten Ansatz vorschlagen, der nahe an dem liegt, was Sie wollen. Ich habe auch versucht, eine Einstellung für Ihre Daten mithilfe der binomialFamilie vorzunehmen, aber es gibt einige Probleme mit Werten zwischen 0 und 1. In diesem Fall benötigen Sie eine zusätzliche Variable, um die jeweiligen Proportionen zu bestimmen. Der Code in den folgenden Zeilen verwendet eine nichtlineare Näherung, um Ihre Ausgabe zu skizzieren.

Zunächst die Daten:

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 nichtlinearen kleinsten Quadraten müssen Sie Anfangswerte für die Suche nach idealen Parametern definieren. Wir verwenden den nächsten Code mit Basisfunktion nls(), um diese Anfangswerte zu erhalten:

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

Mit diesem ersten Ansatz von Parametern skizzieren wir das Diagramm mit geom_smooth(). Wir verwenden wieder nls(), um die richtigen Parameter zu finden:

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

Die Ausgabe: