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

