Calculando a tendência de um ângulo quando ele cruza 360 -> 0
Tenho uma variável que mede um ângulo que descreve a posição relativa de dois objetos (ou seja, pode variar de 0 a 359) e gostaria de quantificar como isso mudou ao longo do tempo.
Por exemplo, aqui temos a posição relativa dos dois itens variando 1 grau por ano:
year <- seq(1981, 2020)
angle <- c(seq(341, 359), seq(0, 20))

No entanto, tomar a inclinação aqui não faz sentido por causa do "cruzamento" que ocorre no ano 2000. Tenho várias amostras diferentes, algumas com esse problema e outras não. Não sei a priori quais amostras terão esse problema, nem quando o crossover ocorrerá, então não posso simplesmente aplicar algum tipo de compensação (ou seja, adicionar 360 aos últimos 20 anos).
Existe uma maneira aceita de calcular tendências angulares, considerando o fato de que 0 = 360?
Respostas
Pense no ângulo$y$a qualquer momento$t$como o acúmulo de pequenas mudanças no ângulo. Simbolicamente, quando$f(t)$é a taxa de variação do ângulo no tempo$t$e$t_0$é o início das observações,
$$y(t) = y(t_0) + \int_{t_0}^t f(t)\,\mathrm{d}t.$$
seu problema é esse$y(t)$foi gravado modulo$360$graus - talvez com algum erro$\epsilon(t).$Ou seja, você observou apenas os valores
$$y^{*}(t) = y(t) + \epsilon(t) \mod 360.$$
Você pode, no entanto, reconstruir$y(t) + \epsilon(t)$desde que você tenha observações suficientemente frequentes. Por tempos sucessivos$t \lt s,$perceber
$$y^{*}(s) - y^{*}(t) = y(s) - y(t) + \epsilon(s) - \epsilon(t) \mod 360 = \int_t^s f(t)\,\mathrm{d}t + \delta$$
Onde$\delta$é igual à contribuição dos erros$\epsilon(s)-\epsilon(t)$ mais, talvez, algum múltiplo inteiro de$360$sempre que houver uma quebra angular entre$y^{*}(t)$e$y^{*}(s).$Agora, desde que o tamanho do erro total$|\epsilon(s)-\epsilon(t)|$é menos do que$180$graus e desde que o ângulo não gire mais de uma vez, podemos descobrir se ocorreu uma quebra: se$|\epsilon(s)-\epsilon(t)| \gt 180,$adicionar ou subtrair$360$graus de$\delta$colocá-lo no intervalo de$-180$para$+180$graus.
Embora não possamos observar esses erros diretamente, se estivermos amostrando com frequência suficiente para fazer os incrementos$y(t_i) - y(t_{i-1})$bastante pequeno, simplesmente aplicamos esse ajuste às diferenças observadas. Desta forma,
Em qualquer momento$|y^{*}(s)-y^{*}(t)| \gt 180,$adicionar ou subtrair$360$graus de$\delta$colocá-lo no intervalo de$-180$para$+180$graus.
Equivalentemente, calcule as diferenças módulo$180$mas expressá-los no intervalo de$-180$para$+180$graus em vez de (como é convencional) a faixa de$0$para$360.$
Vamos chamar o valor ajustado$\delta^{*}(t,s),$de modo a
$$y^{*}(s) - y^{*}(t) = \int_t^s f(t)\,\mathrm{d}t + \delta(t,s)^{*}.$$
Isso é igualdade, não igualdade módulo$360.$ Podemos agora remover o efeito de gravar os ângulos módulo$360$somando essas diferenças ajustadas. Quando as observações são feitas às vezes$t_0 \lt t_1\lt \cdots \lt t_n,$temos
$$\begin{aligned} y^{*}(t_i) &= y^{*}(t_0) + \left[y^{*}(t_1) - y^{*}(t_0)\right] + \cdots + \left[y^{*}(t_i) - y^{*}(t_{i-1})\right] \\ &=y(t_0) + \int_{t_0}^{t_i} f(t)\,\mathrm{d}t + \delta(t_0,t_1)^{*} + \delta(t_1,t_2)^{*} + \cdots + \delta(t_{i-1},t_i)^{*} \\ &= y(t_i) + \left[\epsilon(t_i) - \epsilon(t_0)\right]. \end{aligned}$$
O problema com módulo de computação$360$se foi: agora você pode usar qualquer procedimento que desejar para modelar a resposta$y^{*}(t).$
Aqui está uma ilustração com um conjunto de dados bastante difícil. Os dados foram gerados de acordo com o modelo$y(t) = 30t \mod 360$e observado anualmente de 1980 a 2020 com iid Erro normalmente distribuído do desvio padrão$60$graus (uma grande quantidade).

A tendência é quase imperceptível nos dados brutos, mas o algoritmo de ajuste de ângulo os alinhou visivelmente. Podemos ajustar um modelo de mínimos quadrados aos dados ajustados, por exemplo, produzindo este resultado:

A escala vertical expandida para os dados brutos mostra detalhes do ajuste e seus desvios. Aliás, neste exemplo, a estimativa da inclinação é$28.0 \pm 0.74$graus, não notavelmente diferente do verdadeiro valor de$30$graus (o valor-p para esta comparação é$1.1\%$).
Terminarei observando que quando o desvio padrão dos erros$\epsilon(t)$é grande (maior que$180/2/\sqrt{2} \approx 64$graus, aproximadamente), às vezes o ajuste angular ficará incorreto. Isso aparecerá nos resíduos do modelo como uma mudança repentina em um valor em torno de 360 graus. Assim, uma análise rotineira dos resíduos do modelo pode detectar tais problemas, permitindo modificar os valores ajustados para um melhor ajuste. Os detalhes disso dependerão do seu modelo e do procedimento de adaptação.
Este R
código criou as figuras. Em "ajustar os ângulos" mostra como o ajuste do ângulo pode ser calculado de forma eficiente.
#
# Specify the data-generation process.
#
year <- 1980:2020 # Dates to use
beta <- 30 # Annual rate of change
sigma <- 60 # Error S.D.
#
# Generate the data.
#
set.seed(17)
angle <- (year * beta + rnorm(length(year), 0, sigma)) %% 360
X <- data.frame(year, angle)
#
# Adjust the angles.
#
X$`total angle` <- with(X, {
d <- (diff(angle) + 180) %% 360 - 180
cumsum(c(angle[1], d))
})
#
# Fit a model to the adjusted angles.
#
fit <- lm(`total angle` ~ year, X)
#
# Analyze the fit.
#
b <- coefficients(fit)
y.hat <- predict(fit)
#--Compute dates the fit must wrap around from 360 to 0:
y.breaks <- seq(floor(min(y.hat) / 360)*360, max(y.hat), by=360)
year.breaks <- (y.breaks - b[1]) / b[2]
#--Make the plots:
u <- ceiling(max(X$`total angle`)/360)
par(mfcol=c(1,2))
#--The fits:
plot(X$year, X$angle, pch=19, ylim=c(0, 360), yaxp=c(0, 360, 4),
col="gray", ylab="Angle (degrees)", xlab="Year",
main="Raw Data and Fit")
for (x in year.breaks)
abline(c(-x * b[2], b[2]), col="Red", lwd=2)
plot(X$year, X$`total angle`, ylim=c(0,u*360), yaxp=c(0, u*360, u),
xlab="Year", ylab="Total angle",
main="Adjusted Data and Fit")
abline(fit, col="Red", lwd=2)
#--The raw data:
plot(X$year, X$angle, ylim=c(0,u*360), yaxp=c(0, u*360, u),
pch=19, col="gray", ylab="Angle (degrees)", xlab="Year",
main="Raw Data")
plot(X$year, X$`total angle`, ylim=c(0,u*360),
yaxp=c(0, u*360, u),
xlab="Year", ylab="Total angle",
main="Adjusted Data")
par(mfcol=c(1,1))