Obtener los valores iniciales correctos para un modelo nls en R

Estoy tratando de ajustar un modelo de ley de potencia simple a un conjunto de datos que es el siguiente:

mydf:

rev     weeks
17906.4 1
5303.72 2
2700.58 3
1696.77 4
947.53  5
362.03  6

El objetivo es pasar la línea eléctrica y usarla para predecir revvlaues para las próximas semanas. Un montón de investigación me ha llevado a la nlsfunción, que implementé de la siguiente manera.

newMod <- nls(rev ~ a*weeks^b, data=modeldf, start = list(a=1,b=1))
predict(newMod, newdata = data.frame(weeks=c(1,2,3,4,5,6,7,8,9,10)))

Si bien esto funciona para un lmmodelo, recibo un singular gradienterror, que entiendo que tiene que ver con mis valores iniciales ay b. Intenté diferentes valores, incluso llegando a trazar esto en Excel, pasar un solitario, obtener una ecuación, luego usar los valores de la ecuación, pero aún obtengo el error. Miré un montón de respuestas como esta y probé la segunda respuesta (no pude entender la primera), pero no obtuve ningún resultado.

Realmente podría usar algo de ayuda aquí sobre cómo encontrar los valores iniciales correctos. O, alternativamente, qué otra función puedo usar en lugar de nls.

En caso de que quiera recrear mydfcon facilidad:

mydf <- data.frame(rev=c(17906.4, 5303.72, 2700.58 ,1696.77 ,947.53 ,362.03), weeks=c(1,2,3,4,5,6)) 

Este es un problema común con los modelos de mínimos cuadrados no lineales; Si sus valores iniciales están muy lejos de lo óptimo, el algoritmo puede no converger, incluso si se comporta bien cerca del óptimo.

Si se inicia mediante la adopción de los registros de ambas partes y ajustar un modelo lineal, se obtiene estimaciones de y como la pendiente y la intersección (9,947 y -2,011) (edit: eso es logaritmo natural)$\log(a)$$b$

Si los usa para guiar los valores iniciales de y todo parece funcionar bien:$a$$b$

 newMod <- nls(rev ~ a*weeks^b, data=mydf, start = list(a=exp(9.947),b=-2.011))
 predict(newMod, newdata = data.frame(weeks=c(1,2,3,4,5,6,7,8,9,10)))
 [1] 17919.2138  5280.7001  2584.0109  1556.1951  1050.1230   761.4947   580.3091   458.6027
 [9]   372.6231   309.4658

Tratar

 newMod <- nls(rev ~ a*weeks^b, data=mydf, startlist(a=17919.2127344,b=-1.76270557120))

Me han pedido que amplíe un poco esta respuesta. Este problema es tan simple que me sorprende que nls falle en eso. Sin embargo, el verdadero problema es con todo el enfoque R y la filosofía del ajuste de modelos no lineales. En el mundo real, uno escalaría x para situarse entre -1 y 1 e y e y para situarse entre 0 y 1 (y = ax ^ b). Eso probablemente sería suficiente para lograr que nls converja. Por supuesto, como señala Glen, puede ajustar el modelo log-lineal correspondiente. Eso se basa en el hecho de que existe una transformación simple que linealiza el modelo. Ese no suele ser el caso. El problema con las rutinas R como nls es que no ofrecen soporte para volver a parametrizar el modelo. En este caso, la reparametrización es simple, solo reescala / recientere x e y. Sin embargo, una vez que se ajuste al modelo, el usuario tendrá diferentes parámetros ayb de los originales. Si bien es simple calcular los originales a partir de estos, la otra dificultad es que no es tan simple en general obtener las desviaciones estándar estimadas para estas estimaciones de parámetros. Esto se realiza mediante el método delta que involucra el hessiano de la probabilidad logarítmica y algunas derivadas. El software de estimación de parámetros no lineales debería proporcionar estos cálculos automáticamente, de modo que la reparametrización del modelo sea fácilmente compatible. Otra cosa que el software debería soportar es la noción de fases. Puede pensar en ajustar primero el modelo con la versión de Glen como la fase 1. El modelo "real" se ajusta en la etapa 2. la otra dificultad es que no es tan simple en general obtener las desviaciones estándar estimadas para estas estimaciones de parámetros. Esto se realiza mediante el método delta que involucra el hessiano de la probabilidad logarítmica y algunas derivadas. El software de estimación de parámetros no lineales debería proporcionar estos cálculos automáticamente, de modo que la reparametrización del modelo sea fácilmente compatible. Otra cosa que el software debería soportar es la noción de fases. Puede pensar en ajustar primero el modelo con la versión de Glen como la fase 1. El modelo "real" se ajusta en la etapa 2. la otra dificultad es que no es tan simple en general obtener las desviaciones estándar estimadas para estas estimaciones de parámetros. Esto se realiza mediante el método delta que involucra el hessiano de la probabilidad logarítmica y algunas derivadas. El software de estimación de parámetros no lineales debería proporcionar estos cálculos automáticamente, de modo que la reparametrización del modelo sea fácilmente compatible. Otra cosa que el software debería soportar es la noción de fases. Puede pensar en ajustar primero el modelo con la versión de Glen como la fase 1. El modelo "real" se ajusta en la etapa 2. El software de estimación de parámetros no lineales debería proporcionar estos cálculos automáticamente, de modo que la reparametrización del modelo sea fácilmente compatible. Otra cosa que el software debería soportar es la noción de fases. Puede pensar en ajustar primero el modelo con la versión de Glen como la fase 1. El modelo "real" se ajusta en la etapa 2. El software de estimación de parámetros no lineales debería proporcionar estos cálculos automáticamente, de modo que la reparametrización del modelo sea fácilmente compatible. Otra cosa que el software debería soportar es la noción de fases. Puede pensar en ajustar primero el modelo con la versión de Glen como la fase 1. El modelo "real" se ajusta en la etapa 2.

Encajo su modelo con AD Model Builder que admite fases de forma natural. En la primera fase solo se estimó a. Esto lleva a tu modelo al estadio de béisbol. En la segunda fase, se estima que ayb obtienen la solución. AD Model Builder calcula automáticamente las desviaciones estándar para cualquier función de los parámetros del modelo a través del método delta, de modo que fomenta la reparametrización estable del modelo.

El algoritmo Levenberg-Marquardt puede ayudar:

modeldf <- data.frame(rev=c(17906.4, 5303.72, 2700.58 ,1696.77 ,947.53 ,362.03), weeks=c(1,2,3,4,5,6))

require(minpack.lm)
fit <- nlsLM(rev ~ a*weeks^b, data=modeldf, start = list(a=1,b=1))

require(broom)
fit_data <- augment(fit)

plot(.fitted~rev, data=fit_data)

En mi experiencia, una buena manera de encontrar valores iniciales para los parámetros de los modelos NLR es usar un algoritmo evolutivo. De una población inicial (100) de estimaciones aleatorias (padres) en un espacio de búsqueda, elija las mejores 20 (descendencia) y utilícelas para ayudar a definir una búsqueda en una población subsiguiente. Repita hasta la convergencia. No hay necesidad de gradientes o arpillera, solo evaluaciones SSE. Si no eres demasiado codicioso, esto a menudo funciona. El problema que las personas suelen tener es que están utilizando una búsqueda local (Newton-Raphson) para realizar el trabajo de una búsqueda global. Como siempre, se trata de utilizar la herramienta correcta para el trabajo en cuestión. Tiene más sentido usar una búsqueda global de EA para encontrar valores iniciales para la búsqueda local de Newton, y luego dejar que esto se reduzca al mínimo. Pero, como con todas las cosas, el diablo está en los detalles.