Argel Gastélum Arellánez
2012-Jan-23 21:55 UTC
[R-es] Ajuste no lineal con nls no converge
Hola compañeros de la lista.
Estoy tratando de hacer el ejercicio de realizar en R un ajuste no
lineal publicado en un artículo, originalmente hecho en SAS con el
método Gauss-Newton. Se trata de encontrar las constantes Km, Kic, Kiu y
k3 del modelo de inhibición mixta de Michaelis-Menten en su forma
integrada. El artículo proporciona los datos crudos utilizados para el
análisis (Los pongo al final de este mensaje). Mi duda es que no obtengo
el mismo resultado que el que indica dicho artículo, estoy usando
"nls"
con la opción algorithm = "default" y "port". El resultado
que obtengo
con "default" es el siguiente:
AJUSTE.default <- nls(t ~
((Km+Km*P0/Kic+Km*S0/Kic)*log((S0-(Pt-P0))/S0)+(-Km/Kic+1+P0/Kiu+S0/Kiu)*((S0-(Pt-P0))-S0)+(-1/(2*Kiu))*((S0-(Pt-P0))^2-S0^2))*(-1/(E*k3)),
DATOS, start = list (Km = 50, Kic = 6, Kiu = 6, k3 = 0.02), control =
nls.control(warnOnly = TRUE), algorithm = "default")
AJUSTE.default
Nonlinear regression model
model: t ~ ((Km + Km * P0/Kic + Km * S0/Kic) * log((S0 - (Pt - P0))/S0)
+ (-Km/Kic + 1 + P0/Kiu + S0/Kiu) * ((S0 - (Pt - P0)) - S0) + (-1/(2 *
Kiu)) * ((S0 - (Pt - P0))^2 - S0^2)) * (-1/(E * k3))
data: DATOS
Km Kic Kiu k3
4.373e+01 6.808e+00 -5.972e+14 1.825e-02
residual sum-of-squares: 159.1
Number of iterations till stop: 7
Achieved convergence tolerance: 1.397
Reason stopped: singular gradient
y al usar algorithm = "port" obtengo:
AJUSTE.port <- nls(t ~
((Km+Km*P0/Kic+Km*S0/Kic)*log((S0-(Pt-P0))/S0)+(-Km/Kic+1+P0/Kiu+S0/Kiu)*((S0-(Pt-P0))-S0)+(-1/(2*Kiu))*((S0-(Pt-P0))^2-S0^2))*(-1/(E*k3)),
DATOS, start = list (Km = 50, Kic = 6, Kiu = 6, k3 = 0.02), control =
nls.control(warnOnly = TRUE), algorithm = "port")
AJUSTE.port
Nonlinear regression model
model: t ~ ((Km + Km * P0/Kic + Km * S0/Kic) * log((S0 - (Pt - P0))/S0)
+ (-Km/Kic + 1 + P0/Kiu + S0/Kiu) * ((S0 - (Pt - P0)) - S0) + (-1/(2 *
Kiu)) * ((S0 - (Pt - P0))^2 - S0^2)) * (-1/(E * k3))
data: DATOS
Km Kic Kiu k3
4.426e+01 6.928e+00 1.137e+06 1.839e-02
residual sum-of-squares: 157.1
Algorithm "port", convergence message: singular convergence (7)
El resultado que indica el artículo es: Km = 46.1, Kic = 9.2, Kiu = 1.2
x 10^64, k3 = 0.0183. Como verán, las constantes Km, Kic y k3 obtenidas
en R son muy cercanas a las reportadas por SAS, mientras que el valor de
Kiu es diferente. Al final de cuentas la interpretación es la misma, la
constante Kiu resulta demasiado extrema (negativa con "default" y muy
grande con "port") para tener sentido, pero no entiendo por qué no
converge a un resultado en R cuando en SAS sí se obtiene (de acuerdo con
el artículo). ¿Qué argumento de nls se podría modificar para llegar a un
resultado similar al que se reporta con SAS?
Agradezco de antemano la ayuda brindada.
Saludos.
--
Argel.
Los datos usados son los siguientes ("Bezerra y Dias, 2007. Utilization
of Integrated Michaelis-Menten Equation to Determine Kinetic Constants.
Biochemistry and molecular biology education. Vol. 35, No. 2, pp.
145?150"):
# Datos del artículo: Bezerra y Días (2007).
# Se trata de evaluar los parámetros cinéticos de una fosfatasa
alcalina, sobre 4-nitrofenil fosfato como sustrato.
# Las reacciones se llevaron a cabo a 37 °C, en volúmenes de 2.75 mL,
con 9.5 µg de enzima y 8 concentraciones de sustrato, en amortiguador
Tris/HCl 0.1 M, pH 9.0.
# La cantidad de producto se determinó espectrofotométricamente a 405 nM.
# Unidades de las variables:
# Tiempo (t): s
# Concentración de producto (Pt): µM
# Concentración inicial de sustrato (S0): µM
# Concentración inicial de producto (P0): µM
# Enzima (E): µg/2.75 mL
#
"t" "Pt" "S0" "P0" "E"
0 0 25 0 9.5
10 0.566 25 0 9.5
20 1.154 25 0 9.5
30 1.704 25 0 9.5
40 2.281 25 0 9.5
50 2.709 25 0 9.5
60 3.124 25 0 9.5
0 0 34.7 0 9.5
10 0.634 34.7 0 9.5
20 1.296 34.7 0 9.5
30 2.01 34.7 0 9.5
40 2.639 34.7 0 9.5
50 3.257 34.7 0 9.5
60 3.844 34.7 0 9.5
0 0 45 0 9.5
10 0.966 45 0 9.5
20 1.749 45 0 9.5
30 2.558 45 0 9.5
40 3.275 45 0 9.5
50 4.014 45 0 9.5
60 4.611 45 0 9.5
0 0 90 0 9.5
10 0.952 90 0 9.5
20 1.953 90 0 9.5
30 2.939 90 0 9.5
40 3.92 90 0 9.5
50 4.825 90 0 9.5
60 5.766 90 0 9.5
0 0 200 0 9.5
10 1.261 200 0 9.5
20 2.608 200 0 9.5
30 3.906 200 0 9.5
40 5.135 200 0 9.5
50 6.423 200 0 9.5
60 7.612 200 0 9.5
0 0 500 0 9.5
10 1.363 500 0 9.5
20 2.928 500 0 9.5
30 4.37 500 0 9.5
40 6.015 500 0 9.5
50 7.78 500 0 9.5
60 9.363 500 0 9.5
0 0 1000 0 9.5
10 1.588 1000 0 9.5
20 3.125 1000 0 9.5
30 4.772 1000 0 9.5
40 6.463 1000 0 9.5
50 7.77 1000 0 9.5
60 9.907 1000 0 9.5
0 0 2000 0 9.5
10 2.044 2000 0 9.5
20 3.816 2000 0 9.5
30 5.346 2000 0 9.5
40 7.141 2000 0 9.5
50 8.803 2000 0 9.5
60 10.441 2000 0 9.5
Hola Argel,
Muchas gracias por los datos y el codigo ;-)
Despues de modificar las opciones nls.control(), obtengo:
# -----------
# datos
# -----------
DATOS <- structure(list(t = c(0L, 10L, 20L, 30L, 40L, 50L, 60L, 0L, 10L,
20L, 30L, 40L, 50L, 60L, 0L, 10L, 20L, 30L, 40L, 50L, 60L, 0L,
10L, 20L, 30L, 40L, 50L, 60L, 0L, 10L, 20L, 30L, 40L, 50L, 60L,
0L, 10L, 20L, 30L, 40L, 50L, 60L, 0L, 10L, 20L, 30L, 40L, 50L,
60L, 0L, 10L, 20L, 30L, 40L, 50L, 60L), Pt = c(0, 0.566, 1.154,
1.704, 2.281, 2.709, 3.124, 0, 0.634, 1.296, 2.01, 2.639, 3.257,
3.844, 0, 0.966, 1.749, 2.558, 3.275, 4.014, 4.611, 0, 0.952,
1.953, 2.939, 3.92, 4.825, 5.766, 0, 1.261, 2.608, 3.906, 5.135,
6.423, 7.612, 0, 1.363, 2.928, 4.37, 6.015, 7.78, 9.363, 0, 1.588,
3.125, 4.772, 6.463, 7.77, 9.907, 0, 2.044, 3.816, 5.346, 7.141,
8.803, 10.441), S0 = c(25, 25, 25, 25, 25, 25, 25, 34.7, 34.7,
34.7, 34.7, 34.7, 34.7, 34.7, 45, 45, 45, 45, 45, 45, 45, 90,
90, 90, 90, 90, 90, 90, 200, 200, 200, 200, 200, 200, 200, 500,
500, 500, 500, 500, 500, 500, 1000, 1000, 1000, 1000, 1000, 1000,
1000, 2000, 2000, 2000, 2000, 2000, 2000, 2000), P0 = c(0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L), E = c(9.5, 9.5, 9.5, 9.5, 9.5, 9.5,
9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5,
9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5,
9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5,
9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5)), .Names =
c("t",
"Pt", "S0", "P0", "E"), class =
"data.frame", row.names = c(NA,
-56L))
# -----------
# modelos
# -----------
# default
AJUSTE.default <- nls(t ~
((Km+Km*P0/Kic+Km*S0/Kic)*log((S0-(Pt-P0))/S0)+(-Km/Kic+1+P0/Kiu+S0/Kiu)*((S0-(Pt-P0))-S0)+(-1/(2*Kiu))*((S0-(Pt-P0))^2-S0^2))*(-1/(E*k3)),
DATOS, start = list (Km = 50, Kic = 6, Kiu = 6, k3 = 0.02), control
nls.control(warnOnly = TRUE, maxiter = 50000, tol = 1e-05), algorithm
"default")
#Warning message:
#In nls(t ~ ((Km + Km * P0/Kic + Km * S0/Kic) * log((S0 - (Pt - P0))/S0) +
:
# singular gradient
summary(AJUSTE.default)
#Error in chol2inv(object$m$Rmat()) :
# element (4, 4) is zero, so the inverse cannot be computed
# port
AJUSTE.port <- nls(t ~
((Km+Km*P0/Kic+Km*S0/Kic)*log((S0-(Pt-P0))/S0)+(-Km/Kic+1+P0/Kiu+S0/Kiu)*((S0-(Pt-P0))-S0)+(-1/(2*Kiu))*((S0-(Pt-P0))^2-S0^2))*(-1/(E*k3)),
DATOS, start = list (Km = 50, Kic = 6, Kiu = 6, k3 = 0.02), control
nls.control(warnOnly = TRUE, maxiter = 50000, tol = 1e-05), algorithm
"port")
#Warning message:
#In nls(t ~ ((Km + Km * P0/Kic + Km * S0/Kic) * log((S0 - (Pt - P0))/S0) +
:
# Convergence failure: singular convergence (7)
summary(AJUSTE.port)
#Formula: t ~ ((Km + Km * P0/Kic + Km * S0/Kic) * log((S0 - (Pt - P0))/S0)
+
# (-Km/Kic + 1 + P0/Kiu + S0/Kiu) * ((S0 - (Pt - P0)) - S0) +
# (-1/(2 * Kiu)) * ((S0 - (Pt - P0))^2 - S0^2)) * (-1/(E *
# k3))
#Parameters:
# Estimate Std. Error t value Pr(>|t|)
#Km 4.426e+01 3.704e+00 11.950 <2e-16 ***
#Kic 6.928e+00 2.117e+00 3.273 0.0019 **
#Kiu 1.137e+06 1.220e+10 0.000 0.9999
#k3 1.839e-02 7.075e-04 25.991 <2e-16 ***
#---
#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Residual standard error: 1.738 on 52 degrees of freedom
#
#Algorithm "port", convergence message: singular convergence (7)
Como veras, aunque el modelo "default" definitivamente no converge,
aunque
el modelo "port" si. Lo que no entiendo es por que la diferencia
entre los
resultados que envias y lo que obtuve (ejecute el codigo tal y como lo
enviaste y no obtuve lo mismo).
Te sugeriria revisar un poco las expresiones a la derecha de "t ~ ";
seguramente puede haber un termino que aparece varias veces.
Esta es mi sessionInfo()
> sessionInfo()
R version 2.14.0 Patched (2011-11-12 r57642)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
loaded via a namespace (and not attached):
[1] tools_2.14.0
Saludos,
Jorge.-
2012/1/23 Argel Gastélum Arellánez <>
> Hola compañeros de la lista.
>
> Estoy tratando de hacer el ejercicio de realizar en R un ajuste no lineal
> publicado en un artículo, originalmente hecho en SAS con el método
> Gauss-Newton. Se trata de encontrar las constantes Km, Kic, Kiu y k3 del
> modelo de inhibición mixta de Michaelis-Menten en su forma integrada. El
> artículo proporciona los datos crudos utilizados para el análisis (Los
> pongo al final de este mensaje). Mi duda es que no obtengo el mismo
> resultado que el que indica dicho artículo, estoy usando "nls"
con la
> opción algorithm = "default" y "port". El resultado que
obtengo con
> "default" es el siguiente:
>
> AJUSTE.default <- nls(t ~ ((Km+Km*P0/Kic+Km*S0/Kic)*log(**
> (S0-(Pt-P0))/S0)+(-Km/Kic+1+**P0/Kiu+S0/Kiu)*((S0-(Pt-P0))-**
> S0)+(-1/(2*Kiu))*((S0-(Pt-P0))**^2-S0^2))*(-1/(E*k3)), DATOS, start >
list (Km = 50, Kic = 6, Kiu = 6, k3 = 0.02), control = nls.control(warnOnly
> = TRUE), algorithm = "default")
>
> AJUSTE.default
>
> Nonlinear regression model
> model: t ~ ((Km + Km * P0/Kic + Km * S0/Kic) * log((S0 - (Pt - P0))/S0) +
> (-Km/Kic + 1 + P0/Kiu + S0/Kiu) * ((S0 - (Pt - P0)) - S0) + (-1/(2 * Kiu))
> * ((S0 - (Pt - P0))^2 - S0^2)) * (-1/(E * k3))
> data: DATOS
> Km Kic Kiu k3
> 4.373e+01 6.808e+00 -5.972e+14 1.825e-02
> residual sum-of-squares: 159.1
>
> Number of iterations till stop: 7
> Achieved convergence tolerance: 1.397
> Reason stopped: singular gradient
>
>
>
> y al usar algorithm = "port" obtengo:
>
> AJUSTE.port <- nls(t ~ ((Km+Km*P0/Kic+Km*S0/Kic)*log(**
> (S0-(Pt-P0))/S0)+(-Km/Kic+1+**P0/Kiu+S0/Kiu)*((S0-(Pt-P0))-**
> S0)+(-1/(2*Kiu))*((S0-(Pt-P0))**^2-S0^2))*(-1/(E*k3)), DATOS, start >
list (Km = 50, Kic = 6, Kiu = 6, k3 = 0.02), control = nls.control(warnOnly
> = TRUE), algorithm = "port")
>
> AJUSTE.port
>
> Nonlinear regression model
> model: t ~ ((Km + Km * P0/Kic + Km * S0/Kic) * log((S0 - (Pt - P0))/S0) +
> (-Km/Kic + 1 + P0/Kiu + S0/Kiu) * ((S0 - (Pt - P0)) - S0) + (-1/(2 * Kiu))
> * ((S0 - (Pt - P0))^2 - S0^2)) * (-1/(E * k3))
> data: DATOS
> Km Kic Kiu k3
> 4.426e+01 6.928e+00 1.137e+06 1.839e-02
> residual sum-of-squares: 157.1
>
> Algorithm "port", convergence message: singular convergence (7)
>
>
>
> El resultado que indica el artículo es: Km = 46.1, Kic = 9.2, Kiu = 1.2 x
> 10^64, k3 = 0.0183. Como verán, las constantes Km, Kic y k3 obtenidas en R
> son muy cercanas a las reportadas por SAS, mientras que el valor de Kiu es
> diferente. Al final de cuentas la interpretación es la misma, la constante
> Kiu resulta demasiado extrema (negativa con "default" y muy
grande con
> "port") para tener sentido, pero no entiendo por qué no converge
a un
> resultado en R cuando en SAS sí se obtiene (de acuerdo con el artículo).
> ¿Qué argumento de nls se podría modificar para llegar a un resultado
> similar al que se reporta con SAS?
>
> Agradezco de antemano la ayuda brindada.
>
> Saludos.
>
> --
> Argel.
>
>
>
>
>
>
>
> Los datos usados son los siguientes ("Bezerra y Dias, 2007.
Utilization of
> Integrated Michaelis-Menten Equation to Determine Kinetic Constants.
> Biochemistry and molecular biology education. Vol. 35, No. 2, pp.
145–150"):
>
> # Datos del artículo: Bezerra y Días (2007).
> # Se trata de evaluar los parámetros cinéticos de una fosfatasa alcalina,
> sobre 4-nitrofenil fosfato como sustrato.
> # Las reacciones se llevaron a cabo a 37 °C, en volúmenes de 2.75 mL, con
> 9.5 µg de enzima y 8 concentraciones de sustrato, en amortiguador Tris/HCl
> 0.1 M, pH 9.0.
> # La cantidad de producto se determinó espectrofotométricamente a 405 nM.
> # Unidades de las variables:
> # Tiempo (t): s
> # Concentración de producto (Pt): µM
> # Concentración inicial de sustrato (S0): µM
> # Concentración inicial de producto (P0): µM
> # Enzima (E): µg/2.75 mL
> #
> "t" "Pt" "S0" "P0" "E"
> 0 0 25 0 9.5
> 10 0.566 25 0 9.5
> 20 1.154 25 0 9.5
> 30 1.704 25 0 9.5
> 40 2.281 25 0 9.5
> 50 2.709 25 0 9.5
> 60 3.124 25 0 9.5
> 0 0 34.7 0 9.5
> 10 0.634 34.7 0 9.5
> 20 1.296 34.7 0 9.5
> 30 2.01 34.7 0 9.5
> 40 2.639 34.7 0 9.5
> 50 3.257 34.7 0 9.5
> 60 3.844 34.7 0 9.5
> 0 0 45 0 9.5
> 10 0.966 45 0 9.5
> 20 1.749 45 0 9.5
> 30 2.558 45 0 9.5
> 40 3.275 45 0 9.5
> 50 4.014 45 0 9.5
> 60 4.611 45 0 9.5
> 0 0 90 0 9.5
> 10 0.952 90 0 9.5
> 20 1.953 90 0 9.5
> 30 2.939 90 0 9.5
> 40 3.92 90 0 9.5
> 50 4.825 90 0 9.5
> 60 5.766 90 0 9.5
> 0 0 200 0 9.5
> 10 1.261 200 0 9.5
> 20 2.608 200 0 9.5
> 30 3.906 200 0 9.5
> 40 5.135 200 0 9.5
> 50 6.423 200 0 9.5
> 60 7.612 200 0 9.5
> 0 0 500 0 9.5
> 10 1.363 500 0 9.5
> 20 2.928 500 0 9.5
> 30 4.37 500 0 9.5
> 40 6.015 500 0 9.5
> 50 7.78 500 0 9.5
> 60 9.363 500 0 9.5
> 0 0 1000 0 9.5
> 10 1.588 1000 0 9.5
> 20 3.125 1000 0 9.5
> 30 4.772 1000 0 9.5
> 40 6.463 1000 0 9.5
> 50 7.77 1000 0 9.5
> 60 9.907 1000 0 9.5
> 0 0 2000 0 9.5
> 10 2.044 2000 0 9.5
> 20 3.816 2000 0 9.5
> 30 5.346 2000 0 9.5
> 40 7.141 2000 0 9.5
> 50 8.803 2000 0 9.5
> 60 10.441 2000 0 9.5
>
> ______________________________**_________________
> R-help-es mailing list
> R-help-es@r-project.org
>
https://stat.ethz.ch/mailman/**listinfo/r-help-es<https://stat.ethz.ch/mailman/listinfo/r-help-es>
>
[[alternative HTML version deleted]]
Argel Gastélum Arellánez
2012-Jan-24 00:26 UTC
[R-es] Ajuste no lineal con nls no converge
Hola Olivier, muchas gracias por tus comentarios. Pongo algunas notas más abajo. El 23/01/12 16:43, Olivier Nuñez escribió:> Argel, > > Kiu = 1.2 x 10^64 es un valor bastante grande para una constante de > inhibición!Efectivamente, Kiu = 1.2 x 10^64 es demasiado grande para una constante de inhibición. El problema no es de interpretación del resultado, es más del por qué "nls" no converge hacia los mismos valores que ellos reportan para el mismo set de datos crudos, es más, sólo no converge. Es precisamente por lo que los autores del artículo ajustan a los mismos datos el modelo que considera sólo la inhibición competitiva.> Y sospecho que el algoritmo de estos autores tampoco converge con > estos datos.Los autores mencionan que en SAS encuentran convergencia con este modelo mixto hacia los valores que mencioné anteriormente. Yo realmente desconozco SAS y no tengo acceso a él como para comprobarlo. Pensé que con R y "nls" debería llegar a converger también puesto que con algorithm = "default" se trataría también del método Gauss-Newton que ellos usaron en SAS.> De hecho si pones trace=TRUE, verás que este parámetro se va al > infinito y el problema parece más bien de identificación del modelo: >Ya lo había hecho antes de consultar a la lista y había notado que Kiu se va incrementando hacia valores muy grandes y sin mucho sentido (el comportamiento esperado, ya que el modelo mixto no se ajusta bien a los datos, sino que se ajusta mejor a un modelo con sólo inhibición competitiva), pero como los autores mencionan convergencia creí que en R también pasaría así. Lo que me interesa por el momento es más bien poder obtener los mismos valores ajustados que ellos reportan para los parámetros, sabiendo de antemano que no son los que mejor se ajustan (explico esto con un poco más de detalle más abajo).> > AJUSTE.default <- nls(t ~ > ((Km+Km*P0/Kic+Km*S0/Kic)*log((S0-(Pt-P0))/S0)+(-Km/Kic+1+P0/Kiu+S0/Kiu)*((S0-(Pt-P0))-S0)+(-1/(2*Kiu))*((S0-(Pt-P0))^2-S0^2))*(-1/(E*k3)), > DATOS, start = list (Km = 50, Kic = 6, Kiu = 6, k3 = 0.02), algorithm > = "default",trace=TRUE) > 9689.914 : 50.00 6.00 6.00 0.02 > 6670.874 : 34.93095937 5.79629581 14.16884101 0.01610926 > 1896.282 : 37.29205443 5.77977412 32.30063332 0.01671887 > 840.2903 : 37.3846129 5.7802131 88.8218065 0.0167428 > 536.3691 : 37.38476020 5.78018175 361.40233750 0.01674285 > 463.0871 : 3.738474e+01 5.780188e+00 3.765039e+03 1.674285e-02 > 455.831 : 3.738453e+01 5.780068e+00 3.377163e+05 1.674283e-02 > 455.8272 : 3.738363e+01 5.780040e+00 3.325798e+08 1.674261e-02 > Error en nls(t ~ ((Km + Km * P0/Kic + Km * S0/Kic) * log((S0 - (Pt - > P0))/S0) + : > singular gradient > > Un saludo. OlivierComento un poco más sobre el artículo. El objetivo en general es ajustar los siguientes modelos de Michaelis-Menten: Sin inhibición, con inhibición competitiva, con inhibición no competitiva, con inhibición acompetitiva y con inhibición mixta (todos en sus formas integradas), y de acuerdo con los resultados de los valores de los parámetros, del número de parámetros y de las sumas de cuadrados decidir cuál es el que mejor se ajusta a los datos. En R, "nls" ha reproducido con exactitud los parámetros estimados por los autores para los modelos sin inhibición y con inhibición competitiva, pero no para los que consideran inhibición acompetitiva e inhibición mixta (donde "nls" reporta no convergencia y los autores reportan convergencia con valores de Kiu muy grandes, del orden de 10^38 a 10^64). Los autores mencionan que ellos tampoco encontraron convergencia con el modelo de inhibición no competitiva. Sin embargo, al usar en R y "nls" las opciones warnOnly = TRUE, "trrosace" y "port" observo que aunque no reporta convergencia los valores de los parámetros ajustados son muy similares a lo reportado por los autores, excepto para las constantes de inhibición que resultan demasiado grandes, donde "nls" reporta no convergencia y valores finales del orden de 10^5 para Kiu. Espero haber aclarado un poco la falta de información preliminar. De nuevo muchas gracias por tus comentarios. Saludos. -- Argel.
Argel Gastélum Arellánez
2012-Jan-24 00:52 UTC
[R-es] Ajuste no lineal con nls no converge
Hola Jorge, muchas gracias por tus comentarios. Pongo algunas notas más abajo. El 23/01/12 17:19, Jorge I Velez escribió:> Hola Argel, > > Muchas gracias por los datos y el codigo ;-) > > Despues de modificar las opciones nls.control(), obtengo: > > # ----------- > # datos > # ----------- > DATOS <- structure(list(t = c(0L, 10L, 20L, 30L, 40L, 50L, 60L, 0L, 10L, > 20L, 30L, 40L, 50L, 60L, 0L, 10L, 20L, 30L, 40L, 50L, 60L, 0L, > 10L, 20L, 30L, 40L, 50L, 60L, 0L, 10L, 20L, 30L, 40L, 50L, 60L, > 0L, 10L, 20L, 30L, 40L, 50L, 60L, 0L, 10L, 20L, 30L, 40L, 50L, > 60L, 0L, 10L, 20L, 30L, 40L, 50L, 60L), Pt = c(0, 0.566, 1.154, > 1.704, 2.281, 2.709, 3.124, 0, 0.634, 1.296, 2.01, 2.639, 3.257, > 3.844, 0, 0.966, 1.749, 2.558, 3.275, 4.014, 4.611, 0, 0.952, > 1.953, 2.939, 3.92, 4.825, 5.766, 0, 1.261, 2.608, 3.906, 5.135, > 6.423, 7.612, 0, 1.363, 2.928, 4.37, 6.015, 7.78, 9.363, 0, 1.588, > 3.125, 4.772, 6.463, 7.77, 9.907, 0, 2.044, 3.816, 5.346, 7.141, > 8.803, 10.441), S0 = c(25, 25, 25, 25, 25, 25, 25, 34.7, 34.7, > 34.7, 34.7, 34.7, 34.7, 34.7, 45, 45, 45, 45, 45, 45, 45, 90, > 90, 90, 90, 90, 90, 90, 200, 200, 200, 200, 200, 200, 200, 500, > 500, 500, 500, 500, 500, 500, 1000, 1000, 1000, 1000, 1000, 1000, > 1000, 2000, 2000, 2000, 2000, 2000, 2000, 2000), P0 = c(0L, 0L, > 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, > 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, > 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, > 0L, 0L, 0L, 0L, 0L, 0L), E = c(9.5, 9.5, 9.5, 9.5, 9.5, 9.5, > 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, > 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, > 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, > 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5)), .Names = c("t", > "Pt", "S0", "P0", "E"), class = "data.frame", row.names = c(NA, > -56L)) > > # ----------- > # modelos > # ----------- > # default > AJUSTE.default <- nls(t ~ > ((Km+Km*P0/Kic+Km*S0/Kic)*log((S0-(Pt-P0))/S0)+(-Km/Kic+1+P0/Kiu+S0/Kiu)*((S0-(Pt-P0))-S0)+(-1/(2*Kiu))*((S0-(Pt-P0))^2-S0^2))*(-1/(E*k3)), > DATOS, start = list (Km = 50, Kic = 6, Kiu = 6, k3 = 0.02), control = > nls.control(warnOnly = TRUE, maxiter = 50000, tol = 1e-05), algorithm > = "default") > #Warning message: > #In nls(t ~ ((Km + Km * P0/Kic + Km * S0/Kic) * log((S0 - (Pt - > P0))/S0) + : > # singular gradient > > summary(AJUSTE.default) > #Error in chol2inv(object$m$Rmat()) : > # element (4, 4) is zero, so the inverse cannot be computed > > # port > AJUSTE.port <- nls(t ~ > ((Km+Km*P0/Kic+Km*S0/Kic)*log((S0-(Pt-P0))/S0)+(-Km/Kic+1+P0/Kiu+S0/Kiu)*((S0-(Pt-P0))-S0)+(-1/(2*Kiu))*((S0-(Pt-P0))^2-S0^2))*(-1/(E*k3)), > DATOS, start = list (Km = 50, Kic = 6, Kiu = 6, k3 = 0.02), control = > nls.control(warnOnly = TRUE, maxiter = 50000, tol = 1e-05), algorithm > = "port") > #Warning message: > #In nls(t ~ ((Km + Km * P0/Kic + Km * S0/Kic) * log((S0 - (Pt - > P0))/S0) + : > # Convergence failure: singular convergence (7) > summary(AJUSTE.port) > > #Formula: t ~ ((Km + Km * P0/Kic + Km * S0/Kic) * log((S0 - (Pt - > P0))/S0) + > # (-Km/Kic + 1 + P0/Kiu + S0/Kiu) * ((S0 - (Pt - P0)) - S0) + > # (-1/(2 * Kiu)) * ((S0 - (Pt - P0))^2 - S0^2)) * (-1/(E * > # k3)) > > #Parameters: > # Estimate Std. Error t value Pr(>|t|) > #Km 4.426e+01 3.704e+00 11.950 <2e-16 *** > #Kic 6.928e+00 2.117e+00 3.273 0.0019 ** > #Kiu 1.137e+06 1.220e+10 0.000 0.9999 > #k3 1.839e-02 7.075e-04 25.991 <2e-16 *** > #--- > #Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > # > #Residual standard error: 1.738 on 52 degrees of freedom > # > #Algorithm "port", convergence message: singular convergence (7) > > Como veras, aunque el modelo "default" definitivamente no converge, > aunque el modelo "port" si. Lo que no entiendo es por que la > diferencia entre los resultados que envias y lo que obtuve (ejecute el > codigo tal y como lo enviaste y no obtuve lo mismo).Ejecuté el código que me envías modificado y obtengo exactamente el mismo resultado que encuentras. Sin embargo, creo que tampoco el algoritmo "port" converge, a menos que yo esté interpretando mal el mensaje último: " Convergence failure: singular convergence (7) "> Te sugeriria revisar un poco las expresiones a la derecha de "t ~ "; > seguramente puede haber un termino que aparece varias veces. >Los modelos están capturados exactamente como los reportan los autores del artículo en SAS. Mi intención por ahora es tratar de obtener el mismo resultado que ellos con los mismos modelos pero en R y con nls. Comenté un poco más sobre este artículo en mi respuesta a Olivier.> Saludos, > Jorge.-De nuevo muchas gracias por tu ayuda. Saludos. -- Argel. [[alternative HTML version deleted]]
Hola,
He probado con otra variante, visto que utilizando nls y los algoritmos de
aproximación (default y port) están dando resultados todavía lejanos a los
del artículo, además de los problemas de convergencia o de singularidad.
La aproximación que he utilizado es la siguiente:
1. He utilizado varios algoritmos de los que proporciona el paquete
nls2. En este paquete la aproximación es más de "fuerza bruta".
2. De los varios algortimos que ofrece el paquete: fuerza-bruta,
random-search, plinear. Los resultados con fuerza-bruta y random-search son
equivalentes. He desechado el de plinear por no ajustar correctamente.
3. Para el algoritmo "brute-force", se puede definir un grid de
evaluación que permite al algoritmo evaluar diferentes opciones, es algo
más lento pero no excesivo.
4. En este grid, he jugado con varios extremos para el valor de
"Kiu"
que es el que se va de madre. He probado con valores desde 0 a 1000, hasta
quedarme con el último caso en el que lo he iterado entre 10^11 y 10^20.
5. Además de la librería "nls2", he utilizado también la librería
"nlstools" que ofrece una función "overview()" que
proporciona el intervalo
de confianza de cada unos de los parámetros del modelo (intervalos de
confianza del 2.5% y del 97.5%).
El código de todo esto es el siguiente:
################
#### Fuerza Bruta
# Resultados paper SAS: Km = 46.1, Kic = 9.2, Kiu = 1.2 x 10^64, k3 = 0.0183
library(nls2)
library(nlstools)
mi.me<-t ~
((Km+Km*P0/Kic+Km*S0/Kic)*log((S0-(Pt-P0))/S0)+(-Km/Kic+1+P0/Kiu+S0/Kiu)*((S0-(Pt-P0))-S0)+(-1/(2*Kiu))*((S0-(Pt-P0))^2-S0^2))*(-1/(E*k3))
gr.id<-expand.grid(
Km=seq(30,80,len=10),
Kic=seq(6,12,len=10),
Kiu=seq(10^11,10^20,len=10),
k3=seq(0.01,0.04,len=10)
)
### Brute-Force
aj.bf <- nls2(
mi.me,
data=d.a,
start = gr.id,
algorithm = "brute-force"
)
summary(aj.bf)
overview(aj.bf)
################
Y el resultado de overview(aj.bf) es el siguiente (tan sólo incluyo el
intervalo de confianza):
------
Asymptotic confidence interval:
2.5% 97.5%
Km 4.451537e+01 5.992907e+01
Kic 1.481808e+00 1.185153e+01
Kiu -1.039715e+17 1.039717e+17
k3 1.921825e-02 2.078175e-02
------
Como detalles y preguntas adicionales:
- No he probado con varlores superiones de Kiu. Bueno, sí que he probado
con valores de 10^60 y 10^61 pero el modelo no converge. Supongo que ya
aparecen problemas de precisión...
- En el paquete nlstools hay ejemplos de ajuste de curvas de
Michaelis-Menten dentro del punto de su help "compet_mich". Puede
resultar
ilustrativo ver el tipo de ajustes que se utilizan.
- Tambén en el paquete nlstools, hay uhtilidades adicinales
(nlsConfRegions o nlsResiduals) que añaden detalles al análisis.
- De los datos que has aportado, P0 siempre es cero (0). Aunque son
datos empíricos, ¿ser siempre cero no podría ayudar a simplificar la
fórmula del modelo?. Vaya eliminar P0 de la fórmula y así quitar un
parámetro de ajuste.
- Otro tema que me genera ciertas dudas es el de que los datos
experimentales de "t" tienen la pinta de diferentes experimentos
donde "t"
siempre se mueve en un rango (0 a 60). Y esto se repite 8 veces. ¿No puede
ayudar este patrón a ajustar un primer modelo con cada una de las celdas?...
Por último otra duda que puede ayudar. En el artículo donde comentan esta
aproximación vía código SAS, ¿se incluyen los intervalos de confianza para
cada parámetro?. A lo mejor los resultados que ya obtenemos están en esos
intervalos (lo digo más que nada por Kiu).
Saludos,
Carlos Ortega
www.qualityexcellence.es
El 24 de enero de 2012 01:52, Argel Gastélum Arellánez <
argel.gastelum@gmail.com> escribió:
> Hola Jorge, muchas gracias por tus comentarios. Pongo algunas notas
> más abajo.
>
> El 23/01/12 17:19, Jorge I Velez escribió:
> > Hola Argel,
> >
> > Muchas gracias por los datos y el codigo ;-)
> >
> > Despues de modificar las opciones nls.control(), obtengo:
> >
> > # -----------
> > # datos
> > # -----------
> > DATOS <- structure(list(t = c(0L, 10L, 20L, 30L, 40L, 50L, 60L, 0L,
10L,
> > 20L, 30L, 40L, 50L, 60L, 0L, 10L, 20L, 30L, 40L, 50L, 60L, 0L,
> > 10L, 20L, 30L, 40L, 50L, 60L, 0L, 10L, 20L, 30L, 40L, 50L, 60L,
> > 0L, 10L, 20L, 30L, 40L, 50L, 60L, 0L, 10L, 20L, 30L, 40L, 50L,
> > 60L, 0L, 10L, 20L, 30L, 40L, 50L, 60L), Pt = c(0, 0.566, 1.154,
> > 1.704, 2.281, 2.709, 3.124, 0, 0.634, 1.296, 2.01, 2.639, 3.257,
> > 3.844, 0, 0.966, 1.749, 2.558, 3.275, 4.014, 4.611, 0, 0.952,
> > 1.953, 2.939, 3.92, 4.825, 5.766, 0, 1.261, 2.608, 3.906, 5.135,
> > 6.423, 7.612, 0, 1.363, 2.928, 4.37, 6.015, 7.78, 9.363, 0, 1.588,
> > 3.125, 4.772, 6.463, 7.77, 9.907, 0, 2.044, 3.816, 5.346, 7.141,
> > 8.803, 10.441), S0 = c(25, 25, 25, 25, 25, 25, 25, 34.7, 34.7,
> > 34.7, 34.7, 34.7, 34.7, 34.7, 45, 45, 45, 45, 45, 45, 45, 90,
> > 90, 90, 90, 90, 90, 90, 200, 200, 200, 200, 200, 200, 200, 500,
> > 500, 500, 500, 500, 500, 500, 1000, 1000, 1000, 1000, 1000, 1000,
> > 1000, 2000, 2000, 2000, 2000, 2000, 2000, 2000), P0 = c(0L, 0L,
> > 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
> > 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
> > 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
> > 0L, 0L, 0L, 0L, 0L, 0L), E = c(9.5, 9.5, 9.5, 9.5, 9.5, 9.5,
> > 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5,
> > 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5,
> > 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5,
> > 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5)), .Names =
c("t",
> > "Pt", "S0", "P0", "E"), class
= "data.frame", row.names = c(NA,
> > -56L))
> >
> > # -----------
> > # modelos
> > # -----------
> > # default
> > AJUSTE.default <- nls(t ~
> >
>
((Km+Km*P0/Kic+Km*S0/Kic)*log((S0-(Pt-P0))/S0)+(-Km/Kic+1+P0/Kiu+S0/Kiu)*((S0-(Pt-P0))-S0)+(-1/(2*Kiu))*((S0-(Pt-P0))^2-S0^2))*(-1/(E*k3)),
> > DATOS, start = list (Km = 50, Kic = 6, Kiu = 6, k3 = 0.02), control
> > nls.control(warnOnly = TRUE, maxiter = 50000, tol = 1e-05), algorithm
> > = "default")
> > #Warning message:
> > #In nls(t ~ ((Km + Km * P0/Kic + Km * S0/Kic) * log((S0 - (Pt -
> > P0))/S0) + :
> > # singular gradient
> >
> > summary(AJUSTE.default)
> > #Error in chol2inv(object$m$Rmat()) :
> > # element (4, 4) is zero, so the inverse cannot be computed
> >
> > # port
> > AJUSTE.port <- nls(t ~
> >
>
((Km+Km*P0/Kic+Km*S0/Kic)*log((S0-(Pt-P0))/S0)+(-Km/Kic+1+P0/Kiu+S0/Kiu)*((S0-(Pt-P0))-S0)+(-1/(2*Kiu))*((S0-(Pt-P0))^2-S0^2))*(-1/(E*k3)),
> > DATOS, start = list (Km = 50, Kic = 6, Kiu = 6, k3 = 0.02), control
> > nls.control(warnOnly = TRUE, maxiter = 50000, tol = 1e-05), algorithm
> > = "port")
> > #Warning message:
> > #In nls(t ~ ((Km + Km * P0/Kic + Km * S0/Kic) * log((S0 - (Pt -
> > P0))/S0) + :
> > # Convergence failure: singular convergence (7)
> > summary(AJUSTE.port)
> >
> > #Formula: t ~ ((Km + Km * P0/Kic + Km * S0/Kic) * log((S0 - (Pt -
> > P0))/S0) +
> > # (-Km/Kic + 1 + P0/Kiu + S0/Kiu) * ((S0 - (Pt - P0)) - S0) +
> > # (-1/(2 * Kiu)) * ((S0 - (Pt - P0))^2 - S0^2)) * (-1/(E *
> > # k3))
> >
> > #Parameters:
> > # Estimate Std. Error t value Pr(>|t|)
> > #Km 4.426e+01 3.704e+00 11.950 <2e-16 ***
> > #Kic 6.928e+00 2.117e+00 3.273 0.0019 **
> > #Kiu 1.137e+06 1.220e+10 0.000 0.9999
> > #k3 1.839e-02 7.075e-04 25.991 <2e-16 ***
> > #---
> > #Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> > #
> > #Residual standard error: 1.738 on 52 degrees of freedom
> > #
> > #Algorithm "port", convergence message: singular convergence
(7)
> >
> > Como veras, aunque el modelo "default" definitivamente no
converge,
> > aunque el modelo "port" si. Lo que no entiendo es por que
la
> > diferencia entre los resultados que envias y lo que obtuve (ejecute el
> > codigo tal y como lo enviaste y no obtuve lo mismo).
>
> Ejecuté el código que me envías modificado y obtengo exactamente el
> mismo resultado que encuentras. Sin embargo, creo que tampoco el
> algoritmo "port" converge, a menos que yo esté interpretando mal
el
> mensaje último: " Convergence failure: singular convergence (7) "
>
> > Te sugeriria revisar un poco las expresiones a la derecha de "t ~
";
> > seguramente puede haber un termino que aparece varias veces.
> >
>
> Los modelos están capturados exactamente como los reportan los
> autores del artículo en SAS. Mi intención por ahora es tratar de obtener
> el mismo resultado que ellos con los mismos modelos pero en R y con nls.
> Comenté un poco más sobre este artículo en mi respuesta a Olivier.
>
> > Saludos,
> > Jorge.-
>
> De nuevo muchas gracias por tu ayuda.
>
> Saludos.
>
> --
> Argel.
>
>
> [[alternative HTML version deleted]]
>
>
> _______________________________________________
> R-help-es mailing list
> R-help-es@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es
>
>
[[alternative HTML version deleted]]
Argel Gastélum Arellánez
2012-Jan-24 19:41 UTC
[R-es] Ajuste no lineal con nls no converge
Hola Carlos, muchas gracias por tu ayuda. Me parece muy interesante la forma de aplicar nls2. Pongo algunas notas abajo. El 24/01/12 09:23, Carlos Ortega escribió:> Hola, > > He probado con otra variante, visto que utilizando nls y los > algoritmos de aproximación (default y port) están dando resultados > todavía lejanos a los del artículo, además de los problemas de > convergencia o de singularidad. > antEn realidad el parámetro que sale bastante alejado de los resultados del artículo es Kiu, los otros dan resultados muy similares. Puesto que uno de los objetivos es comprobar que el modelo mixto no es el que mejor se ajusta a los datos... ¿Sería válido quedarse con el resultado final que se obtiene con nls y "port", aunque el algoritmo señale que existe "Convergence failure: singular convergence (7)"?, los valores que obtiene "port" son (±SE): Km = 44.3 (±3.7), Kic = 6.9 (±2.1), Kiu = 1.1 x 10^6 (±1.2 x 10^10) y k3 = 0.0184 (±7.1 x 10^-4), SSE = 157.1 mientras que en el artículo reportan con SAS que (no menciona si es ±SE, SD o CI): Km = 46.1 (±1.6), Kic = 9.2 (±0.9), Kiu = 1.2 x 10^64 (±ND) y k3 = 0.0183 (±0.0002), SSE = 229.1> La aproximación que he utilizado es la siguiente: > > 1. He utilizado varios algoritmos de los que proporciona el paquete > nls2. En este paquete la aproximación es más de "fuerza bruta". > 2. De los varios algortimos que ofrece el paquete: fuerza-bruta, > random-search, plinear. Los resultados con fuerza-bruta y > random-search son equivalentes. He desechado el de plinear por no > ajustar correctamente. > 3. Para el algoritmo "brute-force", se puede definir un grid de > evaluación que permite al algoritmo evaluar diferentes opciones, > es algo más lento pero no excesivo. > 4. En este grid, he jugado con varios extremos para el valor de "Kiu" > que es el que se va de madre. He probado con valores desde 0 a > 1000, hasta quedarme con el último caso en el que lo he iterado > entre 10^11 y 10^20. > 5. Además de la librería "nls2", he utilizado también la librería > "nlstools" que ofrece una función "overview()" que proporciona el > intervalo de confianza de cada unos de los parámetros del modelo > (intervalos de confianza del 2.5% y del 97.5%). > > El código de todo esto es el siguiente: > > ################ > #### Fuerza Bruta > # Resultados paper SAS: Km = 46.1, Kic = 9.2, Kiu = 1.2 x 10^64, k3 = > 0.0183 > library(nls2) > library(nlstools) > > mi.me <http://mi.me><-t ~ > ((Km+Km*P0/Kic+Km*S0/Kic)*log((S0-(Pt-P0))/S0)+(-Km/Kic+1+P0/Kiu+S0/Kiu)*((S0-(Pt-P0))-S0)+(-1/(2*Kiu))*((S0-(Pt-P0))^2-S0^2))*(-1/(E*k3)) > > gr.id <http://gr.id><-expand.grid( > Km=seq(30,80,len=10), > Kic=seq(6,12,len=10), > Kiu=seq(10^11,10^20,len=10), > k3=seq(0.01,0.04,len=10) > ) > ### Brute-Force > aj.bf <http://aj.bf> <- nls2( > mi.me <http://mi.me>, > data=d.a, > start = gr.id <http://gr.id>, > algorithm = "brute-force" > ) > > summary(aj.bf <http://aj.bf>) > overview(aj.bf <http://aj.bf>) > ################ > > Y el resultado de overview(aj.bf <http://aj.bf>) es el siguiente (tan > sólo incluyo el intervalo de confianza): > > ------ > Asymptotic confidence interval: > 2.5% 97.5% > Km 4.451537e+01 5.992907e+01 > Kic 1.481808e+00 1.185153e+01 > Kiu -1.039715e+17 1.039717e+17 > k3 1.921825e-02 2.078175e-02 >El valor de Km obtenidos por fuerza bruta resulta un poco más alejado del artículo que el último obtenido usando nls con "port". Aunque considerando los intervalos de confianza no habría mucha diferencia. Estimate Std. Error t value Pr(>|t|) Km 5.222e+01 3.841e+00 13.60<2e-16 *** Kic 6.667e+00 2.584e+00 2.58 0.0127 * Kiu 1.000e+11 5.181e+16 0.00 1.0000 k3 2.000e-02 3.896e-04 51.34<2e-16 ***> Como detalles y preguntas adicionales: > > * No he probado con valores superiones de Kiu. Bueno, sí que he > probado con valores de 10^60 y 10^61 pero el modelo no converge. > Supongo que ya aparecen problemas de precisión... >Y al final el problema es muy similar al encontrado con nls, que no se llega al valor reportado de Kiu obtenido con SAS, aunque aquí la ventaja es que por fuerza bruta con nls2 sí converge.> * En el paquete nlstools hay ejemplos de ajuste de curvas de > Michaelis-Menten dentro del punto de su help "compet_mich". Puede > resultar ilustrativo ver el tipo de ajustes que se utilizan. >Ya había revisado la ayuda de nlstools. Ahí ajustan los modelos de michaelis menten sin inhibición, con inhibición competitiva y con inhibición no competitiva (en sus formas no integradas), pero a datos de velocidad inicial de reacción vs concentración de sustrato, la cual es la manera clásica (después del método de Lineweaver-Burk) de obtener los valores de Km, Vmáx, Kic y Kiu por medio de regresión no lineal. De hecho, esta ayuda fue muy ilustrativa y me permitió entender en general cómo usar nls. En el artículo también lo hacen así, calculando las velocidades iniciales de reacción en ausencia y presencia del inhibidor desde el tiempo cero (que es el mismo producto de la reacción), y al ajustar los modelos no integrados de Michaelis Menten usando nls con "port" se obtiene exactamente los valores reportados con SAS en el documento para el mismo tipo de ajuste.> * Tambén en el paquete nlstools, hay uhtilidades adicinales > (nlsConfRegions o nlsResiduals) que añaden detalles al análisis. > * De los datos que has aportado, P0 siempre es cero (0). Aunque son > datos empíricos, ¿ser siempre cero no podría ayudar a simplificar > la fórmula del modelo?. Vaya eliminar P0 de la fórmula y así > quitar un parámetro de ajuste. >Aquí lo que pasó fue que sólo envié los datos donde P0 es siempre cero pues el ajuste que me daba problemas era con este set de datos, pero existe otro set de datos donde P0 = 5. El total de datos (tanto con P0 = 0 como P0 = 5) se utiliza para ajustar todos los modelos (integrados y no integrados) y obtener los parámetros en presencia del inhibidor desde el tiempo cero de la reacción. Esto debido a que otro de los objetivos del artículo fue demostrar que el ajuste de los modelos integrados de Michaelis Menten es más preciso para obtener los valores de los parámetros que al ajustar los modelos no integrados (incluso en ausencia de una concentración inicial del inhibidor), puesto que no se tiene la impresición adicional de medir la velocidad inicial en presencia del inhibidor desde el inicio de la reacción.> * Otro tema que me genera ciertas dudas es el de que los datos > experimentales de "t" tienen la pinta de diferentes experimentos > donde "t" siempre se mueve en un rango (0 a 60). Y esto se repite > 8 veces. ¿No puede ayudar este patrón a ajustar un primer modelo > con cada una de las celdas?... >Tienes razón, t (tiempo en segundos) indica que hay ocho experimentos, monitoreados de 0 a 60 s, donde se determina la concentración de producto (Pt) liberado al medio de reacción, a una concentración de producto (S0). P0 es la concentración de producto añadido en el tiempo cero y E la concentración de enzima. Existe otro set de datos donde P0 = 5. Esto es lo que los autores denominan "curvas de progreso de reacción" (Pt vs t). La velocidad inicial de reacción para cada concentración de sustrato se puede calcular por regresión lineal de las diferentes curvas Pt vs t (una por cada valor de S0), para de esta forma ajustar los modelos no integrados de Michaelis-Menten a las curvas Vinicial vs S0, como lo hacen en la ayuda de nlstools. Sin embargo, el objetivo de los autores fue demostrar que ajustar los modelo integrados de Michaelis-Menten directamente a las curvas de progreso Pt vs t es más preciso que ajustar los modelos no integrados a las curvas Vinicial vs S0, lo que parecen lograr de acuerdo con los resultados mostrados en el documento.> Por último otra duda que puede ayudar. En el artículo donde comentan > esta aproximación vía código SAS, ¿se incluyen los intervalos de > confianza para cada parámetro?. A lo mejor los resultados que ya > obtenemos están en esos intervalos (lo digo más que nada por Kiu). >Más arriba puse los intervalos de error mostrados en el artículo, pero no indican si son error estándar, desviación estándar o intervalos de confianza. Al parecer los resultados obtenidos con R y nls caen dentro de estos intervalos, excepto para los valores de Kiu. Cabe notar que al ajustar los modelos integrados al total de los datos (con P0 = 0 y P0 = 5 simultáneamente) se obtienen casi exactamente los valores reportados por los autores con SAS, sin que nls tenga problemas para encontrar convergencia, excepto para el modelo mixto, en el que el valor de Kiu vuelve a ser cada vez más grande. Es más, los autores reportan en este caso que el modelo con inhibición no competitiva no converge, mientras que en R con nls sí se encuentra convergencia.> Saludos, > Carlos Ortega > www.qualityexcellence.es <http://www.qualityexcellence.es>De nuevo muchas gracias por tu ayuda Carlos. Saludos. -- Argel. [[alternative HTML version deleted]]
Hola Argel, Como referencia adicional, por si puede dar algún resultado o ayudar a enfocar el problema de una manera alternativa, hay un libro "Non Linear Regression with R" que también tiene un paquete asociado en R (nlrwr): http://www.amazon.es/Nonlinear-Regression-R-Use/dp/0387096159/ref=sr_1_sc_1?ie=UTF8&qid=1327441411&sr=8-1-spell http://cran.at.r-project.org/web/packages/nlrwr/index.html En este libro también se estudia el caso de esta ecuación y se contempla lo que has comentado de P0 con diferentes valores al que se le puede seguir aplicando la función nls() pero especificando que P0 indica realmente un factor. Extender el conjunto de datos con el caso que comentas de P0=5 podría ayudar al modelo. De todas formas llegado un punto de bloqueo también puedes escribir a los autores del paper y preguntarles por detalles de cómo codificaron su modelo, si han llegado o pueden utilizar R para replicar resultados. El impacto de su paper ya está conseguido, vaya supongo que te darían acceso a lo que han hecho... Saludos, Carlos Ortega www.qualityexcellence.es El 24 de enero de 2012 20:41, Argel Gastélum Arellánez < argel.gastelum@gmail.com> escribió:> Hola Carlos, muchas gracias por tu ayuda. Me parece muy interesante > la forma de aplicar nls2. Pongo algunas notas abajo. > > El 24/01/12 09:23, Carlos Ortega escribió: > > Hola, > > > > He probado con otra variante, visto que utilizando nls y los > > algoritmos de aproximación (default y port) están dando resultados > > todavía lejanos a los del artículo, además de los problemas de > > convergencia o de singularidad. > > ant > > En realidad el parámetro que sale bastante alejado de los > resultados del artículo es Kiu, los otros dan resultados muy similares. > Puesto que uno de los objetivos es comprobar que el modelo mixto no es > el que mejor se ajusta a los datos... ¿Sería válido quedarse con el > resultado final que se obtiene con nls y "port", aunque el algoritmo > señale que existe "Convergence failure: singular convergence (7)"?, los > valores que obtiene "port" son (±SE): > > Km = 44.3 (±3.7), Kic = 6.9 (±2.1), Kiu = 1.1 x 10^6 (±1.2 x 10^10) y k3 > = 0.0184 (±7.1 x 10^-4), SSE = 157.1 > > mientras que en el artículo reportan con SAS que (no menciona si es ±SE, > SD o CI): > > Km = 46.1 (±1.6), Kic = 9.2 (±0.9), Kiu = 1.2 x 10^64 (±ND) y k3 > 0.0183 (±0.0002), SSE = 229.1 > > > La aproximación que he utilizado es la siguiente: > > > > 1. He utilizado varios algoritmos de los que proporciona el paquete > > nls2. En este paquete la aproximación es más de "fuerza bruta". > > 2. De los varios algortimos que ofrece el paquete: fuerza-bruta, > > random-search, plinear. Los resultados con fuerza-bruta y > > random-search son equivalentes. He desechado el de plinear por no > > ajustar correctamente. > > 3. Para el algoritmo "brute-force", se puede definir un grid de > > evaluación que permite al algoritmo evaluar diferentes opciones, > > es algo más lento pero no excesivo. > > 4. En este grid, he jugado con varios extremos para el valor de "Kiu" > > que es el que se va de madre. He probado con valores desde 0 a > > 1000, hasta quedarme con el último caso en el que lo he iterado > > entre 10^11 y 10^20. > > 5. Además de la librería "nls2", he utilizado también la librería > > "nlstools" que ofrece una función "overview()" que proporciona el > > intervalo de confianza de cada unos de los parámetros del modelo > > (intervalos de confianza del 2.5% y del 97.5%). > > > > El código de todo esto es el siguiente: > > > > ################ > > #### Fuerza Bruta > > # Resultados paper SAS: Km = 46.1, Kic = 9.2, Kiu = 1.2 x 10^64, k3 > > 0.0183 > > library(nls2) > > library(nlstools) > > > > mi.me <http://mi.me><-t ~ > > > ((Km+Km*P0/Kic+Km*S0/Kic)*log((S0-(Pt-P0))/S0)+(-Km/Kic+1+P0/Kiu+S0/Kiu)*((S0-(Pt-P0))-S0)+(-1/(2*Kiu))*((S0-(Pt-P0))^2-S0^2))*(-1/(E*k3)) > > > > gr.id <http://gr.id><-expand.grid( > > Km=seq(30,80,len=10), > > Kic=seq(6,12,len=10), > > Kiu=seq(10^11,10^20,len=10), > > k3=seq(0.01,0.04,len=10) > > ) > > ### Brute-Force > > aj.bf <http://aj.bf> <- nls2( > > mi.me <http://mi.me>, > > data=d.a, > > start = gr.id <http://gr.id>, > > algorithm = "brute-force" > > ) > > > > summary(aj.bf <http://aj.bf>) > > overview(aj.bf <http://aj.bf>) > > ################ > > > > Y el resultado de overview(aj.bf <http://aj.bf>) es el siguiente (tan > > sólo incluyo el intervalo de confianza): > > > > ------ > > Asymptotic confidence interval: > > 2.5% 97.5% > > Km 4.451537e+01 5.992907e+01 > > Kic 1.481808e+00 1.185153e+01 > > Kiu -1.039715e+17 1.039717e+17 > > k3 1.921825e-02 2.078175e-02 > > > El valor de Km obtenidos por fuerza bruta resulta un poco más alejado > del artículo que el último obtenido usando nls con "port". Aunque > considerando los intervalos de confianza no habría mucha diferencia. > > Estimate Std. Error t value Pr(>|t|) > Km 5.222e+01 3.841e+00 13.60<2e-16 *** > Kic 6.667e+00 2.584e+00 2.58 0.0127 * > Kiu 1.000e+11 5.181e+16 0.00 1.0000 > k3 2.000e-02 3.896e-04 51.34<2e-16 *** > > > > > Como detalles y preguntas adicionales: > > > > * No he probado con valores superiones de Kiu. Bueno, sí que he > > probado con valores de 10^60 y 10^61 pero el modelo no converge. > > Supongo que ya aparecen problemas de precisión... > > > > Y al final el problema es muy similar al encontrado con nls, que no > se llega al valor reportado de Kiu obtenido con SAS, aunque aquí la > ventaja es que por fuerza bruta con nls2 sí converge. > > > * En el paquete nlstools hay ejemplos de ajuste de curvas de > > Michaelis-Menten dentro del punto de su help "compet_mich". Puede > > resultar ilustrativo ver el tipo de ajustes que se utilizan. > > > > Ya había revisado la ayuda de nlstools. Ahí ajustan los modelos de > michaelis menten sin inhibición, con inhibición competitiva y con > inhibición no competitiva (en sus formas no integradas), pero a datos de > velocidad inicial de reacción vs concentración de sustrato, la cual es > la manera clásica (después del método de Lineweaver-Burk) de obtener los > valores de Km, Vmáx, Kic y Kiu por medio de regresión no lineal. De > hecho, esta ayuda fue muy ilustrativa y me permitió entender en general > cómo usar nls. > > En el artículo también lo hacen así, calculando las velocidades > iniciales de reacción en ausencia y presencia del inhibidor desde el > tiempo cero (que es el mismo producto de la reacción), y al ajustar los > modelos no integrados de Michaelis Menten usando nls con "port" se > obtiene exactamente los valores reportados con SAS en el documento para > el mismo tipo de ajuste. > > > * Tambén en el paquete nlstools, hay uhtilidades adicinales > > (nlsConfRegions o nlsResiduals) que añaden detalles al análisis. > > * De los datos que has aportado, P0 siempre es cero (0). Aunque son > > datos empíricos, ¿ser siempre cero no podría ayudar a simplificar > > la fórmula del modelo?. Vaya eliminar P0 de la fórmula y así > > quitar un parámetro de ajuste. > > > > Aquí lo que pasó fue que sólo envié los datos donde P0 es siempre > cero pues el ajuste que me daba problemas era con este set de datos, > pero existe otro set de datos donde P0 = 5. El total de datos (tanto con > P0 = 0 como P0 = 5) se utiliza para ajustar todos los modelos > (integrados y no integrados) y obtener los parámetros en presencia del > inhibidor desde el tiempo cero de la reacción. > > Esto debido a que otro de los objetivos del artículo fue demostrar > que el ajuste de los modelos integrados de Michaelis Menten es más > preciso para obtener los valores de los parámetros que al ajustar los > modelos no integrados (incluso en ausencia de una concentración inicial > del inhibidor), puesto que no se tiene la impresición adicional de medir > la velocidad inicial en presencia del inhibidor desde el inicio de la > reacción. > > > * Otro tema que me genera ciertas dudas es el de que los datos > > experimentales de "t" tienen la pinta de diferentes experimentos > > donde "t" siempre se mueve en un rango (0 a 60). Y esto se repite > > 8 veces. ¿No puede ayudar este patrón a ajustar un primer modelo > > con cada una de las celdas?... > > > > Tienes razón, t (tiempo en segundos) indica que hay ocho > experimentos, monitoreados de 0 a 60 s, donde se determina la > concentración de producto (Pt) liberado al medio de reacción, a una > concentración de producto (S0). P0 es la concentración de producto > añadido en el tiempo cero y E la concentración de enzima. Existe otro > set de datos donde P0 = 5. Esto es lo que los autores denominan "curvas > de progreso de reacción" (Pt vs t). > > La velocidad inicial de reacción para cada concentración de > sustrato se puede calcular por regresión lineal de las diferentes curvas > Pt vs t (una por cada valor de S0), para de esta forma ajustar los > modelos no integrados de Michaelis-Menten a las curvas Vinicial vs S0, > como lo hacen en la ayuda de nlstools. > > Sin embargo, el objetivo de los autores fue demostrar que ajustar > los modelo integrados de Michaelis-Menten directamente a las curvas de > progreso Pt vs t es más preciso que ajustar los modelos no integrados a > las curvas Vinicial vs S0, lo que parecen lograr de acuerdo con los > resultados mostrados en el documento. > > > Por último otra duda que puede ayudar. En el artículo donde comentan > > esta aproximación vía código SAS, ¿se incluyen los intervalos de > > confianza para cada parámetro?. A lo mejor los resultados que ya > > obtenemos están en esos intervalos (lo digo más que nada por Kiu). > > > > Más arriba puse los intervalos de error mostrados en el artículo, > pero no indican si son error estándar, desviación estándar o intervalos > de confianza. Al parecer los resultados obtenidos con R y nls caen > dentro de estos intervalos, excepto para los valores de Kiu. > > Cabe notar que al ajustar los modelos integrados al total de los > datos (con P0 = 0 y P0 = 5 simultáneamente) se obtienen casi exactamente > los valores reportados por los autores con SAS, sin que nls tenga > problemas para encontrar convergencia, excepto para el modelo mixto, en > el que el valor de Kiu vuelve a ser cada vez más grande. Es más, los > autores reportan en este caso que el modelo con inhibición no > competitiva no converge, mientras que en R con nls sí se encuentra > convergencia. > > > Saludos, > > Carlos Ortega > > www.qualityexcellence.es <http://www.qualityexcellence.es> > > De nuevo muchas gracias por tu ayuda Carlos. > > Saludos. > > -- > Argel. > > [[alternative HTML version deleted]] > > > _______________________________________________ > R-help-es mailing list > R-help-es@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-help-es > >[[alternative HTML version deleted]]
Argel Gastélum Arellánez
2012-Jan-24 22:58 UTC
[R-es] Ajuste no lineal con nls no converge
Hola Carlos, El 24/01/12 15:54, Carlos Ortega escribió:> Hola Argel, > > Como referencia adicional, por si puede dar algún resultado o ayudar a > enfocar el problema de una manera alternativa, hay un libro "Non > Linear Regression with R" que también tiene un paquete asociado en R > (nlrwr): > > http://www.amazon.es/Nonlinear-Regression-R-Use/dp/0387096159/ref=sr_1_sc_1?ie=UTF8&qid=1327441411&sr=8-1-spell > <http://www.amazon.es/Nonlinear-Regression-R-Use/dp/0387096159/ref=sr_1_sc_1?ie=UTF8&qid=1327441411&sr=8-1-spell> > > http://cran.at.r-project.org/web/packages/nlrwr/index.html > > En este libro también se estudia el caso de esta ecuación y se > contempla lo que has comentado de P0 con diferentes valores al que se > le puede seguir aplicando la función nls() pero especificando que P0 > indica realmente un factor. Extender el conjunto de datos con el caso > que comentas de P0=5 podría ayudar al modelo.Acabo de conseguir el libro, pero apenas lo voy comenzando a leer. Estudiaré con detalle la sección que comentas. Me queda la siguiente duda al tratar de ajustar el modelo mixto con nls, ¿sería válido quedarse con el resultado final que se obtiene con nls y "port", aunque el algoritmo señale que existe "Convergence failure: singular convergence (7)"? Lo anterior considerando que los valores de Km, Kic y k3 obtenidos con ambos sistemas de análisis son muy similares (excepto por Kiu que en SAS es muy grande y que en nls no converge pero aumenta con tendencia hacia valores muy grandes) la interpretación sería la misma que la de los autores, donde la constante Kiu que alcanza valores tan grandes hace que realmente carezca de sentido incluirla para modelar la reacción enzimática en estudio... quedando sólo las constantes Km, Kic y k3 en el modelo que ahora sería de inhibición competitiva.> De todas formas llegado un punto de bloqueo también puedes escribir a > los autores del paper y preguntarles por detalles de cómo codificaron > su modelo, si han llegado o pueden utilizar R para replicar > resultados. El impacto de su paper ya está conseguido, vaya supongo > que te darían acceso a lo que han hecho... >Buena idea, espero obtener respuesta y poder complementar todo lo que hasta ahora se ha discutido en este hilo.> Saludos, > Carlos Ortega > www.qualityexcellence.es <http://www.qualityexcellence.es>De nuevo, muchas gracias por tus sugerencias. Saludos. -- Argel. [[alternative HTML version deleted]]
Hola, Con el algoritmo "port", ví en la ayuda de nls lo siguiente: The algorithm = "port" code appears unfinished, and does not even check that the starting value is within the bounds. Use with caution, especially where bounds are supplied. Por este mensaje que no ayuda a tener mucha confianza, y por que ya lo habíais estado explorando vosotros es por lo que me decidí a probar otras cosas. Aunque tengas un mensaje de no-convergencia, como tienes una solución, es de suponer que no será óptima. Si decides utilizarla tendrás que adjunar a tu resultado (al de cada coeficiente) el error (intervalo de confianza), la bondad del ajuste, algún análisis de residuos por si pudieras aplicar alguna transformación, etc. Todo esto lo consigues con funciones que estan en los paquetes nlstools y en nlrwr. No he visto tampoco si es posible una representación grafica de la solución ajustada sobre los datos experimentales, que vaya suele aparecer de las primeras evidencias para ver la bondad del ajuste. Saludos, Carlos Ortega www.qualityexcellence.es El 24 de enero de 2012 23:58, Argel Gastélum Arellánez < argel.gastelum@gmail.com> escribió:> Hola Carlos, > > El 24/01/12 15:54, Carlos Ortega escribió: > > Hola Argel, > > > > Como referencia adicional, por si puede dar algún resultado o ayudar a > > enfocar el problema de una manera alternativa, hay un libro "Non > > Linear Regression with R" que también tiene un paquete asociado en R > > (nlrwr): > > > > > http://www.amazon.es/Nonlinear-Regression-R-Use/dp/0387096159/ref=sr_1_sc_1?ie=UTF8&qid=1327441411&sr=8-1-spell > > < > http://www.amazon.es/Nonlinear-Regression-R-Use/dp/0387096159/ref=sr_1_sc_1?ie=UTF8&qid=1327441411&sr=8-1-spell > > > > > > http://cran.at.r-project.org/web/packages/nlrwr/index.html > > > > En este libro también se estudia el caso de esta ecuación y se > > contempla lo que has comentado de P0 con diferentes valores al que se > > le puede seguir aplicando la función nls() pero especificando que P0 > > indica realmente un factor. Extender el conjunto de datos con el caso > > que comentas de P0=5 podría ayudar al modelo. > > Acabo de conseguir el libro, pero apenas lo voy comenzando a leer. > Estudiaré con detalle la sección que comentas. > > Me queda la siguiente duda al tratar de ajustar el modelo mixto con > nls, ¿sería válido quedarse con el resultado final que se obtiene con > nls y "port", aunque el algoritmo señale que existe "Convergence > failure: singular convergence (7)"? > > Lo anterior considerando que los valores de Km, Kic y k3 obtenidos > con ambos sistemas de análisis son muy similares (excepto por Kiu que en > SAS es muy grande y que en nls no converge pero aumenta con tendencia > hacia valores muy grandes) la interpretación sería la misma que la de > los autores, donde la constante Kiu que alcanza valores tan grandes hace > que realmente carezca de sentido incluirla para modelar la reacción > enzimática en estudio... quedando sólo las constantes Km, Kic y k3 en el > modelo que ahora sería de inhibición competitiva. > > > De todas formas llegado un punto de bloqueo también puedes escribir a > > los autores del paper y preguntarles por detalles de cómo codificaron > > su modelo, si han llegado o pueden utilizar R para replicar > > resultados. El impacto de su paper ya está conseguido, vaya supongo > > que te darían acceso a lo que han hecho... > > > > Buena idea, espero obtener respuesta y poder complementar todo lo > que hasta ahora se ha discutido en este hilo. > > > Saludos, > > Carlos Ortega > > www.qualityexcellence.es <http://www.qualityexcellence.es> > > De nuevo, muchas gracias por tus sugerencias. > > Saludos. > > -- > Argel. > > [[alternative HTML version deleted]] > > > _______________________________________________ > R-help-es mailing list > R-help-es@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-help-es > >[[alternative HTML version deleted]]
Argel Gastélum Arellánez
2012-Jan-25 17:50 UTC
[R-es] Ajuste no lineal con nls no converge
Hola Carlos, qué tal. Muchas gracias por tus comentarios. Pongo algunas notas más abajo. El 25/01/12 07:55, Carlos Ortega escribió:> Hola, > > Con el algoritmo "port", ví en la ayuda de nls lo siguiente: > > The|algorithm = "port"|code appears unfinished, and does not even > check that the starting value is within the bounds. Use with caution, > especially where bounds are supplied.Había visto esta aclaración en la ayuda, pero como me daba los mismos resultados que "default" (excepto para cuando no había convergencia) no me había precupado demasiado por eso, aunque lo tenía en mente por si me salía algo raro en la estimación de los parámetros. Creo entender que el problema principalmente es con los límites asignados a los parámetros, que por el momento no estoy asignando, aunque tal vez sí sea conveniente puesto que todos los parámetros serían mayores a cero.> Por este mensaje que no ayuda a tener mucha confianza, y por que ya lo > habíais estado explorando vosotros es por lo que me decidí a probar > otras cosas. > > Aunque tengas un mensaje de no-convergencia, como tienes una solución, > es de suponer que no será óptima. Si decides utilizarla tendrás que > adjunar a tu resultado (al de cada coeficiente) el error (intervalo de > confianza), la bondad del ajuste, algún análisis de residuos por si > pudieras aplicar alguna transformación, etc. Todo esto lo consigues > con funciones que estan en los paquetes nlstools y en nlrwr. >Aquí pasa algo curioso, para el modelo mixto evaluado en SAS por los autores del artículo, enuentran convergencia con un valor de SSE = 229.1, mientras que nls con "port" no encuentra convergencia (debido a la tendencia a crecer en el valor de Kiu), sin embargo, el valor final que se supone no óptimo tiene un valor de SSE = 157.1 ¿Esto indicaría que probablemente el ajuste obtenido en R (aún sin encontrar convergencia) es mejor que el obtenido en SAS?> No he visto tampoco si es posible una representación grafica de la > solución ajustada sobre los datos experimentales, que vaya suele > aparecer de las primeras evidencias para ver la bondad del ajuste. > > Saludos, > Carlos Ortega > www.qualityexcellence.es <http://www.qualityexcellence.es>De nuevo muchas gracias por tu ayuda. Saludos. -- Argel. [[alternative HTML version deleted]]
Argel Gastélum Arellánez
2012-Jan-27 18:04 UTC
[R-es] Ajuste no lineal con nls no converge
Hola Carlos, muchas gracias por seguir revisando el tema. Pongo algunas notas abajo. El 27/01/12 04:24, Carlos Ortega escribió:> Hola Argel, > > He continuado mirando este tema, pero buscando la función que hay que > ajustar podía expresarse analíticamente de otra manera. También he > mirado de qué forma queda expresada la función si se considera que P0 > es cero y que Kiu tiende a infinito. >Aquí no creo conveniente considerar que siempre P0 = 0 para simplificar el modelo, porque en realidad puede tomar valores iniciales. No lo puse en los datos que envié inicialmente porque el ajuste que no me daba convergencia era con este set de datos donde P0 = 0 (en el que los autores reportan que sí converge), pero el set de datos completo incluye una serie de 7 experimentos adicionales donde P0 = 5, y el análisis completo involucra el ajuste con los dos grupos de datos.> Pero en el camino me he encontrado con un detalle que no sé si es un > error o es algo que no lo estoy yo haciendo bien. > En la expresión: > > ((Km+Km*P0/Kic+Km*S0/Kic)*log( > (S0-(Pt-P0))/S0)+(-Km/Kic+1+P0/Kiu+S0/Kiu)*((S0-(Pt-P0))-S0)+(-1/(2*Kiu))*((S0-(Pt-P0))^2-S0^2))*(-1/(E*k3)) > > Mira el término, que aunque tiene muchos paréntesis: > > ((S0-(Pt-P0))-S0) > > Esto desarrollado es: > ( (S0 - Pt + P0) - S0) = ( S0 - Pt + P0 - S0 ) = (P0 - Pt) .... ¿Se > cancelan los S0? >Tienes razón, se cancelan los S0 y queda sólo (P0-Pt). Así reportan los autores que ingresaron el modelo en SAS. Ya lo había notado y lo cambié para ingresarlo en nls dejando sólo (P0-Pt) y el resultado es el mismo, no converge, aunque los valores finales de "port" aunque son similares a los reportados por los autores con SAS, al parecer ajustan mejor, de acuerdo con las sumas de cuadrados que comentaba en el mensaje anterior.> No he probado a ajustar la función si esto se confirma, por si hay > cambios en el ajuste nls... > > Saludos, > Carlos Ortega > www.qualityexcellence.es <http://www.qualityexcellence.es>Saludos. -- Argel. [[alternative HTML version deleted]]
Hola Angel, Podrias por favor compartir la referencia de la que hablas, i.e., el paper en el que usan SAS? Puede ser un link. Aunque he trabajado relativamente poco con modelos nls, te sugeriria utilizar los parametros estimados como valores iniciales en el proceso de optimizacion; esta estrategia me ha funcionado en varias ocasiones. El ultimo ejemplo de ?nls ilustra dicha estrategia. Un saludo, Jorge.- 2012/1/27 Argel Gastélum Arellánez <>> Hola Carlos, muchas gracias por seguir revisando el tema. Pongo > algunas notas abajo. > > El 27/01/12 04:24, Carlos Ortega escribió: > > Hola Argel, > > > > He continuado mirando este tema, pero buscando la función que hay que > > ajustar podía expresarse analíticamente de otra manera. También he > > mirado de qué forma queda expresada la función si se considera que P0 > > es cero y que Kiu tiende a infinito. > > > > Aquí no creo conveniente considerar que siempre P0 = 0 para > simplificar el modelo, porque en realidad puede tomar valores iniciales. > No lo puse en los datos que envié inicialmente porque el ajuste que no > me daba convergencia era con este set de datos donde P0 = 0 (en el que > los autores reportan que sí converge), pero el set de datos completo > incluye una serie de 7 experimentos adicionales donde P0 = 5, y el > análisis completo involucra el ajuste con los dos grupos de datos. > > > Pero en el camino me he encontrado con un detalle que no sé si es un > > error o es algo que no lo estoy yo haciendo bien. > > En la expresión: > > > > ((Km+Km*P0/Kic+Km*S0/Kic)*log( > > > (S0-(Pt-P0))/S0)+(-Km/Kic+1+P0/Kiu+S0/Kiu)*((S0-(Pt-P0))-S0)+(-1/(2*Kiu))*((S0-(Pt-P0))^2-S0^2))*(-1/(E*k3)) > > > > Mira el término, que aunque tiene muchos paréntesis: > > > > ((S0-(Pt-P0))-S0) > > > > Esto desarrollado es: > > ( (S0 - Pt + P0) - S0) = ( S0 - Pt + P0 - S0 ) = (P0 - Pt) .... ¿Se > > cancelan los S0? > > > > Tienes razón, se cancelan los S0 y queda sólo (P0-Pt). Así reportan > los autores que ingresaron el modelo en SAS. Ya lo había notado y lo > cambié para ingresarlo en nls dejando sólo (P0-Pt) y el resultado es el > mismo, no converge, aunque los valores finales de "port" aunque son > similares a los reportados por los autores con SAS, al parecer ajustan > mejor, de acuerdo con las sumas de cuadrados que comentaba en el mensaje > anterior. > > > No he probado a ajustar la función si esto se confirma, por si hay > > cambios en el ajuste nls... > > > > Saludos, > > Carlos Ortega > > www.qualityexcellence.es <http://www.qualityexcellence.es> > > Saludos. > > -- > Argel. > > [[alternative HTML version deleted]] > > > _______________________________________________ > R-help-es mailing list > R-help-es@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-help-es > >[[alternative HTML version deleted]]
Argel Gastélum Arellánez
2012-Jan-27 22:15 UTC
[R-es] Ajuste no lineal con nls no converge
Hola Jorge, qué tal... El 27/01/12 14:17, Jorge I Velez escribió:> Hola Angel, > > Podrias por favor compartir la referencia de la que hablas, i.e., el > paper en el que usan SAS? Puede ser un link. >Claro que sí, la referencia es la siguiente: "Bezerra y Dias, 2007. Utilization of Integrated Michaelis-Menten Equation to Determine Kinetic Constants. Biochemistry and molecular biology education. Vol. 35, No. 2, pp. 145--150". El artículo se puede desacargar de aquí: http://onlinelibrary.wiley.com/doi/10.1002/bmb.32/pdf> Aunque he trabajado relativamente poco con modelos nls, te sugeriria > utilizar los parametros estimados como valores iniciales en el proceso > de optimizacion; esta estrategia me ha funcionado en varias ocasiones. > El ultimo ejemplo de ?nls ilustra dicha estrategia. >También lo había intentado ya y sigue sin converger.> Un saludo, > Jorge.-Saludos. -- Argel. [[alternative HTML version deleted]]