¿Por qué nls() me da errores de "matriz de gradiente singular en las estimaciones iniciales de los parámetros"?

Question

Tengo algunos datos básicos sobre la reducción de emisiones y el coste por coche:

q24 <- read.table(text = "reductions  cost.per.car
    50  45
    55  55
    60  62
    65  70
    70  80
    75  90
    80  100
    85  200
    90  375
    95  600
    ",header = TRUE, sep = "")

Sé que se trata de una función exponencial, por lo que espero poder encontrar un modelo que se ajuste a:

    model <- nls(cost.per.car ~ a * exp(b * reductions) + c, 
         data = q24, 
         start = list(a=1, b=1, c=0))

pero me da un error:

Error in nlsModel(formula, mf, start, wts) : 
  singular gradient matrix at initial parameter estimates

He leído a través de un montón de preguntas sobre el error que estoy viendo y estoy deduciendo que el problema es probablemente que necesito mejor/diferente start valores (el initial parameter estimates tiene un poco más de sentido) pero no estoy seguro, dados los datos que tengo, de cómo haría para estimar mejores parámetros.

Answer 1

Automáticamente Encontrar buenos valores de partida para un modelo no lineal es un arte. (Es relativamente fácil en el caso de conjuntos de datos puntuales cuando se pueden trazar los datos y hacer algunas buenas conjeturas visualmente). Un enfoque es linealizar el modelo y utilizar las estimaciones por mínimos cuadrados.

En este caso, el modelo tiene la forma

$$\mathbb{E}(Y) = a \exp(b x) + c$$

para parámetros desconocidos $a,b,c$ . La presencia del exponencial nos anima a utilizar logaritmos, pero la adición de $c$ hace que sea difícil hacerlo. Sin embargo, observe que si $a$ es positivo, entonces $c$ será menor que el menor valor esperado de $Y$ --y por lo tanto podría ser un poco menos que el más pequeño observado valor de $Y$ . (Si $a$ podría ser negativo también tendrá que considerar un valor de $c$ que es un poco mayor que el mayor valor observado de $Y$ .)

Entonces, ocupémonos de $c$ utilizando como estimación inicial $c_0$ algo así como la mitad del mínimo de las observaciones $y_i$ . El modelo se puede reescribir ahora sin ese espinoso término aditivo como

$$\mathbb{E}(Y) - c_0 \approx a \exp(b x).$$

Que podemos tomar el registro de:

$$\log(\mathbb{E}(Y) - c_0) \approx \log(a) + b x.$$

Es una aproximación lineal al modelo. Ambos $\log(a)$ y $b$ se puede estimar con mínimos cuadrados.

Aquí está el código revisado:

c.0 <- min(q24$cost.per.car) * 0.5
model.0 <- lm(log(cost.per.car - c.0) ~ reductions, data=q24)
start <- list(a=exp(coef(model.0)[1]), b=coef(model.0)[2], c=c.0)
model <- nls(cost.per.car ~ a * exp(b * reductions) + c, data = q24, start = start)

Su salida (para los datos del ejemplo) es

Nonlinear regression model
  model: cost.per.car ~ a * exp(b * reductions) + c
   data: q24
        a         b         c 
 0.003289  0.126805 48.487386 
 residual sum-of-squares: 2243

Number of iterations to convergence: 38 
Achieved convergence tolerance: 1.374e-06

La convergencia parece buena. Vamos a trazarla:

plot(q24)
p <- coef(model)
curve(p["a"] * exp(p["b"] * x) + p["c"], lwd=2, col="Red", add=TRUE)

Ha funcionado bien.

Al automatizar esto, podría realizar algunos análisis rápidos de los residuos, como comparar sus extremos con la dispersión en el ( $y$ ). También podría necesitar un código análogo para tratar la posibilidad $a\lt 0$ Lo dejo como un ejercicio.

---

Otro método para estimar los valores iniciales se basa en la comprensión de su significado, que puede basarse en la experiencia, la teoría física, etc. Un ejemplo extendido de un ajuste no lineal (moderadamente difícil) cuyos valores iniciales pueden determinarse de esta manera se describe en mi respuesta en https://stats.stackexchange.com/a/15769 .

Análisis visual de un gráfico de dispersión (para determinar las estimaciones iniciales de los parámetros) se describe e ilustra en https://stats.stackexchange.com/a/32832 .

En algunas circunstancias, se realiza una secuencia de ajustes no lineales en la que cabe esperar que las soluciones cambien lentamente. En ese caso, suele ser conveniente (y rápido) utilizar las soluciones anteriores como estimaciones iniciales para las siguientes . Recuerdo haber utilizado esta técnica (sin comentarios) en https://stats.stackexchange.com/a/63169 .

Answer 2

Esta biblioteca fue capaz de resolver mi problema con nls's singular gradient : http://www.r-bloggers.com/a-better-nls/ Un ejemplo:

library(minpack.lm)
nlsLM(function, start=list(variable=2,variable2=12))

Answer 3

Así que... Creo que leí mal esto como una función exponencial. Todo lo que necesitaba era poly()

model <- lm(cost.per.car ~ poly(reductions, 3), data=q24)
new.data <- data.frame(reductions = c(91,92,93,94))
predict(model, new.data)

plot(q24)
lines(q24$reductions, predict(model, list(reductions = q24$reductions)))

O bien, utilizando lattice :

xyplot(cost.per.car ~ reductions, data = q24,
       panel = function(x, y) {
         panel.xyplot(x, y)
         panel.lines(x, predict(model,list(reductions = x) ))
       }, 
       xlab = "Reductions", 
       ylab = "Cost per car")