Réponse à la dose - Ajustement de la courbe globale à l'aide de R

Question

J'ai les données de réponse à la dose suivantes et je souhaite tracer le modèle de réponse à la dose et la courbe d'ajustement global. [xdata = concentration de médicament; ydata (0-5) = valeurs de réponse à différentes concentrations du médicament]. J'ai tracé la courbe Std sans problème.

Std Courbe en forme de données:

df <- data.frame(xdata = c(1000.00,300.00,100.00,30.00,10.00,3.00,1.00,0.30, 
          0.10,0.03,0.01,0.00), 
       ydata = c(91.8,95.3,100,123,203,620,1210,1520,1510,1520,1590, 
          1620)) 

nls.fit <- nls(ydata ~ (ymax*xdata/(ec50 + xdata)) + Ns*xdata + ymin, data=df, 
       start=list(ymax=1624.75, ymin = 91.85, ec50 = 3, Ns = 0.2045514)) 

Courbe de réponse dose en forme de données:

df <- data.frame(
     xdata = c(10000,5000,2500,1250,625,312.5,156.25,78.125,39.063,19.531,9.766,4.883, 
       2.441,1.221,0.610,0.305,0.153,0.076,0.038,0.019,0.010,0.005), 
     ydata1 = c(97.147, 98.438, 96.471, 73.669, 60.942, 45.106, 1.260, 18.336, 9.951, 2.060, 
        0.192, 0.492, -0.310, 0.591, 0.789, 0.075, 0.474, 0.278, 0.399, 0.217, 1.021, -1.263), 
     ydata2 = c(116.127, 124.104, 110.091, 111.819, 118.274, 78.069, 52.807, 40.182, 26.862, 
        15.464, 6.865, 3.385, 10.621, 0.299, 0.883, 0.717, 1.283, 0.555, 0.454, 1.192, 0.155, 1.245), 
     ydata3 = c(108.410, 127.637, 96.471, 124.903, 136.536, 104.696, 74.890, 50.699, 47.494, 23.866, 
        20.057, 10.434, 2.831, 2.261, 1.085, 0.399, 1.284, 0.045, 0.376, -0.157, 1.158, 0.281), 
     ydata4 = c(107.281, 118.274, 99.051, 99.493, 104.019, 99.582, 87.462, 75.322, 47.393, 42.459, 
        8.311, 23.155, 3.268, 5.494, 2.097, 2.757, 1.438, 0.655, 0.782, 1.128, 1.323, 0.645), 
     ydata0 = c(109.455, 104.989, 101.665, 101.205, 108.410, 101.573, 119.375, 101.757, 65.660, 35.672, 
        31.613, 12.323, 25.515, 17.283, 7.170, 2.771, 2.655, 0.491, 0.290, 0.535, 0.298, 0.106)) 

Quand j'ai essayé d'obtenir les paramètres d'ajustement en utilisant le script R ci-dessous, je reçois l'erreur suivante:

Erreur dans nls (yda TA1 ~ BAS + (TOP - BAS)/(1 + 10^((logEC50 - xdata) *:
gradient singulier

nls.fit1 <- nls(ydata1 ~ BOTTOM + (TOP-BOTTOM)/(1+10**((logEC50-xdata)*hillSlope)), data=df, 
       start=list(TOP = max(df$ydata1), BOTTOM = min(df$ydata1),hillSlope = 1.0, logEC50 = 4.310345e-08)) 

nls.fit2 <- nls(ydata2 ~ BOTTOM + (TOP-BOTTOM)/(1+10**((logEC50-xdata)*hillSlope)), data=df, 
       start=list(TOP = max(df$ydata2), BOTTOM = min(df$ydata2),hillSlope = 1.0, logEC50 = 4.310345e-08)) 

nls.fit3 <- nls(ydata3 ~ BOTTOM + (TOP-BOTTOM)/(1+10**((logEC50-xdata)*hillSlope)), data=df, 
       start=list(TOP = max(df$ydata3), BOTTOM = min(df$ydata3),hillSlope = 1.0, logEC50 = 4.310345e-08)) 

nls.fit4 <- nls(ydata4 ~ BOTTOM + (TOP-BOTTOM)/(1+10**((logEC50-xdata)*hillSlope)), data=df, 
       start=list(TOP = max(df$ydata4), BOTTOM = min(df$ydata4),hillSlope = 1.0, logEC50 = 4.310345e-08)) 

nls.fit5 <- nls(ydata0 ~ BOTTOM + (TOP-BOTTOM)/(1+10**((logEC50-xdata)*hillSlope)), data=df, 
       start=list(TOP = max(df$ydata0), BOTTOM = min(df$ydata0),hillSlope = 1.0, logEC50 = 4.310345e-08)) 

S'il vous plaît me conseiller sur la façon de résoudre ce problème

Answer 1

Première note que le rapport de la plus grande valeur de xdata au plus petit est de 2 millions donc nous voulons probablement utiliser log(xdata) à la place de xdata.

Maintenant, en faisant cette modification, nous obtenons le modèle log-logistic LL2.4 de 4 paramètres du drc package mais avec un paramétrage légèrement différent que dans la question. En supposant que vous êtes d'accord avec ces changements, nous pouvons adapter le premier modèle comme suit. Voir ?LL2.4 pour les détails du paramétrage et voir les exemples pertinents au bas de ?ryegrass. Ici df est le df montré dans la question - le modèle LL2.4 lui-même rend la transformation log(xdata).

library(drc) 

fm1 <- drm(ydata1 ~ xdata, data = df, fct = LL2.4()) 
fm1 
plot(fm1) 

Ici nous ajustons tous les 5 modèles et visuellement nous voyons de l'intrigue à la fin que les ajustements sont assez bons.

library(drc) 

fun <- function(yname) { 
    fo <- as.formula(paste(yname, "~ xdata")) 
    fit <- do.call("drm", list(fo, data = quote(df), fct = quote(LL2.4()))) 
    plot(fit) 
    fit 
} 

par(mfrow = c(3, 2)) 
L <- Map(fun, names(df)[-1]) 
par(mfrow = c(1, 1)) 

sapply(L, coef) 

donnant:

     ydata1 ydata2 ydata3 ydata4 ydata0 
    b:(Intercept) -1.37395 -1.1411 -1.1337 -1.0633 -1.6525 
    c:(Intercept) 0.70388 1.9364 1.5800 1.3751 5.7010 
    d:(Intercept) 101.02741 122.0825 120.8042 108.2420 107.9106 
    e:(Intercept) 6.17225 5.0686 4.3215 3.7139 3.2813 

et les ajustements graphiques suivants (cliquez sur l'image pour l'élargir):