Modelos para Análisis de Bioensayos
Álvaro Alonso Fernández, Department of Life Sciences (Ecology), University of Alcalá, Spain
17 de noviembre de 2017
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 (por ejemplo ED10 o ED90).
Manual de drc. Cuidado que apartir de 2016 hay cambios sustanciales, y algunas de las funciones del paper ya no están disponibles.
Instalamos el paquete drc y sus dependencias install.packages("drc")
MODELOS
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) (Función 1)
Función 1
Existe el mismo modelo con solo tres parámetros, siendo el límite inferior del anterior igual a cero.
En caso de que exista hormesis se debe emplear el modelo de Brain-Cousens (Brain and Cousens 1989).
Una vez seleccionado el modelo y realizado (varias especies, varios compuestos, etc.), es interesante comparar parámetros entre las diferentes curvas. Para ello utilizan un modelo simple que se ajusta y se utiliza el modelo delta (van der Vaart 1998, en el capítulo 3). Con ello se consiguen calcular los errores estandard para las funciones de los parámetros.
Para el caso de los herbicidas se utiliza la EFFECTIVE DOSAGE (ED) y el SELECTIVITY INDEX (SI). Por ejemplo la ED50 sería la mitad de la máxima respuesta alcanzada (que sería el 100%). El SI es la ratio entre ED de dos curvas distintas.
SOBRE EL PAQUETE drc
La principal función es multdrc.
El la Tabla 1 del paper se especifican las diferentes funciones para los diferentes modelos.
Una vez realizado el modelo se pueden realizar varias cosas, por ejemplo un ANOVA entre dos modelos (anova), dibujar los modelos (plot), o plotdrc con las observaciones, todas las curvas juntas, etc.
TRABAJAMOS CON EL PAQUETE drc
Empezamos cargando el paquete:
library(drc)## Warning: package 'drc' was built under R version 3.4.3## Loading required package: MASS## Warning: package 'MASS' was built under R version 3.4.3##
## '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, getInitialTRABAJAMOS CON EL PAQUETE drc: EJEMPLO 1 DATOS FINNEY71
Se emplean datos que estan en el paquete para el ejemplo finney71. For each of six concentration of an insecticid the number of insects affected (out of the number of insects) was recorded. A data frame with 6 observations on the following 3 variables:
dose: a numeric vector
total: a numeric vector
affected: a numeric vector
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 0dose<-(finney71$dose)
affected<-(finney71$affected)
plot(dose, affected)
Ahora hacemos el modelo a través de la función drm. El argumento fct indica el tipo de modelo, en este caso de cuatro parámetros logístico. Para ello hacemos lo siguiente sin poner ninguna restricción de valor máximo ni mínimo:
modelo1<-drm(affected ~ dose, data=finney71, fct=LL.4())
plot(modelo1)
summary(modelo1)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) -3.123049 0.761812 -4.0995 0.054669 .
## c:(Intercept) -0.044097 2.310337 -0.0191 0.986505
## d:(Intercept) 49.988650 5.241182 9.5377 0.010815 *
## e:(Intercept) 4.975754 0.475453 10.4653 0.009007 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 2.344228 (2 degrees of freedom)Ahora lo repetimos limitando el valor máximo a 50 y el mínimo a cero, y dando el nombre a cada una de las variables del modelo. Si se pone NA no se hace ninguna limitación:
modelo2<-drm(affected ~ dose, data=finney71, fct=L.4(fixed = c(NA, 0, 50, NA), names = c("slope", "lower limit", "upper limit", "EC50")))
plot(modelo2)
summary(modelo2)##
## Model fitted: Logistic (ED50 as parameter) (2 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## slope:(Intercept) -0.630776 0.085379 -7.388 0.00179 **
## EC50:(Intercept) 5.242461 0.223590 23.447 1.961e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 2.520149 (4 degrees of freedom)Ahora se puede realizar un anova para comparar dos modelos entre si, en este caso al ser los mismos datos los modelos son iguales (no significativos):
anova(modelo1, modelo2, details=TRUE,test=NULL)##
## 1st model
## fct: L.4(fixed = c(NA, 0, 50, NA), names = c("slope", "lower limit",
## 2nd model
## fct: LL.4()## ANOVA table
##
## ModelDf RSS Df F value p value
## 2nd model 4 25.405
## 1st model 2 10.991 2 1.3114 0.4326Se puede aplicar el test de Neill para ver la falta de ajuste del modelo (tanto con concentraciones replicadas como no replicadas), se utiliza el argumento neill.test:
neill.test(modelo1)## Grouping used
##
## 1 2 3 4 5
## 1 1 1 2 1## Neill's lack-of-fit test
##
## F value p value
## 37.3867 0.1032neill.test(modelo2)## Grouping used
##
## 1 2 3
## 2 3 1## Neill's lack-of-fit test
##
## F value p value
## 1.3998 0.3220Para conocer los intervalos de confianza se utiliza confint:
confint(modelo1)## 2.5 % 97.5 %
## b:(Intercept) -6.400861 0.1547641
## c:(Intercept) -9.984676 9.8964821
## d:(Intercept) 27.437664 72.5396360
## e:(Intercept) 2.930043 7.0214648confint(modelo2)## 2.5 % 97.5 %
## slope:(Intercept) -0.8678246 -0.393727
## EC50:(Intercept) 4.6216749 5.863248TRABAJAMOS CON EL PAQUETE drc: EJEMPLO 2 DATOS DAPHNIDS
Cargamos lo datos del ejemplo, en este caso tenemos dos tiempos (24 y 48h), dosis, respuesta (inmóviles) y total de individuos. Se cuentan los mismos individuos en los dos periodos:
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 48hplot(daphnids$no~daphnids$dose, daphnids$time)
Realizamos el modelo dividiendo los afectados respecto al total frente a la concentración no/total~dose y en función del time. En este caso se opta por un modelo logístico de dos parámetros LL.2(), type = "binomial", con valor máximo de 1 y mínimo de o:
modelo3 <- drm(no/total~dose, time, weights = total,
data = daphnids, fct = LL.2(), type = "binomial")
modelFit(modelo3)## Goodness-of-fit test
##
## Df Chisq value p value
##
## DRC model 12 13.873 0.3089plot(modelo3, xlab="Log Dose",ylab="Response (0-1)")
summary(modelo3)##
## 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 ' ' 1Para conocer los intervalos de confianza se utiliza confint:
confint(modelo3)## 2.5 % 97.5 %
## b:24h -1.609647 -0.7380275
## b:48h -2.396943 -1.3024081
## e:24h 3062.857245 7205.2096381
## e:48h 1141.062395 1877.0683753REFERENCES
Brain P, Cousens R (1989). “An Equation to Describe Dose Responses where there is Stimulation of Growth at Low Dose” Weed Research, 29, 93-96.
van der Vaart AW (1998). Asymptotic Statistics. Cambridge University Press, Cambridge.
