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]]

