[R-es] Ajuste no lineal con nls no converge

Argel Gastélum Arellánez argel.gastelum en gmail.com
Lun Ene 23 22:55:14 CET 2012


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



Más información sobre la lista de distribución R-help-es