$\displaystyle=-\frac{8\nu_{ab}\sigma(\tau-1)\sqrt{\tau}\left(10\sigma^{2}-2\sigma\tau+3\tau^{2}\right)}{25\sqrt{\pi}(\sigma+\tau)^{7/2}},$
(101h)

where $\tau=T_{a}/T_{b}$, $\sigma=m_{a}/m_{b}$ are the temperature and mass
ratios of the colliding species, respectively. Similarly, the non-vanishing
terms  $\Delta\mathcal{C}_{ab,p^{\prime}j^{\prime}}^{Fpj}$ in the gyro-moment
expansion of the field component $\Delta\mathcal{C}_{ab}^{Fpj}$ (see Eq. (74))
are

$\displaystyle\Delta\mathcal{C}_{ab30}^{F10}$
$\displaystyle=-\frac{4\sqrt{\frac{6}{\pi}}\sigma^{3/2}\tau\left(-\sqrt{(\sigma+1)(\sigma+\tau)}+\sigma+1\right)}{5(\sigma+\tau)^{5/2}},$
(102a) $\displaystyle\Delta\mathcal{C}_{ab11}^{F10}$
$\displaystyle=\frac{8\sigma\sqrt{\sigma\tau}\left(-\sqrt{(\sigma+1)\tau^{3}}+\sigma\left(\sqrt{\tau(\sigma+\tau)}-\sqrt{(\sigma+1)\tau}\right)+\sqrt{\tau(\sigma+\tau)}\right)}{5\sqrt{\pi}(\sigma+\tau)^{3}},$
(102b) $\displaystyle\Delta\mathcal{C}_{ab10}^{F30}$
$\displaystyle=\frac{4\sqrt{\frac{6}{\pi}}\sqrt{\sigma}\tau\left(\sqrt{(\sigma+1)\tau^{3}}+\tau\left(-\sqrt{\sigma+\tau}\right)+\sigma\left(-3\tau\sqrt{\sigma+\tau}+\sqrt{(\sigma+1)\tau}+2\sqrt{\sigma+\tau}\right)\right)}{5(\sigma+\tau)^{3}},$
(102c) $\displaystyle\Delta\mathcal{C}_{ab30}^{F30}$
$\displaystyle=\frac{36\sigma^{3/2}\tau\left(\sigma\left(5\tau-\sqrt{\tau}-2\right)-\left(\sqrt{\tau}-3\right)\tau\right)}{25\sqrt{\pi}(\sigma+\tau)^{7/2}},$
(102d) $\displaystyle\Delta\mathcal{C}_{ab11}^{F30}$
$\displaystyle=\frac{12\sqrt{\frac{6}{\pi}}\sigma^{3/2}\tau\left(\sigma\left(-5\tau+\sqrt{\tau}+2\right)+\left(\sqrt{\tau}-3\right)\tau\right)}{25(\sigma+\tau)^{7/2}},$
(102e) $\displaystyle\Delta\mathcal{C}_{ab10}^{F11}$
$\displaystyle=-\frac{8\sqrt{\sigma\tau}\left(\sigma\left(-3\sqrt{\tau^{3}(\sigma+\tau)}+\sqrt{\sigma+1}\tau+2\sqrt{\tau(\sigma+\tau)}\right)-\sqrt{\tau^{3}(\sigma+\tau)}+\sqrt{\sigma+1}\tau^{2}\right)}{5\sqrt{\pi}(\sigma+\tau)^{3}},$
(102f) $\displaystyle\Delta\mathcal{C}_{ab30}^{F11}$
$\displaystyle=\Delta\mathcal{C}_{ab11}^{F30},$ (102g)
$\displaystyle\Delta\mathcal{C}_{ab11}^{F11}$
$\displaystyle=\frac{24\sigma^{3/2}\tau\left(\sigma\left(5\tau-\sqrt{\tau}-2\right)-\left(\sqrt{\tau}-3\right)\tau\right)}{25\sqrt{\pi}(\sigma+\tau)^{7/2}}.$
(102h)

We notice that, from the conservation of particle in Eq. (2a), it follows
$\Delta\mathcal{C}_{ab,p^{\prime}j^{\prime}}^{T00}=\Delta\mathcal{C}_{ab,p^{\prime}j^{\prime}}^{F00}$.
Finally, we give the non-vanishing terms in the case of like-species
collision. Because the test component of the OS operator reduces to the
Coulomb collision operator when $a=b$, $\Delta
C_{aap^{\prime}j^{\prime}}^{Tpj}=0$, while the non-vanishing terms of the
field component $\Delta C_{aap^{\prime}j^{\prime}}^{Fpj}$ are given by

$\displaystyle\Delta C_{aa30}^{F30}$
$\displaystyle=\frac{9}{25}\sqrt{\frac{2}{\pi}},$ (103a) $\displaystyle\Delta
C_{aa11}^{F30}$ $\displaystyle=-\frac{6}{25}\sqrt{\frac{3}{\pi}},$ (103b)
$\displaystyle\Delta C_{aa11}^{F11}$
$\displaystyle=\frac{6}{{25}}\sqrt{\frac{2}{\pi}}.$ (103c)

We also provide the non-vanishing lowest order terms in the gyro-moment
expansion of the OS collision operator, that are given by

$\displaystyle\mathcal{C}^{ST10}_{ab10}$
$\displaystyle=-\frac{8\nu_{ab}(\sigma+1)}{3\sqrt{\pi}}\left(\frac{\tau}{\sigma+\tau}\right)^{3/2},$
(104a) $\displaystyle\mathcal{C}^{ST10}_{ab30}$
$\displaystyle=\frac{4\nu_{ab}}{5}\sqrt{\frac{6}{\pi}}\frac{\sqrt{\sigma+1}\tau^{2}}{(\sigma+\tau)^{2}},$
(104b) $\displaystyle\mathcal{C}^{ST10}_{ab11}$
$\displaystyle=-\frac{8\nu_{ab}\sqrt{\sigma+1}\tau^{2}}{5\sqrt{\pi}(\sigma+\tau)^{2}},$
(104c) $\displaystyle\mathcal{C}^{ST20}_{ab20}$
$\displaystyle=-\frac{16\nu_{ab}\sqrt{\tau}\left(5\sigma^{2}(\tau+2)+21\sigma\tau+6\tau^{2}\right)}{45\sqrt{\pi}(\sigma+\tau)^{5/2}},$
(104d) $\displaystyle\mathcal{C}^{ST20}_{ab01}$
$\displaystyle=-\frac{16\nu_{ab}\sqrt{\frac{2}{\pi}}\sqrt{\tau}\left(-5\sigma^{2}(\tau-1)+3\sigma\tau+3\tau^{2}\right)}{45(\sigma+\tau)^{5/2}},$
(104e) $\displaystyle\mathcal{C}^{ST01}_{ab01}$
$\displaystyle=-\frac{16\nu_{ab}\sqrt{\tau}\left(5\sigma^{2}+2(5\sigma+9)\sigma\tau+3\tau^{2}\right)}{45\sqrt{\pi}(\sigma+\tau)^{5/2}},$
(104f) $\displaystyle\mathcal{C}^{ST30}_{ab30}$
$\displaystyle=-\frac{4\nu_{ab}\sqrt{\tau}\left(70\sigma^{2}+56\sigma\tau+31\tau^{2}\right)}{35\sqrt{\pi}(\sigma+\tau)^{5/2}},$
(104g) $\displaystyle\mathcal{C}^{ST30}_{ab11}$
$\displaystyle=-\frac{4\nu_{ab}\sqrt{\frac{2}{3\pi}}\tau^{3/2}(28\sigma+\tau)}{35(\sigma+\tau)^{5/2}},$
(104h) $\displaystyle\mathcal{C}^{ST11}_{ab11}$
$\displaystyle=-\frac{8\nu_{ab}\sqrt{\tau}\left(105\sigma^{2}+98\sigma\tau+47\tau^{2}\right)}{105\sqrt{\pi}(\sigma+\tau)^{5/2}},$
(104i)

with $\mathcal{C}^{STpj}_{ablk}=\mathcal{C}^{STlk}_{abpj}$, and

$\displaystyle\mathcal{C}^{SF10}_{ab10}$
$\displaystyle=\frac{8\nu_{ab}\sigma(\sigma+1)\tau}{3\sqrt{\pi}\sqrt{\sigma(\sigma+\tau)^{3}}},$
(105a) $\displaystyle\mathcal{C}^{SF10}_{ab30}$
$\displaystyle=-\frac{4\nu_{ab}\sqrt{\frac{6}{\pi}}\sqrt{\sigma^{3}(\sigma+1)}\tau}{5(\sigma+\tau)^{2}},$
(105b) $\displaystyle\mathcal{C}^{SF10}_{ab11}$
$\displaystyle=\frac{8\nu_{ab}\sqrt{\sigma^{3}(\sigma+1)}\tau}{5\sqrt{\pi}(\sigma+\tau)^{2}},$
(105c) $\displaystyle\mathcal{C}^{SF20}_{ab20}$
$\displaystyle=\frac{16\nu_{ab}\sigma(\sigma+1)\tau^{3/2}}{9\sqrt{\pi}(\sigma+\tau)^{5/2}},$
(105d) $\displaystyle\mathcal{C}^{SF20}_{ab01}$
$\displaystyle=-\frac{16\nu_{ab}\sqrt{\frac{2}{\pi}}\sigma(\sigma+1)\tau^{3/2}}{9(\sigma+\tau)^{5/2}},$
(105e) $\displaystyle\mathcal{C}^{SF01}_{ab20}$
$\displaystyle=-\frac{16\nu_{ab}\sqrt{\frac{2}{\pi}}\sigma(\sigma+1)\tau^{3/2}}{9(\sigma+\tau)^{5/2}},$
(105f) $\displaystyle\mathcal{C}^{SF01}_{ab01}$
$\displaystyle=\frac{32\nu_{ab}\sigma(\sigma+1)\tau^{3/2}}{9\sqrt{\pi}(\sigma+\tau)^{5/2}},$
(105g) $\displaystyle\mathcal{C}^{SF30}_{ab10}$
$\displaystyle=-\frac{4\nu_{ab}\sqrt{\frac{6}{\pi}}\sqrt{\sigma(\sigma+1)\tau^{3}}}{5(\sigma+\tau)^{2}},$
(105h) $\displaystyle\mathcal{C}^{SF30}_{ab30}$
$\displaystyle=\frac{36\nu_{ab}(\sigma\tau)^{3/2}}{25\sqrt{\pi}(\sigma+\tau)^{5/2}},$
(105i) $\displaystyle\mathcal{C}^{SF30}_{ab11}$
$\displaystyle=-\frac{12\nu_{ab}\sqrt{\frac{6}{\pi}}(\sigma\tau)^{3/2}}{25(\sigma+\tau)^{5/2}},$
(105j) $\displaystyle\mathcal{C}^{SF11}_{ab10}$
$\displaystyle=\frac{8\nu_{ab}\sqrt{\sigma(\sigma+1)\tau^{3}}}{5\sqrt{\pi}(\sigma+\tau)^{2}},$
(105k) $\displaystyle\mathcal{C}^{SF11}_{ab30}$
$\displaystyle=-\frac{12\nu_{ab}\sqrt{\frac{6}{\pi}}(\sigma\tau)^{3/2}}{25(\sigma+\tau)^{5/2}},$
(105l) $\displaystyle\mathcal{C}^{SF11}_{ab11}$
$\displaystyle=\frac{24\nu_{ab}(\sigma\tau)^{3/2}}{25\sqrt{\pi}(\sigma+\tau)^{5/2}}.$
(105m)

for the test and field components, respectively.  In the case of like-species
collision, the non-vanishing terms are, given that

$\displaystyle\mathcal{C}^{S20}_{aa20}$
$\displaystyle=-\frac{64}{45}\sqrt{\frac{2}{\pi}},$ (106a)
$\displaystyle\mathcal{C}^{S20}_{aa01}$
$\displaystyle=-\frac{64}{45}\sqrt{\frac{1}{\pi}},$ (106b)
$\displaystyle\mathcal{C}^{S01}_{aa01}$
$\displaystyle=-\frac{32}{45}\sqrt{\frac{2}{\pi}},$ (106c)
$\displaystyle\mathcal{C}^{S30}_{aa30}$
$\displaystyle=-\frac{361}{175}\sqrt{\frac{2}{\pi}},$ (106d)
$\displaystyle\mathcal{C}^{S30}_{aa11}$
$\displaystyle=-\frac{208}{175}\sqrt{\frac{1}{3\pi}},$ (106e)
$\displaystyle\mathcal{C}^{S11}_{aa11}$
$\displaystyle=-\frac{1187}{525}\sqrt{\frac{2}{\pi}}.$ (106f)

Finally and for completeness, we report the non-vanishing coefficients of the
DK Coulomb collision operator of the test $\mathcal{C}_{ab}^{LT}$ and the
field $\mathcal{C}_{ab}^{LF}$ components respectively, i.e.

$\displaystyle\mathcal{C}_{ab10}^{LT10}$
$\displaystyle=-\frac{8\nu_{ab}(\sigma+1)}{3\sqrt{\pi}}\left(\frac{\tau}{\sigma+\tau}\right)^{3/2},$
(107a) $\displaystyle\mathcal{C}_{ab30}^{LT10}$
$\displaystyle=\frac{4}{5}\sqrt{\frac{6}{\pi}}(\sigma+1)\left(\frac{\tau}{\sigma+\tau}\right)^{5/2},$
(107b) $\displaystyle\mathcal{C}_{ab11}^{LT10}$
$\displaystyle=-\frac{8\nu_{ab}(\sigma+1)}{{5\sqrt{\pi}}\left(\frac{\tau}{\sigma+\tau}\right)^{5/2}},$
(107c) $\displaystyle\mathcal{C}_{ab00}^{LT20}$
$\displaystyle=-\frac{8\nu_{ab}\sqrt{\frac{2}{\pi}}(\tau-1)\sqrt{\frac{\tau}{\sigma}}}{3\left(\frac{\sigma+\tau}{\sigma}\right)^{3/2}},$
(107d) $\displaystyle\mathcal{C}_{ab20}^{LT20}$
$\displaystyle=-\frac{\nu_{ab}8\left(\frac{\tau}{\sigma}\right)^{3/2}(\sigma(10\sigma+\tau+13)+4\tau)}{15\sqrt{\pi}\sigma\left(\frac{\sigma+\tau}{\sigma}\right)^{5/2}},$
(107e) $\displaystyle\mathcal{C}_{ab01}^{LT20}$
$\displaystyle=\frac{8\nu_{ab}\sqrt{\frac{2}{\pi}}\left(\frac{\tau}{\sigma}\right)^{3/2}(-3\sigma\tau+\sigma-2\tau)}{15\sigma\left(\frac{\sigma+\tau}{\sigma}\right)^{5/2}},$
(107f) $\displaystyle\mathcal{C}_{ab00}^{LT01}$
$\displaystyle=\frac{16\nu_{ab}(\tau-1)\sqrt{\frac{\tau}{\sigma}}}{3\sqrt{\pi}\left(\frac{\sigma+\tau}{\sigma}\right)^{3/2}},$
(107g) $\displaystyle\mathcal{C}_{ab20}^{LT01}$
$\displaystyle=\frac{8\sqrt{\frac{2}{\pi}}\left(\frac{\tau}{\sigma}\right)^{3/2}(-3\sigma\tau+\sigma-2\tau)}{15\sigma\left(\frac{\sigma+\tau}{\sigma}\right)^{5/2}},$
(107h) $\displaystyle\mathcal{C}_{ab01}^{LT01}$
$\displaystyle=-\frac{16\left(\frac{\tau}{\sigma}\right)^{3/2}(\sigma(5\sigma-\tau+7)+\tau)}{15\sqrt{\pi}\sigma\left(\frac{\sigma+\tau}{\sigma}\right)^{5/2}},$
(107i) $\displaystyle\mathcal{C}_{ab10}^{LT30}$
$\displaystyle=-\frac{4\nu_{ab}}{5\sigma^{2}}\sqrt{\frac{2}{3\pi}}\sqrt{\frac{\tau}{\sigma}}\left(\frac{\sigma}{\sigma+\tau}\right)^{5/2}\left(10\sigma^{2}(\tau-1)+\sigma(\tau-4)\tau-3\tau^{2}\right),$
(107j) $\displaystyle\mathcal{C}_{ab30}^{LT30}$
$\displaystyle=-\frac{4\nu_{ab}\left(\frac{\tau}{\sigma+\tau}\right)^{3/2}\left(14\sigma^{2}(\tau+8)+70\sigma^{3}+\sigma\tau(19\tau+68)+31\tau^{2}\right)}{35\sqrt{\pi}(\sigma+\tau)^{2}},$
(107k) $\displaystyle\mathcal{C}_{ab11}^{LT30}$
$\displaystyle=\frac{4\nu_{ab}\sqrt{\frac{2}{3\pi}}\left(\frac{\tau}{\sigma+\tau}\right)^{3/2}\left(\sigma^{2}(14-42\tau)+\sigma\tau(3\tau-32)-\tau^{2}\right)}{35(\sigma+\tau)^{2}},$
(107l) $\displaystyle\mathcal{C}_{ab10}^{LT11}$
$\displaystyle=\frac{8\nu_{ab}\sqrt{\frac{\tau}{\sigma}}\left(10\sigma^{2}(\tau-1)+\sigma(\tau-4)\tau-3\tau^{2}\right)}{15\sqrt{\pi}\sigma^{2}\left(\frac{\sigma+\tau}{\sigma}\right)^{5/2}},$
(107m) $\displaystyle\mathcal{C}_{ab30}^{LT11}$
$\displaystyle=\frac{4\nu_{ab}\sqrt{\frac{2}{3\pi}}\left(\frac{\tau}{\sigma+\tau}\right)^{3/2}\left(\sigma^{2}(14-42\tau)+\sigma\tau(3\tau-32)-\tau^{2}\right)}{35(\sigma+\tau)^{2}},$
(107n) $\displaystyle\mathcal{C}_{ab11}^{LT11}$
$\displaystyle=-\frac{8\nu_{ab}\left(\frac{\tau}{\sigma+\tau}\right)^{3/2}\left(7(15\sigma+23)\sigma^{2}+(27\sigma+47)\tau^{2}+2(21\sigma+59)\sigma\tau\right)}{105\sqrt{\pi}(\sigma+\tau)^{2}},$
(107o)

and

$\displaystyle\mathcal{C}_{ab10}^{LF10}$
$\displaystyle=\frac{8\nu_{ab}(\sigma+1)\tau}{3\sqrt{\pi}\sigma}\left(\frac{\sigma}{\sigma+\tau}\right)^{3/2},$
(108a) $\displaystyle\mathcal{C}_{ab30}^{LF10}$
$\displaystyle=-\frac{4\nu_{ab}\sqrt{\frac{6}{\pi}}(\sigma+1)\tau}{5\sigma}\left(\frac{\sigma}{\sigma+\tau}\right)^{5/2},$
(108b) $\displaystyle\mathcal{C}_{ab11}^{LF10}$
$\displaystyle=\frac{8\nu_{ab}(\sigma+1)\tau}{5\sqrt{\pi}\sigma}\left(\frac{\sigma}{\sigma+\tau}\right)^{5/2},$
(108c) $\displaystyle\mathcal{C}_{ab00}^{LF20}$
$\displaystyle=-\frac{8\nu_{ab}\sqrt{\frac{2}{\pi}}(\tau-1)\sqrt{\frac{\tau}{\sigma}}}{3\left(\frac{\sigma+\tau}{\sigma}\right)^{3/2}},$
(108d) $\displaystyle\mathcal{C}_{ab20}^{LF20}$
$\displaystyle=\frac{8\nu_{ab}\sqrt{\frac{\tau}{\sigma}}(\sigma(3\tau-1)+2\tau)}{5\sqrt{\pi}\sigma\left(\frac{\sigma+\tau}{\sigma}\right)^{5/2}},$
(108e) $\displaystyle\mathcal{C}_{ab01}^{LF20}$
$\displaystyle=\frac{8\nu_{ab}\sqrt{\frac{2}{\pi}}\sqrt{\frac{\tau}{\sigma}}(-3\sigma\tau+\sigma-2\tau)}{15\sigma\left(\frac{\sigma+\tau}{\sigma}\right)^{5/2}},$
(108f) $\displaystyle\mathcal{C}_{ab00}^{LF01}$
$\displaystyle=\frac{16(\tau-1)\sqrt{\frac{\tau}{\sigma}}}{3\sqrt{\pi}\left(\frac{\sigma+\tau}{\sigma}\right)^{3/2}},$
(108g) $\displaystyle\mathcal{C}_{ab20}^{LF01}$
$\displaystyle=\frac{8\sqrt{\frac{2}{\pi}}\sqrt{\frac{\tau}{\sigma}}(-3\sigma\tau+\sigma-2\tau)}{15\sigma\left(\frac{\sigma+\tau}{\sigma}\right)^{5/2}},$
(108h) $\displaystyle\mathcal{C}_{ab01}^{LF01}$
$\displaystyle=\frac{32\sqrt{\frac{\tau}{\sigma}}(\sigma(3\tau-1)+2\tau)}{15\sqrt{\pi}\sigma\left(\frac{\sigma+\tau}{\sigma}\right)^{5/2}},$
(108i) $\displaystyle\mathcal{C}_{ab10}^{LF30}$
$\displaystyle=-\frac{4\nu_{ab}\sqrt{\frac{6}{\pi}}\tau(\sigma(3\tau-2)+\tau)}{5\sigma^{2}\left(\frac{\sigma+\tau}{\sigma}\right)^{5/2}},$
(108j) $\displaystyle\mathcal{C}_{ab30}^{LF30}$
$\displaystyle=\frac{12\nu_{ab}\tau(\sigma(5\tau-2)+3\tau)}{7\sqrt{\pi}(\sigma+\tau)^{2}\left(\frac{\sigma+\tau}{\sigma}\right)^{3/2}},$
(108k) $\displaystyle\mathcal{C}_{ab11}^{LF30}$
$\displaystyle=-\frac{12\nu_{ab}\sqrt{\frac{6}{\pi}}\tau(\sigma(5\tau-2)+3\tau)}{35(\sigma+\tau)^{2}\left(\frac{\sigma+\tau}{\sigma}\right)^{3/2}},$
(108l) $\displaystyle\mathcal{C}_{ab10}^{LF11}$
$\displaystyle=\frac{8\nu_{ab}\tau(\sigma(3\tau-2)+\tau)}{5\sqrt{\pi}\sigma^{2}\left(\frac{\sigma+\tau}{\sigma}\right)^{5/2}},$
(108m) $\displaystyle\mathcal{C}_{ab30}^{LF11}$
$\displaystyle=-\frac{12\nu_{ab}\sqrt{\frac{6}{\pi}}\tau(\sigma(5\tau-2)+3\tau)}{35(\sigma+\tau)^{2}\left(\frac{\sigma+\tau}{\sigma}\right)^{3/2}},$
(108n) $\displaystyle\mathcal{C}_{ab11}^{LF11}$
$\displaystyle=\frac{48\nu_{ab}\tau(\sigma(5\tau-2)+3\tau)}{35\sqrt{\pi}(\sigma+\tau)^{2}\left(\frac{\sigma+\tau}{\sigma}\right)^{3/2}},$
(108o)

respectively.  Finally, the non-vanishing terms for like-species collisions
are

$\displaystyle\mathcal{C}_{ab20}^{L20}$
$\displaystyle=-\frac{16}{15}\sqrt{\frac{2}{\pi}},$ (109a)
$\displaystyle\mathcal{C}_{ab01}^{L20}$
$\displaystyle=-\frac{16}{15}\sqrt{\frac{1}{\pi}},$ (109b)
$\displaystyle\mathcal{C}_{ab01}^{L01}$
$\displaystyle=-\frac{8}{15}\sqrt{\frac{2}{\pi}},$ (109c)
$\displaystyle\mathcal{C}_{ab30}^{L30}$
$\displaystyle=-\frac{8}{5}\sqrt{\frac{2}{\pi}},$ (109d)
$\displaystyle\mathcal{C}_{ab11}^{L30}$
$\displaystyle=-\frac{8}{5}\sqrt{\frac{1}{3\pi}},$ (109e)
$\displaystyle\mathcal{C}_{ab11}^{L11}$
$\displaystyle=-\frac{28}{15}\sqrt{\frac{2}{\pi}},$ (109f)

The lowest order gyro-moments of the IS, OS and Coulomb reported above can be
used to obtain reduced-fluid models to study the plasma dynamics in the
Pfirsch-Schlüter regime.

## References

  * Belli and Candy (2017) E. A. Belli and J. Candy, “Implications of advanced collision operators for gyrokinetic simulation,” Plasma Physics and Controlled Fusion 59, 045005 (2017).
  * Barnes  _et al._ (2009) M. Barnes, I. G. Abel, W. Dorland, D. R. Ernst, G. W. Hammett, P. Ricci, B. N. Rogers, A. A. Schekochihin, and T. Tatsuno, “Linearized model Fokker–Planck collision operators for gyrokinetic simulations. ii. numerical implementation and tests,” Physics of Plasmas 16, 072107 (2009).
  * Manas  _et al._ (2015) P. Manas, Y. Camenen, S. Benkadda, W. A. Hornsby, and A. G. Peeters, “Enhanced stabilisation of trapped electron modes by collisional energy scattering in tokamaks,” Physics of Plasmas 22, 062302 (2015).
  * Pan, Ernst, and Crandall (2020) Q. Pan, D. R. Ernst, and P. Crandall, “First implementation of gyrokinetic exact linearized Landau collision operator and comparison with models,” Physics of Plasmas 27, 042307 (2020).
  * Pan, Ernst, and Hatch (2021) Q. Pan, D. R. Ernst, and D. R. Hatch, “Importance of gyrokinetic exact fokker-planck collisions in fusion plasma turbulence,” Physical Review E 103, L051202 (2021).
  * Frei  _et al._ (2021) B. J. Frei, J. Ball, A. C. D. Hoffmann, R. Jorge, P. Ricci, and L. Stenger, “Development of advanced linearized gyrokinetic collision operators using a moment approach,” Journal of Plasma Physics 87, 905870501 (2021).
  * Belli and Candy (2011) E. A. Belli and J. Candy, “Full linearized fokker–planck collisions in neoclassical transport simulations,” Plasma physics and controlled fusion 54, 015015 (2011).
  * Rosenbluth, MacDonald, and Judd (1957) M. N. Rosenbluth, W. M. MacDonald, and D. L. Judd, “Fokker-Planck equation for an inverse-square force,” Physical Review 107, 1 (1957).
  * Dougherty (1964) J. P. Dougherty, “Model Fokker-Planck equation for a plasma and its solution,” The Physics of Fluids 7, 1788 (1964).
  * Hirshman and Sigmar (1976) S. P. Hirshman and D. J. Sigmar, “Approximate Fokker-Planck collision operator for transport theory applications,” The Physics of Fluids 19, 1532 (1976).
  * Abel  _et al._ (2008) I. G. Abel, M. Barnes, S. C. Cowley, W. Dorland, and A. A. Schekochihin, “Linearized model Fokker-Planck collision operators for gyrokinetic simulations. i. theory,” Physics of Plasmas 15, 122509 (2008).
  * Francisquez  _et al._ (2021) M. Francisquez, J. Juno, A. Hakim, G. W. Hammett, and D. R. Ernst, “Improved multispecies dougherty collisions,” arXiv preprint arXiv:2109.10381 (2021).
  * Sugama, Watanabe, and Nunami (2009) H. Sugama, T.-H. Watanabe, and M. Nunami, “Linearized model collision operators for multiple ion species plasmas and gyrokinetic entropy balance equations,” Physics of Plasmas 16, 112503 (2009).
  * Nunami  _et al._ (2015) M. Nunami, M. Nakata, T.-H. Watanabe, and H. Sugama, “Development of linearized collision operator for multiple ion species in gyrokinetic flux-tube simulations,” Plasma and Fusion Research 10, 1403058 (2015).
  * Nakata  _et al._ (2015) M. Nakata, M. Nunami, T.-H. Watanabe, and H. Sugama, “Improved collision operator for plasma kinetic simulations with multi-species ions and electrons,” Computer Physics Communications 197, 61 (2015).
  * Pan, Ernst, and Hatch. (2021) Q. Pan, D. R. Ernst, and D. R. Hatch., “Importance of gyrokinetic exact Fokker-Planck collisions in fusion plasma turbulence,” Physical Review E 103, L051202 (2021).
  * Frei, Hoffmann, and Ricci (2022) B. J. Frei, A. C. D. Hoffmann, and P. Ricci, “Local gyrokinetic collisional theory of the ion-temperature gradient mode,” Journal of Plasma Physics 88, 905880304 (2022).
  * Sugama  _et al._ (2019) H. Sugama, S. Matsuoka, S. Satake, M. Nunami, and T.-H. Watanabe, “Improved linearized model collision operator for the highly collisional regime,” Physics of Plasmas 26, 102108 (2019).
  * Hirshman and Sigmar (1981) S. P. Hirshman and D. J. Sigmar, “Neoclassical transport of impurities in tokamak plasmas,” Nuclear Fusion 21, 1079 (1981).
  * Honda (2014) M. Honda, “Impact of higher-order flows in the moment equations on pfirsch-schlüter friction coefficients,” Physics of Plasmas 21, 092508 (2014).
  * Matsuoka, Sugama, and Idomura (2021) S. Matsuoka, H. Sugama, and Y. Idomura, “Neoclassical transport simulations with an improved model collision operator,” Physics of Plasmas 28, 064501 (2021).
  * Frei, Jorge, and Ricci (2020) B. J. Frei, R. Jorge, and P. Ricci, “A gyrokinetic model for the plasma periphery of tokamak devices,” Journal of Plasma Physics 86, 905860205 (2020).
  * Jenko  _et al._ (2000) F. Jenko, W. Dorland, M. Kotschenreuther, and B. N. Rogers, “Electron temperature gradient driven turbulence,” Physics of plasmas 7, 1904 (2000).
  * Casson  _et al._ (2014) F. J. Casson, C. Angioni, E. A. Belli, R. Bilato, P. Mantica, T. Odstrcil, T. Pütterich, M. Valisa, L. Garzotti, C. Giroud, _et al._ , “Theoretical description of heavy impurity transport and its application to the modelling of tungsten in jet and asdex upgrade,” Plasma Physics and Controlled Fusion 57, 014031 (2014).
  * Landreman and Ernst (2013) M. Landreman and D. R. Ernst, “New velocity-space discretization for continuum kinetic calculations and fokker–planck collisions,” Journal of Computational Physics 243, 130 (2013).
  * Helander and Sigmar (2005) P. Helander and D. J. Sigmar,  _Collisional transport in magnetized plasmas_ , Vol. 4 (Cambridge University Press, 2005).
  * Jorge, Frei, and Ricci (2019) R. Jorge, B. J. Frei, and P. Ricci, “Nonlinear gyrokinetic coulomb collision operator,” Journal of Plasma Physics 85, 905850604 (2019).
  * Snider (2017) R. F. Snider,  _Irreducible Cartesian Tensors_, De Gruyter Studies in Mathematical Physics No. 1 (De Gruyter, 2017).
  * Brizard and Hahm (2007) A. J. Brizard and T. Hahm, “Foundations of nonlinear gyrokinetic theory,” Reviews of modern physics 79, 421 (2007).
  * Jorge, Ricci, and Loureiro (2017) R. Jorge, P. Ricci, and N. F. Loureiro, “A drift-kinetic analytical model for scrape-off layer plasma dynamics at arbitrary collisionality,” Journal of Plasma Physics 83, 905830606 (2017).
  * Gradshteyn and Ryzhik (2014) I. S. Gradshteyn and I. M. Ryzhik,  _Table of integrals, series, and products_ (Academic press, 2014).
  * Fulton  _et al._ (2014) D. P. Fulton, Z. Lin, I. Holod, and Y. Xiao, “Microturbulence in diii-d tokamak pedestal. i. electrostatic instabilities,” Physics of Plasmas 21, 042110 (2014).
  * Kotschenreuther  _et al._ (2017) M. Kotschenreuther, D. R. Hatch, S. Mahajan, P. Valanju, L. Zheng, and X. Liu, “Pedestal transport in h-mode plasmas for fusion gain,” Nuclear Fusion 57, 064001 (2017).
  * Pueschel  _et al._ (2019) M. J. Pueschel, D. R. Hatch, D. R. Ernst, W. Guttenfelder, P. W. Terry, J. Citrin, and J. W. Connor, “On microinstabilities and turbulence in steep-gradient regions of fusion devices,” Plasma Physics and Controlled Fusion 61, 034002 (2019).
  * Xie and Li (2016) H.-S. Xie and B. Li, “Global theory to understand toroidal drift waves in steep gradient,” Physics of Plasmas 23, 082513 (2016).
  * Han  _et al._ (2017) M. Han, Z.-X. Wang, J. Dong, and H. Du, “Multiple ion temperature gradient driven modes in transport barriers,” Nuclear Fusion 57, 046019 (2017).
  * Thomas  _et al._ (2006) D. M. Thomas, A. W. Leonard, T. H. Osborne, R. J. Groebner, W. P. West, and K. H. Burrell, “The effect of plasma collisionality on pedestal current density formation in DIII-D,” Plasma physics and controlled fusion 48, A183 (2006).
  * Lapillonne  _et al._ (2009) X. Lapillonne, S. Brunner, T. Dannert, S. Jolliet, A. Marinoni, L. Villard, T. Görler, F. Jenko, and F. Merz, “Clarifications to the limitations of the s-$\alpha$ equilibrium model for gyrokinetic computations of turbulence,” Physics of Plasmas 16, 032308 (2009).
  * Coppi, Migliuolo, and Pu (1990) B. Coppi, S. Migliuolo, and Y.-K. Pu, “Candidate mode for electron thermal energy transport in multi-kev plasmas,” Physics of Fluids B: Plasma Physics 2, 2322 (1990).
  * Ernst  _et al._ (2005) D. R. Ernst, K. Zeller, N. Basse, L. Lin, M. Porkolab, W. Dorland, and A. Long, “New developments in trapped electron mode turbulence,”  _APS Division of Plasma Physics Meeting Abstracts_ , 50, 235 (2005).
  * Wang  _et al._ (2012) E. Wang, X. Xu, J. Candy, R. J. Groebner, P. B. Snyder, Y. Chen, S. E. Parker, W. Wan, G. Lu, and J. Q. Dong, “Linear gyrokinetic analysis of a diii-d h-mode pedestal near the ideal ballooning threshold,” Nuclear Fusion 52, 103015 (2012).
  * Diamond  _et al._ (2005) P. H. Diamond, S. I. Itoh, K. Itoh, and T. S. Hahm, “Zonal flows in plasma-a review,” Plasma Physics and Controlled Fusion 47, R35 (2005).
  * Rosenbluth and Hinton (1998) M. N. Rosenbluth and F. L. Hinton, “Poloidal flow driven by ion-temperature-gradient turbulence in tokamaks,” Physical Review Letters 80, 724 (1998).
  * Hinton and Rosenbluth (1999) F. L. Hinton and M. N. Rosenbluth, “Dynamics of axisymmetric (E × B) and poloidal flows in tokamaks,” Plasma Physics and Controlled Fusion 41, A653 (1999).
  * Xiao, Catto, and Molvig (2007) Y. Xiao, P. J. Catto, and K. Molvig, “Collisional damping for ion temperature gradient mode driven zonal flow,” Physics of Plasmas 14, 032302 (2007).
  * Spitzer J. and Härm (1953) L. Spitzer J. and R. Härm, “Transport phenomena in a completely ionized gas,” Physical Review 89, 977 (1953).
  * Giacomin  _et al._ (2022) M. Giacomin, A. Pau, P. Ricci, O. Sauter, T. Eich, J. Contributors, _et al._ , “First-principles density limit scaling in tokamaks based on edge turbulent transport and implications for iter,” Physical Review Letters 128, 185003 (2022).
  * Li and Ernst (2011) B. Li and D. R. Ernst, “Gyrokinetic fokker-planck collision operator,” Phys. Rev. Lett. 106, 195002 (2011).
  * Crandall  _et al._ (2020) P. Crandall, D. Jarema, H. Doerk, Q. Pan, G. Merlo, T. Görler, A. B. Navarro, D. Told, M. Maurer, and F. Jenko, “Multi-species collisions for delta-f gyrokinetic simulations: Implementation and verification with gene,” Computer Physics Communications 255, 107360 (2020).


# Doubly Robust-Based Generalized Estimating Equations for the Analysis of
Longitudinal Ordinal Missing Data

José Luiz P. da Silva, Enrico A. Colosimo, Fábio N. Demarqui
Departament of Statistics
Universidade Federal de Minas Gerais

###### Abstract

Generalized Estimation Equations (GEE) are a well-known method for the
analysis of non-Gaussian longitudinal data. This method has computational
simplicity and marginal parameter interpretation. However, in the presence of
missing data, it is only valid under the strong assumption of missing
completely at random (MCAR). Some corrections can be done when the missing
data mechanism is missing at random (MAR): inverse probability weighting
(WGEE) and multiple imputation (MIGEE). In order to obtain consistent
estimates, it is necessary the correct specification of the weight model for
WGEE or the imputation model for the MIGEE. A recent method combining ideas of
these two approaches has doubly robust property. For consistency, it requires
only the weight or the imputation model to be correct. In this work it is
assumed a proportional odds model and it is proposed a doubly robust estimator
for the analysis of ordinal longitudinal data with intermittently missing
response and covariate under the MAR mechanism. Simulation results revealed
better performance of the proposed method compared to WGEE and MIGEE. The
method is applied to a data set related to Analgesia Pain in Childbirth study.

Keywords: Missing at random; Multiple imputation; Proportional Odds Model;
Weighted GEE.

##  1 Introduction

The Generalized Estimating Equation (GEE) method (Liang & Zeger,, 1986) is one
of the most popular approaches for the analysis of non-Gaussian correlated
data. Its main advantage resides in the fact that one is only required to
specify correctly the mean structure of the response for the parameter
estimator to be consistent and asymptotically normal. In its basic formulation
the association parameters among repeated measures were taken as nuisance
parameters. GEE method has computational simplicity (there is no need of
dealing with complex, and in some cases, intractable likelihoods) and it
further allows populational-averaged interpretation of the parameter of
interest.

It is very common for sets of longitudinal data to be incomplete, in the sense
that not all planned observations are actually observed. This problem is
pervasive in longitudinal data because nonresponse can occur any time from the
beginning of the study. Two patterns of missing data can be observed for the
response: (1) dropout, when a subject leaves the study prematurely for reasons
beyond the control of the investigator, leading to a monotone pattern of
nonresponse, or (2) intermittent nonresponse, in which a subject returns to
the study after some occasions of nonresponse. Covariates may also be missing,
leading to limitations in data analysis. In the presence of missing data,
three issues are of main concern: (1) potential serious bias due to systematic
differences between the observed data and the missing data, (2) complications
in data handling and statistical inferences, and (3) loss of efficiency.
Therefore, in order to make valid inferences it is fundamental to know the
missing data mechanism generating the nonresponse and how to handle it.

Little & Rubin, (1987) provided a formal framework for dealing with missing
data by defining the commonly adopted taxonomy of missing data mechanisms. A
nonresponse process is said to be _missing completely at random_ (MCAR) if
missingness is independent of both unobserved and observed data, and _missing
at random_ (MAR) if, conditionally on the observed data, the missingness is
independent of the unobserved data. When the nonresponse process depends on
unobserved quantities it is said to be _missing not at random_ (MNAR).

When data are incomplete, GEE suffers from its frequentist nature and is, in
its basic form, valid only under MCAR (Liang & Zeger,, 1986). The first effort
to make GEE applicable to the more realistic MAR scenario was Multiple
Imputation (MIGEE), proposed by Little & Rubin, (1987), in which the missing
portions of data are multiply imputed taking into account the uncertainty
associated with the predicted values. The completed data sets are analyzed by
standard methods for complete data, and estimates are combined into a final
analysis. Multiple imputation is detailed in the books by Schafer, (1997),
Little & Rubin, (2002) and Carpenter & Kenward, (2013). Later, Robins et al.,
(1995) proposed the Weighted Generalized Estimating Equations (WGEE), which
consists in weighting each observation by the inverse of the probability of
the data being observed. This method produces consistent estimates provided
the weight model is correctly specified.

Doubly robust estimators (DRGEE) arise as a third generalization of ordinary
GEE to deal with data subject to MAR mechanism. DR methods have received
increasingly attention in the literature in the last decade (see Carpenter et
al., (2006), Bang & Robins, (2005), Tsiatis, (2006), Seaman & Copas, (2009),
Chen & Zhou, (2011)). The main idea is to supplement the WGEE with a
predictive model for the missing quantities conditional on the observed ones.
Doubly Robust method requires only the dropout or the conditional model to be
correctly specified in order to provide consistent estimates. In the analysis
of longitudinal binary data, DR methods have been applied by Seaman & Copas,
(2009), Birhanu et al., (2011), for missing responses and by Chen & Zhou,
(2011) for intermittently missing response and a single missing covariate.

Literature of GEE for missing data is comparatively scarce for longitudinal
ordinal response. In Toledano & Gatsonis, (1999), the authors used a weighted
GEE method to accommodate arbitrary patterns of a MCAR missing response and
missingness in a key covariate subject to a MAR mechanism. A recent paper from
Donneau et al., (2014a) compared through a simulation study two multiple
imputation methods (multivariate normal imputation and ordinal imputation
regression) for longitudinal ordinal data subject to dropout. In another paper
the same authors compared joint modeling and fully conditional specification
approaches for non-monotone missingness (Donneau et al.,, 2014b). The above
mentioned papers used single robust versions of GEE and they have treated only
a missing MAR baseline covariate or missing MAR response. Thus the use of a DR
GEE method for ordinal data with simultaneously intermittently missing
response and missing covariate has been in need of further development.

This work was motivated by the Analgesia in Childbirth study which was
conducted in Minas Gerais state, Brazil. The main objective of that study was
to compare two techniques of analgesia for labor pain in $49$ patients. The
response, pain intensity, was subjectively assessed by each patient, and
various clinical covariates were taken until delivery. Response and a
particular covariate (consumption of oxytocin) were missing for some patients
and the MAR mechanism seems to be a reasonable assumption for this data.

In the current paper, it is proposed a doubly robust approach for the analysis
of longitudinal ordinal data with intermittently missing response and
covariate that is MAR. The proposed methodology, to our knowledge new to the
GEE literature of missing ordinal data, can be used for handling arbitrary
patterns of missing data in the ordinal response and missingness in a key
covariate, as those frequently arising in medical studies.

The paper is organized as follows. In Section 2 are defined the notation for
GEE with fully observed data and missing data mechanisms. Section 3 outlines
WGEE and MIGEE approaches. The proposed methodology is established in Section
4. A simulation study is presented in Section 5, in which the finite-sample
biases and standard errors are compared for the standard GEE, MIGEE, WGEE and
doubly robust versions. Data arising from the Analgesia en Childbirth study
are analyzed in Section 6. Paper ends with a discussion and future directions
in Section 7.

##  2 GEE for Complete Data and Missing Data Assumptions

In this section it is introduced generalized estimating equations for the
analysis of fully observed ordinal data. Section 2.1 establishes the model and
notation for longitudinal ordinal data. Section 2.2 presents a series of
assumptions related to mechanism causing data to be missing and necessary to
be considered in order to build valid estimators.

###  2.1 GEE for Longitudinal Ordinal Response

Let $O_{it}\in\left\\{1,2,\ldots,J\right\\}$ be the ordinal response for
subject $i\ (i=1,\ldots,n)$ at time $t\ (t=1,\ldots,T_{i},\ T_{i}\leq T)$. As
the response has $J$ levels it can be defined $Y_{itj}=I(O_{it}=j)$ for
$j=1,\ldots,J$, where $I(A)$ denotes the indicator function. $Y_{itj}$ is
converted into the equivalent $(J-1)$-variate vector
$\boldsymbol{Y}_{it}=(Y_{it1},\ldots,Y_{it(J-1)})^{T}$ and let
$\boldsymbol{Y}_{i}=(Y_{i1}^{T},\ldots,Y_{iT_{i}}^{T})^{T}$ the stacked
response vector. When $J=2$ the response is binary and $\boldsymbol{Y}_{it}$
is a scalar. Let $\boldsymbol{X}_{i}=(X_{i1}^{T},\ldots,X_{iT_{i}}^{T})^{T}$
denotes the $T_{i}\times 1$ covariate vector that may be missing for the
$i$-th subject, and
$\boldsymbol{Z}_{i}=(\boldsymbol{Z}_{i1}^{T},\ldots,\boldsymbol{Z}_{iT_{i}}^{T})^{T}$
the $T_{i}\times q$ matrix of explanatory variables that are always observed.

The marginal distribution of $\boldsymbol{Y}_{it}$ is assumed to be
multinomial (with sample size $\sum_{j=1}^{J}Y_{itj}=1$), that is

$f(\boldsymbol{Y}_{it}|\boldsymbol{X}_{it},\boldsymbol{Z}_{it},\boldsymbol{\beta})=\prod_{j=1}^{J}\mu_{itj}^{y_{itj}},$
(1)

where
$\mu_{itj}=\mu_{itj}(\boldsymbol{\beta})=E(Y_{itj}|\boldsymbol{X}_{i},\boldsymbol{Z}_{i},\boldsymbol{\beta})=Pr(O_{it}=j|\boldsymbol{X}_{i},\boldsymbol{Z}_{i},\boldsymbol{\beta})$,
is the probability of response $j$ at time $t$ and $\boldsymbol{\beta}$ is a
$p\times 1$ vector of parameters. Two common choices for modeling $\mu_{itj}$
are the cumulative logit and probit models. In this work it is assumed a
cumulative logit link, that is,

$\mbox{logit}\left[Pr(O_{it}\leq
j|X_{it},Z_{it})\right]=\beta_{0j}+X_{it}\beta_{x}+\boldsymbol{Z}_{it}^{T}\boldsymbol{\beta}_{z},\
\ \ j=1,\ldots,J-1.$ (2)

Formulation in (2) implies a proportional odds model (McCullagh,, 1980). In
such model the interpretation of $\boldsymbol{\beta}$ is the same regardless
of the number of categories (i.e., it is invariant to combination of
categories). A desired feature is that the exponential of the parameters is
interpreted as an odds ratio (Agresti,, 2013).

Main interest is to make inferences related to the regression parameters
$\boldsymbol{\beta}=(\beta_{01},\ldots,\beta_{0,J-1},\beta_{x},\boldsymbol{\beta}_{z}^{T})^{T}$
associated with the $(J-1)\times 1$ marginal probability vectors

$E(Y_{it}|\boldsymbol{X}_{i},\boldsymbol{Z}_{i})=\boldsymbol{\mu}_{it}(\boldsymbol{\beta})=(\mu_{it1},\ldots,\mu_{it(J-1)})^{T}.$

$\boldsymbol{\mu}_{it}$ is grouped to form a vector
$E(\boldsymbol{Y}_{i}|\boldsymbol{X}_{i},\boldsymbol{Z}_{i})=\boldsymbol{\mu}_{i}=(\boldsymbol{\mu}_{i1}^{T},\ldots,\boldsymbol{\mu}_{iT_{i}}^{T})^{T}$
with the same dimension of $\boldsymbol{Y}_{i}$.

In order to estimate $\boldsymbol{\beta}$ generalized estimation equations are
used (Liang & Zeger, (1986); Lipsitz et al., (1994)), which takes the form

$\boldsymbol{U}(\boldsymbol{\beta})=\sum_{i=1}^{n}\boldsymbol{U}_{i}(\boldsymbol{\beta})=\sum_{i=1}^{n}\boldsymbol{D}_{i}\boldsymbol{V}_{i}^{-1}(\boldsymbol{Y}_{i}-\boldsymbol{\mu}_{i})=\boldsymbol{0},$
(3)

where
$\boldsymbol{D}_{i}=\frac{\partial\boldsymbol{\mu}_{i}}{\partial\boldsymbol{\beta}^{T}}$
and
$\boldsymbol{V}_{i}=\boldsymbol{V}_{i}(\boldsymbol{\beta},\boldsymbol{\alpha})$
is a $T_{i}(J-1)\times T_{i}(J-1)$ “working covariance” matrix usually
decomposed into the form
$\boldsymbol{V}_{i}(\boldsymbol{\beta},\boldsymbol{\alpha})=\boldsymbol{F}_{i}^{1/2}(\boldsymbol{\beta})\boldsymbol{C}_{i}(\boldsymbol{\alpha})\boldsymbol{F}_{i}^{1/2}(\boldsymbol{\beta})$,
where $\boldsymbol{F}_{i}$ is a matrix containing the marginal variances,
$\boldsymbol{F}_{it}$, given by

$\boldsymbol{F}_{it}=\mbox{diag}\left[\mu_{it1}(1-\mu_{it1}),\ldots,\mu_{it,J-1}(1-\mu_{it,J-1})\right],$

and $\boldsymbol{C}_{i}$ is equal to the marginal correlation matrix. The
$(J-1)\times(J-1)$ diagonal blocks of
$\boldsymbol{C}_{i}(\boldsymbol{\alpha})$ are
$\boldsymbol{F}_{it}^{-1/2}\boldsymbol{V}_{it}\boldsymbol{F}_{it}^{-1/2}$,
with
$\boldsymbol{V}_{it}=\mbox{diag}(\boldsymbol{\mu}_{it})-\boldsymbol{\mu}_{it}\boldsymbol{\mu}_{it}^{T}$;
and the $(J-1)\times(J-1)$ off-diagonal blocks of
$\boldsymbol{C}_{i}(\boldsymbol{\alpha})$ are
$\boldsymbol{\alpha}_{itt^{\prime}}$, which represents the correlation between
$\boldsymbol{Y}_{it}$ and $\boldsymbol{Y}_{it^{\prime}}$, $t\neq t^{\prime}$
(Lipsitz et al.,, 1994).

Under mild regularity conditions and correct specification of the marginal
mean model in $(\ref{um_mediaor})$, Liang and Zeger (1986) proved that the
estimator $\boldsymbol{\hat{\beta}}$, obtained by solving $(\ref{um_gee})$, is
consistent and $\sqrt{n}(\boldsymbol{\hat{\beta}}-\boldsymbol{\beta})$
converges in distribution to a $p$-variate normal distribution with mean
$\boldsymbol{0}$ and covariance matrix

$\boldsymbol{\mbox{V}}_{\boldsymbol{\beta}}=\lim_{n\rightarrow\infty}n\boldsymbol{\Sigma}_{0}^{-1}\boldsymbol{\Sigma}_{1}\boldsymbol{\Sigma}_{0}^{-1},$
(4)

where
$\boldsymbol{\Sigma}_{0}=\sum_{i=1}^{n}\boldsymbol{D}_{i}\boldsymbol{V}_{i}^{-1}\boldsymbol{D}_{i}^{T}$,
and
$\boldsymbol{\Sigma}_{1}=\sum_{i=1}^{n}\boldsymbol{D}_{i}\boldsymbol{V}_{i}^{-1}\mbox{Cov}(\boldsymbol{Y}_{i})\boldsymbol{V}_{i}^{-1}\boldsymbol{D}_{i}^{T}$.
In practice, the “sandwich” covariance matrix
$\boldsymbol{\mbox{V}}_{\boldsymbol{\beta}}$ in $(\ref{um_covbeta})$ is
calculated by ignoring the limit and replacing
$(\boldsymbol{\beta},\boldsymbol{\alpha})$ and
$\mbox{Cov}(\boldsymbol{Y}_{i})$ by
$(\boldsymbol{\hat{\beta}},\boldsymbol{\hat{\alpha}})$ and
$(\boldsymbol{Y}_{i}-\boldsymbol{\hat{\mu}}_{i})(\boldsymbol{Y}_{i}-\boldsymbol{\hat{\mu}}_{i})^{T}$,
respectively (Touloumis et al.,, 2013).

###  2.2 Missing Data Framework

For each occasion $t$ it can be defined $R_{it}=0$ if $O_{it}$ and $X_{it}$
are missing, $R_{it}=1$ if $O_{it}$ is missing and $X_{it}$ is observed,
$R_{it}=2$ if $O_{it}$ is observed and $X_{it}$ is missing, and $R_{it}=3$ if
$O_{it}$ and $X_{it}$ are both observed. Let
$\boldsymbol{R}_{i}=(R_{i1},\ldots,R_{iT_{i}})^{T}$, and
$\boldsymbol{\bar{R}}_{it}=(R_{i1},\ldots,R_{i,t-1})$.

By specifying conditional models of the form
$Pr(R_{it}=r_{it}|\boldsymbol{\bar{R}}_{it},\boldsymbol{O}_{i},\boldsymbol{X}_{i},\boldsymbol{Z}_{i})$
it can be obtained
$Pr(\boldsymbol{R}_{i}=\boldsymbol{r}_{i}|\boldsymbol{O}_{i},\boldsymbol{X}_{i},\boldsymbol{Z}_{i})$
through
$\prod_{i=2}^{T_{i}}Pr(R_{it}=r_{it}|\boldsymbol{\bar{R}}_{it},\boldsymbol{O}_{i},\boldsymbol{X}_{i},\boldsymbol{Z}_{i})Pr(R_{i1}=r_{i1}|\boldsymbol{O}_{i},\boldsymbol{X}_{i},\boldsymbol{Z}_{i})$.
Let
$\lambda_{itk}=Pr(R_{it}=k|\boldsymbol{\bar{R}}_{it},\boldsymbol{O}_{i},\boldsymbol{X}_{i},\boldsymbol{Z}_{i})$,
for $k=0,1,2,3$. This general formulation encompasses MCAR, MAR and MNAR
mechanisms. In particular, the MAR mechanism requires

$Pr(\boldsymbol{R}_{i}=\boldsymbol{r}_{i}|\boldsymbol{O}_{i},\boldsymbol{X}_{i},\boldsymbol{Z}_{i})=Pr(\boldsymbol{R}_{i}=\boldsymbol{r}_{i}|\boldsymbol{O}_{i}^{o},\boldsymbol{X}_{i}^{o},\boldsymbol{Z}_{it}),$
(5)

where $\boldsymbol{O}_{i}^{o}$ and $\boldsymbol{X}_{i}^{o}$ denotes the
observed components of $\boldsymbol{O}_{i}$ and $\boldsymbol{X}_{i}$,
respectively. Because $\boldsymbol{R}_{i}$ is modeled through a product of
conditional models, it is natural to make the following further assumption
(Chen & Zhou,, 2011)

$Pr(R_{it}=r_{it}|\boldsymbol{\bar{R}}_{it},\boldsymbol{O}_{i},\boldsymbol{X}_{i},\boldsymbol{Z}_{i})=Pr(R_{it}=r_{it}|\boldsymbol{\bar{R}}_{it},\boldsymbol{\bar{O}}_{it}^{o},\boldsymbol{\bar{X}}_{it}^{o},\boldsymbol{\bar{Z}}_{it}),$
(6)

for each time $t$, where $\boldsymbol{\bar{O}}_{it}^{o}$ and
$\boldsymbol{\bar{X}}_{it}^{o}$ are the histories of observed responses and
covariates up to time $t-1$.

Let
$\pi_{it}=Pr(R_{it}=3|\boldsymbol{O}_{i},\boldsymbol{X}_{i},\boldsymbol{Z}_{i})$
be the marginal probability of observing both $\boldsymbol{O}_{i}$ and
$\boldsymbol{X}_{i}$ at time $t$, given the entire vectors of responses and
covariates. Then, $\pi_{it}$ is expressed by

$\pi_{it}=\sum_{r_{i1},\ldots,r_{i,t-1}}Pr(R_{it}=3,R_{i,t-1}=r_{i,t-1},\ldots,R_{i1}=r_{i1}|\boldsymbol{O}_{i},\boldsymbol{X}_{i},\boldsymbol{Z}_{i}).$

This marginal probability can be expressed in terms of the conditional
probabilities $\lambda_{itk}$’s. Throughout the paper it is required the so-
called _positivity assumption_ , that is, $\pi_{it}$ must be bounded away from
zero. This condition is needed in order to guarantee the existence of
$\sqrt{n}$-consistent estimators of $\boldsymbol{\beta}$ (Robins et al.,,
1995).

##  3 Available Approaches for Missing Data

Multiple imputation and weighted generalized estimation equations are two
commonly used methods available for missing data under MAR mechanism. These
methods are presented in Sections 3.1 and 3.2, respectively. They serve as the
basis for the construction of the doubly robust estimator, presented in
Section 4.

###  3.1 Multiple Imputation Generalized Estimating Equations

A imputation model commonly used to handle intermittently missing response and
covariate, is imputation using chained equations (van Buuren et al., (1999),
van Buuren, (2007)), which is more commonly referred to as full conditional
specification (FCS). This approach specifies conditional distributions for
each incomplete variable, conditional on all others variables in the
imputation model. Starting from an initial imputation, FCS draws imputations
by iterating over the conditional densities.

Denote by $\boldsymbol{\tilde{\beta}}_{m}$ and $\boldsymbol{\tilde{U}}_{m}$,
respectively, the estimate of $\boldsymbol{\beta}$ and its covariance matrix
from the GEE analysis of the $m$-th completed data set, $(m=1,\ldots,M)$.
Following Rubin, (1987), the combined point estimate for the parameter of
interest $\boldsymbol{\beta}$ from the MI is simply the average of the M
complete-data point estimates

$\boldsymbol{\hat{\beta}}_{MI}=\frac{1}{M}\sum_{m=1}^{M}\boldsymbol{\tilde{\beta}}_{m},$

and an estimate of the covariance matrix of $\boldsymbol{\hat{\beta}}_{MI}$ is
given by

$\boldsymbol{\widehat{U}}_{MI}=\boldsymbol{\tilde{W}}+\left(\frac{M+1}{M}\right)\boldsymbol{\tilde{B}},$

where

$\boldsymbol{\tilde{W}}=\frac{1}{M}\sum_{m=1}^{M}\boldsymbol{\tilde{U}}_{m}\ \
\ \mbox{and}\ \ \
\boldsymbol{\tilde{B}}=\frac{1}{M-1}\sum_{m=1}^{M}(\boldsymbol{\tilde{\beta}}_{m}-\boldsymbol{\hat{\beta}_{MI}})(\boldsymbol{\tilde{\beta}}_{m}-\boldsymbol{\hat{\beta}_{MI}})^{\prime}.$

###  3.2 Weighted Generalized Estimating Equations

Robins et al., (1995) proposed a class of weighted estimating equations to
allow for MAR mechanism. In binary longitudinal data, Chen & Zhou, (2011)
extended the method to accommodate arbitrary patterns of missing response and
missing covariate. Their method was adapted here for longitudinal ordinal
responses.

Define a weight matrix
$\boldsymbol{\Delta}_{i}=\left[\delta_{itt^{\prime}}\right]_{T_{i}(J_{i}-1)\times
T_{i}(J_{i}-1)},\ t=1,\ldots,T_{i},t^{\prime}=1,\ldots,T_{i},$ where
$\delta_{itt^{\prime}}=\left\\{I(R_{it}=1,R_{it^{\prime}}=3)+I(R_{it}=3,R_{it^{\prime}}=3)\right\\}/\pi_{itt^{\prime}}$
for $t\neq t^{\prime}$, $\delta_{itt}=I(R_{it}=3)/\pi_{it}$, and
$\pi_{itt^{\prime}}=Pr(R_{it}=1,R_{it^{\prime}}=3|\boldsymbol{O}_{i},\boldsymbol{X}_{i},\boldsymbol{Z}_{i})+Pr(R_{it}=3,R_{it^{\prime}}=3|\boldsymbol{O}_{i},\boldsymbol{X}_{i},\boldsymbol{Z}_{i})$.
Let
$\boldsymbol{M}_{i}=\boldsymbol{F}_{i}^{-1/2}(\boldsymbol{C}_{i}^{-1}\boldsymbol{\cdot}\boldsymbol{\Delta}_{i})\boldsymbol{F}_{i}^{-1/2}$
where $\boldsymbol{A\cdot B}=\left[a_{it}\cdot b_{it}\right]$ denotes the
Hadamard product of matrix $\boldsymbol{A}=\left[a_{it}\right]$ and
$\boldsymbol{B}=\left[b_{it}\right]$.

The weighted generalized estimating equations (WGEE) for $\boldsymbol{\beta}$
are given by

$\boldsymbol{U}(\boldsymbol{\beta},\boldsymbol{\psi})=\sum_{i=1}^{n}\boldsymbol{U}_{i}(\boldsymbol{\beta},\boldsymbol{\psi})=\boldsymbol{0},$
(7)

where
$\boldsymbol{U}_{i}(\boldsymbol{\beta},\boldsymbol{\psi})=\boldsymbol{D}_{i}\boldsymbol{M}_{i}(\boldsymbol{Y}_{i}-\boldsymbol{\mu}_{i})$.
A consistent estimate for $\boldsymbol{\beta}$ can be obtained by solving (7),
under the correct specification of the missing data model.

To model $\lambda_{itk}$ it is adopted a politomic logistic regression, with
$\lambda_{it0}$ as the reference category, that is

$\log\left(\frac{\lambda_{itk}}{\lambda_{it0}}\right)=\boldsymbol{u}_{itk}^{T}\psi_{k},\
\ \ k=1,2,3,$ (8)

where the covariates $\boldsymbol{u}_{itk}$ are some function of
$\left\\{\boldsymbol{\bar{R}}_{it},\boldsymbol{\bar{O}}_{it}^{o},\boldsymbol{\bar{X}}_{it}^{o},\boldsymbol{\bar{Z}}_{it}\right\\}$.

##  4 Doubly Robust GEE for Longitudinal Ordinal Data

Some authors (e.g., Scharfstein et al., (1999), Tsiatis, (2006)) noted that
adding a term of expectation zero, say $\phi(\cdot)$, to the inverse
probability weighted estimators would still result in consistent estimates
under a MAR mechanism. The solutions of these augmented estimating equations
give rise to the so-called _doubly robust_ estimators.

Chen & Zhou, (2011) showed that the optimal $\phi_{opt}$ for missing response
and covariate is given by
$\phi_{opt}=E_{(\boldsymbol{Y}_{i}^{m},\boldsymbol{X}_{i}^{m}|\boldsymbol{Y}_{i}^{o},\boldsymbol{X}_{i}^{o},\boldsymbol{Z}_{i},\boldsymbol{R}_{i})}\left\\{\boldsymbol{D}_{i}\boldsymbol{N}_{i}(\boldsymbol{Y}_{i}-\boldsymbol{\mu}_{i})\right\\}$,
with
$\boldsymbol{N}_{i}=\boldsymbol{F}_{i}^{-1/2}\left\\{\boldsymbol{C}_{i}^{-1}\boldsymbol{\cdot}(\boldsymbol{11}^{T}-\boldsymbol{\Delta}_{i})\right\\}\boldsymbol{F}_{i}^{-1/2}$,
where $\boldsymbol{1}$ is a vector of 1’s of length $T_{i}(J-1)$, and
$\boldsymbol{Y}_{i}^{m}$ and $\boldsymbol{X}_{i}^{m}$ denote the missing
components of $\boldsymbol{Y}_{i}$ and $\boldsymbol{X}_{i}$, respectively.

An improved estimate for $\boldsymbol{\beta}$ can then be obtained by solving
the estimating equations

$\boldsymbol{S}_{1}(\boldsymbol{\theta})=\sum_{i=1}^{n}\boldsymbol{S}_{1i}(\boldsymbol{\theta})=\sum_{i=1}^{n}\left[\boldsymbol{D}_{i}\boldsymbol{M}_{i}(\boldsymbol{Y}_{i}-\boldsymbol{\mu}_{i})+E_{(\boldsymbol{Y}_{i}^{m},\boldsymbol{X}_{i}^{m}|\boldsymbol{Y}_{i}^{o},\boldsymbol{X}_{i}^{o},\boldsymbol{Z}_{i},\boldsymbol{R}_{i})}\left\\{\boldsymbol{D}_{i}\boldsymbol{N}_{i}(\boldsymbol{Y}_{i}-\boldsymbol{\mu}_{i})\right\\}\right]=\boldsymbol{0}.$
(9)

The estimator for $\boldsymbol{\beta}$ in (9) is doubly-robust in the sense
that it is consistent if _at least one_ of the missing data model or the
covariate model is correctly specified.

Applications of doubly robust estimators in longitudinal settings include Bang
& Robins, (2005), Seaman & Copas, (2009), Chen & Zhou, (2011) and Birhanu et
al., (2011). Those developments focus mainly on binary response and have not
been, to our knowledge, investigated with ordinal longitudinal data. In this
work it is considered a longitudinal response measured on a ordinal scale.

The referred expectation in the second part of (9) is over the conditional
distribution of
$(\boldsymbol{Y}_{i}^{m},\boldsymbol{X}_{i}^{m}|\boldsymbol{Y}_{i}^{o},\boldsymbol{X}_{i}^{o},\boldsymbol{Z}_{i},\boldsymbol{R}_{i})$,
which can be written as

$\displaystyle
P(\boldsymbol{Y}_{i}^{m}=\boldsymbol{y}_{i}^{m},\boldsymbol{X}_{i}^{m}=\boldsymbol{x}_{i}^{m}|\boldsymbol{Y}_{i}^{o},\boldsymbol{X}_{i}^{o},\boldsymbol{Z}_{i},\boldsymbol{R}_{i};\boldsymbol{\beta}^{*},\boldsymbol{\gamma})$
$\displaystyle=$ $\displaystyle
P(\boldsymbol{Y}_{i}^{m}=\boldsymbol{y}_{i}^{m},\boldsymbol{X}_{i}^{m}=\boldsymbol{x}_{i}^{m}|\boldsymbol{Y}_{i}^{o},\boldsymbol{X}_{i}^{o},\boldsymbol{Z}_{i};\boldsymbol{\beta}^{*},\boldsymbol{\gamma})$
$\displaystyle=$ $\displaystyle
P(\boldsymbol{Y}_{i}^{m}=\boldsymbol{y}_{i}^{m}|\boldsymbol{Y}_{i}^{o},\boldsymbol{X}_{i}=\boldsymbol{x}_{i},\boldsymbol{Z}_{i};\boldsymbol{\beta}^{*})$
$\displaystyle\times
P(\boldsymbol{X}_{i}^{m}=\boldsymbol{x}_{i}^{m}|\boldsymbol{Y}_{i}^{o},\boldsymbol{X}_{i}^{o},\boldsymbol{Z}_{i};\boldsymbol{\gamma}).$

The multivariate distribution
$P(\boldsymbol{Y}_{i}^{m}=\boldsymbol{y}_{i}^{m}|\boldsymbol{Y}_{i}^{o},\boldsymbol{X}_{i}=\boldsymbol{x}_{i},\boldsymbol{Z}_{i};\boldsymbol{\beta}^{*})$
is expressed through a product of univariate ordinal models.

###  4.1 Estimation for the Nuisance Parameters

The method of maximum likelihood is employed to estimate $\boldsymbol{\psi}$.
The log likelihood of the politomic logistic model for $\boldsymbol{\psi}$ has
the form

$l(\boldsymbol{\boldsymbol{\psi}})=\sum_{i=1}^{n}l_{i}(\boldsymbol{\psi})=\sum_{i=1}^{n}\sum_{t=1}^{T_{i}}\sum_{k=0}^{3}I(R_{it}=k)\log(\lambda_{itk}),$

with corresponding score function given by

$\displaystyle\boldsymbol{S}_{2}(\boldsymbol{\psi})=\sum_{i=1}^{n}\boldsymbol{S}_{2i}(\boldsymbol{\psi})=\sum_{i=1}^{n}\sum_{t=1}^{T_{i}}\sum_{k=0}^{3}\frac{I(R_{it}=k)}{\lambda_{itk}}\frac{\partial\lambda_{itk}}{\partial\boldsymbol{\psi}^{T}}.$

The maximum likelihood estimator $\boldsymbol{\hat{\psi}}$, is obtained by
solving $\boldsymbol{S}_{2}(\boldsymbol{\psi})=\boldsymbol{0}$.

For the missing covariate model, the observed likelihood function for
$\boldsymbol{\gamma}$ is

$L_{3}(\boldsymbol{\gamma})=\prod_{i=1}^{n}\int
Pr(\boldsymbol{X}_{i}|\boldsymbol{Z}_{i},\boldsymbol{Y}_{i})d\boldsymbol{X}_{i}^{m},$

with score function
$\boldsymbol{S}_{3}(\boldsymbol{\gamma})=\sum_{i=1}^{n}\boldsymbol{S}_{3i}(\boldsymbol{\gamma})$,
where $\boldsymbol{S}_{3i}(\boldsymbol{\gamma})=\partial\log\int
Pr(\boldsymbol{X}_{i}|\boldsymbol{Z}_{i},\boldsymbol{Y}_{i}))d\boldsymbol{X}_{i}^{m}/\partial\boldsymbol{\gamma}^{T}$.
Similarly, a consistent estimator of $\boldsymbol{\gamma}$ can be obtained by
solving $\boldsymbol{S}_{3}(\boldsymbol{\gamma})=\boldsymbol{0}$.

###  4.2 Estimation and Inference for the Doubly Robust Method

Let’s denote the vector of all parameters as
$\boldsymbol{\theta}=(\boldsymbol{\beta}^{T},\boldsymbol{\psi}^{T},\boldsymbol{\gamma}^{T})^{T}$.
Our primary interest lies is in estimating $\boldsymbol{\beta}$. Such task can
be accomplished by plugging in the estimates $\boldsymbol{\hat{\psi}}$ and
$\boldsymbol{\hat{\gamma}}$ in (9) and solving the estimating equations for
$\boldsymbol{\beta}$, that is,

$\boldsymbol{S}_{1}(\boldsymbol{\beta},\boldsymbol{\hat{\psi}},\boldsymbol{\hat{\gamma}})=\sum_{i=1}^{n}\boldsymbol{S}_{1i}(\boldsymbol{\beta},\boldsymbol{\hat{\psi}},\boldsymbol{\hat{\gamma}})=0.$
(10)

Second term of $\boldsymbol{S}_{1i}$ can be written, for $\boldsymbol{X}$
discrete, as

$E_{(\boldsymbol{Y}_{i}^{m},\boldsymbol{X}_{i}^{m}|\boldsymbol{Y}_{i}^{o},\boldsymbol{X}_{i}^{o},\boldsymbol{Z}_{i},\boldsymbol{R}_{i})}\left\\{\boldsymbol{D}_{i}\boldsymbol{N}_{i}(\boldsymbol{Y}_{i}-\boldsymbol{\mu}_{i})\right\\}=\sum_{(\boldsymbol{y}_{i}^{m},\boldsymbol{x}_{i}^{m})}w_{ixy}\left\\{\boldsymbol{D}_{i}\boldsymbol{N}_{i}(\boldsymbol{Y}_{i}-\boldsymbol{\mu}_{i})\right\\},$

where the weight $w_{ixy}$ is given by

$\displaystyle w_{ixy}$ $\displaystyle=$ $\displaystyle
P(\boldsymbol{Y}_{i}^{m}=\boldsymbol{y}_{i}^{m}|\boldsymbol{Y}_{i}^{o},\boldsymbol{X}_{i}=\boldsymbol{x}_{i},\boldsymbol{Z}_{i};\boldsymbol{\beta}^{*})\times
P(\boldsymbol{X}_{i}^{m}=\boldsymbol{x}_{i}^{m}|\boldsymbol{Y}_{i}^{o},\boldsymbol{X}_{i}^{o},\boldsymbol{Z}_{i};\boldsymbol{\hat{\gamma}}).$

In the case of $\boldsymbol{X}$ continuous, the second term in
$\boldsymbol{S}_{1i}$ takes the form

$E_{(\boldsymbol{Y}_{i}^{m},\boldsymbol{X}_{i}^{m}|\boldsymbol{Y}_{i}^{o},\boldsymbol{X}_{i}^{o},\boldsymbol{Z}_{i},\boldsymbol{R}_{i})}\left\\{\boldsymbol{D}_{i}\boldsymbol{N}_{i}(\boldsymbol{Y}_{i}-\boldsymbol{\mu}_{i})\right\\}=\int_{(\boldsymbol{Y}_{i}^{m},\boldsymbol{X}_{i}^{m})}w_{ixy}\left\\{\boldsymbol{D}_{i}\boldsymbol{N}_{i}(\boldsymbol{Y}_{i}-\boldsymbol{\mu}_{i})\right\\}d\boldsymbol{Y}_{i}^{m}\boldsymbol{X}_{i}^{m},$

with conditional probability

$\displaystyle w_{ixy}$ $\displaystyle=$ $\displaystyle
P(\boldsymbol{Y}_{i}^{m}=\boldsymbol{y}_{i}^{m}|\boldsymbol{Y}_{i}^{o},\boldsymbol{X}_{i}=\boldsymbol{x}_{i},\boldsymbol{Z}_{i};\boldsymbol{\beta}^{*})\times
P(\boldsymbol{X}_{i}^{m}=\boldsymbol{x}_{i}^{m}|\boldsymbol{Y}_{i}^{o},\boldsymbol{X}_{i}^{o},\boldsymbol{Z}_{i};\boldsymbol{\hat{\gamma}}).$

The expectation in (9) can be cumbersome, depending on the missing data
pattern. In such case, instead of using numerical integration, a Monte Carlo
method can be applied to approximate the corresponding integral.

In this work it is assumed independence working correlation and it is adopted
a sandwich standard error as given in Appendix.

##  5 Simulation Study

A small simulation study taking into account different sample sizes was
conducted in order to quantify the bias and precision under misspecification
of the predictive models. It is considered a study with $T_{i}=T=3$ repeated
ordinal measures (with three categories) and two covariates (quantitative and
qualitative). The true marginal model is

$logit\ Pr(O_{it}\leq
j|X_{it},Z_{it})=\beta_{0j}+\beta_{1}X_{it}+\beta_{2}Z_{it},\ \ \ j=1,2.$ (11)

where $Z_{it}$ is normal with unity variance and mean $(0,0.5,1)$ for
$t=1,2,3$.

The binary covariate $X_{it}$ may be missing at some time points and is
generated according to

$logit\
Pr(X_{it}=1|\bar{X}_{it},Z_{it})=\gamma_{0}+\gamma_{1}X_{i,t-1}+\gamma_{2}Z_{it}.$
(12)

It is assumed $\beta_{01}=-0.4$, $\beta_{02}=1.2$, $\beta_{1}=-0.5$,
$\beta_{2}=0.5$, $\gamma_{0}=log(1)$, $\gamma_{1}=2$ and $\gamma_{2}=2$. The
correlated ordinal responses were generated according to the algorithm
proposed by Touloumis (Anestis Touloumis,, 2013) with constant correlation
between the latent vectors set equal to $\rho=0.9$.

As independent estimating equations were fitted, $R_{it}$ can be defined as
the indicator of observing both $O_{it}$ and $X_{it}$, and it was taken

$log\left(\frac{Pr(R_{it}=1)}{Pr(R_{it}=0)}\right)=\psi_{0t}+\psi_{1}I(R_{i,t-1}=1)+\psi_{2}O_{i,t-1}^{*}+\psi_{3}X_{i,t-1}^{*}+\psi_{4}Z_{it},\
\ t=2,3,$ (13)

where $O_{i,t-1}^{*}=O_{i,t-1}$, if $O_{i,t-1}$ is observed and $0$ otherwise,
and $X_{i,t-1}^{*}=X_{i,t-1}$if $X_{i,t-1}$ is observed and $0$ otherwise. The
true values are taken as $\psi_{02}=6.6$, $\psi_{03}=6$, $\psi_{1}=2$,
$\psi_{2}=-2$, $\psi_{3}=-2$ and $\psi_{4}=2$. It was observed about $24\%$ of
missing observations under this setup.

For comparison purposes, it was considered ordinary GEE for the complete and
available data, respectively, weighted GEE (WGEE), multiple imputation (MIGEE)
by chained equations (van Buuren & Groothuis-Oudshoorn,, 2011) with $M=10$,
and the proposed doubly robust version (DRGEE). In order to investigate
robustness of these methods, the predicted models were also misspecified by
omitting the covariate $X_{t-1}$ from the covariate model (12) or the missing
data model (13).

Results are summarized in Table 1. In each of the $S=1000$ Monte Carlo
replications it was obtained the relative bias percentage for each parameter,
defined as $100\times(\hat{\beta}-\beta)/\beta$, its standard deviation
obtained through the sandwich estimator, and the coverage probability as a
nominal 95%.

Table 1: Relative bias percentage, standard deviation and empirical coverage for $1000$ simulations of incomplete covariate and response data. | Empirical Bias |  | Standard Deviation |  | Empirical Coverage
---|---|---|---|---|---
| $\beta_{01}$ | $\beta_{02}$ | $\beta_{1}$ | $\beta_{2}$ |  | $\beta_{01}$ | $\beta_{02}$ | $\beta_{1}$ | $\beta_{2}$ |  | $\beta_{01}$ | $\beta_{02}$ | $\beta_{1}$ | $\beta_{2}$
| $n=50$
Complete | -4.60 | 7.89 | 15.07 | 7.38 |  | 0.379 | 0.410 | 0.458 | 0.179 |  | 0.94 | 0.94 | 0.97 | 0.95
Available | -19.81 | 15.37 | -8.36 | -11.66 |  | 0.386 | 0.420 | 0.473 | 0.190 |  | 0.94 | 0.92 | 0.97 | 0.94
WGEE($r^{+}$) | -12.33 | 12.02 | 19.96 | 1.74 |  | 0.406 | 0.440 | 0.513 | 0.204 |  | 0.94 | 0.94 | 0.97 | 0.95
WGEE($r^{-}$) | 1.50 | 7.43 | 1.31 | 0.51 |  | 0.413 | 0.450 | 0.508 | 0.201 |  | 0.95 | 0.93 | 0.96 | 0.95
MIGEE($x^{+}$) | 22.17 | -3.44 | -38.26 | -4.73 |  | 0.386 | 0.417 | 0.480 | 0.186 |  | 0.96 | 0.96 | 0.97 | 0.96
MIGEE($x^{-}$) | 15.81 | -0.13 | -26.18 | -2.54 |  | 0.387 | 0.417 | 0.481 | 0.187 |  | 0.95 | 0.96 | 0.97 | 0.96
DRGEE($x^{+},r^{+}$) | -8.86 | 9.67 | 22.20 | 8.99 |  | 0.410 | 0.442 | 0.513 | 0.199 |  | 0.95 | 0.94 | 0.97 | 0.95
DRGEE($x^{-},r^{+}$) | -7.77 | 9.44 | 19.33 | 7.91 |  | 0.471 | 0.507 | 0.579 | 0.204 |  | 0.95 | 0.94 | 0.97 | 0.95
DRGEE($x^{+},r^{-}$) | -7.95 | 10.05 | 20.58 | 8.26 |  | 0.433 | 0.473 | 0.540 | 0.197 |  | 0.95 | 0.93 | 0.96 | 0.95
DRGEE($x^{-},r^{-}$) | 2.86 | 6.42 | -0.28 | 2.99 |  | 0.414 | 0.450 | 0.529 | 0.200 |  | 0.95 | 0.93 | 0.97 | 0.95
| $n=150$
Complete | -1.36 | 1.53 | 5.13 | 1.99 |  | 0.220 | 0.236 | 0.264 | 0.103 |  | 0.94 | 0.94 | 0.96 | 0.95
Available | -17.50 | 8.94 | -18.41 | -17.15 |  | 0.224 | 0.242 | 0.273 | 0.109 |  | 0.91 | 0.92 | 0.94 | 0.91
WGEE($r^{+}$) | -5.84 | 3.41 | 13.47 | 1.36 |  | 0.247 | 0.267 | 0.317 | 0.126 |  | 0.93 | 0.93 | 0.95 | 0.94
WGEE($r^{-}$) | 10.20 | -2.04 | -12.91 | -2.59 |  | 0.253 | 0.279 | 0.307 | 0.120 |  | 0.94 | 0.93 | 0.93 | 0.94
MIGEE($x^{+}$) | 8.22 | -2.57 | -14.36 | -2.40 |  | 0.228 | 0.244 | 0.283 | 0.108 |  | 0.94 | 0.94 | 0.96 | 0.95
MIGEE($x^{-}$) | 9.08 | -1.97 | -16.33 | -3.60 |  | 0.226 | 0.242 | 0.280 | 0.108 |  | 0.93 | 0.93 | 0.95 | 0.94
DRGEE($x^{+},r^{+}$) | -4.21 | 2.64 | 11.90 | 4.12 |  | 0.249 | 0.269 | 0.317 | 0.117 |  | 0.94 | 0.94 | 0.96 | 0.95
DRGEE($x^{-},r^{+}$) | -5.02 | 3.07 | 12.26 | 3.63 |  | 0.320 | 0.345 | 0.396 | 0.119 |  | 0.94 | 0.94 | 0.95 | 0.95
DRGEE($x^{+},r^{-}$) | -1.86 | 1.99 | 7.26 | 2.95 |  | 0.249 | 0.275 | 0.313 | 0.113 |  | 0.94 | 0.94 | 0.95 | 0.94
DRGEE($x^{-},r^{-}$) | 9.88 | -1.86 | -15.60 | -2.86 |  | 0.246 | 0.269 | 0.317 | 0.116 |  | 0.94 | 0.93 | 0.94 | 0.94
| $n=300$
Complete | 1.42 | 0.19 | -0.04 | 0.85 |  | 0.156 | 0.167 | 0.187 | 0.072 |  | 0.96 | 0.94 | 0.94 | 0.95
Available | -14.82 | 7.67 | -23.24 | -18.14 |  | 0.159 | 0.171 | 0.194 | 0.077 |  | 0.93 | 0.91 | 0.92 | 0.87
WGEE($r^{+}$) | -0.24 | 0.75 | 5.61 | 1.23 |  | 0.182 | 0.198 | 0.235 | 0.095 |  | 0.95 | 0.94 | 0.94 | 0.93
WGEE($r^{-}$) | 15.64 | -4.50 | -22.17 | -3.89 |  | 0.184 | 0.205 | 0.222 | 0.088 |  | 0.95 | 0.94 | 0.92 | 0.94
MIGEE($x^{+}$) | 6.23 | -1.86 | -9.95 | -1.44 |  | 0.162 | 0.174 | 0.201 | 0.076 |  | 0.95 | 0.94 | 0.94 | 0.95
MIGEE($x^{-}$) | 10.27 | -2.70 | -18.22 | -3.89 |  | 0.160 | 0.172 | 0.200 | 0.076 |  | 0.94 | 0.92 | 0.92 | 0.94
DRGEE($x^{+},r^{+}$) | -0.06 | 0.83 | 2.82 | 1.28 |  | 0.185 | 0.202 | 0.242 | 0.089 |  | 0.96 | 0.95 | 0.95 | 0.94
DRGEE($x^{-},r^{+}$) | -0.71 | 1.08 | 3.84 | 1.43 |  | 0.194 | 0.213 | 0.251 | 0.088 |  | 0.96 | 0.95 | 0.95 | 0.94
DRGEE($x^{+},r^{-}$) | 2.53 | 0.03 | -2.11 | 0.12 |  | 0.178 | 0.198 | 0.224 | 0.082 |  | 0.96 | 0.96 | 0.94 | 0.94
DRGEE($x^{-},r^{-}$) | 13.86 | -3.69 | -24.29 | -5.51 |  | 0.176 | 0.194 | 0.228 | 0.084 |  | 0.95 | 0.94 | 0.93 | 0.94
| $n=600$
Complete | -1.01 | 0.57 | 1.93 | 0.82 |  | 0.110 | 0.118 | 0.132 | 0.051 |  | 0.95 | 0.93 | 0.94 | 0.94
Available | -16.78 | 7.85 | -22.20 | -18.43 |  | 0.112 | 0.121 | 0.137 | 0.054 |  | 0.91 | 0.88 | 0.92 | 0.78
WGEE($r^{+}$) | -1.00 | 0.28 | 5.43 | 1.25 |  | 0.132 | 0.144 | 0.170 | 0.070 |  | 0.95 | 0.94 | 0.95 | 0.95
WGEE($r^{-}$) | 14.87 | -4.99 | -21.53 | -3.03 |  | 0.133 | 0.149 | 0.158 | 0.064 |  | 0.94 | 0.94 | 0.92 | 0.93
MIGEE($x^{+}$) | 1.76 | -0.66 | -3.63 | -0.39 |  | 0.115 | 0.123 | 0.143 | 0.054 |  | 0.95 | 0.94 | 0.95 | 0.94
MIGEE($x^{-}$) | 6.96 | -2.01 | -14.18 | -3.38 |  | 0.113 | 0.122 | 0.141 | 0.054 |  | 0.94 | 0.93 | 0.93 | 0.93
DRGEE($x^{+},r^{+}$) | -1.49 | 0.59 | 3.19 | 0.98 |  | 0.132 | 0.144 | 0.174 | 0.062 |  | 0.96 | 0.94 | 0.94 | 0.94
DRGEE($x^{-},r^{+}$) | -1.39 | 0.56 | 2.96 | 0.91 |  | 0.136 | 0.149 | 0.178 | 0.062 |  | 0.96 | 0.94 | 0.96 | 0.95
DRGEE($x^{+},r^{-}$) | 0.42 | 0.08 | -0.12 | 0.47 |  | 0.126 | 0.140 | 0.159 | 0.058 |  | 0.96 | 0.95 | 0.95 | 0.93
DRGEE($x^{-},r^{-}$) | 12.24 | -3.82 | -23.18 | -5.35 |  | 0.124 | 0.138 | 0.161 | 0.060 |  | 0.93 | 0.93 | 0.92 | 0.93
“+" indicates correctly specified model and “-" indicates misspecified model omitting the $X_{t}$ predictor |

Specially, MAR missingness impact over the response and the covariate is
observed for all the regression parameters, the largest relative bias occur in
binary covariate $X$. This comes in addition to the natural increase of
parameter uncertainty. Bias in the intercept coefficients imply incorrect
predicted probabilities for the levels of the response, whereas bias for
parameter estimates associated with the regression covariates may erroneously
attenuate or highlight an effect, thus leading to misinterpretations related
to the importance of a given predictor on the longitudinal dynamics of the
ordinal response.

It can be observed that for small sample size even the GEE with for complete
data presents a certain degree of bias. Increasing the sample size allows to
clarify the performance distinctions among the compared methods. WGEE and
MIGEE methods are valid when the model for the weight or the imputation model,
respectively, are correctly specified. In this case it is noted that both
methods give good results for large sample sizes, the main distinction between
them being due to the greater variability of the estimates for the weighted
estimator.

Doubly robust method requires the simultaneous specification of two predictive
model. When at least one of them is correctly specified the resulting
estimator is still consistent. Estimates are, on average, closer to those
obtained with fully observed data compared to WGEE or MIGEE. This behavior is
systematic and it can be observed to all parameters. By increasing the sample
size the estimates from DRGEE present empirical bias in general smaller than
their single robust competitors. This is specially true for the parameter
associated with the incomplete binary variable $X$. Regarding the uncertainty
of parameter estimates, it is noted that the variability in DRGEE is greater
than multiple imputation, but of the same order as of the weighted method.
Further, the efficiency of the doubly robust estimates appears relatively more
sensitive to misspecification of the weight model than the covariate model.
Empirical coverage rates were acceptable for correctly specified WGEE and
MIGEE as well as for DRGEE when at least one of predictive models are
correctly specified.

Figure 1 shows boxplots of the percentage relative bias for the methods
expected to be valid. The boxs represent, respectively, GEE for complete and
available data, WGEE($r^{+}$), MIGEE($x^{+}$), DRGEE($x^{+},r^{+}$),
DRGEE($x^{-},r^{+}$) and DRGEE($x^{+},r^{-}$). As the degree of bias is
different in the parameter estimates, for easy of visualization, they are
represented in different scales. It can be seen that the estimates with larger
variability were those associated with the first intercept and the covariate
with incomplete values. Proposed method presents median relative bias close to
zero and very similar to those observed with complete data for all sample
sizes. Variability of the DR estimators is slightly larger than MI but of the
same order as the weighted estimator. As expected, for all methods it can be
noticed a decrease in the relative bias with increasing sample size,
reflecting their theoretical asymptotic consistency.

Figure 1: Boxplot of the relative bias for parameter estimates under correctly
specified models

Figure 2 allows the comparison of methods incorrectly specified. The boxs
represent, respectively, GEE for complete and available data, WGEE($r^{-}$),
MIGEE($x^{-}$), and DRGEE($x^{-},r^{-}$).

Figure 2: Boxplot of the relative bias for the estimates with incorrectly
specified models

When a key covariate is omitted from the weight model and/or the imputation
one, it is expected all methods to be biased. This is specially true for the
first intercept and for the incomplete binary covariate. Bias for
MIGEE($x^{-}$) seemed to get smaller than WGEE($r^{-}$) as sample sizes
increases. Bias for DRGEE($x^{-},r^{-}$) was comparable to WGEE($r^{-}$) and
slightly higher than multiple imputation for large sample sizes. In terms of
variability of the estimates, the same pattern is observed as with the
correctly specified models. That is, the multiple imputation is more
efficient, followed by the doubly robust estimator and the weighed estimator.

##  6 Data Analysis: Analgesia in Childbirth

This study was conducted in Minas Gerais state, Brazil, in order to compare
two techniques of analgesia for labor pain. There were $49$ patients who were
monitored during their entire labor period until childbirth. Pain intensity
was subjectively assessed by the patient, and measurements of blood pressure,
maternal heart rate, consumption of oxytocin, sedation level, signs of
respiratory depression, apnea, and other variables were recorded. One of the
techniques used was epidural analgesia (the gold standard), which is a local
anesthetic. The other one, whose efficiency was to be compared to the gold
standard, involved continuous intravenous infusion of remifentanil.

The response of interest is the intensity of pain as measured by a Visual
Analog Scale (VAS) ($1$: tolerable and mild pain; $2$: moderate pain that
causes discomfort; $3$: intense and unbearable pain). Three measurements
($0,60,90$ minutes) were selected for data analysis. Predictor variables
considered were treatment GROUP ($0$: peridural; $1$: remifentanil), AGE (in
years), DU (uterine dilatation), and OXYT (consumption of oxytocin). The OXYT
is a time-varying ordinal covariate, coded as $1$, if no consumption, $2$, if
consumption equals to $10$ or $30$, and $3$, if consumption equals or above
$45$. Thes eother covariates were chosen after a previous exploratory
analysis.

The response and oxytocin consumption was missing for $9$ patients at the time
$60$ and for $18$ patients at time $90$. Missing is due to childbirth happened
before $60$ or $90$ minutes. Therefore a MAR mechanism seems to be a
reasonable assumption for this data set. The others covariates in the analysis
were fully observed.

For the ordinal response it was used the following proportional odds model

$logit\ Pr(PAIN_{itj}\leq
j|\boldsymbol{u_{it}})=\beta_{0j}+\boldsymbol{u}_{it}^{T}\boldsymbol{\beta},\
\ \ j=1,2,\ \ t=1,2,3,$ (14)

where $\boldsymbol{u}_{it}$ is the covariate vector at time $t$, and it is
formed by TIME, GROUP, AGE, DU and OXIT. In response model, TIME predictor was
expressed in hours rather than minutes.

When using WGEE or DRGGE it is necessary to correctly model $\pi_{i}$ in order
to obtain consistent estimates of $\boldsymbol{\beta}$. For the missing data
process $R_{it}$ was defined as the indicator of observing both $PAIN_{it}$
and $OXIT_{it}$, and take the following form

$log\left(\frac{Pr(R_{it}=1)}{Pr(R_{it}=0)}\right)=\psi_{0t}+\boldsymbol{w}_{it}^{T}\boldsymbol{\psi},\
\ t=2,3,$ (15)

where $\boldsymbol{w}_{it}$ includes GROUP, AGE, DU, histories of OXYT and
PAIN, and the previous indicator of missing data.

The distribution of the missing covariate OXYT also needs to be specified in a
predictive model. With this aim, it was assumed a proportional odds model of
the form

$logito\ Pr(OXYT_{itj}\leq
j|\boldsymbol{v_{it}})=\gamma_{0j}+\boldsymbol{v}_{it}^{T}\boldsymbol{\gamma},\
\ \ j=1,2,\ \ t=2,3$ (16)

where $\boldsymbol{v}_{it}$ includes main effects for GROUP, AGE, and DU. It
can be observed that the estimate of $\boldsymbol{\gamma}$ is not of interest,
however it is necessary to model the missing mechanism related to covariate as
close as possible to true in order to obtain valid estimates of
$\boldsymbol{\beta}$. The same is true for the missing data process. All
predictors in this model process were maintained since an overspecification is
better than a underspecification.

Results from four methods are shown in Table 2. The first one is the usual GEE
method using the available data; the second is the weighted method (WGEE)
using model (15) for the weights; the third is the multiple imputation by
chained equation (MIGEE) in the _R_ package _mice_ ; and the fourth, labeled
DRGEE, is the proposed doubly robust method using (15) and (16) for the weight
and the covariate models, respectively. It was used an independent working
correlation.

Table 2: Regression Parameters for the Analgesia in Birth Data | Available |  | WGEE |  | MIGEE |  | DRGEE
---|---|---|---|---|---|---|---
Parameter | Est | SE | P |  | Est | SE | P |  | Est | SE | P |  | Est | SE | P
INTERCEPT1 | 1.796 | 1.308 | 0.170 |  | 1.827 | 1.267 | 0.149 |  | 1.599 | 1.237 | 0.196 |  | 1.649 | 1.315 | 0.210
INTERCEPT2 | 3.302 | 1.314 | 0.012 |  | 3.261 | 1.310 | 0.013 |  | 3.105 | 1.242 | 0.012 |  | 3.069 | 1.351 | 0.023
TIME | -0.182 | 0.286 | 0.524 |  | -0.157 | 0.315 | 0.619 |  | -0.110 | 0.273 | 0.687 |  | -0.192 | 0.293 | 0.511
GROUP | -1.221 | 0.445 | 0.006 |  | -1.244 | 0.439 | 0.005 |  | -1.056 | 0.409 | 0.010 |  | -1.008 | 0.390 | 0.010
AGE | -0.066 | 0.035 | 0.057 |  | -0.072 | 0.033 | 0.027 |  | -0.062 | 0.031 | 0.046 |  | -0.071 | 0.035 | 0.045
DU | -0.362 | 0.158 | 0.022 |  | -0.354 | 0.168 | 0.035 |  | -0.372 | 0.156 | 0.017 |  | -0.340 | 0.166 | 0.040
OXYT(=2) | 0.727 | 0.554 | 0.189 |  | 0.712 | 0.554 | 0.199 |  | 0.954 | 0.542 | 0.078 |  | 0.821 | 0.531 | 0.122
OXYT(=3) | 1.285 | 0.510 | 0.012 |  | 1.431 | 0.503 | 0.004 |  | 1.410 | 0.488 | 0.004 |  | 1.468 | 0.475 | 0.002

TIME effect is non significant for all the four methods. All methods provide
the same conclusion for effects of GROUP. The negative effect for GROUP means
that the chance of women feel mild pain is lower among the group receiving the
remifentanil compared to the peridural group (the estimated odds is
$e^{-1.008}=0.365$ in the doubly robust method). All methods also agree with
respect to the effect of DU. That is, for each increase of $1$ cm of uterine
dilation the chance of the parturient feel mild pain decreases (in the DRGEE
it is $e^{-0.340}=0.712$, for example). It can be noticed that $p$-value for
AGE effect goes from a non-significant of $0.057$ in the standard GEE to a
significant one in DRGEE, as well as for the other two missing data
approaches. The conclusion is that older women have lower chance of
experiencing mild pain than young women. For the OXIT covariate all methods
reached the same conclusion.

##  7 Discussion

When longitudinal ordinal data are of interest, doubly robust estimator is a
nice alternative. Doubly robust method combines ideas from weighting and
imputation and has been applied elsewhere for estimation of means, causal
inference, and in the longitudinal setting for binary response data (Bang &
Robins, (2005), Carpenter et al., (2006), Seaman & Copas, (2009), Chen & Zhou,
(2011), Li et al., (2013)). However, as far as we know, it has not been
investigated for the longitudinal ordinal case. A doubly robust estimator is
attractive in the sense that it needs only the correct specification of at
least one of the models, but not necessarily both. Simulation results have
indicated that, when at least the covariate model or missing data model is
correct, the doubly robust estimators are consistent and present small-sample
bias comparable to single robust alternatives MIGEE or WGEE. The proposed
method presented good coverage probabilities, as well as its competitors but
with a slight larger variance than multiple imputation. Simulation results
also indicated that the bias of doubly robust estimators when both the
covariate model and the missing data model are incorrect was the same
magnitude of misspecified WGEE or MIGEE. We hope that, in practical
applications, none of the predictive models would be grossly misspecified and
then the proposed estimator would have a great potential of reducing the bias
if the MAR assumption is correct.

When the assumed independent working correlation structure differs from the
true underlying structure, there is no price to pay in terms of the
consistency and asymptotic normality of $\boldsymbol{\beta}$, but such a poor
choice may result in loss of efficiency (Molenberghs & Verbeke,, 2005).
However, modeling of the association structure in the presence of missing data
remains a challenge, specially with longitudinal ordinal data, because there
is no direct way of modeling the association parameters. Future research
involves the investigation of the impact of other association structures in
the doubly robust estimates.

In the proposed doubly robust estimator marginal means were modeled by
cumulative logits. This implies a proportional odds model that in some cases
may not be valid. Another future possible extension of the proposed model is,
therefore, to allow non-proportional odds for a subset of the explanatory
variables (Peterson & Harrell Jr,, 1990).

## References

  * Agresti, (2013) Agresti, Alan. 2013.  Categorical data analysis.  John Wiley & Sons.
  * Anestis Touloumis, (2013) Anestis Touloumis. 2013.  SimCorMultRes: Simulating Correlated Multinomial Responses.  Version 1.1.
  * Bang & Robins, (2005) Bang, Heejung, & Robins, James M. 2005.  Doubly robust estimation in missing data and causal inference models.  Biometrics, 61(4), 962–973.
  * Birhanu et al., (2011) Birhanu, Teshome, Molenberghs, Geert, Sotto, Cristina, & Kenward, Michael G. 2011\.  Doubly robust and multiple-imputation-based generalized estimating equations.  Journal of biopharmaceutical statistics, 21(2), 202–225.
  * Carpenter & Kenward, (2013) Carpenter, James R., & Kenward, Michael G. 2013.  Multiple Imputation and its Application.  Wiley & Sons.
  * Carpenter et al., (2006) Carpenter, James R, Kenward, Michael G, & Vansteelandt, Stijn. 2006.  A comparison of multiple imputation and doubly robust estimation for analyses with missing data.  Journal of the Royal Statistical Society: Series A (Statistics in Society), 169(3), 571–584.
  * Chen & Zhou, (2011) Chen, Baojiang, & Zhou, Xiao-Hua. 2011.  Doubly robust estimates for binary longitudinal data analysis with missing response and missing covariates.  Biometrics, 67(3), 830–842.
  * Donneau et al., (2014a) Donneau, AF, Mauer, M, Molenberghs, Geert, & Albert, A. 2014a.  A simulation study comparing multiple imputation methods for incomplete longitudinal ordinal data.  Communications in Statistics-Simulation and Computation.
  * Donneau et al., (2014b) Donneau, Anne-Françoise, Mauer, Murielle, Lambert, Philippe, Molenberghs, Geert, & Albert, Adelin. 2014b.  Simulation-based study comparing multiple imputation methods for non-monotone missing ordinal data in longitudinal settings.  Journal of Biopharmaceutical Statistics.
  * Li et al., (2013) Li, Lingling, Shen, Changyu, Li, Xiaochun, & Robins, James M. 2013.  On weighting approaches for missing data.  Statistical Methods in Medical Research, 22(1), 14–30.
  * Liang & Zeger, (1986) Liang, Kung-Yee, & Zeger, Scott L. 1986.  Longitudinal data analysis using generalized linear models.  Biometrika, 73(1), 13–22.
  * Lipsitz et al., (1994) Lipsitz, Stuart R, Kim, Kyungmann, & Zhao, Lueping. 1994.  Analysis of repeated categorical data using generalized estimating equations.  Statistics in medicine, 13(11), 1149–1163.
  * Little & Rubin, (1987) Little, Roderick JA, & Rubin, Donald B. 1987.  Statistical analysis with missing data.  Vol. 539.  Wiley New York.
  * Little & Rubin, (2002) Little, Roderick JA, & Rubin, Donald B. 2002.  Statistical analysis with missing data.
  * McCullagh, (1980) McCullagh, Peter. 1980.  Regression models for ordinal data.  Journal of the royal statistical society. Series B (Methodological), 109–142.
  * Molenberghs & Verbeke, (2005) Molenberghs, Geert, & Verbeke, Geert. 2005.  Models for discrete longitudinal data.  Springer.
  * Peterson & Harrell Jr, (1990) Peterson, Bercedis, & Harrell Jr, Frank E. 1990.  Partial proportional odds models for ordinal response variables.  Applied Statistics, 205–217.
  * Pierce, (1982) Pierce, Donald A. 1982.  The asymptotic effect of substituting estimators for parameters in certain types of statistics.  The Annals of Statistics, 475–478.
  * Robins et al., (1995) Robins, James M, Rotnitzky, Andrea, & Zhao, Lue Ping. 1995.  Analysis of semiparametric regression models for repeated outcomes in the presence of missing data.  Journal of the American Statistical Association, 90(429), 106–121.
  * Rubin, (1987) Rubin, Donald B. 1987.  Multiple Imputation for Nonresponse in Surveys.  John Wiley.
  * Schafer, (1997) Schafer, Joseph L. 1997.  Analysis of incomplete multivariate data.  CRC press.
  * Scharfstein et al., (1999) Scharfstein, Daniel O, Rotnitzky, Andrea, & Robins, James M. 1999.  Adjusting for nonignorable drop-out using semiparametric nonresponse models.  Journal of the American Statistical Association, 94(448), 1096–1120.
  * Seaman & Copas, (2009) Seaman, Shaun, & Copas, Andrew. 2009.  Doubly robust generalized estimating equations for longitudinal data.  Statistics in medicine, 28(6), 937–955.
  * Toledano & Gatsonis, (1999) Toledano, Alicia Y, & Gatsonis, Constantine. 1999.  Generalized estimating equations for ordinal categorical data: arbitrary patterns of missing responses and missingness in a key covariate.  Biometrics, 55(2), 488–496.
  * Touloumis et al., (2013) Touloumis, Anestis, Agresti, Alan, & Kateri, Maria. 2013.  GEE for Multinomial Responses Using a Local Odds Ratios Parameterization.  Biometrics, 1–8.
  * Tsiatis, (2006) Tsiatis, Anastasios A. 2006.  Semiparametric theory and missing data.  Springer.
  * van Buuren et al., (1999) van Buuren, S., Boshuizen, H. C., & Knook, D. L. 1999.  Multiple Imputation of Missing Blood Pressure Covariates in Survival Analysis.  Statistics in Medicine, 18, 681–694.
  * van Buuren, (2007) van Buuren, Stef. 2007.  Multiple imputation of discrete and continuous data by fully conditional specification.  Statistical methods in medical research, 16(3), 219–242.
  * van Buuren & Groothuis-Oudshoorn, (2011) van Buuren, Stef, & Groothuis-Oudshoorn, Karin. 2011.  mice: Multivariate Imputation by Chained Equations in R.  Journal of Statistical Software, 45(3), 1–67.

##  Appendix A Appendix

###  A.1 Asymptotic Variance

To state the asymptotic properties of $\boldsymbol{\hat{\beta}}$, let
$\boldsymbol{S}_{1i}(\boldsymbol{\beta},\boldsymbol{\psi},\boldsymbol{\gamma})$
be the individual’s contribution to the estimating equations for
$\boldsymbol{\beta}$,
$\boldsymbol{S}_{2i}(\boldsymbol{\psi})$ be the individual’s contribution to
the estimating equations for $\boldsymbol{\psi}$, and
$\boldsymbol{S}_{3i}(\boldsymbol{\gamma})$ be the individual’s contribution to
the estimating equations for $\boldsymbol{\gamma}$.

Define
$\boldsymbol{\Gamma}(\boldsymbol{\beta},\boldsymbol{\psi},\boldsymbol{\gamma})=E\left\\{\partial\boldsymbol{S}_{1i}(\boldsymbol{\beta},\boldsymbol{\psi},\boldsymbol{\gamma})/\partial\boldsymbol{\beta}^{T}\right\\}$,
$\boldsymbol{I}_{12}(\boldsymbol{\beta},\boldsymbol{\psi},\boldsymbol{\gamma})=E\left\\{\partial\boldsymbol{S}_{1i}(\boldsymbol{\beta},\boldsymbol{\psi},\boldsymbol{\gamma})/\partial\boldsymbol{\psi}^{T}\right\\}$,
$\boldsymbol{I}_{13}(\boldsymbol{\beta},\boldsymbol{\psi},\boldsymbol{\gamma})=E\left\\{\partial\boldsymbol{S}_{1i}(\boldsymbol{\beta},\boldsymbol{\psi},\boldsymbol{\gamma})/\partial\boldsymbol{\gamma}^{T}\right\\}$,
$\boldsymbol{I}_{2}(\boldsymbol{\psi})=E\left\\{\partial\boldsymbol{S}_{2i}(\boldsymbol{\psi})/\partial\boldsymbol{\psi}^{T}\right\\}$,
$\boldsymbol{I}_{3}(\boldsymbol{\gamma})=E\left\\{\partial\boldsymbol{S}_{3i}(\boldsymbol{\gamma})/\partial\boldsymbol{\gamma}^{T}\right\\}$,
and
$\boldsymbol{Q}_{i}(\boldsymbol{\beta},\boldsymbol{\psi},\boldsymbol{\gamma})=\boldsymbol{S}_{1i}(\boldsymbol{\beta},\boldsymbol{\psi},\boldsymbol{\gamma})-\boldsymbol{I}_{12}(\boldsymbol{\beta},\boldsymbol{\psi},\boldsymbol{\gamma})\boldsymbol{I}_{2}^{-1}(\boldsymbol{\psi})\boldsymbol{S}_{2i}(\boldsymbol{\psi})-\boldsymbol{I}_{13}(\boldsymbol{\beta},\boldsymbol{\psi},\boldsymbol{\gamma})\boldsymbol{I}_{3}^{-1}(\boldsymbol{\gamma})\boldsymbol{S}_{3i}(\boldsymbol{\gamma})$.

###### Theorem 1

If either the missing data model or the covariate model is correctly
specified, then

$n^{1/2}(\boldsymbol{\hat{\beta}}-\boldsymbol{\beta_{0}})\longrightarrow
N(\boldsymbol{0},\boldsymbol{\Gamma}^{-1}(\boldsymbol{\beta}_{0},\boldsymbol{\psi}_{0},\boldsymbol{\gamma}_{0})\boldsymbol{\Sigma}\left\\{\boldsymbol{\Gamma}^{-1}(\boldsymbol{\beta}_{0},\boldsymbol{\psi}_{0},\boldsymbol{\gamma}_{0})\right\\}^{T}),$
(17)

where $\boldsymbol{\beta}_{0}$ is the true value of $\boldsymbol{\beta}$,
$\boldsymbol{\psi}_{0}$ and $\boldsymbol{\gamma}_{0}$ are the probability
limits of $\boldsymbol{\hat{\psi}}$ and $\boldsymbol{\hat{\gamma}}$, and
$\boldsymbol{\Sigma}=E\left\\{\boldsymbol{Q}_{i}(\boldsymbol{\beta}_{0},\boldsymbol{\psi}_{0},\boldsymbol{\gamma}_{0})\boldsymbol{Q}_{i}^{T}(\boldsymbol{\beta}_{0},\boldsymbol{\psi}_{0},\boldsymbol{\gamma}_{0})\right\\}$.

Inferences for $\boldsymbol{\beta}$ follows by replacing the unknown
quantities in (17) by its consistent estimators. We make use of “generalized
information equality” (Pierce,, 1982) that
$E\left\\{\partial\boldsymbol{S}_{1i}(\boldsymbol{\beta},\boldsymbol{\psi},\boldsymbol{\gamma})/\partial\boldsymbol{\psi}^{T}\right\\}=-E\left\\{\boldsymbol{S}_{1i}(\boldsymbol{\beta},\boldsymbol{\psi},\boldsymbol{\gamma})\boldsymbol{S}_{2i}^{T}(\boldsymbol{\psi})\right\\}$,
and
$E\left\\{\partial\boldsymbol{S}_{1i}(\boldsymbol{\beta},\boldsymbol{\psi},\boldsymbol{\gamma})/\partial\boldsymbol{\gamma}^{T}\right\\}=-E\left\\{\boldsymbol{S}_{1i}(\boldsymbol{\beta},\boldsymbol{\psi},\boldsymbol{\gamma})\boldsymbol{S}_{3i}^{T}(\boldsymbol{\gamma})\right\\}$.
Similarly (Robins et al.,, 1995),
$E\left\\{\partial\boldsymbol{S}_{2i}(\boldsymbol{\psi})/\partial\boldsymbol{\psi}^{T}\right\\}=-Var\left\\{\boldsymbol{S}_{2i}(\boldsymbol{\psi})\right\\}$,
and
$E\left\\{\partial\boldsymbol{S}_{3i}(\boldsymbol{\gamma})/\partial\boldsymbol{\gamma}^{T}\right\\}=-Var\left\\{\boldsymbol{S}_{3i}(\boldsymbol{\gamma})\right\\}$.

The matrix $\boldsymbol{\Gamma}$ is replaced by
$\boldsymbol{\hat{\Gamma}}=n^{-1}\sum_{i=1}^{n}\left\\{\partial\boldsymbol{S}_{1i}(\boldsymbol{\hat{\theta}})/\partial\boldsymbol{\beta}^{T}\right\\}$,
and $\boldsymbol{\Sigma}$ by
$\boldsymbol{\hat{\Sigma}}=n^{-1}\sum_{i=1}^{n}\left\\{\boldsymbol{\hat{Q}}_{i}\boldsymbol{\hat{Q}}_{i}^{T}\right\\}$,
$\boldsymbol{\hat{Q}}_{i}=\boldsymbol{S}_{1i}(\boldsymbol{\hat{\theta}})-\boldsymbol{\hat{I}}_{12}(\boldsymbol{\hat{\theta}})\boldsymbol{\hat{I}}_{2}^{-1}(\boldsymbol{\hat{\psi}})\boldsymbol{S}_{2i}(\boldsymbol{\hat{\psi}})-\boldsymbol{\hat{I}}_{13}(\boldsymbol{\hat{\theta}})\boldsymbol{\hat{I}}_{3}^{-1}(\boldsymbol{\hat{\gamma}})\boldsymbol{S}_{3i}(\boldsymbol{\hat{\gamma}})$,
$\boldsymbol{\hat{I}}_{12}(\boldsymbol{\hat{\theta}})=n^{-1}\sum_{i=1}^{n}\left\\{\partial\boldsymbol{S}_{1i}(\boldsymbol{\hat{\theta}})/\partial\boldsymbol{\psi}^{T}\right\\}$,
$\boldsymbol{\hat{I}}_{13}(\boldsymbol{\hat{\theta}})=n^{-1}\sum_{i=1}^{n}\left\\{\partial\boldsymbol{S}_{1i}(\boldsymbol{\hat{\theta}})/\partial\boldsymbol{\gamma}^{T}\right\\}$,
$\boldsymbol{\hat{I}}_{2}(\boldsymbol{\hat{\psi}})=n^{-1}\sum_{i=1}^{n}\left\\{\partial\boldsymbol{S}_{2i}(\boldsymbol{\hat{\psi}})/\partial\boldsymbol{\psi}^{T}\right\\}$,
$\boldsymbol{\hat{I}}_{3}(\boldsymbol{\hat{\gamma}})=n^{-1}\sum_{i=1}^{n}\left\\{\partial\boldsymbol{S}_{3i}(\boldsymbol{\hat{\gamma}})/\partial\boldsymbol{\gamma}^{T}\right\\}$.

The proof is similar to Chen & Zhou, (2011) and is omitted here.