Jak dopasować regresję, taką jak in R?

9

Mam pewne dane szeregów czasowych, w których mierzoną zmienną są dyskretne dodatnie liczby całkowite (liczby). Chcę sprawdzić, czy z czasem (lub nie) występuje trend wzrostowy. Zmienna niezależna (x) jest w zakresie 0-500, a zmienna zależna (y) jest w zakresie 0-8.

Myślałem, że odpowiem na to, dopasowując regresję formy y = floor(a*x + b)za pomocą zwykłych najmniejszych kwadratów (OLS).

Jak mógłbym to zrobić za pomocą R (lub Python)? Czy istnieje już dla niego pakiet, czy może lepiej napisać własny algorytm?

PS: Wiem, że to nie jest idealna technika, ale muszę przeprowadzić stosunkowo prostą analizę, którą właściwie potrafię zrozumieć - moje tło to biologia, a nie matematyka. Wiem, że naruszam założenia dotyczące błędu mierzonej zmiennej i niezależności pomiarów w czasie.

r regression python

— afaulconbridge
źródło

5

Chociaż matematycznie naturalne jest wypróbowanie regresji tej formy, kryje się za nią błąd statystyczny: termin błędu będzie teraz silnie skorelowany z przewidywaną wartością. To dość silne naruszenie założeń OLS. Zamiast tego użyj techniki liczenia, jak sugeruje odpowiedź Grega Snowa. (Z radością jednak głosowałem za tym pytaniem, ponieważ odzwierciedla ono prawdziwą myśl i spryt. Dziękuję, że

— zadałeś

11

Możesz dopasować model, który podajesz za pomocą funkcji nls(nieliniowej najmniejszych kwadratów) R, ale jak powiedziałeś, to naruszy wiele założeń i prawdopodobnie nie będzie miało większego sensu (mówisz, że przewidywany wynik jest losowy na etapie funkcja, a nie wartości całkowite wokół płynnie rosnącego związku).

Najczęstszym sposobem dopasowania danych zliczania jest regresja Poissona przy użyciu glmfunkcji in R. Pierwszy przykład na stronie pomocy to regresja Poissona, chociaż jeśli nie znasz się na statystykach, najlepiej skonsultować się ze statystykami, aby upewnić się że robisz wszystko poprawnie.

Jeśli wartość 8 jest absolutnym maksimum (niemożliwym do zobaczenia wyższego wyniku, nie tylko to widziałeś), możesz rozważyć regresję logistyczną proporcjonalności szans, istnieje kilka narzędzi, aby to zrobić w pakietach R, ale możesz naprawdę powinieneś zaangażować statystyk, jeśli chcesz to zrobić.

— Greg Snow
źródło

„mówisz, że przewidywany wynik jest losowy wokół funkcji krokowej, a nie liczb całkowitych wokół płynnie rosnącego związku” - tego nie wziąłem pod uwagę. W końcu poszedłem z regresją Poissona przez glm. To nie jest idealny wybór, ale „wystarczająco dobry” do tego, czego potrzebowałem.

— afaulconbridge

10

$\def\lf{\lfloor}\def\rf{\rfloor}\def\pnorm{\mathrm{pnorm}}$ Oczywiste jest, że sugestia Grega jest pierwszą rzeczą do wypróbowania: regresja Poissona jest naturalnym modelem w wielu konkretnych przypadkach sytuacje.

Jednak model, który sugerujesz, może wystąpić na przykład, gdy obserwujesz zaokrąglone dane: iid normalne błędy .

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

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

ϵ_{i}

$\epsilon_i$

Myślę, że to interesujące, aby zobaczyć, co można z tym zrobić. Oznaczam przez cdf standardowej zmiennej normalnej. Jeśli , to przy użyciu znanych notacji komputerowych. $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*}$

Obserwujesz punkty danych . Prawdopodobieństwo dziennika jest podane przez Nie jest to identyczne z najmniejszymi kwadratami. Możesz spróbować zmaksymalizować to za pomocą metody numerycznej. Oto ilustracja w 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");

zaokrąglony model liniowy

Na czerwono i niebiesko linie znalezione przez numeryczną maksymalizację tego prawdopodobieństwa i odpowiednio najmniejszych kwadratów. Zielone schody to dla znalezione z maksymalnego prawdopodobieństwa ... to sugeruje, że możesz użyć najmniejszych kwadratów, do tłumaczenia o 0,5, i uzyskać mniej więcej ten sam wynik; lub te najmniejsze kwadraty dobrze pasują do modelu gdzie jest najbliższą liczbą całkowitą. Zaokrąglone dane są tak często spotykane, że jestem pewien, że jest to znane i zostało gruntownie zbadane ... $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
źródło

4

+1 Uwielbiam tę technikę i rzeczywiście przedłożyłem artykuł na ten temat do dziennika analizy ryzyka kilka lat temu. (Niektórzy analitycy ryzyka są bardzo zainteresowani danymi o wartościach interwałowych). Został on odrzucony jako „zbyt matematyczny” dla ich odbiorców. :-(. Jedna wskazówka: przy stosowaniu metod numerycznych zawsze dobrym pomysłem jest podanie dobrych wartości początkowych dla rozwiązania. Rozważ zastosowanie OLS do surowych danych w celu uzyskania tych wartości, a następnie „

— wypoleruj

Tak, to dobra sugestia. Właściwie w takim przypadku wybieram wartości zdalne, aby podkreślić, że „to działa”, ale w praktyce Twoja sugestia byłaby jedynym rozwiązaniem, aby uniknąć rozpoczynania od bardzo płaskiego regionu, w zależności od danych ...

— Elvis