Estimador de pendiente robusto a valores repetidos

Question

Tengo un conjunto de datos hidrológicos que contiene muchos valores repetidos (en mi caso, 0s), y quiero entender si ha habido una tendencia a través del tiempo. He aquí un ejemplo idealizado del aspecto de los datos:

# make sample data with linear increase
year <- seq(1, 51)
value <- seq(0, 500, 10)

# replace just over 50% of values with 0s
value[seq(1,51,2)] <- 0

plot(year, value)

Normalmente, utilizaría el método no paramétrico Prueba de Mann-Kendall para determinar si existe un cambio significativo a lo largo del tiempo, y la Estimador de la pendiente de Theil-Sen para determinar la pendiente de ese cambio.

En este caso, la prueba de Mann-Kendall rechaza la hipótesis nula de ausencia de cambios a lo largo del tiempo y tiene una tau positiva, interpretada como un aumento estadísticamente significativo a lo largo del tiempo. Sin embargo, el estimador de Theil-Sen arroja un valor de 0:

manken <- rkt::rkt(year, value)
manken$tau  # Kendall tau = 0.235
manken$sl   # p-value = 0.009
manken$B    # Theil-Sen estimator = 0

Creo que esto ocurre porque el estimador Theil-Sen devuelve la mediana de la pendiente de todos los pares. Por lo tanto, siempre que haya >50% de puntos que tengan el mismo valor, la pendiente mediana (creo) siempre será igual a 0.

Pregunta: ¿Existen estimadores estadísticos de la pendiente más adecuados para datos con muchos valores repetidos?

Por si sirve de algo, un modelo lineal ( lm(value ~ year) ) también arroja una pendiente positiva significativa (p=0,001). Y lo mismo ocurre con la muy científica "prueba del globo ocular".

Edición 8/3/2020: Para mayor contexto, mi eje y real es el número de días con caudal cero al año. Así pues, este problema surge cuando hay un arroyo que fluye durante todo el año durante >50% de los años, pero que se seca durante parte o la totalidad de algunos años. Me gustaría saber si la "sequedad" (frecuencia/duración) está aumentando.

Answer 1

Para mayor contexto, mi eje y real es el número de días con cero flujo por año.

Sus datos simulados no contienen incertidumbre y, por tanto, no son muy útiles. Además, ¿cómo puede tener un año más de 365 días? Voy a simular su dependiente como una variable de recuento, es decir, con una distribución de Poisson.

No soy un experto en modelos de recuento y no he visto sus datos reales, por lo que otras distribuciones (como la binomial negativa) podrían ser mejores para su modelo. Incluso podría ser necesario utilizar una distribución con un límite superior (si tiene valores cercanos a 365 días).

Utilizaré una distribución binomial para simular sus valores cero. De nuevo, otras distribuciones podrían representar mejor tus datos.

set.seed(42)

year <- seq(1, 51)
value <- rpois(length(year), lambda = exp(year * 0.07 + 1))

# replace about 50% of values with 0s
value[as.logical(rbinom(length(value), 1, 0.5))] <- 0
mean(value == 0)
#[1] 0.5686275

DF <- data.frame(year, value)

plot(value ~ year, data = DF)

Ahora podemos ajustar un modelo de obstáculos. Un modelo de obstáculos combina dos modelos. El primero modela si los valores son cero o distintos de cero. El segundo modela los valores distintos de cero. Ambos modelos son modelos lineales generalizados .

library(pscl)
fit <- hurdle(value ~ year, dist = "poisson", zero.dist = "binomial", data = DF)

summary(fit)
#Call:
#hurdle(formula = value ~ year, dist = "poisson", zero.dist = "binomial")
#
#Pearson residuals:
#    Min      1Q  Median      3Q     Max 
#-0.7026 -0.6698 -0.6171  1.1072  2.0128 
#
#Count model coefficients (truncated poisson with log link):
#            Estimate Std. Error z value Pr(>|z|)    
#(Intercept) 1.458876   0.151480   9.631   <2e-16 ***
#year        0.058816   0.003802  15.471   <2e-16 ***
#Zero hurdle model coefficients (binomial with logit link):
#             Estimate Std. Error z value Pr(>|z|)
#(Intercept) -0.881742   0.618411  -1.426    0.154
#year         0.003785   0.020518   0.184    0.854
#---
#Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
#
#Number of iterations in BFGS optimization: 11 
#Log-likelihood: -73.36 on 4 Df

Como ves, la salida nos dice que la probabilidad de que un valor sea cero es independiente del año (así lo hemos simulado). El modelo de recuento muestra una intercepción y una pendiente fuertemente significativas (nótese el enlace logarítmico). Vamos a trazar las predicciones del modelo de recuento:

curve(predict(fit, type = "count", newdata = data.frame(year = x)), add = TRUE, col = "red")

Creo que los modelos de obstáculos podrían ayudarte, pero tendrías que investigar un poco más qué supuestos serían sensatos respecto a las distribuciones y las funciones de enlace. Por supuesto, para ello sería útil tener conocimientos mecanicistas sobre por qué se producen valores distintos de cero y qué podría causar el aumento con el tiempo. Sería útil disponer de predictores adicionales.

Me gustaría saber si la "sequedad" (frecuencia/duración) está aumentando.

El modelo cero te diría si la frecuencia de años con caudal cero depende del tiempo (en los datos simulados no es así). El modelo de recuento indicaría si el número de días con caudal cero en años secos ("gravedad" de la sequía) depende del tiempo.

Tenga en cuenta que los modelos hurdle son para datos con inflación cero, suponen que intervienen dos "procesos". Uno controla si un valor es distinto de cero, el otro la magnitud de los valores distintos de cero. Sus datos simulados apoyan esta suposición. Sus datos reales podrían no estar inflados a cero.

Answer 2

El estimador Theil-Sen es robusto frente a una distribución de errores sesgada o con colas gruesas, pero sigue suponiendo que el modelo sigue siendo lineal y tiene una pendiente única y bien definida. Es decir, asume que la especificación sigue siendo de la forma

$$ y = \beta_1 x + \beta_0 + \epsilon $$

Pero en lugar de asumir $\epsilon \sim \mathcal{N}(0, \sigma^2)$ no impone ninguna restricción a la distribución de $\epsilon$ excepto quizás la media 0. Según Wikipedia, es robusto hasta el 29% de los puntos que se cambian.

Sin embargo, la forma en que generó los datos (denominada modelo mixto ) no supone que haya una única pendiente, sino que hay dos clases, cada una con una pendiente diferente, que se han mezclado aleatoriamente.

Kendall $\rm{T}$ sigue funcionando, porque está claro que la mezcla no es un orden aleatorio, pero la estimación de la pendiente no es aplicable porque se han violado gravemente los supuestos.

Existe una forma estándar de ajustar modelos a los datos generados por modelos de mezclas: los modelos de variables latentes y los modelos de Algoritmo de maximización de expectativas . Sus datos podrían llamarse una "mezcla de regresiones". Hay un paquete de R que puede manejar este caso llamado flexmix . Así es como yo usaría ese paquete para ajustar tus datos falsos. El k=2 le está diciendo que hay dos clases, que sabemos que a priori.

# generate 51 equally spaced points along a line
year <- seq(1, 51)
value <- seq(0, 500, 10)

# add a little bit of noise to prevent likelihood underflow
value <- value + rnorm(n=51, mean=0, sd=1)

# replace just over 50% of values with 0s
value[seq(1,51,2)] <- 0

#install.packages("flexmix")
library(flexmix)

1model <- flexmix(value ~ year, k=2)
summary(model)

plot(year, value, col = clusters(model), pch=19)
abline(parameters(model)[1:2, 1], col = "black", lty=2)
abline(parameters(model)[1:2, 2], col = "red", lty=2)

Call:
flexmix(formula = value ~ year, k = 2)

       prior size post>0 ratio
Comp.1   0.5   26     26 1.000
Comp.2   0.5   25     26 0.962

'log Lik.' -118.8863 (df=7)
AIC: 251.7726   BIC: 265.2954

El algoritmo EM funciona adivinando a qué clase pertenece cada punto. Empieza suponiendo que cada punto tiene una probabilidad del 50% de pertenecer a cada clase. A continuación, ajusta un ponderado modelo de regresión para cada clase. A continuación, basándose en los dos modelos de regresión ajustados, vuelve atrás y actualiza las probabilidades de estar en cada clase para cada punto. Por ejemplo, si inicialmente se suponía que un punto tenía las mismas probabilidades de pertenecer a una u otra clase, pero después de la primera iteración queda muy cerca de la línea de regresión de la clase 1 y muy lejos de la línea de regresión de la clase 2, sus probabilidades se actualizarán al 80% para la clase 1 y al 20% para la clase 2. Este proceso se repite hasta que se alcanza la convergencia. En ese momento, tenemos una estimación bastante buena de la clase de la que procede cada punto, y dos líneas de regresión distintas; debido a la ponderación, podemos imaginar que cada línea se ajustó sólo a los puntos que probablemente pertenecen a la misma clase.

El algoritmo EM es bueno pero no perfecto. El hiperparámetro k debe elegirse con mucho cuidado. Aunque se garantiza que la probabilidad aumente con cada iteración, el algoritmo a veces puede ser inestable y converger a soluciones diferentes si se ajusta a una submuestra aleatoria diferente de los datos. En algunos casos, la verosimilitud puede llegar hasta el infinito; esto ocurre con el conjunto de datos falsos, ¡porque todos los datos se encuentran en una línea perfectamente recta! (Añadir un poco de ruido aleatorio soluciona ese problema, que de todas formas es muy poco probable que ocurra en los datos del mundo real). Sin embargo, si se cumplen los supuestos, puede ser una técnica muy potente.