1
$\begingroup$

I have a variable that measures an angle describing the relative position of two objects (i.e., can range from 0-359), and I would like to quantify how this has changed through time.

For example, here we have the relative position of the two items changing by 1 degree per year:

year <- seq(1981, 2020)
angle <- c(seq(341, 359), seq(0, 20))

a scatterplot showing the angle for each year

However, taking the slope here is meaningless because of the "crossover" that occurs in the year 2000. I have a number of different samples, and some have this issue and some don't. I don't know a priori which samples will have this issue, nor when the crossover occurs, so I can't just apply some sort of offset (i.e., add 360 to the last 20 years).

Is there an accepted way to calculate angular trends, accounting for the fact that 0 = 360?

$\endgroup$
0
1
$\begingroup$

Think of the angle $y$ at any time $t$ as the accumulation of small changes in the angle. Symbolically, when $f(t)$ is the rate of change of the angle at time $t$ and $t_0$ is the beginning of the observations,

$$y(t) = y(t_0) + \int_{t_0}^t f(t)\,\mathrm{d}t.$$

Your problem is that $y(t)$ has been recorded modulo $360$ degrees -- perhaps with some error $\epsilon(t).$ That is, you have observed only the values

$$y^{*}(t) = y(t) + \epsilon(t) \mod 360.$$

You can, however, reconstruct $y(t) + \epsilon(t)$ provided you have sufficiently frequent observations. For successive times $t \lt s,$ notice

$$y^{*}(s) - y^{*}(t) = y(s) - y(t) + \epsilon(s) - \epsilon(t) \mod 360 = \int_t^s f(t)\,\mathrm{d}t + \delta$$

where $\delta$ equals the contribution of the errors $\epsilon(s)-\epsilon(t)$ plus, perhaps, some integral multiple of $360$ whenever there has been an angular break between $y^{*}(t)$ and $y^{*}(s).$ Now, provided the size of the total error $|\epsilon(s)-\epsilon(t)|$ is less than $180$ degrees and provided the angle didn't go around more than once, we can figure out whether a break occurred: if $|\epsilon(s)-\epsilon(t)| \gt 180,$ add or subtract $360$ degrees from $\delta$ to place it into the interval from $-180$ to $+180$ degrees.

Although we can't observe these errors directly, if we are sampling frequently enough to make the increments $y(t_i) - y(t_{i-1})$ fairly small, we simply apply this adjustment to the observed differences. Thus,

Whenever $|y^{*}(s)-y^{*}(t)| \gt 180,$ add or subtract $360$ degrees from $\delta$ to place it into the interval from $-180$ to $+180$ degrees.

Equivalently, compute the differences modulo $180$ but express them in the range from $-180$ to $+180$ degrees rather than (as is conventional) the range from $0$ to $360.$

Let's call the adjusted value $\delta^{*}(t,s),$ so that

$$y^{*}(s) - y^{*}(t) = \int_t^s f(t)\,\mathrm{d}t + \delta(t,s)^{*}.$$

This is equality, not equality modulo $360.$ We may now remove the effect of recording the angles modulo $360$ by summing these adjusted differences. When the observations are made at times $t_0 \lt t_1\lt \cdots \lt t_n,$ we have

$$\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}$$

The problem with computation modulo $360$ is gone: you may now use any procedure you like to model the response $y^{*}(t).$


Here is an illustration with a fairly difficult dataset. The data were generated according to the model $y(t) = 30t \mod 360$ and observed annually from 1980 through 2020 with iid Normally distributed error of standard deviation $60$ degrees (a large amount).

Figure showing raw and adjusted data

The trend is barely discernible in the raw data, but the angle adjustment algorithm has visibly aligned them. We may fit a least squares model to the adjusted data, for instance, producing this result:

Figure showing the fits

The expanded vertical scale for the raw data shows details of the fit and their deviations from it. Incidentally, in this example the estimate of the slope is $28.0 \pm 0.74$ degrees, not remarkably different from the true value of $30$ degrees (the p-value for this comparison is $1.1\%$).

I will end by remarking that when the standard deviation of the errors $\epsilon(t)$ is large (greater than $180/2/\sqrt{2} \approx 64$ degrees, roughly), sometimes the angular adjustment will be incorrect. This will show up in the model residuals as a sudden change by a value around 360 degrees. Thus, a routine analysis of the model residuals can detect such problems, enabling you to modify the adjusted values for a better fit. The details of this will depend on your model and the fitting procedure.


This R code created the figures. At "adjust the angles" it shows how the angle adjustment can be computed efficiently.

#
# 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))
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.