¿Cómo ajustar una regresión como en R?

9

Tengo algunos datos de series temporales en los que la variable medida es enteros positivos discretos (recuentos). Quiero probar si hay una tendencia al alza con el tiempo (o no). La variable independiente (x) está en el rango 0-500 y la variable dependiente (y) está en el rango 0-8.

Pensé que respondía esto ajustando una regresión de la forma y = floor(a*x + b)usando mínimos cuadrados ordinarios (MCO).

¿Cómo haría para hacer esto usando R (o Python)? ¿Existe un paquete para él o es mejor que escriba mi propio algoritmo?

PD: Sé que esta no es la técnica ideal, pero necesito hacer un análisis relativamente simple que realmente pueda entender: mi experiencia es la biología, no las matemáticas. Sé que estoy violando supuestos sobre el error en la variable medida y la independencia de las mediciones con el tiempo.

r regression python afaulconbridge
fuente

55

Aunque es matemáticamente natural intentar una regresión de esta forma, detrás de él se esconde un error estadístico: el término de error ahora estará fuertemente correlacionado con el valor predicho. Esa es una violación bastante fuerte de los supuestos de OLS. En su lugar, utilice una técnica basada en el conteo como lo sugiere la respuesta de Greg Snow. (Con mucho gusto upvoted esta pregunta, sin embargo, debido a que refleja un pensamiento real y la inteligencia Gracias por preguntar aquí.!)

whuber

11

Podrías ajustar el modelo en el que declaras usando la función nls(mínimos cuadrados no lineales) R, pero como dijiste, eso violará muchos de los supuestos y probablemente aún no tenga mucho sentido (estás diciendo que el resultado predicho es aleatorio en torno a un paso función, no valores enteros alrededor de una relación que aumenta sin problemas).

La forma más común de ajustar los datos de conteo es usar la regresión de Poisson usando la glmfunción R, el primer ejemplo en la página de ayuda es una regresión de Poisson, aunque si no está tan familiarizado con las estadísticas, sería mejor consultar con un estadístico para asegurarse que estás haciendo las cosas correctamente

Si el valor de 8 es un máximo absoluto (imposible ver un conteo más alto, no solo eso es lo que vio), entonces puede considerar la regresión logística de probabilidades proporcionales, hay un par de herramientas para hacer esto en paquetes R, pero usted Realmente debería involucrar a un estadístico si desea hacer esto.

Greg Snow
fuente

"Está diciendo que el resultado predicho es aleatorio en torno a una función de paso, no valores enteros en torno a una relación que aumenta sin problemas" --- Eso es algo que no había considerado. Al final, fui con la regresión de Poisson por glm. No es la elección perfecta, pero "lo suficientemente bueno" para lo que necesitaba.

afaulconbridge

10

$\def\lf{\lfloor}\def\rf{\rfloor}\def\pnorm{\mathrm{pnorm}}$ Está claro que la sugerencia de Greg es lo primero que hay que intentar: la regresión de Poisson es el modelo natural en muchos, muchos concretos situaciones

Sin embargo, el modelo que está sugiriendo puede ocurrir, por ejemplo, cuando observa datos redondeados: con errores normales de id .

Y_{i} = ⌊ a x_{i} + b + ϵ_{i} ⌋,

$Y_i = \lf ax_i + b + \epsilon_i \rf,$

ϵ_{i}

$\epsilon_i$

Creo que es interesante ver qué se puede hacer con él. Denoto por el cdf de la variable normal estándar. Si , entonces usando notaciones familiares de computadora. $F$ $\epsilon \sim \mathcal N(0,\sigma^2)$

\begin{aligned} P (⌊ a x + b + ϵ ⌋ = k) & = F (\frac{k - b + 1 - a x}{σ}) - F (\frac{k - b - a x}{σ}) \\ = p n o r m (k + 1 - a x - b, s d = σ) - p n o r m (k - a x - b, s d = σ), \end{aligned}

$\begin{align*} \mathbb P\left(\lf ax + b + \epsilon \rf = k\right) &= F\left({k-b+1-ax\over \sigma}\right) - F\left({k-b-ax\over \sigma}\right)\\ &= \pnorm(k+1-ax-b,sd=\sigma) - \pnorm(k-ax-b,sd=\sigma),\end{align*}$

puntos de datos . La probabilidad de registro viene dada por Esto no es idéntico a los mínimos cuadrados. Puede intentar maximizar esto con un método numérico. Aquí hay una ilustración en R: $(x_i,y_i)$

ℓ (a, b, σ) = \sum_{i} \log (F (\frac{y_{i} - b + 1 - a x_{i}}{σ}) - F (\frac{y_{i} - b - a x_{i}}{σ})) .

$\ell(a,b,\sigma) = \sum_i \log\left( F\left({y_i-b+1-ax_i\over \sigma}\right) - F\left({y_i-b-ax_i\over \sigma}\right) \right).$

log_lik <- function(a,b,s,x,y)
  sum(log(pnorm(y+1-a*x-b, sd=s) - pnorm(y-a*x-b, sd=s)));

x <- 0:20
y <- floor(x+3+rnorm(length(x), sd=3))
plot(x,y, pch=19)
optim(c(1,1,1), function(p) -log_lik(p[1], p[2], p[3], x, y)) -> r
abline(r$par[2], r$par[1], lty=2, col="red")
t <- seq(0,20,by=0.01)
lines(t, floor( r$par[1]*t+r$par[2]), col="green")

lm(y~x) -> r1
abline(r1, lty=2, col="blue");

modelo lineal redondeado

En rojo y azul, las líneas encuentran por maximización numérica de esta probabilidad y mínimos cuadrados, respectivamente. La escalera verde es para encontrada desde la probabilidad máxima ... esto sugiere que podría usar mínimos cuadrados, hasta una traducción de por 0.5, y obtener aproximadamente el mismo resultado; o, esos mínimos cuadrados se ajustan bien al modelo donde es el entero más cercano. Los datos redondeados se cumplen con tanta frecuencia que estoy seguro de que esto se conoce y se ha estudiado ampliamente ... $ax+b$ $\lf ax +b\rf$ $a,b$ $b$

Y_{i} = [a x_{i} + b + ϵ_{i}],

$Y_i = [ a x_i + b +\epsilon_i],$

[x] = ⌊ x + 0.5 ⌋

$[x] = \lf x + 0.5 \rf$

Elvis
fuente

44

+1 Me encanta esta técnica y en realidad envié un artículo sobre ella a una revista de análisis de riesgos hace unos años. (Algunos analistas de riesgo están bastante interesados en los datos con valor de intervalo). Se rechazó por ser "demasiado matemático" para su audiencia. :-(. Un consejo: cuando se utilizan métodos numéricos, siempre es una buena idea proporcionar buenos valores iniciales para la solución. Considere aplicar OLS a los datos sin procesar para obtener esos valores, luego "pulirlos" con el optimizador numérico.

whuber

Sí, esta es una buena sugerencia. De hecho, en ese caso elijo valores remotos para enfatizar que "funciona", pero en la práctica su sugerencia sería la única solución para evitar comenzar desde una región muy plana, dependiendo de los datos ...

Elvis

¿Cómo ajustar una regresión como en R?

Respuestas: