[R-es] ajuste a curva e interpolación

[Argel Gastélum Arellánez]

El 26/01/11 16:58, David escribió:
> Hola Carlos, probaré con nls
>
> En efecto, se trata de un conjunto de datos de absorción obtenidos en el espectrofotómetro. En estos análisis, el kit trae unas muestras de las cuales se conoce la concentración del producto a investigar, lo que permite construir una curva estandar. Con esta curva, y a partir de los datos del espectrofotómetro, se trata de calcular la concentración del producto en muestras desconocidas.
>
> Matemáticamente no sé por qué el fabricante recomienda esta curva (sigmoidal de 4 parámetros, no lineal), pero tampoco es lo que realmente me interesa. Lo que necesito saber es si con R puedo realizar estos calculos
>
> 1.- construir la curva estandar con datos XY de absorción/concentración
> 2.- interpolar los datos de absorción(Y) para obtener los datos de concentración(X)
>
> No sé si era esto a lo que te referías. Muchas gracias.
>
> David
>
> --- El mié, 26/1/11, Carlos Ortega<coforfe en gmail.com>  escribió:
>
> De: Carlos Ortega<coforfe en gmail.com>
> Asunto: Re: [R-es] ajuste a curva e interpolación
> Para: "David"<dmenor1 en yahoo.es>
> CC: R-help-es en r-project.org
> Fecha: miércoles, 26 de enero, 2011 18:27
>
> Hola David,
> De tu correo, tan sólo he entendido que te haría falta una solución en R que te permitiera ajustar un conjunto de puntos. Sí eso existe, y como dices que utilizas una curva sigmoidad de 4 parámetros, es no-lineal, la solución en R se llama "nls()".
>
> A partir de este punto, el que sea una curva sigmoidal, el que se pueda ajustar con 4 o más/menos parámetros, dependerá de si puedes/quieres darnos más detalles: típicos resultados de un ELISA, incluso un conjunto de datos de muestra, o una referencia válida a la que podamos acceder para ver si esta solución es la válida.
>
> En general, si te quieres beneficiar de la ayuda que desde esta lista se te pueda prestar, el planteamiento que hagas de tu problema específico tienes que abstraerlo a ideas en lo posible de fácil digestión. De otro modo, tan sólo podrá ayudarte, en el mejor de los casos, un colega de campo/problema. 
>
> Saludos,Carlos Ortega
> www.qualityexcellence.eswww.datanalytics.com/blog
>
> 2011/1/26 David<dmenor1 en yahoo.es>
>
> Hola, quisiera saber qué paquete de R utilizar para calcular los datos de concentración de un ELISA, a partir de los datos del espectrofotómetro y los valores para construir una curva estandar
>
>
>
> el fabricante sugiere realizar el ajuste con una curva sigmoidal de 4 parámetros.
>
>
>
> Gracias

     Hola David, qué tal. De acuerdo con tu descripción, creo que lo que 
quieres es hacer algo como el código que pego abajo, seguramente se 
podrá adaptar a las necesidades de tu trabajo. Espero te sirva.

     Saludos.

--
     Argel.



#--------------------------------------------------------
# AJUSTE DEL MODELO LOGÍSTICO A CRECIMIENTO BACTERIANO.

# DATOS:

# Tiempo (h).

T <- c(0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 30)

# Densidad óptica.

DO <- c(0.0530, 0.0557, 0.0785, 0.1360, 0.1906, 0.2465, 0.2536, 0.4107, 
0.4606, 0.4498, 0.4294, 0.4119, 0.4285, 0.4093)

# Estimación inicial de parámetros (para no estar adivinando):

K.inicial <- 1.01*max(DO)

K.inicial

C.inicial <- (K.inicial/DO[1])-1

C.inicial

N.para.MU.inicial <- DO[length(DO)]

N.para.MU.inicial

T.para.MU.inicial <- T[length(T)]

T.para.MU.inicial

MU.inicial <- 
(log((-1+(K.inicial/N.para.MU.inicial))/C.inicial))/(-T.para.MU.inicial)

MU.inicial

# Crear el marco de datos:

Datos <- data.frame(x=T,  y=DO)

# Realizar el ajuste no lineal de la curva logística a los datos:

fit <-  nls(y  ~  K/(1+C*(exp(-MU*T))),  Datos, start = list (K = 
K.inicial, C = C.inicial, MU = MU.inicial))

fit

summary(fit)

# Suma de cuadrados de la regresión (o suma de cuadrados de los residuos):

SC.reg <- fit$m$deviance()

SC.reg

# Suma de cuadrados totales:

SC.tot <- sum((DO[1:length(DO)]-mean(DO))^2)

SC.tot

# Error estándar de la regresión (o desviación estándar de la regresión):

Sy.x <- summary(fit)$sigma

Sy.x

# Coeficiente de determinación (asumiéndolo sólo como una aproximación 
lineal):

R.cuadrada <- 1-(SC.reg/SC.tot)

R.cuadrada

# Estimar los intervalos de confianza al 95 % (usar la librería "nlrwr" 
(nonlinear regression with R)):

library(nlrwr)

IC95 <- confint2(fit)

IC95

# Graficar los resultados:

Coeficientes <- coef(fit)

Coeficientes

K <- (as.vector(Coeficientes))[1]

C <- (as.vector(Coeficientes))[2]

MU <- (as.vector(Coeficientes))[3]

K

C

MU

Logistica <- function(x){K/(1+C*(exp(-MU*x)))}
# Esta función también se puede usar para hacer predicciones puntuales, 
por ejemplo:

Logistica(x=10)

plot(x=T, y=DO)

points(x=c(0:30), Logistica(c(0:30)), type="l")

# ¿Bandas de confianza y de predicción? Esto al parecer aplica la 
aproximación lineal al problema de obtener las bandas de confianza

se.fit.confianza <- sqrt(apply(fit$m$gradient(),1,function(x) 
sum(vcov(fit)*outer(x,x))))

se.fit.prediccion <- sqrt((se.fit.confianza^2)+(summary(fit)$sigma)^2)

matplot(
     T,
     cbind(
         predict(fit)+outer(se.fit.confianza, qnorm(
         c(.5, .025,.975))),
         (predict(fit)+outer(se.fit.prediccion, qnorm(
         c(.5, .025,.975))))[,-1]
     ),
     type="l",
     lty=c(1,2,2,3,3),
     col=c("black", "blue", "blue", "brown", "brown"),
     xlab="Tiempo (h)",
     ylab="Densidad Óptica a 600 nm",
     main="Ajuste No Lineal del Modelo Logístico a Crecimiento Bacteriano"
)

points(T, DO)

legend(
     "bottomright",
     legend=c(
         "Datos",
         "Regresión No Lineal",
         paste(
             "DO =", round(K, 4), "/ ( 1 +", round(C, 4),
             "e^( -", round(MU, 4), "T ))"
         ),
         "Banda de Confianza",
         "Banda de Predicción"
     ),
     col=c("black", "black", NA, "blue", "brown"),
     lty=c(NA, 1, NA, 2, 3),
     pch=c(21, NA, NA, NA, NA)
)

# Para graficar sólo la banda de confianza:

matplot(
     T,
     predict(fit)+outer(se.fit.confianza,qnorm(c(.5, .025,.975))),
     type="l",
     lty=c(1,2,2),
     col=c("black", "blue", "blue"),
     xlab="Tiempo (h)",
     ylab="Densidad Óptica a 600 nm"
)

points(T, DO)

# Para graficar sólo la banda de predicción:

matplot(
     T,
     predict(fit)+outer(se.fit.prediccion,qnorm(c(.5, .025,.975))),
     type="l",
     lty=c(1,2,2),
     col=c("black", "brown", "brown"),
     xlab="Tiempo (h)",
     ylab="Densidad Óptica a 600 nm"
)

points(T, DO)

#--------------------------------------------------------