Como ajustar uma regressão como

9

Eu tenho alguns dados de séries temporais em que a variável medida é números inteiros positivos discretos (contagens). Quero testar se há uma tendência ascendente ao longo do tempo (ou não). A variável independente (x) está no intervalo de 0 a 500 e a variável dependente (y) está no intervalo de 0 a 8.

Eu pensei que respondesse isso ajustando uma regressão da forma y = floor(a*x + b)usando mínimos quadrados ordinários (OLS).

Como eu faria isso usando R (ou Python)? Existe um pacote existente para ele ou é melhor escrever meu próprio algoritmo?

PS: Eu sei que essa não é a técnica ideal, mas preciso fazer uma análise relativamente simples que eu possa realmente entender - minha formação é biologia, não matemática. Sei que estou violando suposições sobre erro na variável medida e independência das medidas ao longo do tempo.

r regression python afaulconbridge
fonte

5

Embora seja matematicamente natural tentar uma regressão dessa forma, por trás dela há um erro estatístico: o termo do erro agora estará fortemente correlacionado com o valor previsto. Essa é uma violação bastante forte das suposições do OLS. Em vez disso, use uma técnica baseada em contagem, conforme sugerido pela resposta de Greg Snow. (De bom grado votei esta questão, no entanto, porque ela reflete algum pensamento e esperteza reais. Obrigado por perguntar aqui!)

whuber

11

Você pode ajustar o modelo que declara usando a função nls(mínimos quadrados não lineares) R, mas como você disse que violará muitas das suposições e ainda provavelmente não fará muito sentido (você está dizendo que o resultado previsto é aleatório em torno de uma etapa função, não valores inteiros em torno de um relacionamento que aumenta suavemente).

A maneira mais comum de ajustar os dados de contagem é usar a regressão de Poisson usando a glmfunção in R, o primeiro exemplo na página de ajuda é uma regressão de Poisson, embora se você não estiver familiarizado com as estatísticas, seria melhor consultar um estatístico para garantir que você está fazendo as coisas corretamente.

Se o valor de 8 for um máximo absoluto (impossível de ver uma contagem mais alta, não é apenas isso que você viu), você pode considerar a regressão logística de probabilidades proporcionais, existem algumas ferramentas para fazer isso em pacotes R, mas você realmente deve envolver um estatístico se você quiser fazer isso.

Greg Snow
fonte

"você está dizendo que o resultado previsto é aleatório em torno de uma função step, não valores inteiros em um relacionamento que aumenta suavemente" --- Isso é algo que eu não havia considerado. No final, eu fui com a regressão de Poisson por glm. Não é a escolha perfeita, mas "suficientemente boa" para o que eu precisava.

precisa saber é o seguinte

10

$\def\lf{\lfloor}\def\rf{\rfloor}\def\pnorm{\mathrm{pnorm}}$ É claro que a sugestão de Greg é a primeira coisa a tentar: a regressão de Poisson é o modelo natural em muitos aspectos concretos. situações.

No entanto, o modelo que você está sugerindo pode ocorrer, por exemplo, quando você observa dados arredondados: com erros normais do iid .

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

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

ϵ_{i}

$\epsilon_i$

Eu acho que é interessante dar uma olhada no que pode ser feito com isso. Denoto por o cdf da variável normal padrão. Se , então usando notações familiares do computador. $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*}$

Você observa pontos de dados . A probabilidade do log é dada por Isso não é idêntico aos mínimos quadrados. Você pode tentar maximizar isso com um método numérico. Aqui está uma ilustração em 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 linear arredondado

Em vermelho e azul, as linhas encontradas por maximização numérica dessa probabilidade e mínimos quadrados, respectivamente. A escada verde é para encontrada a partir da probabilidade máxima ... isso sugere que você poderia usar mínimos quadrados, até uma tradução de por 0,5 e obter aproximadamente o mesmo resultado; ou, esses mínimos quadrados se encaixam bem no modelo onde é o número inteiro mais próximo. Os dados arredondados são tão frequentemente encontrados que tenho certeza de que isso é conhecido e que foi estudado extensivamente ... $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
fonte

4

+1 Adoro essa técnica e, na verdade, enviei um artigo para uma revista de análise de risco alguns anos atrás. (Alguns analistas de risco estão bastante interessados em dados com intervalo). Foi rejeitado como sendo "matemático demais" para o público. :-(. Uma dica: ao usar métodos numéricos, é sempre uma boa idéia fornecer bons valores iniciais para a solução. Considere aplicar OLS aos dados brutos para obter esses valores e, em seguida, "aperfeiçoá-los" com o otimizador numérico.

whuber

Sim, esta é uma boa sugestão. De fato, nesse caso, escolho valores remotos para enfatizar que "funciona", mas, na prática, sua sugestão seria a única solução para evitar começar de uma região muito plana, dependendo dos dados ...

Elvis

Como ajustar uma regressão como

Respostas: