Ajuster une courbe sigmoïdale aux points avec ggplot

Question

J'ai une base de données simple pour les mesures de réponse d'un traitement médicamenteux à différentes doses:

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)

Je peux tracer ces points en utilisant le code suivant et obtenir les bases de la visualisation que je voudrais:

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 prochaine chose que je voudrais faire avec ces données est d'ajouter une courbe sigmoïdale au graphique, qui correspond aux points tracés pour chaque médicament. Par la suite, je voudrais calculer la CE50 pour cette courbe. Je me rends compte que je n'ai peut-être pas toute la plage de la courbe sigmoïdale dans mes données, mais j'espère obtenir la meilleure estimation possible avec ce que j'ai. En outre, le point final pour le médicament_1 ne suit pas la tendance attendue d'une courbe sigmoïdale, mais ce n'est en fait pas inattendu car les solutions dans lesquelles se trouve le médicament peuvent inhiber les réponses à des concentrations élevées (chaque médicament est dans une solution différente). Je voudrais exclure ce point des données.

Je suis coincé à l'étape d'ajustement d'une courbe sigmoïdale à mes données. J'ai examiné d'autres solutions pour ajuster les courbes sigmoïdales aux données, mais aucune ne semble fonctionner.

Un article qui est très proche de mon problème est le suivant: (sigmoïde) ajustement de la courbe glm en r

Sur cette base, j'ai essayé:

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

Cela donne l'erreur suivante et semble par défaut pour tracer des lignes droites:

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

J'ai également essayé la solution de ce lien: Ajuster une courbe sigmoïdale à ces données oxy-Hb

Dans ce cas, j'obtiens l'erreur suivante:

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

et aucune ligne n'est ajoutée au tracé.

J'ai essayé de rechercher ces deux erreurs mais je n'arrive pas à trouver une raison qui ait du sens avec mes données.

Toute aide serait très appréciée!

Answer 1

Comme je l'ai dit dans un commentaire, je n'utiliserais que geom_smooth()pour un problème très simple; dès que je rencontre des problèmes, j'utilise à la nlsplace.

Ma réponse est très similaire à celle de @ Duck, avec les différences suivantes:

• Je montre des ajustements non pondérés et pondérés (à variance inverse).
• Pour que les ajustements pondérés fonctionnent, j'ai dû utiliser le nls2package, qui fournit un algorithme légèrement plus robuste
• J'utilise SSlogis()pour obtenir la sélection automatique des paramètres initiaux (auto-démarrant)
• Je fais toutes les prédictions en dehors de ggplot2, puis je les alimentegeom_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)
)

Je crois que la CE50 est équivalente au xmidparamètre ... notez les grandes différences entre les estimations pondérées et non pondérées ...

Answer 2

Je suggérerais la prochaine approche qui est proche de ce que vous voulez. J'ai également essayé avec un paramètre pour vos données en utilisant la binomialfamille, mais il y a quelques problèmes concernant les valeurs entre 0 et 1. Dans ce cas, vous auriez besoin d'une variable supplémentaire afin de déterminer les proportions respectives. Le code dans les lignes suivantes utilise une approximation non linéaire afin d'esquisser votre sortie.

Au départ, les données:

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

Dans les moindres carrés non linéaires, vous devez définir des valeurs initiales pour la recherche de paramètres idéaux. Nous utilisons le code suivant avec la fonction de base nls()pour obtenir ces valeurs initiales:

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

Avec cette approche initiale des paramètres, nous esquissons le graphique en utilisant geom_smooth(). Nous utilisons nls()à nouveau pour trouver les bons paramètres:

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

Le résultat: