Cálculo EC50
Concentración-respuesta: EC50 con paquete drc
junio 01, 2022
Autor: Álvaro Alonso Fernández
Departamento de Ciencias de la Vida
Universidad de Alcalá (España)
1 Introducción
Vamos a utilizar diferentes modelos concentración-respuesta para poder calcular un parámetro muy utilizado en ecotoxicología la EC (Concentración efectiva). Para ello utilizaremos el paquete drc. En el tutorial utilizaremos ejemplos que trae el propio paquete.
Paper de referencia: Ritz C and Streibig JC (2005) Bioassay Analysis using R. Journal of Statistical Software 12: 1-22. Hay una nueva versión en 2016 de este paquete.
Describe como ajustar varios modelos no-lineales de regresión para curvas de dosis-respuesta. Para ello describe el paquete drc. Entre las varias funciones que incluye permite la comparación de parámetros entre sí (por ejemplo EC10 o EC90).
Manual de drc. Cuidado que desde 2016 hay cambios sustanciales y algunas de las funciones del paper ya no están disponibles.
Según los autores los dos modelos más ampliamente utilizados son EL MODELO LOGÍSTICO (LOGISTIC MODEL) y EL MODELO GOMPERTZ.
Para el modelo logístico hay cuatro parámetros: e (ED50), d (upper limit), c (lower limit), b (relative slope around e):
2 Paquete drc
Instalamos el paquete drc y sus dependencias install.packages("drc")
library(drc)## Loading required package: MASS##
## 'drc' has been loaded.## Please cite R and 'drc' if used for a publication,## for references type 'citation()' and 'citation('drc')'.##
## Attaching package: 'drc'## The following objects are masked from 'package:stats':
##
## gaussian, getInitiallibrary(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.1 v dplyr 1.0.6
## v tidyr 1.2.0 v stringr 1.4.0
## v readr 1.4.0 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x dplyr::select() masks MASS::select()3 Un primer modelo
Vamos a hacer un primer modelo sin demasiado cuidado a la hora de seleccionar el tipo de modelo. Vamos simplemnete a ver como funcionan los principales elementos del paquete drc.
3.1 Nuestros datos
Cargamos nuestros datos. Son datos que muestran el crecimiento de raices sometidos a un herbicida. Para más detalles sobre la naturaleza de los datos se puede consultar:
str(ryegrass)## 'data.frame': 24 obs. of 2 variables:
## $ rootl: num 7.58 8 8.33 7.25 7.38 ...
## $ conc : num 0 0 0 0 0 0 0.94 0.94 0.94 1.88 ...Hacemos el primer modelo, lo llamamos con el original nombre de modelo1
modelo1<- drm(rootl~conc, data=ryegrass, fct=LL.4(names = c("Slope", "Lower Limit", "Upper Limit", "ED50")))El modelo es logístico con cuatro parámetros LL.4 (b, c, d, y e). A cada uno de ellos le damos su nombre correspondiente.
3.2 Representación gráfica y resumen
plot(modelo1, type="all")#al poner "all" representamos todas las observaciones que tenemos en nuestros datos
Le pedimos a R con summary un resumen del modelo que acabamos de construir:
summary(modelo1)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## Slope:(Intercept) 2.98222 0.46506 6.4125 2.960e-06 ***
## Lower Limit:(Intercept) 0.48141 0.21219 2.2688 0.03451 *
## Upper Limit:(Intercept) 7.79296 0.18857 41.3272 < 2.2e-16 ***
## ED50:(Intercept) 3.05795 0.18573 16.4644 4.268e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 0.5196256 (20 degrees of freedom)Aquí ya vemos nuestro primer resultado: la dosis efectiva 50, mejor dicho nuestra concentración efectiva 50 (EC50), que sería 3.05. Para más detalles consulta la referencia antes citada.
3.3 Calculamos ECs
Ahora vamos a calcular la EC10, la EC50 y la EC99 mediante la función ED:
ED(modelo1, c(10,50,99), interval="delta")##
## Estimated effective doses
##
## Estimate Std. Error Lower Upper
## e:1:10 1.46371 0.18677 1.07411 1.85330
## e:1:50 3.05795 0.18573 2.67053 3.44538
## e:1:99 14.27608 3.57468 6.81943 21.73273Con “delta” calculamos los límites al 95% para nuestra CE. La EC50 es 3.05 con límte superior de 3.44 e inferior de 2.67, es decir 3.05(2.67-3.44).
3.4 Algunos consejos
Existen varios tipos de modelos para los datos ecotoxicológicos, que cambian ligeramente unos de otros. Básicamente, el logístico (log logistic) asume simetría respecto a un punto intermedio, los modelos Weibull asumen asimetría respecto a un punto intermedio, etc. Nuestros datos pueden ser continuos (por ejemplo crecimiento, o porcentaje de mortalidad) o binomiales (muerto-vivo, etc.), es decir hay infinidad de posibilidades. Respecto a los parámetros del modelo, los hay de 2, 3, 4 y 5 parámetros. Por ejemplo, si nuetros datos solo pueden tomar un máximo de 1 y un mínimo de cero, entonces fijaríamos el “lower” y el “upper” a 0 y 1, respectivamente. En caso contrario no indicariamos nada. Normalmente cuantos menos parámetros pongamos libres mejor es el modelo resultante.
Resumen de modelos en Ritz C and Streibig JC (2005) Bioassay Analysis using R. Journal of Statistical Software 12: 1-22
4 Hagamos más modelos
Ahora utilizaremos una nueva base de datos del paquete drc. El nombre es finney71 y para poder utilizarla debemos leer nuestros datos y ver su estructura con str:
data(finney71)
finney71## dose total affected
## 1 10.2 50 44
## 2 7.7 49 42
## 3 5.1 46 24
## 4 3.8 48 16
## 5 2.6 50 6
## 6 0.0 49 0str(finney71)## 'data.frame': 6 obs. of 3 variables:
## $ dose : num 10.2 7.7 5.1 3.8 2.6 0
## $ total : int 50 49 46 48 50 49
## $ affected: int 44 42 24 16 6 0Son datos en los que se muestra para cada concentración de tóxico empleada dose el número total de organismos utilizados total y los que han mostrado respuesta a esa concentración affected. Ojo que el número de animales totales en cada concentración no es el mismo.
4.1 Representación gráfica de los datos brutos
Una figura sencilla con el paquete ggplot2 para representar nuestros datos brutos:
library(ggplot2)
gr1<-ggplot(finney71, aes(dose, affected))+ geom_point()
gr1
Parece que si hay una concentración-respuesta. Son los datos brutos (total de afectados en cada concentración, pero el total de animales cambia en cada concentración).
Ahora calculamos el porcentaje de afectados. Para ello añadimos una columna de la siguiente forma utilizando el paquete tidyverse:
library(tidyverse)
datos<-finney71 %>%
mutate(porcentaje = (finney71$affected*100)/finney71$total)
str(datos)## 'data.frame': 6 obs. of 4 variables:
## $ dose : num 10.2 7.7 5.1 3.8 2.6 0
## $ total : int 50 49 46 48 50 49
## $ affected : int 44 42 24 16 6 0
## $ porcentaje: num 88 85.7 52.2 33.3 12 ...datos## dose total affected porcentaje
## 1 10.2 50 44 88.00000
## 2 7.7 49 42 85.71429
## 3 5.1 46 24 52.17391
## 4 3.8 48 16 33.33333
## 5 2.6 50 6 12.00000
## 6 0.0 49 0 0.00000Ya tenemos nuestra columna con la variable respuesta de porcentaje de afectados por el tóxico porcentaje.
4.2 Cuatro modelos
Hagamos cuatro modelos diferentes model1``model2``model3y model4 con el mismo conjunto de datos:
#De cuatro parámetros sin fijar ninguno
model1<-drm(porcentaje ~ dose, data=datos, fct=LL.4())
#Fijamos el lower:
model2<-drm(porcentaje ~ dose, data=datos, fct=L.4(fixed = c(NA, 0, NA, NA), names = c("slope", "lower limit", "upper limit", "EC50")))
#Fijamos lower y upper:
model3<-drm(porcentaje ~ dose, data=datos, fct=L.4(fixed = c(NA, 0, 100, NA), names = c("slope", "lower limit", "upper limit", "EC50")))
#Con tres parámetros
model4<- drm(porcentaje ~ dose, data=datos,fct=LL.3(fixed=c(NA, 100, NA),names = c("Slope", "Upper Limit", "EC50")))4.3 Representación gráfica de los cuatro modelos
Vamos a poner los cuatro modelos juntos para ver las pequeñas diferencias en su representación debido al diferente ajuste:
par(mfrow=c(2,2))
plot(model1,main="Model1 fct=LL.4())",ylim = c(0, 100),col="red")
plot(model2,main="Model2 (fixed = c(NA, 0, NA, NA)",ylim = c(0, 100),col="darkgreen")
plot(model3,main="Model3 (fixed = c(NA, 0, 100, NA)",ylim = c(0, 100),col="blue")
plot(model4,main="Model4 LL.3",ylim = c(0, 100))
4.4 Datos de los modelos
Con la función summary vamos a ver la estimación de los diferentes parámetros del modelo:
summary(model1)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) -3.40895 0.69891 -4.8775 0.039557 *
## c:(Intercept) 0.37135 4.01511 0.0925 0.934740
## d:(Intercept) 96.76762 7.15493 13.5246 0.005423 **
## e:(Intercept) 4.69983 0.32252 14.5722 0.004676 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 4.113979 (2 degrees of freedom)summary(model2)##
## Model fitted: Logistic (ED50 as parameter) (3 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## slope:(Intercept) -0.862373 0.094352 -9.1399 0.002768 **
## upper limit:(Intercept) 89.703984 2.689573 33.3525 5.925e-05 ***
## EC50:(Intercept) 4.590425 0.147723 31.0746 7.322e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 2.920287 (3 degrees of freedom)summary(model3)##
## Model fitted: Logistic (ED50 as parameter) (2 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## slope:(Intercept) -0.66504 0.10365 -6.4161 0.003033 **
## EC50:(Intercept) 5.04656 0.23424 21.5446 2.745e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 5.472125 (4 degrees of freedom)summary(model4)##
## Model fitted: Log-logistic (ED50 as parameter) with lower limit at 0 (2 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## Slope:(Intercept) -3.18011 0.24621 -12.916 0.0002072 ***
## EC50:(Intercept) 4.81614 0.11904 40.458 2.23e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 3.036288 (4 degrees of freedom)Los resultados son bastante parecidos. En el primer modelo la EC50 es la letra e.
4.5 Comparación estadística entre los modelos y cálculo de parámetros
Podemos comparar dos modelos entre sí por medio de anova del siguiente modo:
anova(model1, model2, details=TRUE,test=NULL)##
## 1st model
## fct: L.4(fixed = c(NA, 0, NA, NA), names = c("slope", "lower limit",
## 2nd model
## fct: LL.4()## ANOVA table
##
## ModelDf RSS Df F value p value
## 2nd model 3 25.584
## 1st model 2 33.850 1 -0.4884 1.0000anova(model1, model3, details=TRUE,test=NULL)##
## 1st model
## fct: L.4(fixed = c(NA, 0, 100, NA), names = c("slope", "lower limit",
## 2nd model
## fct: LL.4()## ANOVA table
##
## ModelDf RSS Df F value p value
## 2nd model 4 119.78
## 1st model 2 33.85 2 2.5385 0.2826anova(model2, model3, details=TRUE,test=NULL)##
## 1st model
## fct: L.4(fixed = c(NA, 0, 100, NA), names = c("slope", "lower limit",
## 2nd model
## fct: L.4(fixed = c(NA, 0, NA, NA), names = c("slope", "lower limit",## ANOVA table
##
## ModelDf RSS Df F value p value
## 2nd model 4 119.777
## 1st model 3 25.584 1 11.0450 0.0449anova(model2, model4, details=TRUE,test=NULL)##
## 1st model
## fct: LL.3(fixed = c(NA, 100, NA), names = c("Slope", "Upper Limit",
## 2nd model
## fct: L.4(fixed = c(NA, 0, NA, NA), names = c("slope", "lower limit",## ANOVA table
##
## ModelDf RSS Df F value p value
## 2nd model 4 36.876
## 1st model 3 25.584 1 1.3241 0.3333##etc.El modelo 2 y 3 están en el límite de las diferencias, las demás comparaciones muestran que los modelos son iguales.
Ahora calculamos las EC10, EC50 y EC99 con sus intervalos de confianza al 95% con la función ED:
ED(model1, c(10,50,99), interval="delta") ##
## Estimated effective doses
##
## Estimate Std. Error Lower Upper
## e:1:10 2.46694 0.29300 1.20627 3.72761
## e:1:50 4.69983 0.32252 3.31214 6.08752
## e:1:99 18.09228 5.66244 -6.27125 42.45581ED(model2, c(10,50,99), interval="delta") ##
## Estimated effective doses
##
## Estimate Std. Error Lower Upper
## e:1:10 2.04254 0.24139 1.27434 2.81074
## e:1:50 4.59043 0.14772 4.12031 5.06054
## e:1:99 9.91889 0.66932 7.78880 12.04897ED(model3, c(10,50,99), interval="delta") ##
## Estimated effective doses
##
## Estimate Std. Error Lower Upper
## e:1:10 1.74263 0.47813 0.41514 3.07013
## e:1:50 5.04656 0.23424 4.39621 5.69690
## e:1:99 11.95614 1.18568 8.66416 15.24812ED(model4, c(10,50,99), interval="delta")##
## Estimated effective doses
##
## Estimate Std. Error Lower Upper
## e:1:10 2.41342 0.14247 2.01787 2.80897
## e:1:50 4.81614 0.11904 4.48563 5.14665
## e:1:99 20.42846 2.33813 13.93679 26.920144.6 Grado de ajuste de los modelos: Neill test y Akaike
Para comprobar el grado de ajuste de cada uno de nuestros modelos lo que hacemos es aplicar la función neill.test. Si no es significativo el ajuste del modelo es bueno:
neill.test(model1)## Grouping used
##
## 1 2 3 4 5
## 1 1 1 2 1## Neill's lack-of-fit test
##
## F value p value
## 47.8382 0.0914neill.test(model2)## Grouping used
##
## 1 2 3 4
## 2 1 2 1## Neill's lack-of-fit test
##
## F value p value
## 1.2721 0.3765neill.test(model3)## Grouping used
##
## 1 2 3
## 2 3 1## Neill's lack-of-fit test
##
## F value p value
## 1.9821 0.2539neill.test(model4)## Grouping used
##
## 1 2 3
## 2 3 1## Neill's lack-of-fit test
##
## F value p value
## 0.0628 0.8183En todos ellos encontramos un buen grado de ajuste.
Podemos utilizar el Akaike’s information criterion para comprobar cuál sería el mejor modelo-ajuste. Cuanto más bajo sea el valor mejor es el modelo. Para ello se emplea la función mselect:
mselect(model1, fctList = list(W1.3(fixed=c(NA, 100, NA)),W1.4(), W2.3(fixed=c(NA, 100, NA)), W2.4()),linreg=TRUE) ## logLik IC Lack of fit Res var
## Cubic -10.28180 30.56359 NA 5.408655
## W2.4 -11.21353 32.42707 NA 7.378568
## LL.4 -13.70414 37.40828 NA 16.924822
## W1.3 -16.02586 38.05173 NA 18.348405
## W2.3 -16.56935 39.13870 NA 21.992566
## W1.4 -15.23664 40.47328 NA 28.208280
## Lin -21.13045 48.26090 NA 100.591901
## Quad -20.89667 49.79334 NA 124.067525En este caso se compara el modelo 1 con otros. El IC más bajo es el mejor, pero con diferencias de hasta 10 puntos entre los modelos se consideran bastante parecidos, es decir el LL4 parece adecuado.
Repetimos con model2
mselect(model2, fctList = list(W1.3(fixed=c(NA, 100, NA)),W1.4(), W2.3(fixed=c(NA, 100, NA)), W2.4()),linreg=TRUE)## logLik IC Lack of fit Res var
## Cubic -10.28180 30.56359 NA 5.408655
## W2.4 -11.21353 32.42707 NA 7.378568
## L.4 -12.86428 33.72856 NA 8.528076
## W1.3 -16.02586 38.05173 NA 18.348405
## W2.3 -16.56935 39.13870 NA 21.992566
## W1.4 -15.23664 40.47328 NA 28.208280
## Lin -21.13045 48.26090 NA 100.591901
## Quad -20.89667 49.79334 NA 124.067525En caso de no decidirnos podemos hacer un promedio de modelos a través de la función maED:
maED(model2,
list(W2.4(),LL.4(), W1.4()),c(10, 50, 99),
interval="kang")## ED10 ED50 ED99 Weight
## L.4 2.042543 4.590425 9.918885 0.32153068
## W2.4 2.249853 4.551633 9.242265 0.61636610
## LL.4 2.466943 4.699830 18.092282 0.05107197
## W1.4 2.655020 5.139162 52.766471 0.01103126## Estimate Lower Upper
## e:1:10 2.198754 1.342679 3.054828
## e:1:50 4.578156 3.889834 5.266478
## e:1:99 10.391934 4.102836 16.681032Vemos que en general los resultados son bastante parecidos entre los diferentes modelos. EL valor promedio de EC50 es de 4.5
5 Hagamos modelos con dos tiempos
Vamos a aplicar el paquete a un conjunto de datos que incluyen el factor tiempo. Los datos son proporcionados por el paquete y se denominan daphnids:
data(daphnids)
daphnids## dose no total time
## 1 105.00 0 20 24h
## 2 400.07 1 20 24h
## 3 600.10 3 20 24h
## 4 1199.20 3 20 24h
## 5 1999.33 5 20 24h
## 6 3198.52 6 20 24h
## 7 5596.91 7 20 24h
## 8 9595.57 17 20 24h
## 9 105.00 0 20 48h
## 10 400.07 0 20 48h
## 11 600.10 6 20 48h
## 12 1199.20 8 20 48h
## 13 1999.33 11 20 48h
## 14 3198.52 16 20 48h
## 15 5596.91 18 20 48h
## 16 9595.57 20 20 48hSon datos tomados a 24 y 48 horas de exposición. Los vamos a representar gráficamente para ver el aspecto que tienen:
par(mfrow=c(1,1))
plot(daphnids$no~daphnids$dose, daphnids$time)
Hacemos un modelo con el número de afectados entre el total respecto a la dosis y teniendo en cuenta el tiempo. Optamos por binomial (0-1). norepresenta el número de animales afectados por el tóxico y total el número total de animales utilizados en cada concentración de tóxico (dose). Hay dos tiempos (24 y 48h) (time)
modeldaph <- drm(no/total~dose, time, weights = total,
data = daphnids, fct = LL.2(), type = "binomial")
modelFit(modeldaph)## Goodness-of-fit test
##
## Df Chisq value p value
##
## DRC model 12 13.873 0.3089plot(modeldaph, xlab="Log Dose",ylab="Response (0-1)")
Y ahora los parámetros:
summary(modeldaph)##
## Model fitted: Log-logistic (ED50 as parameter) with lower limit at 0 and upper limit at 1 (2 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:24h -1.17384 0.22236 -5.2791 1.298e-07 ***
## b:48h -1.84968 0.27922 -6.6244 3.488e-11 ***
## e:24h 5134.03344 1056.74197 4.8584 1.184e-06 ***
## e:48h 1509.06539 187.76008 8.0372 9.037e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ED(modeldaph, c(10,50,99), interval="delta")##
## Estimated effective doses
##
## Estimate Std. Error Lower Upper
## e:24h:10 789.83 243.66 312.25 1267.40
## e:24h:50 5134.03 1056.74 3062.86 7205.21
## e:24h:99 257367.77 222101.90 -177943.95 692679.50
## e:48h:10 460.06 103.37 257.46 662.65
## e:48h:50 1509.07 187.76 1141.06 1877.07
## e:48h:99 18097.52 7013.17 4351.95 31843.08Las posibilidades del paquete drc son muchas y es conveniente una lectura adecuada de los pros y contras de cada tipo de modelo. Normalmente las diferencias entre modelos no suelen ser muy grandes, especialmente si tenemos una buena replicacion y suficientes puntos de observación (en este caso 4-5 concentraciones o más de tóxico).
Espero que sea de utilidad.
6 CRÉDITOS
Departamento de Ciencias de la Vida
Universidad de Alcalá (España)
Comentarios
Publicar un comentario