[R-es] Ajuste no lineal con nls no converge
Argel Gastélum Arellánez
argel.gastelum en gmail.com
Mar Ene 24 01:26:09 CET 2012
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. Olivier
Comento 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.
Más información sobre la lista de distribución R-help-es