# regression with circular response variable

I am analyzing some data where the variables of interests are circular (angles). I use R and the circular package.

In my dataset, every observation consists in a 2D Euclidean vector representing a movement on a plane (the length is constant, only the angle varies), together with a measure of a response to that movement, expressed as displacement in cartesian coordinates (i.e. the displacement along the x and y axis). I would like to fit a model where the initial movement (either the x-y cartesian representation or only the angle, since the length is constant) is the response variable. Ultimately, I would like to use the fitted model to estimate the initial movement vector in a related dataset where I only know the response to the movement.

My measurement of the response to the movement are affected by noise, independent over x and y dimensions, with possibly different variances and different additive biases. To represent graphically the problem, what I want to do is estimate the direction of the original vectors (gray thick arrows in the figure below) starting from the noisy measurements (the thin black arrows).

I have tried fitting different models:

1. one model with both predictor and dependent variable being circular (the angle of the initial movement vector, and the angle computed from the (x,y) displacement)
2. one multivariate, linear model, with the cartesian (x,y) components of the initial movement as dependent variables, and the measured x,y displacement as linear predictors;
3. one model with circular dependent variable, and linear predictor (this last one with no success)

First, I wasn't able to fit the 3rd model for some reason that I don't understand. I report here a reproducible example

require(circular)
rho <- rep(8,500)

# add measurement noise (different for x and y)
mX <- as.numeric(rnorm(500,0,2) + rho * cos(theta))
mY <- as.numeric(rnorm(500,-1,1) + rho * sin(theta))
linearPred <- cbind(mX,mY)

# fit model
mcl <- lm.circular(y=theta, x=linearPred,init=c(1,1), type="c-l",verbose=T)

Here is the output:

> mcl <- lm.circular(y=theta, x=linearPred,init=c(1,1), type="c-l",verbose=T)
Iteration  1 :    Log-Likelihood =  2.392981
Iteration  2 :    Log-Likelihood =  1.082084
Iteration  3 :    Log-Likelihood =  0.8013503
Iteration  4 :    Log-Likelihood =  NA
Error in while (diff > tol) { :
valore mancante dove è richiesto TRUE/FALSE

(the last line says: missing value where a TRUE/FALSE was required). Can anyone shed light on this? I don't understand where this error comes from.

Second question, which model would you suggest to use among the first two? Here is the code that I used for the two models

# fit sin(theta) and cos(theta) with a multivariate linear model
mmv <- lm(cbind(sin(theta),cos(theta)) ~ mX+mY)

# "circular-circular" model
mcc <- lm.circular(y=theta, x=angularPred, type="c-c")

I can compute the angle from the multivariate fitted values as atan2(mmv$fitted[,1],mmv$fitted[,2]). Both seems to perform similarly in terms of mean angular error. To compare the two I computed the correlation between the predicted angles and the initial angle theta (Jammalamadaka - Sarma correlation coefficient), and the multivariate model seems to perform slightly better:

fitted.mmv <- circular(atan2(mmv$fitted[,1],mmv$fitted[,2]),units="radians",type="angles",zero=0,rotation="counter")
> cor.circular(fitted.mmv,theta) # multivariate
[1] 0.9851422
> cor.circular(mcc\$fitted,theta) # "circular-circular"
[1] 0.7262862

However, the distribution of residuals of the multivariate model shows a strange pattern (figure below).

Is this pattern in the residuals a problem? Can anyone gives some advice on this? Which model would you use in this case? Is there any other possible approach that you would recommend? Any advice is appreciated, thanks!