sintaxis para el modelo nls con punto de ruptura

Question

Así que el saga continúa...

Así que estoy tratando de ajustar el modelo $$Runoff = \begin{cases} \beta_0 + \beta_1Pcp & \text{if $(Ant + Pcp) < Thold$;}\\ \beta_2 + \beta_3Pcp & \text{if $(Ant + Pcp) \geq Thold$;} \end{cases}$$ Donde $Pcp$ y $Ant$ son variables observadas, y $Thold$ y se estiman los demás parámetros. No creo que pueda hacer esto con el segmentado ya que el punto de ruptura no está en las variables predictoras, pero creo que debería poder hacerlo de forma bastante sencilla como una regresión no lineal por mínimos cuadrados. No soy muy bueno con la sintaxis nls, así que cualquier ayuda sería apreciada

Answer 1

En general, este es un problema tan desagradable que no deberíamos aplicar optimizadores automáticos como nls en R . Sin embargo, es fácilmente solucionable al observar que el modelo es lineal, y cumple los supuestos de la estimación por mínimos cuadrados ordinarios (MCO), condicionada al valor de Thresh . Por lo tanto, se pueden buscar soluciones de forma fiable variando sistemáticamente Thresh en un rango razonable de valores.

---

Para ilustrar, He simulado algunos datos de esta forma en R utilizando el modelo equivalente

$$Y = \beta_0 + \beta_1 x_1 + \beta_2 I(x_2 \lt \tau) + \beta_3 x_1 I(x_2 \lt \tau) + \varepsilon$$

donde $\beta_0, \ldots, \beta_3$ son los coeficientes (¡pero no son los mismos que los de la pregunta!), $I$ es la función indicadora, $\tau$ es el parámetro del umbral, $x_1$ representa valores de Pcp , $x_2$ representa valores de Pcp + Ant y $\varepsilon$ tiene una distribución normal de media cero y varianza desconocida $\sigma^2$ .

Obsérvese que los parámetros (intercepción, pendiente) de las dos líneas son $(\beta_0, \beta_1)$ cuando $x_2\ge \tau$ y $(\beta_0+\beta_2, \beta_1+\beta_3)$ de lo contrario.

Una estimación de $\tau$ es un valor (que normalmente no será único) que minimiza la suma de los residuos OLS al cuadrado. Esto equivale a la solución de máxima verosimilitud. Para mostrar lo difícil que puede ser el problema de encontrar esta estimación, he representado la suma de los residuos al cuadrado frente a valores de prueba de $\tau$ . El gráfico de la izquierda muestra un ejemplo de este perfil. El valor óptimo de $\tau$ está marcada con una línea vertical de puntos rojos. Su patrón irregular, localmente constante y no diferenciable hace casi imposible que cualquier optimizador de propósito general encuentre el mínimo de forma fiable. (Podría ser fácilmente atrapado en un mínimo local lejos del mínimo global. He aplicado optimize a este problema como una comprobación, y en algunos ejemplos eso es exactamente lo que le ocurrió).

Dada esta estimación de $\tau$ el modelo es lineal y puede ajustarse mediante MCO. Este ajuste se muestra en el gráfico de la derecha. El modelo real consta de dos líneas de pendiente $\pm 1$ cruzando en $(0,1)$ y el umbral verdadero es $-1/3$ . En este conjunto de datos, $x_1$ y $x_2$ estaban casi sin correlación. Una fuerte correlación entre ellos hará que las estimaciones no sean fiables.

Los cuadrados naranjas y los círculos azules distinguen los casos $x_2 \lt \tau$ y $x_2 \ge \tau$ respectivamente. Los dos subconjuntos de datos se ajustan por separado con líneas rectas.

Se puede obtener una solución completa de máxima verosimilitud sustituyendo la suma de los residuos al cuadrado por el logaritmo negativo de la verosimilitud y aplicando los métodos estándar de MLE para obtener regiones de confianza para los parámetros y para comprobar las hipótesis sobre ellos.

---

#
# Create a dataset with a specified model.
#
n <- 80
beta <- c(1,1,0,-2)
threshold <- -1/3
sigma <- 1
x1 <- seq(1-n,n-1,2)/n
set.seed(17)
x2 <- rnorm(n)
i <- x2 < threshold
x <- cbind(1, x1, i, x1*i)
y <- x %*% beta + rnorm(n, sd=sigma)
#
# Display the SSR profile for the threshold.
#
f <- function(threshold) lm(y ~ x1*I(x2 < threshold))
z <- seq(-1,1,0.5/n)*2*sd(x2)                   # Search range
w <- sapply(z, function(z) sum(resid(f(z))^2))  # Sum of squares of residuals

par(mfrow=c(1,2))
plot(z, w, lwd=2, type="l", xlab="Threshold", ylab="SSR",main="Profile")
t.opt <- (z[which.min(w)] + z[length(z)+1 - which.min(rev(w))])/2
abline(v=t.opt, lty=3, lwd=3, col="Red")
#
# Fit the model to the data.
#
fit <- f(t.opt)
#
# Report and display the fit.
#
summary(fit)
plot(x1, y, pch=ifelse(x2 < t.opt, 21, 22), bg=ifelse(x2 < t.opt, "Blue", "Orange"),
     main="Data and Fit")
b <- coef(fit)
abline(b[c(1,2)], col="Orange", lwd=3, lty=1)
abline(b[c(1,2)] + b[c(3,4)], col="Blue", lwd=3, lty=1)