Two-stage Circular-circular Regression with Zero-inflation: Application to Medical Sciences

TL;DR

This paper introduces a novel two-stage Bayesian circular-circular regression model for zero-inflated data, demonstrated with medical case studies, improving analysis accuracy in medical decision-making.

Contribution

It develops the first model to handle zero-inflated circular responses and covariates simultaneously, using Mobius transformation and Bayesian MCMC estimation.

Findings

Simulation shows superior performance over existing methods.

Applied to medical datasets on astigmatism and gait abnormalities.

Enhances decision-making in medical treatments.

Abstract

This paper considers the modeling of zero-inflated circular measurements concerning real case studies from medical sciences. Circular-circular regression models have been discussed in the statistical literature and illustrated with various real-life applications. However, there are no models to deal with zero-inflated response as well as a covariate simultaneously. The Mobius transformation based two-stage circular-circular regression model is proposed, and the Bayesian estimation of the model parameters is suggested using the MCMC algorithm. Simulation results show the superiority of the performance of the proposed method over the existing competitors. The method is applied to analyse real datasets on astigmatism due to cataract surgery and abnormal gait related to orthopaedic impairment. The methodology proposed can assist in efficient decision making during treatment or post-operative care.

Figures

Tables

Table 1. Table 1: Results of the simulation study with 15 % percent 15 15\% zeros in both the response and covariate.

n=50
Parameters	circular mean (c.d.)	Bias	CP	Parameters	mean (s.d.)	Relative Bias	CP
${\theta }_{01}=0$	-0.019(0.193)	-0.019	0.71	$r_1=0.9$	0.907(0.083)	0.008	0.90
${\theta }_{02}=0$	0.019(0.128)	0.019	0.74	$r_2=1.2$	1.177(0.100)	-0.019	0.92
${\theta }_{11}=0$	0.020(0.208)	0.020	0.70	${\rho }_1=0.85$	0.845(0.030)	-0.006	0.95
${\theta }_{12}=0$	-0.018(0.110)	-0.018	0.75	${\rho }_2=0.85$	0.846(0.030)	-0.005	0.95
n=100
Parameters	circular mean (c.d.)	Bias	CP	Parameters	mean (s.d.)	Relative Bias	CP
${\theta }_{01}=0$	-0.003(0.061)	-0.003	0.77	$r_1=0.9$	0.903(0.058)	0.003	0.91
${\theta }_{02}=0$	-0.004 (0.038)	-0.004	0.76	$r_2=1.2$	1.192(0.055)	0.007	0.91
${\theta }_{11}=0$	0.002(0.066)	0.002	0.77	${\rho }_1=0.85$	0.847(0.022)	-0.004	0.92
${\theta }_{12}=0$	0.003(0.031)	0.003	0.76	${\rho }_2=0.85$	0.848(0.020)	-0.002	0.95

Table 2. Table 2: Results of the simulation study with 10 % percent 10 10\% zeros in both the response and covariate.

n=50
Parameters	circular mean (c.d.)	Bias	CP	Parameters	mean (s.d.)	Relative Bias	CP
${\theta }_{01}=0$	0.025(0.198)	0.025	0.71	$r_1=1.2$	1.177(0.127)	-0.019	0.90
${\theta }_{02}=0$	-0.017(0.114)	-0.017	0.77	$r_2=1.2$	1.183(0.087)	-0.014	0.91
${\theta }_{11}=0$	-0.023(0.180)	-0.023	0.71	${\rho }_1=0.85$	0.843(0.030)	-0.008	0.95
${\theta }_{12}=0$	0.014(0.098)	0.014	0.78	${\rho }_2=0.85$	0.844(0.030)	-0.007	0.95
n=100
Parameters	circular mean (c.d.)	Bias	CP	Parameters	mean (s.d.)	Relative Bias	CP
${\theta }_{01}=0$	-0.005(0.044)	-0.005	0.79	$r_1=1.2$	1.189(0.080)	-0.009	0.93
${\theta }_{02}=0$	-0.011 (0.046)	-0.011	0.76	$r_2=1.2$	1.190(0.056)	-0.008	0.94
${\theta }_{11}=0$	0.005(0.039)	0.005	0.78	${\rho }_1=0.85$	0.847(0.021)	-0.004	0.95
${\theta }_{12}=0$	0.012(0.038)	0.012	0.76	${\rho }_2=0.85$	0.847(0.020)	-0.004	0.95

Table 3. Table 3: Results of the simulation study with 10 % percent 10 10\% zeros in response only.

n=50
Parameters	circular mean (c.d.)	Bias	CP	Parameters	mean (s.d.)	Relative Bias	CP
${\theta }_{01}=0$	-0.025(0.127)	-0.025	0.79	$r_1=0.9$	0.888(0.078)	-0.013	0.91
${\theta }_{02}=\frac{\pi }{2}$	1.552(0.087)	-0.019	0.77	$r_2=1.5$	1.448(0.122)	-0.034	0.88
${\theta }_{11}=0$	0.025(0.138)	0.025	0.80	${\rho }_1=0.85$	0.843(0.030)	-0.008	0.95
${\theta }_{12}=0$	0.016(0.062)	0.016	0.79	${\rho }_2=0.85$	0.844(0.031)	-0.007	0.93
n=100
Parameters	circular mean (c.d.)	Bias	CP	Parameters	mean (s.d.)	Relative Bias	CP
${\theta }_{01}=0$	-0.060(0.022)	-0.060	0.86	$r_1=0.9$	0.894(0.034)	-0.007	0.93
${\theta }_{02}=\frac{\pi }{2}$	1.587 (0.030)	0.016	0.78	$r_2=1.5$	1.481(0.082)	-0.013	0.91
${\theta }_{11}=0$	0.059(0.024)	0.059	0.88	${\rho }_1=0.85$	0.846(0.021)	-0.005	0.93
${\theta }_{12}=0$	-0.013(0.021)	-0.013	0.77	${\rho }_2=0.85$	0.847(0.021)	-0.004	0.93

Table 4. Table 4: Results of simulation study with 10 % percent 10 10\% zeros in covariate only.

n=50
Parameters	circular mean (c.d.)	Bias	CP	Parameters	mean (s.d.)	Relative Bias	CP
${\theta }_{01}=\frac{\pi }{2}$	1.543(0.075)	0.028	0.72	$r_1=0.9$	0.897(0.086)	-0.003	0.92
${\theta }_{02}=0$	0.009(0.095)	0.009	0.81	$r_2=1.5$	1.441(0.116)	-0.039	0.90
${\theta }_{11}=0$	0.029(0.165)	0.029	0.75	${\rho }_1=0.85$	0.844(0.031)	-0.007	0.95
${\theta }_{12}=0$	-0.009(0.070)	-0.009	0.81	${\rho }_2=0.85$	0.841(0.031)	-0.011	0.95
n=100
Parameters	circular mean (c.d.)	Bias	CP	Parameters	mean (s.d.)	Relative Bias	CP
${\theta }_{01}=\frac{\pi }{2}$	1.532(0.060)	-0.039	0.76	$r_1=0.9$	0.897(0.056)	-0.037	0.91
${\theta }_{02}=0$	0.008 (0.027)	0.008	0.83	$r_2=1.5$	1.484(0.083)	-0.011	0.91
${\theta }_{11}=0$	0.041(0.061)	0.041	0.76	${\rho }_1=0.85$	0.847(0.021)	-0.004	0.95
${\theta }_{12}=0$	-0.008(0.019)	-0.008	0.83	${\rho }_2=0.85$	0.847(0.021)	-0.004	0.94

Table 5. Table 5: Results of the simulation study without zero-inflation.

n=50
Parameters	circular mean (c.d.)	Bias	CP	Parameters	mean (s.d.)	Relative Bias	CP
${\theta }_{01}=\frac{\pi }{2}$	1.559(0.012)	-0.025	0.95	$r_1=0.3$	0.303(0.049)	0.010	0.94
${\theta }_{02}=\frac{\pi }{2}$	1.572(0.001)	0.009	0.94	$r_2=0.3$	0.318(0.044)	0.060	0.95
${\theta }_{11}=0$	0.030(0.013)	0.030	0.96	${\rho }_1=0.85$	0.842(0.031)	-0.009	0.95
${\theta }_{12}=0$	0.007(0.037)	0.007	0.94	${\rho }_2=0.85$	0.847(0.030)	-0.004	0.94
n=100
Parameters	circular mean (c.d.)	Bias	CP	Parameters	mean (s.d.)	Relative Bias	CP
${\theta }_{01}=\frac{\pi }{2}$	1.567(0.002)	-0.004	0.93	$r_1=0.3$	0.301(0.032)	0.003	0.92
${\theta }_{02}=\frac{\pi }{2}$	1.568 (0.006)	-0.003	0.93	$r_2=0.3$	0.308(0.029)	0.027	0.93
${\theta }_{11}=0$	0.003(0.003)	0.003	0.94	${\rho }_1=0.85$	0.847(0.020)	-0.004	0.95
${\theta }_{12}=0$	0.004(0.016)	0.004	0.91	${\rho }_2=0.85$	0.848(0.019)	-0.002	0.97

Table 6. Table 6: Results of the simulation study with approximately 35 % percent 35 35\% zeros in the response and 50 % percent 50 50\% zeros in the covariate.

Model I
Parameters	circular mean (c.d.)	Bias	CP	Parameters	mean (s.d.)	Relative Bias	CP
${\theta }_{01}=0.070$	0.073(0.004)	0.003	0.56	$r_1=0.9$	0.898(0.060)	-0.002	0.91
${\theta }_{02}=0.070$	0.063(0.001)	-0.007	0.47	$r_2=1.2$	1.186(0.057)	-0.012	0.90
${\theta }_{11}=0$	-0.002(0.011)	-0.002	0.57	${\rho }_1=0.93$	0.928(0.010)	-0.077	0.94
${\theta }_{12}=0$	-0.002(0.008)	-0.002	0.44	${\rho }_2=0.95$	0.949(0.008)	-0.001	0.95
Model II
Parameters	circular mean (c.d.)	Bias	CP	Parameters	mean (s.d.)	Relative Bias	CP
${\theta }_{01}=0.070$	0.092(0.006)	0.022	0.54	$r_1=0.9$	0.894(0.064)	-0.007	0.88
${\theta }_{02}=0.070$	0.036(0.001)	-0.034	0.09	$r_2=1.2$	1.019(0.010)	-0.151	0.08
${\theta }_{11}=0$	-0.046(0.007)	-0.046	0.55	${\rho }_1=0.93$	0.922(0.011)	-0.084	0.92
${\theta }_{12}=0$	-0.029(0.001)	-0.029	0.10	${\rho }_2=0.95$	0.990(0.003)	0.042	0.17
Model III
Parameters	circular mean (c.d.)	Bias	CP	Parameters	mean (s.d.)	Relative Bias	CP
${\theta }_{01}=0.070$	0.088(0.004)	0.018	0.55	$r_1=0.9$	0.899(0.057)	-0.001	0.89
${\theta }_{02}=0.070$	0.037(0.001)	-0.033	0.10	$r_2=1.2$	1.016(0.009)	-0.153	0.05
${\theta }_{11}=0$	-0.020(0.005)	-0.020	0.55	${\rho }_1=0.93$	0.928(0.011)	-0.077	0.96
${\theta }_{12}=0$	-0.031(0.001)	-0.031	0.09	${\rho }_2=0.95$	0.991(0.003)	0.043	0.14

Table 7. Table 7: Results of the simulation study under mis-specified model with p = 0.1 𝑝 0.1 p=0.1 .

Model I
Parameters	circular mean (c.d.)	Bias	CP	Parameters	mean (s.d.)	Relative Bias	CP
${\theta }_{01}=0.070$	0.075(0.003)	0.005	0.51	$r_1=0.9$	0.901(0.064)	0.001	0.90
${\theta }_{02}=0.070$	0.060(0.001)	-0.010	0.48	$r_2=1.2$	1.170(0.053)	-0.025	0.83
${\theta }_{11}=0$	-0.018(0.004)	-0.018	0.51	${\rho }_1=0.93$	0.934(0.010)	0.004	0.92
${\theta }_{12}=0$	-0.002(0.001)	-0.002	0.48	${\rho }_2=0.95$	0.955(0.007)	0.005	0.84
Model II
Parameters	circular mean (c.d.)	Bias	CP	Parameters	mean (s.d.)	Relative Bias	CP
${\theta }_{01}=0.070$	0.061(0.003)	-0.009	0.47	$r_1=0.9$	0.912(0.057)	0.013	0.79
${\theta }_{02}=0.070$	0.035(0.001)	-0.035	0.02	$r_2=1.2$	1.005(0.003)	-0.163	0.01
${\theta }_{11}=0$	-0.033(0.004)	-0.033	0.46	${\rho }_1=0.93$	0.937(0.010)	0.008	0.82
${\theta }_{12}=0$	-0.033(0.001)	-0.033	0.02	${\rho }_2=0.95$	0.998(0.001)	0.051	0.02
Model III
Parameters	circular mean (c.d.)	Bias	CP	Parameters	mean (s.d.)	Relative Bias	CP
${\theta }_{01}=0.070$	0.056(0.003)	-0.014	0.46	$r_1=0.9$	0.911(0.057)	0.011	0.87
${\theta }_{02}=0.070$	0.035(0.001)	-0.035	0.02	$r_2=1.2$	1.005(0.009)	-0.163	0.00
${\theta }_{11}=0$	-0.003(0.004)	-0.003	0.48	${\rho }_1=0.93$	0.935(0.011)	0.005	0.90
${\theta }_{12}=0$	-0.033(0.001)	-0.033	0.02	${\rho }_2=0.95$	0.998(0.003)	0.051	0.03

Table 8. Table 8: Results of the simulation study under mis-specified model with p = 0.2 𝑝 0.2 p=0.2 .

Model I
Parameters	circular mean (c.d.)	Bias	CP	Parameters	mean (s.d.)	Relative Bias	CP
${\theta }_{01}=0.070$	0.065(0.002)	-0.005	0.51	$r_1=0.9$	0.919(0.067)	0.021	0.87
${\theta }_{02}=0.070$	0.056(0.001)	-0.014	0.51	$r_2=1.2$	1.141(0.051)	-0.049	0.72
${\theta }_{11}=0$	-0.016(0.003)	-0.016	0.51	${\rho }_1=0.93$	0.940(0.009)	0.011	0.73
${\theta }_{12}=0$	-0.007(0.001)	-0.007	0.51	${\rho }_2=0.95$	0.961(0.007)	0.012	0.61
Model II
Parameters	circular mean (c.d.)	Bias	CP	Parameters	mean (s.d.)	Relative Bias	CP
${\theta }_{01}=0.070$	0.055(0.001)	-0.015	0.33	$r_1=0.9$	0.950(0.040)	0.056	0.54
${\theta }_{02}=0.070$	0.034(0.001)	-0.036	0.00	$r_2=1.2$	1.002(0.001)	-0.165	0.00
${\theta }_{11}=0$	-0.044(0.001)	-0.044	0.34	${\rho }_1=0.93$	0.963(0.008)	0.035	0.39
${\theta }_{12}=0$	-0.033(0.001)	-0.033	0.00	${\rho }_2=0.95$	0.999(0.001)	0.052	0.00
Model III
Parameters	circular mean (c.d.)	Bias	CP	Parameters	mean (s.d.)	Relative Bias	CP
${\theta }_{01}=0.070$	0.061(0.002)	-0.009	0.52	$r_1=0.9$	0.924(0.063)	0.027	0.84
${\theta }_{02}=0.070$	0.033(0.001)	-0.037	0.00	$r_2=1.2$	1.002(0.001)	-0.165	0.00
${\theta }_{11}=0$	-0.011(0.003)	-0.011	0.50	${\rho }_1=0.93$	0.941(0.009)	0.012	0.70
${\theta }_{12}=0$	-0.032(0.001)	-0.032	0.00	${\rho }_2=0.95$	0.999(0.001)	0.052	0.00

Table 9. Table 9: Computation Time in minutes.

	Model I	Model II	Model III
Table 6	6.75	4.95	5.01
Table 7	6.95	4.82	5.81
Table 8	7.66	5.21	5.95

Equations

$$Y={\beta }_0\frac{X+{\beta }_1}{1+\overline{{\beta }_1}X}ϵ,$$

$$Y={\beta }_0\frac{X+{\beta }_1}{1+\overline{{\beta }_1}X}{ϵ}_1,$$

$$X=b_0\frac{V+b_1}{1+\overline{b_1}V}{ϵ}_2,$$

$${\theta }_Y=\{\begin{array}{ll}0,& \text{mbox}if{\theta }_{Y^{\ast }}\in (-{\delta }_Y,{\delta }_Y),\\ {\theta }_{Y^{\ast }},& \text{mbox}otherwise,\end{array}$$

$${\theta }_X=\{\begin{array}{ll}0,& \text{mbox}if{\theta }_{X^{\ast }}\in (-{\delta }_X,{\delta }_X),\\ {\theta }_{X^{\ast }},& \text{mbox}otherwise,\end{array}$$

$$Y^{\ast }={\beta }_0\frac{X^{\ast }+{\beta }_1}{1+\overline{{\beta }_1}{X}^{\ast }}{ϵ}_1,$$

$$X^{\ast }=b_0\frac{V+b_1}{1+\overline{b_1}V}{ϵ}_2,$$

$$\pi ({\Theta }_1,{\theta }_{Y^{\ast }}|D)\propto \pi ({\Theta }_1)\times \prod _{i=1}^n[f_W({\theta }_{Y_i};{\mu }_{1i},{\rho }_1)\text{mathbbm}1({\theta }_{Y_i}\ne 0)$$

$$+K{f}_{TW}({\theta }_{Y_i^{\ast }};{\mu }_{1i},{\theta }_{X_i^{\ast }},-{\delta }_Y,{\delta }_Y)\text{mathbbm}1({\theta }_{Y_i}=0)],$$

$$\pi ({\Theta }_2,{\theta }_{X^{\ast }}|D)\propto \pi ({\Theta }_2)\times \prod _{i=1}^n[f_W({\theta }_{X_i};{\mu }_{i2},{\rho }_2)\text{mathbbm}1({\theta }_{X_i}\ne 0)$$

$$+K{f}_{TW}({\theta }_{X_i^{\ast }};{\mu }_{i2},{\rho }_2,-{\delta }_X,{\delta }_X)\text{mathbbm}1({\theta }_{X_i}=0)],$$

$$\pi ({\theta }_{X_i^{\ast }}|{\Theta }_1,{\Theta }_2,D)\equiv \{\begin{array}{ll}\text{mathbbm}1({\theta }_{X_i^{\ast }}={\theta }_{X_i}),& \text{mbox}if{\theta }_{X_i}\ne 0\\ f_{TW}({\theta }_{X_i\ast };{\mu }_{2i},{\rho }_2,-{\delta }_X,{\delta }_X),& \text{mbox}otherwise,\end{array}$$

$$\pi ({\theta }_{Y_i^{\ast }}|{\Theta }_1,{\Theta }_2,D)\equiv \{\begin{array}{ll}\text{mathbbm}1({\theta }_{Y_i^{\ast }}={\theta }_{Y_i}),& \text{mbox}if{\theta }_{Y_i}\ne 0\\ f_{TW}({\theta }_{Y_i^{\ast }};{\mu }_{1i},{\rho }_1,-{\delta }_Y,{\delta }_Y),& \text{mbox}otherwise,\end{array}$$

$$\pi ({\Theta }_1|{\Theta }_2,{\theta }_{X^{\ast }},{\theta }_{Y^{\ast }},D)\propto \pi ({\Theta }_1)\prod _{i=1}^n{f}_W({\theta }_{Y_i^{\ast }};{\mu }_{1i},{\rho }_1),$$

$$\pi ({\Theta }_2|{\Theta }_1,{\theta }_{X^{\ast }},{\theta }_{Y^{\ast }},D)\propto \pi ({\Theta }_2)\prod _{i=1}^n{f}_W({\theta }_{X_i^{\ast }};{\mu }_{2i},{\rho }_2),$$

$${\alpha }_{MH}({\nu }_p,{\nu }_p^{\prime })=\min \left\{1,\frac{\pi ({\nu }_p^{\prime }|-)\pi ({\nu }_p^{\prime },{\nu }_p)}{\pi ({\nu }_p|-)\pi ({\nu }_p,{\nu }_p^{\prime })}\right\},$$

$$Y={\beta }_0\frac{X^{(s)}+{\beta }_1}{1+\overline{{\beta }_1}{X}^{(s)}}ϵ,$$

$${\theta }_{X_1}=\{\begin{array}{ll}0,& \text{mbox}if{\theta }_{X_1^{\ast }}\in \left(-{\delta }_{X_1},{\delta }_{X_1}\right),\\ {\theta }_{X_1^{\ast }},& \text{mbox}otherwise,\end{array}$$

$$Y^{\ast }={\beta }_0\frac{X^{(a)}+{\beta }_1}{1+\overline{{\beta }_1}{X}^{(a)}}{ϵ}_1,$$

$$X_1^{\ast }=b_0\frac{W+b_1}{1+\overline{b_1}W}{ϵ}_2,$$

Full text

Keywords: Abnormal gait, Astigmatism, Latent variable, Metropolis-Hastings algorithm, Rose diagram, Truncated wrapped Cauchy.

1 Introduction

In many real-world experiments, measurements are taken in angles, and such variables are commonly labeled as ‘circular variables’ or ‘directional variables’ in the statistics literature (Mardia and Jupp,, 2000). The need for developing statistical methods to study circular regression is important due to its wide application in various fields of science, e.g., orthopaedics (Jammalamadaka and SenGupta,, 2001, p-4), meteorology (Kato et al.,, 2008), ecology (Fisher and Lee,, 1992), and geoscience (Rivest,, 1997). Sarma and Jammalamadaka, (1993) proposed a generic circular-circular regression by joint modeling of the sine and cosine of the angular response on the trigonometric polynomial function of the angular covariate. In this context, some rotational models have been proposed where the predicted mean direction of the response is a fixed rotation of the covariate. See Mackenzie, (1957) and Rivest, (1997) for more details. Bhattacharya and Sengupta, (2009) proposed a Bayesian hierarchical framework for linear–circular regression. Downs and Mardia, (2002) proposed the Möbius transformation based regression link function for circular-circular regression to analyze the relationships between the peak systolic blood pressure times on different days, and wind directions at two different times. Later, Kato et al., (2008) also considered the Möbius transformation based link function by reparameterizing the model of Downs and Mardia, (2002) and applied this model to study the association between the spawning time of fish and the time of low tide. Other applications of the Möbius transformation based circular-circular regression model in opthalmology and genetics are discussed in Jha and Biswas, (2018) and Rueda et al., (2015), respectively. In this article, we consider the problem of modeling circular measurements motivated by the following case studies from orthopaedics and ophthalmology.

1.1 Abnormal Gait Data

The ability to walk safely and efficiently is essential for an independent and productive life. Gait disorders lead to a loss of personal freedom, and increased risk of injuries that result in a marked reduction in the quality of life (Mahlknecht et al.,, 2013). In the last couple of decades, due to the increased awareness of the importance of gait control, researchers across the globe focused their research programs on the clinical and therapeutic aspects of walking and balance (Nutt et al.,, 2011). Orthopaedic impairments of the lower extremities are the most common reasons for non-neurological gait disorders in adults (Mahlknecht et al.,, 2013). Patients with gait disorders caused by orthopaedic impairments are often referred to a physical therapist for treatment. Therapists visually observe the patient’s gait to determine effective treatment protocol or evaluate the effect of a therapeutic intervention. Various kinematic features like sagittal plane angular movements (flexion-extension) at the ankle and knee are extremely important in determining abnormal gait Therapists visually observe the patient’s gait to determine effective treatment protocol or evaluate the effect of a therapeutic intervention. Various kinematic features like sagittal plane angular movements (flexion-extension) at the ankle and knee are extremely important in determining abnormal gait(Jammalamadaka et al.,, 1986). Bio-engineers are interested to analyse the interlimb coordination concerning the flexion-extension to make an efficient design for the body-powered knee-ankle prosthesis (Ingraham et al.,, 2016). The human walking cycle consists of a period when the foot is on the ground, followed by the next phase of the forward movement known as swing. Recent research suggests that the ankle extensors provide necessary kinetic energy for the initiation of swing (Jessop and Pain,, 2016). However, there is a lack of knowledge about the biomechanical process, which leads to the high power output of the ankle extensors (Lipfert et al.,, 2014). In this article, we consider a study conducted at the Bio-Engineering Unit, University of Calcutta, on a total of 29 individuals, including 10 healthy persons and 19 patients with unilateral orthopaedic impairment. Measurements on ankle and knee flexion-extension were recorded using electrogoniometers while subjects were walking at their own pace along an 8-meter long walkway turning around at each end. Previous research found no statistically significant difference between the motions of the right and left limb of a healthy person. Therefore, the angular movements of patient’s healthy limb can be used for comparison with the affected side in the presence of orthopaedic impairments (Roaas and Andersson,, 1982). The principal objective of this study is to assess the recovery of orthopaedic patients and compare their gait with that of a healthy person.

1.2 Cataract Surgery Data

A cataract is a clouding that develops in the natural lens of the eye or its envelope. In due course of time, the lens loses its transparency and leads to partial or total loss of vision. This has been documented to be the most significant cause of bilateral blindness in India (Thulasiraj et al.,, 2003; Murthy et al., 2008a, ; Mohan,, 1989). Cataract surgery is the removal of the opaque natural lens from the eye, and an artificial intra-ocular lens implant is then inserted to restore vision. India is a signatory to the World Health Organization resolution on VISION 2020: The Right to Sight. Efforts from all stakeholders have resulted in an increased number of cataract surgeries performed in India (Murthy et al., 2008b, ). It is well-known that one common side effect of the cataract surgery is that the incision causes unwanted changes to the natural corneal shape causing an astigmatic eye. The refractive error of the astigmatic eye induces several focal points in different directions. For example, the image may be perfectly focused on the retina in the horizontal (sagittal) plane, but not in the vertical (tangential) plane. There are two types of astigmatism based on the axes of the principal meridians - regular or irregular. In this paper, we consider regular astigmatism, where the principal meridians are perpendicular. For regular astigmatism, its axis ranges from [math] to [math]. In general, regular astigmatism can be subdivided into three types: (i) With-the-rule astigmatism (WTR) - the vertical meridian ([math]) is steepest (e.g. an American football lying on its side), (ii) Against-the-rule astigmatism (ATR)- the horizontal meridian ([math]) is steepest (e.g. an American football standing on its end), and (iii) Oblique astigmatism - the steepest curve lies in ($$,$$) or ($$,$$) (Morlet et al.,, 2001; Mimouni et al.,, 2017). The eye affected by WTR astigmatism sees vertical lines more sharply than horizontal lines. The situation is reversed for ATR astigmatism. Oblique astigmatism is worse than WTR and ATR astigmatism because most of the standard objects (e.g. letters) in our surroundings are horizontal or vertical. The objects get distorted horizontally or vertically for WTR and ATR astigmatism, whereas the distortion is severe for oblique astigmatism. The visual distortions for different types of astigmatism are displayed in Figure 1. It is of primary interest for ophthalmologists to study post-operative astigmatism and visual recovery over time (Zheng et al.,, 1997). In this article, we consider a study conducted at Disha Eye Hospital and Research Center, Barrackpore, West Bengal, India, over a period of two years (2008-10). In total, 54 patients were operated, and the axes of astigmatism were measured on the 1st, 7th, and 15th days after the surgery. See Bakshi, (2010) for a detailed description of the study and data description. The main objective of regular monitoring is to study the process of visual recovery and identify patients who require additional care to minimize post-operative trauma. In particular, medical practitioners are interested in foreseeing the improvement of patients based on previous inspections.

1.3 The challenge in analysing Cataract Surgery and Abnormal Gait Data

There are some critical issues involved in the modeling of astigmatism data under consideration. According to medical practitioners, if the axis is closer to [math], [math] or [math], then it is not a matter of serious concern as most of the standard objects in nature are horizontal or vertical. Therefore, it is preferred that the axis is closer to [math], [math] or [math]. In order to make only one preferred direction, the observed angles are multiplied with 4 and then transformed by taking mod [math]. Consequently, the preferred angle reduces to [math] ( =$$$) and hence, the multimodal distribution becomes a distribution having a single mode at [math]. The effect of this transformation on circular distribution is illustrated in Figure [2](#S1.F2). In this particular study, the measurements were taken up to the precision of [math]. Therefore, the transformed variables become zero-inflated due to the high concentration of observations censored in the interval (-$,$). The transformed axes of astigmatism after 7 days and 15 days of the surgery are presented using circular plots in Figure [3](#S1.F3). A similar issue arises related to the measurements on ankle extension in abnormal gait data. The measurements were taken up to the precision of [math] and averaged over 20$strides of each individual. Almost$50%$of the observations on healthy individuals, and$60%$ of the observations on orthopaedically impaired patients, are recorded as [math].

Numerous studies focusing on zero-inflated random variables in the linear setup are available in the literature. Tobin, (1958), Heckman, (1974, 1979) proposed some models in the context of linear regression where the responses are zero-inflated. Lambert, (1992) considered Poisson regression for count data with excess zeros in the response variable. Bhuyan et al., (2019) discussed the case when both the responses and covariates are zero-inflated and proposed estimation methodology under a Bayesian setup. See Min and Agresti, (2002) for a detailed review of the zero-inflated regression models. However, in all of these works, the response and the covariate are linear in nature, where zero-inflation is caused by censorship through a selection mechanism. The case of circular random variables is considerably different from linear ones. The difference mainly arises due to the topology of the circle, where the zero cannot be considered to be located at the boundary of the sample space. In this context, a spike at zero may be interpreted in the same way as a spike at any other angle, since the origin can be fixed arbitrarily without loss of generality. There are only a few works on the analysis of zero-inflated circular data in the literature. The modeling of a circular random variable with point-accumulation was studied in Biswas et al., (2016). In the context of circular-circular regression, Jha and Biswas, (2018) proposed a model with a zero-inflated response variable and discussed the associated inferential issues. However, there is no model available in the literature for the case of a zero-inflated circular covariate. To avoid such difficulty, Jha and Biswas, (2018) analysed a subset of the dataset, considering only non-zero covariate values. It is important to note that the conventional circular-circular regression model and the model proposed by Jha and Biswas, (2018) are not appropriate for handling excess zeros in the covariate, and provide biased results. Moreover, the existing models are not capable of joint modeling of the periodical observations on astigmatism recorded over three different inspections. To model such data, we propose a two-stage circular-circular regression model based on continuous latent variables. In contrast to the assumption of randomly occurring zeros in Jha and Biswas, (2018), we consider a more realistic assumption that zero inflation occurs due to censoring. In Section 2, we discuss the modeling approach and propose an estimation methodology under a Bayesian setup using the Markov chain Monte Carlo (MCMC) algorithm. Some generalisations and a special case of the proposed model are also discussed in the same section. The performance of the proposed method is compared with the existing competitors through simulation in Section 3. Analysis of real datasets on post-operative astigmatism and abnormal gait are presented in Section 4. The key findings are summarised and concluded with some discussions on future research in Section 5.

2 Proposed Model and Methodology

In this section, we first describe the Möbius transformation based circular-circular regression model proposed by Kato et al., (2008), which is a reparameterization of the model proposed by Downs and Mardia, (2002). This reparameterization induces a nice geometry to the regression link function. When the angular error follows wrapped Cauchy distribution, it provides some advantages in terms of distributional properties as the wrapped Cauchy distribution is closed under rotation and Möbius transformation. See Kato et al., (2008) for details. It is important to note that the rotational model proposed by Mackenzie, (1957), where the predicted mean direction of the response is a fixed rotation of the covariate, is a special case of the Möbius transformation based model. Unlike the rotational model, this model is also appropriate when there is a high concentration of observed responses on a section of the unit circle. This unique feature makes the model parameters easily interpretable and provides interesting insights.

Let us represent circular random variables $\theta_{Y}$ and $\theta_{X}$ as complex random variables $Y=e^{i\theta_{Y}}$ and $X=e^{i\theta_{X}}$, respectively taking values on the circumference of a unit circle. The circular-circular regression model of Kato et al., (2008) is represented as:

[TABLE]

where $\beta_{0}$, $\epsilon$ $\in\{z:z\in\mathbbm{C};|z|=1\}$, $\beta_{1}\in\mathbbm{C}$, and the angular error $\arg(\epsilon)$ follows a wrapped Cauchy distribution with mean direction [math]. The regression link function, which is a form of Möbius transformation, is a mapping from a unit circle onto itself. Here, $\beta_{0}$ is the rotation parameter because the multiplication by a unit complex number is an anti-clockwise rotation by the argument of the same unit complex number. In this case, the predicted mean $\mu(\theta_{x})$ given $x$ is obtained by rotating $x_{\beta}$ by $\alpha=\arg(\beta_{0})$, where $x_{\beta}=\frac{x+\beta_{1}}{1+\overline{\beta_{1}}x}$ is the intersection of the unit circle with the line joining $-x$ and $\beta_{1}$ (see Figure 4). For $|\beta_{1}|>1$, the regression link function can be geometrically represented as a straight line connecting $1/\overline{\beta_{1}}$ and $\frac{\beta_{1}}{|\beta_{1}|}\frac{\beta_{1}}{|\beta_{1}|}\overline{x}$. Therefore, it also covers the cases in which the predicted mean direction depends on the conjugate of $x$ (i.e. when there is a reflection of $x$). The intersection of this line with the unit circle is then rotated by $\alpha$ to obtain the predicted mean direction $\mu$. If $|\beta_{1}|$ is closer to 1 and $x$ is uniformly distributed, then $x_{\beta}$ is highly concentrated around $\frac{\beta_{1}}{|\beta_{1}|}$. The distribution of $Y$ becomes independent of $x$ if $|\beta_{1}|=1$. For $|\beta_{1}|=0$, the predicted mean direction is just a rotation of $x$, i.e. $\mu=\beta_{0}x$. See Kato et al., (2008) for more details.

To model the axis of astigmatism recursively based on three consecutive inspections, we extend the circular-circular regression model (1) in a two-stage setup as:

[TABLE]

[TABLE]

where $\beta_{0}$, $b_{0}$, $\epsilon_{1}$, $\epsilon_{2}$ $\in\{z:z\in\mathbbm{C};|z|=1\}$; $\beta_{1}$ ,$b_{1}$ $\in\mathbbm{C}$, $Y=e^{i\theta_{Y}}$, $X=e^{i\theta_{X}}$ and $V=e^{i\theta_{V}}$. We assume $\arg(\epsilon_{i})\sim WC(0,\rho_{i})$ for $i=1,2$, where $WC(\mu,\rho)$ represents the wrapped Cauchy distribution with parameters $\mu\in[0,2\pi)$ and $\rho\in[0,1]$. Also, $\arg(\epsilon_{1})$ and $\arg(\epsilon_{2})$ are assumed to be independently distributed. In the literature of Econometrics, $V$ is known as instrumental variable (Cameron and Trivedi,, 2005, p-97), which may not directly effect the response $Y$ but it can induce changes only through the covariate $X$. This modeling approach allows us to compare the visual recovery of the patients over two consecutive weeks.

Note that the circular-circular regression model, given by equation (1), is inadequate for handling zero-inflated data (Jha and Biswas,, 2018). Similarly, the aforementioned two-stage model will give very poor fit to the data due to point accumulation at zero for both the response and the covariate. In order to model zero-inflated circular data, we first define circular latent variables $\theta_{Y^{*}}$ and $\theta_{X^{*}}$ as

[TABLE]

and

[TABLE]

respectively, where $\delta_{Y}$, $\delta_{X}$ are constants taking values in $[0,\pi)$. Now we propose the two-stage circular-circular regression model based on the aforementioned latent variables as

[TABLE]

[TABLE]

where $\beta_{1}$, $b_{1}$ $\in\mathbbm{C}$, $Y^{*}=e^{i\theta_{Y^{*}}}$, $X^{*}=e^{i\theta_{X^{*}}}$. In the absence of instrumental variable, one can still use (4) and (5) for modeling of zero-inflated response and covariate with $|b_{1}|=1$. Note that the above model reduces to the two-stage circular-circular regression model without zero-inflation, given by (2) and (3), when $\delta_{X}=\delta_{Y}=0$.

2.1 Bayesian Estimation

We propose a Bayesian estimation of the model parameters involved in the two-stage model given by (4) and (5), using the MCMC algorithm based on data augmentation. In the conventional frequentist approach, computational challenges arise in the model fitting due to intractable numerical integration involved in the log-likelihood function. Unlike the usual frequentist set-up, the proposed Bayesian approach provides simultaneous estimation of the model parameters and a natural framework for prediction over unobserved data. Thus, by generating posterior predictive densities, rather than point estimates, we can make probability statements giving greater flexibility in presenting results. For instance, we can discuss findings concerning specific hypotheses or in terms of credible intervals which can offer a more intuitive understanding for the practitioners.

Let us denote the observed data by $D=\{(\theta_{V_{i}},\theta_{X_{i}},\theta_{Y_{i}}),i=1\ldots n\}$, and the parameter vector by $\Theta_{i}=(\theta_{0i},r_{i},\theta_{1i},\rho_{i})$ for $i=1,2$, where $\beta_{0}=e^{i\theta_{01}}$, $\beta_{1}=r_{1}e^{i\theta_{11}}$, $b_{0}=e^{i\theta_{02}}$, $b_{1}=r_{2}e^{i\theta_{12}}$, $r_{1},r_{2}\in[0,\infty)$. Also, denote the density of the wrapped Cauchy distribution and the truncated wrapped Cauchy distribution by $f_{W}(\theta;\mu,\rho)$ and $f_{TW}(\theta;\mu,\rho,-\delta,\delta)=K^{-1}f_{W}(\theta;\mu,\rho)\mathbbm{1}(\theta\in(-\delta,\delta))$, respectively, where $K=\int_{-\delta}^{\delta}f_{W}(\theta;\mu,\rho)d\theta$. The joint posterior density of the model parameters and the latent variable involved in equation (4) is given by

[TABLE]

where $\mu_{1i}=\arg(\beta_{0}\frac{X_{i}^{*}+\beta_{1}}{1+\overline{\beta_{1}}X_{i}*})$, and $\pi(\Theta_{1})$ denotes the prior density of $\Theta_{1}$. Similarly, the joint posterior density of the model parameters and the latent variable involved in equation (5), is given by

[TABLE]

where $\mu_{2i}=\arg(b_{0}\frac{V_{i}+b_{1}}{1+\overline{b_{1}}V_{i}})$, and $\pi(\Theta_{2})$ denotes the prior density of $\Theta_{2}$. The full conditional densities of the latent variables $\theta_{X_{i}^{*}}$ and $\theta_{Y_{i}^{*}}$ have the following closed-form expressions

[TABLE]

and

[TABLE]

respectively. We propose an algorithm for generating samples from the truncated wrapped Cauchy distribution which is discussed in Subsection 2.2. The full conditional densities of the model parameters $\Theta_{1}$ and $\Theta_{2}$ are given by

[TABLE]

[TABLE]

respectively. Note that the full conditionals of $\Theta_{1}$ and $\Theta_{2}$ cannot be expressed in closed form. Therefore, we employ the Metropolis-Hastings algorithm for generating samples from the posterior densities of the parameters and the detailed algorithm is provided in Subsection 2.2. In our context, $\delta_{Y}$ and $\delta_{X}$ are known a priori. However, in many situations, these are not known and one can consider a suitable prior as discussed in Subsection 2.2.

2.2 Sampling Algorithms

2.2.1 Sample Generation from Truncated Wrapped Cauchy Distribution

Let $\theta_{Z}$ be a circular random variable following truncated wrapped Cauchy distribution with pdf $f_{TW}(\theta_{z};\mu,\rho,a,b)$, where $a,b\in[-\pi,\pi)$. The support of $f_{TW}(\cdot)$ is $(a,b)$, if $a<b$, and $(a,\pi)\cup[-\pi,b)$, otherwise. To generate samples from $f_{TW}(\theta_{z};\mu,\rho,a,b)$, for $a<b$, we simulate observations from $f_{W}(\theta_{z};\mu,\rho)$ and accept the observations which lie in $(a,b)$. Similarly, for $a>b$, we accept the observations lying in $(a,\pi)\cup[-\pi,b)$. However, the acceptance rate is very low for small values of $(b-a)1(a<b)+\left[2\pi-(b-a)\right]{1(a>b)}$. For example, the acceptance rate is approximately $0.03\%$ when $a=\pi-0.035$, $b=-\pi+0.035$, $\mu=0$, and $\rho=0.95$. Thus, we propose a novel algorithm for generating samples from the truncated wrapped Cauchy distribution based on the geometry of the Möbius transformation.

If $\theta_{Z}$ is uniformly distributed in $[0,2\pi)$ and $Z=e^{i\theta_{z}}$, then $\arg\left(\frac{\psi-Z}{1-\overline{\psi}Z}\right)$ follows wrapped Cauchy distribution with parameters $\mu=\arg(\psi)$, and $\rho=|\psi|$, where $\psi\in\{c\in\mathbbm{C}:|c|\leq 1\}$ (Kato et al.,, 2008). The geometry of the transformation $\eta(Z)=\frac{\psi-Z}{1-\overline{\psi}Z}$ is presented in Figure 4 with $x=-Z$, $\beta_{1}=\psi$ and $\alpha=0$. Note that, $\eta(Z)$ is the point on the circumference of the unit circle which is situated at the intersection of the line joining $Z$ and $\psi$ and the unit circle. It is easy to see that $\arg\{\eta(Z)\}$ follows wrapped Cauchy distribution with $\mu=\arg(\psi)$ and $\rho=|\psi|$. To simulate an observation $U$ from the truncated wrapped Cauchy distribution, first generate a unit complex number $\xi$ uniformly in the region between $\eta(A)$ and $\eta(B)$ where $A=e^{ia}$ and $B=e^{ib}$. Then, consider the argument of its inverse Möbius transformation $\eta^{-1}(\xi)$ as $U$. This construction is diagrammatically illustrated in Figure 5 and summarized in a simple algorithm below.

Algorithm 1: Sampling from Truncated Wrapped Cauchy Distribution

• Step 1:

Take $A=e^{ia}$, $B=e^{ib}$ and choose a point $c$ in the support of $f_{TW}(\theta_{z};\mu,\rho,a,b)$.

• Step 2:

Generate a random unit complex number $\xi$ uniformly from the arc joining $\eta(A)$ and $\eta(B)$ containing $\eta(C)$, where $C=e^{ic}$.

• Step 3:

Take $U=\arg\{(\eta^{-1}(\xi)\}$.

2.2.2 Metropolis-Hastings Algorithm

For the purpose of Bayesian estimation, we consider the following prior distributions for the parameter vectors $\Theta_{1}$ and $\Theta_{2}$. The joint prior distribution $\pi(\Theta_{i})$ can be expressed as the product of $\pi(\theta_{0i})\propto 1$, $\pi(\theta_{1i})\propto 1$, $\pi(r_{i})\propto e^{-r_{i}^{2}}$ and $\pi(\rho_{i})\propto\rho_{i}^{a_{\rho_{i}}-1}(1-\rho_{i})^{a_{\rho_{i}}-1}$, for $i=1,2$. It can be easily verified that the joint posterior density is proper for $a_{\rho_{i}}>1$. Similar priors have been considered by Ravindran and Ghosh, (2011) in the context of circular-circular and circular-linear regressions. For the purpose of implementing Metropolis-Hastings Algorithm, we consider the proposal distributions for $\theta_{01},\theta_{11},\theta_{02},\theta_{12},\rho_{1},\rho_{2}$ to be uniform, and the proposal distributions for $r_{1},r_{2}$ are chosen to be exponential. As discussed before, $\delta_{Y}$ and $\delta_{X}$ may be unknown and one can employ the following MCMC algorithm with priors $\pi(\delta_{Y})\propto 1$ and $\pi(\delta_{X})\propto 1$, and uniform proposal densities.

**Algorithm 2: Metropolis Hastings Algorithm

**

• Step 1:

Sample $\theta_{X^{*}}$ from the density $\pi(\theta_{X^{*}}|\Theta_{1},\Theta_{2},D)$ given by equation (6), and then sample $\theta_{Y*}$ from the density $\pi(\theta_{Y^{*}}|\Theta_{1},\Theta_{2},D)$ given by equation (7).

• Step 2:

Generate the model parameters sequentially from the corresponding proposal densities and denote it as $\nu_{p}^{\prime}$. Given the previous value of $\nu_{p}$ and the current draw $\nu_{p}^{\prime}$, return $\nu_{p}^{\prime}$ with probability

[TABLE]

where $\pi(d,w)$ denotes the proposal density at $d$ with parameter $w$, and $\pi(\cdot|)$ denotes the full conditional densities given by equations (8) and (9). Otherwise, repeat the previous value $\nu_{p}$.

• Step 3:

Repeat Step 1 and Step 2 until convergence.

2.3 Special case

Under certain restrictions on the model parameters, the geometry of the Möbius transformation based link function provides interesting insights when there is a preferred direction for the response variable. Without loss of generality, we consider [math] to be the preferred direction. Then, we consider a special case of the circular-circular regression model given by equation (1), where $\arg(\beta_{0})=\arg(\beta_{1})=0$ and $\beta_{1}\in[-1,1]$. Under this restriction, the symmetry about the preferred direction for the predicted responses and the corresponding covariate values is maintained as $\arg\{\mu(\theta_{x})\}=-\arg\{\mu(-\theta_{x})\}$. This can be seen directly from the geometry as shown in Figure 4 by taking $\beta_{1}\in[-1,1]$ and $\alpha=0$. In this case, the larger the value of $\beta_{1}$, the more is the tendency of the predicted response to move towards [math]. The constraint $\beta_{1}\in[-1,1]$ ensures that both the predicted mean direction and the corresponding covariate value lie on the same semi-circle having the diameter as the line segment joining [math] and [math]. This implies that the predicted mean direction lies on the shortest arc joining the respective covariate value and [math] when $\beta_{1}>0$. This special case also ensures that $\mu(\theta_{x_{1}})<\mu(\theta_{x_{2}})$, for all $\theta_{x_{1}}<\theta_{x_{2}}$, where $\theta_{x_{1}},\theta_{x_{2}}\in[0,\pi)$. That is, the closer the covariate value is to [math], the closer is the corresponding predicted mean response to [math].

This special case can be reparameterized and presented as the centered regression model proposed by Downs and Mardia, (2002) with both the angular location parameters fixed at 0. The corresponding regression link function can be written as $\arg\{\mu(\theta_{x})\}=2\arctan\{\omega\tan(\frac{\theta_{x}}{2})\}$, where $\omega\in[0,\infty)$. As the angular parameters are fixed at 0, there is a one-one monotonic correspondence between $\omega$ and $\beta_{1}$, given by $\beta_{1}=\frac{1-\omega}{1+\omega}$. Therefore, the parameter $\omega$, if greater than or less than 1, indicates the general tendency of the predicted response to move towards or away from [math] as compared to its covariate value, respectively.

For the two-stage circular-circular regression model given by equations (4) and (5), we consider $\arg(b_{0})=\arg(b_{1})=\arg(\beta_{0})=\arg(\beta_{1})=0$, and $b_{1},\beta_{1}\in[-1,1]$. As compared to the unrestricted model, this special case enforces all the predicted response to either move towards or away from zero as compared to its respective covariate value. Therefore, the two stages of the regression model can be compared in terms of their tendency to move towards [math] on the basis of the posterior distributions of $\beta_{1}$ and $b_{1}$. In contrast to the general model, we consider uniform prior and proposal distributions for $r_{1}$ and $r_{2}$. The choices of prior and proposal distributions for all the other parameters are considered same as discussed in Subsection 2.2.2. When there is a preferred direction, it is easier to compare two regression models based on the estimates obtained from the special case as compared to Model I. This special case is parsimonious and easy to interpret. However, a drawback of this special case as compared to the unrestricted model is that all the predicted mean responses simultaneously move towards or away from the preferred direction. Therefore, this model is not appropriate for prediction when some of the observed responses are closer, and some are farther from the preferred direction compared to their respective covariate values. We consider this special case for data analysis to provide some insights into the cataract surgery data based on the posterior distributions of the model parameters.

2.4 Some Generalisations

The proposed two-stage model, given by (4) and (5), can be extended to the case when there are multiple circular covariates. Jha and Biswas, (2017) proposed a multiple circular-circular regression model (MCR2) based on Möbius transformation which is represented as:

[TABLE]

where $X^{(s)}=\frac{\sum_{j=1}^{k}p_{j}A_{j}X_{j}}{|\sum_{j=1}^{k}p_{j}A_{j}X_{j}|}$, $A_{1}=1$, $A_{j}\in\{z:z\in\mathbb{C};|z|=1\}$ for $j=2,\ldots k,$ $p_{j}\in[0,1]$ for $j=1,\ldots k,$ $\sum_{j=1}^{k}p_{j}=1$, and $X_{1},\ldots,X_{k}$ are covariates. Note that the model (10) reduces to (1) for $k=1$. Now, without loss of generality, we consider that $X_{1}$ is a zero-inflated covariate and define a latent circular variable $X^{*}_{1}$ as:

[TABLE]

where $\theta_{X_{1}^{\ast}}=\arg(X_{1}^{\ast})$. Note that the the response $Y$ is also zero-inflated and we consider the latent response $Y^{\ast}$ as defined in Section 2. Then, the generalised two-stage circular-circular regression model for zero-inflated data with multiple covariates is given as:

[TABLE]

and

[TABLE]

where $X^{(a)}=\frac{p_{1}A_{1}X_{1}^{*}+\sum_{j=2}^{k}p_{j}A_{j}X_{j}}{|p_{1}A_{1}X_{1}^{*}+\sum_{j=2}^{k}p_{j}A_{j}X_{j}|}$, $W=\frac{\sum_{i=1}^{l}q_{i}B_{i}W_{i}}{|\sum_{i=1}^{l}q_{i}B_{i}W_{i}|}$, $B_{1}=1$, $B_{i}\in\{z:z\in\mathbb{C};|z|=1\}$ for $i=2,\ldots l$, $q_{i}\in[0,1]$, $\sum_{i=1}^{l}q_{i}=1$, and $W_{1},\ldots,W_{l}$ are covariates used to regress $X_{1}^{\ast}$. For the purpose of Bayesian estimation, one can consider Dirichlet priors for $(p_{1},\ldots p_{k})$ and $(q_{1},\ldots q_{l})$ with parameter vectors $\textbf{1}_{k\times 1}$ and $\textbf{1}_{l\times 1}$, respectively. The priors for $\arg(A_{j})$, for $j=2,\ldots,k$, and $\arg(B_{i})$ for $i=2,\ldots,l$, can be taken as uniform. The MCMC algorithm mentioned in Subsection 2.2 is readily extended to this generalised model. Further, one can generalise the proposed model for modeling data with multiple points of accumulation due to asymmetric censoring intervals.

3 Simulation Studies

In order to study the performance of the proposed method, we generate data considering different choices of parameter values for two different sample sizes, 50 and 100. We generate $\theta_{V}$ from von Mises distribution with mean 0 and concentration parameter $2$. The five different sets of parameters values of $\Theta_{1}$ and $\Theta_{2}$ are chosen, keeping $\delta_{X}=\delta_{Y}=0.035$ radians ([math]), such that the approximate proportions of zeros in the response and covariate are given by $(0.15,0.15)$, $(0.10,0.10)$, $(0.10,0)$, $(0,0.10)$ and $(0,0)$, respectively. We consider the priors as discussed in the Subsection 2.2 with $a_{\rho_{i}}=2$ for $i=1,2$. We generate 100,000 samples from the posterior distributions of the associated model parameters using the MCMC algorithm and find the posterior mean and standard deviation of the linear parameters based on every 10th iterate discarding the first 40,000 iterations as burn-in. Similarly, we compute the posterior circular mean and circular dispersion for the circular parameters. This is repeated 500 times and the average estimates are reported in Tables 1-5. Note that, the circular dispersion (c.d.) for $n$ circular observations $\phi_{1},\ldots\phi_{n}$ is given by $1-\bar{R}$, where $\bar{R}=\frac{||\sum_{i=1}^{n}z_{i}||}{n}$ and $z_{i}=(\cos\phi_{i},\sin\phi_{i})$ for each $i=1,\ldots,n$. We also report the coverage probability (CP) corresponding to all the parameters. As expected, the standard deviation (s.d.) and c.d. decrease as the sample size increases (See Tables 1-5). In this limited simulation study, our proposed method seems to perform reasonably well.

3.1 Model Comparison and Sensitivity Analysis

For the purpose of comparison, we consider three different models, Model I, Model II and Model III, given by (4) and (5), (2) and (3), and (4) and (3), respectively. Note that Model I accounts for zero-inflation in both the response and the covariate, however, Model II does not account for zero-inflation, and Model III accounts for zero-inflation in the response only. We first compare the performance of parameter estimates associated with Model I, Model II and Model III when the data is generated from Model I with sample size $n=100$ and $\delta_{X}=\delta_{Y}=0.070$ radians ([math]). We consider the priors as discussed in the Subsection 2.2 with $a_{\rho_{i}}=2$ for $i=1,2$. The averages of the estimates over 500 replications are reported in Table 6. It is observed that the performance of most of the estimates based on Model II and Model III are biased compared to those of Model I. The estimates of the parameters associated with the regression of $X$ on $V$ based on Model II and Model III have large bias and very low CP. It is important to note that the estimates associated with the regression of $Y$ on $X$ based on Model III, which incorporates zero-inflation in the response, are biased compared to those of Model I that also incorporates zero-inflation in the covariate. These observations clearly indicate that both Model II and Model III are inadequate to model the underlying relationship in the data when there is significant zero-inflation in both response and covariate.

In order to carry out a sensitivity analysis, we consider a misspecified simulation model, where $\theta_{Y}$ and $\theta_{X}$ are first generated as before and then contaminated with randomly occurring zeros with probability $p$. Note that the aforementioned model reduces to the proposed two-stage circular-circular regression model, given by (4)-(5), for $p=0$. We generate data of sample size $n=100$ for two different choices of $p$ with $\delta=0.070$, and compare the different models with respect to bias and CP based on 500 replications. The results are presented in Table 7 and 8. It is observed that the performance of the estimates for all the models deteriorate with respect to both bias and CP as $p$ increases. The bias and CP for most of the parameters based on Model I are considerably better compared to Model II and III. In particular, the CP of the parameters associated with the regression of $X$ on $V$ based on Model II and Model III are close to zero. Even with the higher value of $p$, the estimates obtained from Model I are robust compared to those of Model II and Model III.

Next, we compare the run-time of different models per dataset in Table 9, under the simulation settings presented in Tables 6-8, using a computer equipped with 32 GB RAM and 2.59 GHz Intel Core(TM) i7-8850H processor. As expected, the computing time increases with the increase in the proportion of zeros for Model I and Model III. We also observed the same pattern for the simulation settings corresponding to Tables 1-5. The difference among the run-times for the three different models is less than 3 minutes, even with large proportion of zeros. As it turns out, Model I seems to have an edge concerning the trade-off between computational time and model performance in presence of zero-inflation.

4 Case Studies

As discussed before, applications of circular-circular regression are commonly found in the domain of orthopaedics and ophthalmology. In this Section, we consider two different case studies from each of these application areas. Significant amount of zero-inflation is observed in both ‘Cataract Surgery Data’ and ‘Abnormal Gait Data’. The results from analysis of these datasets based on the proposed methodology along with the practical implications are discussed below.

4.1 Analysis of Cataract Surgery Data

As discussed in Section 1, we consider a dataset on astigmatism observed at three different inspections from a study on cataract surgery conducted at Disha Eye Hospital and Research Center. In the proposed two-stage setup, we first consider the response ($\theta_{Y}$) and the covariate ($\theta_{X}$) as the axis of astigmatism after 15 days and 7 days of the surgery, respectively. Among the 54 observations, there are $31\%$ and $35\%$ zeros in response and covariate, respectively. In this study, measurements on the axis of astigmatism just after a day of the surgery are also available. Next, we model the axis of astigmatism after 7 days of the surgery ($\theta_{X}$) with covariate as the axis of astigmatism just after a day of surgery ($\theta_{V}$). In order to apply the methodology provided in Section 2, we consider $\delta_{X}=\delta_{Y}=0.035$ radians ([math]) and $a_{\rho_{i}}=2$ for $i=1,2$. This particular choice of $\delta_{X}$ and $\delta_{Y}$ is considered as the original axes of astigmatism are censored in the interval (-$$,$$).

We generate 250,000 samples from the posterior distributions of the associated model parameters using the MCMC algorithm and find the posterior mean/circular mean and s.d./c.d. based on every 10th iterate discarding the first 70,000 iterations as burn-in. The convergence of the chains is monitored graphically and using Geweke’s diagnostic test after taking cosine transformation of the samples for the circular parameters. The mean/circular mean, s.d/c.d. of the parameters and the results for Geweke’s diagnostic test are reported in Appendix A of the Supplementary Material. We also report the $95\%$ highest posterior density (HPD) credible interval. The algorithm for finding the HPD credible interval for circular parameters is provided in Jha, (2017). As expected, the $95\%$ HPD credible interval for both $r_{1}$ and $r_{2}$ does not contain 1. This indicates that $\theta_{Y}$ and $\theta_{X}$ are dependent on $\theta_{X}$ and $\theta_{V}$, respectively.

In order to compare the recovery processes in consecutive weeks, we present spoke plots in Figure 6 and Figure 7 based on the estimates from Model I and its special case, respectively. In the spoke plots, the predicted axes of astigmatism in the outer circle are joined with the corresponding axes of astigmatism a week before in the inner circle. It is important to note that there is a preferred direction for the axis of astigmatism, and it makes practically more sense to consider the equivalence of symmetric points about [math]. It can be seen from the spoke plots based on Model I (See Figure 6), that the predicted responses for symmetric observations around [math] are not symmetric. For example, the predicted mean direction for [math] is rotated towards [math], while for [math], it is rotated away from [math]. However, for the special case, the symmetric points rotate symmetrically towards the preferred direction (See Figure 7). Hence, the equivalence with respect to the preferred direction is maintained. Also, the predicted axes of astigmatism in the second week (astigmatism at day 15 based on day 7) are more attracted towards [math] than those of the first week (astigmatism at day 7 based on day 1) for both the models. This indicates that the recovery in the second week is more than the recovery in the first week. This finding is also supported by the fact that $Pr[\beta_{1}-b_{1}>0|D]=0.96$, based on the estimates obtained from the special case of Model I.

Medical practitioners are more interested in identifying patients whose performances either improve or deteriorate. For this purpose, we consider Model I and compute the posterior predictive distributions of the axis of astigmatism at day 15 and day 7 with three different initial conditions fixed at [math] (normal), [math] (intermediate), and [math] (serious) on day 1. We also compute the posterior predictive distribution of the axis of astigmatism at day 15 based on those three initial conditions on day 7. In Figure 8, we present rose diagrams based on the aforementioned posterior predictive distributions to analyse the recovery of patients over time. In all the cases, it is seen that the predictive distribution is concentrated around [math] if the initial axis of astigmatism is close to [math] (See the top panel of Figure 8). For both intermediate and serious cases, improvement is more prominent in the second week compared to the first week as the posterior predictive distribution shifts more towards [math] in the second week (See middle and bottom panel of Figure 8). For the serious case, the patients who show early sign of recovery improve more at the end of two weeks than those with delayed recovery (See second and third column at the bottom panel in Figure 8).

To compare the predicted values with the observed values, we convert the predicted values in degrees and round off to the nearest integer divisible by 4. In general, the patients not affected by astigmatism after 7 days of the surgery remain unaffected in the near future. In this dataset, 17 patients remain unaffected by astigmatism during the study period and the fitted model identifies all of them correctly. We detect the improvement in 17 out of the 21 patients. However, the deterioration is detected for only 6 out of the 13 patients. Therefore, one can conclude that the patients whose conditions have improved and seem to stay good in the future require less monitoring while the rest of the patients need more frequent monitoring and care.

4.2 Analysis of Abnormal Gait Data

We consider a study on gait control conducted at Bio-Engineering Unit, University of Calcutta. The objective of this study is to analyse the angular movements of the patient’s affected limb with that of a healthy person while keeping the angular movement of the healthy limb fixed at a specified level within the range of natural variation. For this purpose, we consider the response ($\theta_{Y}$) and the covariate ($\theta_{X}$) as the angle of ankle extension of the affected and the unaffected legs, respectively, of the 19 individuals with unilateral orthopaedic impairment. We also jointly model the angle of ankle extension of the non-dominated leg ($\theta_{Y}$) as the response with that of the dominated leg ($\theta_{X}$) as covariate based on the measurements taken from 10 healthy individuals. As mentioned before, more than $50\%$ of the observations are recorded as [math] in both response and covariate. For the purpose of joint modeling, we consider a shared-parameter model replacing $\beta_{0}$ with $\beta_{0}\gamma^{W}$ in (4), where $\gamma\in\{z:z\in\mathbbm{C};|z|=1\}$, and $W$ takes value $1$ for orthopaedically impaired individuals, and [math], otherwise. In order to apply the methodology provided in Section 2, we consider $|b_{1}|=1$, and $\delta_{X}=\delta_{Y}=0.0087$ radians ([math]). We generate 150,000 samples from the posterior distributions of the associated model parameters using the MCMC algorithm and discard the first 50,000 iterations as burn-in. Here also, no convergence issues are found based on Geweke’s diagnostic test. The parameter estimates along with $95\%$ HPD credible interval are reported in Appendix A of the Supplementary Material.

Based on the results, it can be seen from Figure 9 that even the $5\%$ confidence ellipse of the posterior distribution of $(Re(\beta_{1}),Im(\beta_{1}))$ contains $(0,0)$. This indicates that the rotation model is sufficient to explain the regression for this particular dataset under consideration. To study the difference between healthy and orthopaedically impaired person, we plot the posterior density of $\theta_{\gamma}=\arg(\gamma)$ in Figure 10. It is clearly visible that the bulk of the distribution is near 0. This is also supported from Figure 11, and by the fact that $95\%$ HPD credible interval of $\theta_{\gamma}=\arg(\gamma)$ includes [math] (see Appendix A of the Supplementary Material). This indicates that the angle of ankle extension of the affected leg of an orthopaedically impaired patient is not significantly different from that of the non-dominating leg of a healthy person when the angle of extension of the other leg are fixed at the same level. Moreover, the $95\%$ HPD credible interval of $\theta_{01}$ also contains zero, which indicates no significant difference between the angles of ankle extension for both legs for healthy individuals as well as the orthopaedically impaired patients. This is also visible from the posterior predictive distributions of the angle of ankle extension of the affected/non-dominating leg where the modes of the posterior predictive distributions are very close to the respective angle of extension of the healthy/dominating leg (see Figure 12).

5 Discussion

In this paper, a Bayesian methodology has been developed for a circular-circular regression model with point-accumulation in covariate and response, unlike the existing frequentist methods that only model the cases with point-accumulation in the response variable. Circular-circular regression is not well studied from a Bayesian perspective. This paper makes an attempt in that direction and possibly for the first time Bayesian estimation is proposed for the Möbius transformation based circular-circular regression model. The proposed method is applied to analyse two real datasets on post operative astigmatism and abnormal gait due to orthopaedic impairments. Our analysis provides interesting insights and predictions associated with the recovery of patients. Findings from the analysis of cataract surgery data may help to reduce post-operative trauma and improve the patient’s experience after cataract surgery. Similarly, bio-engineers can exploit the analytical findings from abnormal gait data to develop efficient prosthesis.

In general, the methodology is applicable for conventional circular-circular regression with or without point-accumulation in response and/or covariate. Also, one can implement the proposed methodology with other choices of link function and/or choices of angular error distributions. As the latent variables involved in the modeling are continuous, one can easily modify the proposed methodology for a discrete response variable with multiple points of accumulation. Therefore, the scope of the proposed method goes far beyond the particular case studies under consideration. For example, one can model wind direction data where the recorded measurements are among some pre-specified discrete directions (e.g. North, West, East, and South). Although we have not considered missing data in our current analysis, a simple data augmentation technique has to be incorporated into the proposed methodology when the missingness is ignorable. Another possible direction of future research is to develop a bias-reduced estimation methodology for zero-inflated circular-circular regression using the approach proposed by Schwartz and Giles, (2016). In some real-life scenarios, a circular response may depend on both linear and circular covariates. It will be an interesting problem to model such data and develop associated estimation methodology. As the Möbius transformation based circular-circular regression model is not readily extendable to the case with linear covariates, it can be considered as future work.

Acknowledgement

The work of Dr. Prajamitra Bhuyan was supported in part by the Lloyd’s Register Foundation programme on data-centric engineering at the Alan Turing Institute, UK. The authors are thankful to Prof. Anup Dewanji, Dr. Arnab Chakraborty, Prof. Debasis Sengupta, Dr. Jayabrata Biswas, Dr. Sourabh Bhattacharya, and Mr. Sudipta Kundu for many helpful comments and suggestions.