R - cómo obtener el error estándar para un punto de ruptura/parámetro en una regresión por partes

Question

Soy nuevo en esto de hacer preguntas, así que espero que este sea un lugar razonable para esto.

Estoy tratando de encajar una regresión por partes. Espero que mi respuesta (y) aumente desde x0 hasta algún punto a lo largo de x, luego hasta la meseta después de un punto de ruptura; a partir de ahí, pendiente =0. El modelo que me gustaría ajustar tiene 3 parámetros: intercepción, pendiente y punto de ruptura. En el mejor de los casos, compararía este modelo con un modelo sólo de intercepción/nulo usando AIC o algún método similar.

Empecé con el paquete segmentado, con esta pista para fijar la segunda pendiente a 0: https://stat.ethz.ch/pipermail/r-help/2007-July/137625.html Pero después de actualizar un montón de paquetes, todo el asunto (que ya era un poco quisquilloso) dejó de funcionar por completo.

Así que, ahora, estoy usando optim(), pero no estoy seguro de cómo obtener apropiadamente una estimación del error, particularmente alrededor del parámetro del punto de ruptura. Vea el código de abajo para ejemplos de intentos usando bblme y nls (ambos funcionan, pero ambos parecen fallar en mis conjuntos de datos).

Según un colega usuario de SAS, proc nlin parece ser capaz de ajustar este modelo fácilmente y genera estimaciones de error estándar para los parámetros... pero, estoy decidido a no usar SAS.

Este intercambio (y otros) parece indicar que hay preocupaciones generales con respecto a las estimaciones de error para este tipo de modelo, en particular el parámetro del punto de ruptura - y hace referencia a un post potencialmente útil de Venables que lamentablemente ya no está disponible. Análisis del punto de cambio usando R's nls()

(Nota, mis datos están ponderados por una variable que representa el esfuerzo. Cuando está disponible, puede alimentar wt a una declaración de pesos=(), de lo contrario, utilizo un conjunto de datos expandido o peso la suma de los cuadrados en el modelo).

Simular datos segmentados con un tamaño de muestra decreciente

set.seed(15)
a.sim<-0 #  intercept
b.sim<-0.5 #  slope for segment 1
n<-12
brk<-4 # breakpoint
x <- c(1:n)
y <- numeric(n)
y[1:brk]<-a.sim+b.sim*x[1:brk]+rnorm(brk, 0, .2)
y[(brk+1):n] <- (a.sim+ b.sim* brk) + rnorm((n-brk), 0, .2) # second, flat segment
y[n]<-y[n]-.30*y[n] #subtract from last point, as these are often outside of the normal pattern
wt<-c(rep(50, n-4), c(40, 40, 35, 5)) #weight variable 

dat<-as.data.frame(cbind(x, y, wt)) # make dataframe 
dat.expand <- dat[ rep( seq(dim(dat)[1]), dat$wt),]# second dataset with rows repeated based on weight

with(dat, plot(x,y, ylim=c(0, max(y)), pch=16, cex=wt/(10)))### plot, with symbols representing weight variable

Usar optim para resolver

mod<-function(par, data){
  a <- par[1]
  b <-par[2]
  x.star <-par[3]
  yfit<-c(NA,length(data$y))
      small.x<-I(data$x<=x.star)
  yfit[small.x==T]<-(a+b*data$x[small.x==T]) 
      yfit[small.x!=T]<-(a+b*x.star) 
      sum((data$y-yfit)^2)
}
fit1<-optim(par=c(a=.5, b=.5, x.star=3), mod, data=dat, hessian=T)

esto casi siempre funciona bien y proporciona un ajuste razonable. Para el proceso de producción, le doy una gran tabla de valores potenciales de inicio, para encontrar un punto de partida razonable antes de optimizarlo realmente.

bblme para cis

Hablé con un estadístico, que recomendó esto como un posible método para generar intervalos de confianza alrededor del punto de ruptura. Funciona para mis datos en algunos casos, pero normalmente da errores.

library(bbmle)
mod2<-function(a,b,x.star,sig){
yfit<-c(NA,length(y))
small.x<-I(x<=x.star)
yfit[small.x==T]<-(a+b*x[small.x==T]) 
yfit[small.x!=T]<-(a+b*x.star) 
-sum(dnorm(y,yfit,1, log=TRUE))}

fit3<-mle2(mod2, start=list(a=0,b=0.5,x.star=brk), data=dat)
ci<-confint(fit3) 

método nls

Esto básicamente nunca funciona para mi tipo de datos (funciona si simulo conjuntos de datos mucho más grandes) Tengo el error de la matriz de gradiente singular o El error en .... factor de paso 0.000488281 reducido por debajo del 'minFactor' de 0.000976562...

nls.mod<-nls(y ~ ifelse(x <x.star, a+b*x,a+b*x.star), 
data = dat, weights=wt,
start = c(x.star=brk,b=0.4, a=0))

¿JAGS?

Así que, también tengo esto trabajando en JAGS/r2jags, de nuevo con la ayuda de un estadístico. Preferiría no usar este enfoque si hay una forma directa y frecuente de obtener una estimación razonable del error alrededor del punto de ruptura. No estoy buscando la perfección aquí. Sólo algo razonable

Answer 1

Hmm, así que creo que me di cuenta de los problemas que tenía con los segmentos. Tenía que ver con la declaración de peso (no funciona el peso del lm y el modelo segmentado).

Segmentar parece la mejor apuesta para mí. Hace un buen trabajo estimando los puntos de ruptura, incluso con mis cortos conjuntos de datos. No es terriblemente difícil restringir la segunda pendiente a 0, y tiene una prueba incorporada para la importancia del cambio de pendiente. Si alguien tiene ganas de explicar la dificultad de estimar el error alrededor de un punto de ruptura estimado, ¡soy todo oídos!

library(segmented)

con el segundo segmento no limitado a 0

out.lm <- lm(y~x, data=dat)
o<-segmented(out.lm, seg.Z= ~x, weights=wt, psi=10)
with(dat, plot(x,y, ylim=c(0, max(y)), pch=16, cex=wt/(13), main="segmented"))
lines(x=dat$x, y=fitted(o), col="blue")
lines.segmented(o, col="blue")

fijar desde vito para fijar la pendiente a 0

o2<-lm(y~1)
xx<- -x
o3<-segmented(o2,seg.Z=~xx, weights=wt, psi=list(xx=-4))

points(x,fitted(o3),col="green")
ci<-confint(o3, rev.sgn=TRUE)
lines.segmented(o3, col="green", rev.sgn=TRUE, lwd=2)

prueba de cambio de pendiente

davies.test(o3,~xx)

Answer 2

La distribución del estimador del punto de ruptura es complicada y no se pueden utilizar métodos estándar para ello. Afortunadamente, el paquete strucchange implementa pruebas de puntos de ruptura e intervalos de confianza (ver las referencias en strucchange::puntos de ruptura), que se pueden utilizar de manera muy sencilla:

Reescribí un proceso de generación de datos de una manera ligeramente más concisa (también asumiendo un cambio permanente en la intercepción y la pendiente en lugar de los cambios aleatorios específicos del tiempo como hiciste). Ten en cuenta que si tienes sólo 12 valores, no creo que consigas ninguna estimación relevante...

## DGP
set.seed(123)
n <- 100
breakpoint <- 50
x <- rnorm(n)
err <- rnorm(n, sd=0.1)
y <- 1.2+0.9*x +ifelse(1:n >breakpoint, 0.1+0.1*x,0)+err

library(strucchange)
b <- breakpoints(y~x, breaks=1) #selects 47
confint(b) # great, 50 is in the confidence interval!