# Accommodating informative visit times for analysing irregular longitudinal
data: a sensitivity analysis approach with balancing weights estimators

Sean Yiu and Li Su E-mail address<EMAIL_ADDRESS>corresponding
author MRC Biostatistics Unit, School of Clinical Medicine, University of
Cambridge, Cambridge, CB2 0SR, UK

###### Abstract

Irregular visit times in longitudinal studies can jeopardise marginal
regression analyses of longitudinal data by introducing selection bias when
the visit and outcome processes are associated. Inverse intensity weighting is
a useful approach to addressing such selection bias when the _visiting at
random_ assumption is satisfied, i.e., visiting at time $t$ is independent of
the longitudinal outcome at $t$, given the observed covariate and outcome
histories up to $t$. However, the visiting at random assumption is
unverifiable from the observed data, and informative visit times often arise
in practice, e.g., when patients’ visits to clinics are driven by ongoing
disease activities. Therefore, it is necessary to perform sensitivity analyses
for inverse intensity weighted estimators (IIWEs) when the visit times are
likely informative. However, research on such sensitivity analyses is limited
in the literature. In this paper, we propose a new sensitivity analysis
approach to accommodating informative visit times in marginal regression
analysis of irregular longitudinal data. Our sensitivity analysis is anchored
at the visiting at random assumption and can be easily applied to existing
IIWEs using standard software such as the _coxph_ function of the _R_ package
_Survival_. Moreover, we develop novel balancing weights estimators of
regression coefficients by exactly balancing the covariate distributions that
drive the visit and outcome processes to remove the selection bias after
weighting. Simulations show that, under both correct and incorrect model
specifications, our balancing weights estimators perform better than the
existing IIWEs using weights estimated by maximum partial likelihood. We
applied our methods to data from a clinic-based cohort study of psoriatic
arthritis and provide an _R Markdown_ tutorial to demonstrate their
implementation.

_Keywords_ : Balancing weights; Informative observations; Inverse probability
weighting; Marginal model.

##  1 Introduction

###  1.1 Informative visit times

In long-term follow-up studies, irregular longitudinal data arise when
patients are assessed repeatedly over time, but at uncommon and irregularly
spaced time points. One common cause of irregular longitudinal data is
outcome-dependent visit times, where patients with poorer health outcomes have
more frequent visits and shorter gap times between their visits (Pullenayegum
and Lim, 2016). In this case, because observed data are over-represented by
patients with poorer health outcomes, the selection bias due to the
association between the visit and outcome processes needs to be addressed in
marginal regression analyses of the irregular longitudinal data (Pullenayegum
and Lim, 2016).

Several approaches have been proposed to handle irregular longitudinal data;
see a comprehensive review by Pullenayegum and Lim (2016). Among them, inverse
intensity weighting (IIW), which was first proposed by Lin et al. (2004), is a
useful approach when the _visiting at random_ assumption holds, that is, when
visiting at time $t$ is independent of the longitudinal outcome at $t$, given
the observed covariate and outcome histories up to $t$. The aim of IIW is to
remove the selection bias by creating a pseudo-population that is
representative of the study population that could have been observed
continuously from the study baseline until the study end (i.e., the target
population in marginal regression analyses) after weighting. Specifically, the
observed longitudinal data are weighted by the inverse of visit intensities
estimated from a semi-parametric Cox model given the observed history.
Provided that the visiting at random assumption is satisfied and the Cox model
is correctly specified, the inverse intensity weighted estimators (IIWEs) of
the parameters in a marginal regression model of the longitudinal outcome
using generalized estimating equations (GEE) are consistent (Lin et al., 2004;
Bůžková and Lumley, 2009; Pullenayegum and Feldman, 2013).

Unfortunately, the visiting at random assumption is unverifiable from the
observed data and in practice, the visit times are possibly informative for
the outcome process, even after conditioning on the observed history. For
example, patients’ visits to clinics could be partly driven by ongoing disease
activities, in addition to their disease histories. Therefore, it is necessary
to assess the sensitivity of conclusions drawn from the IIWEs to the
violations of the visiting at random assumption. However, there is very
limited research on sensitivity analysis approaches for irregular longitudinal
data with informative visit times in the literature. Smith et al. (2022)
recently proposed a sensitivity analysis approach for an augmented IIWE by
applying exponential tilting to the observed outcome. But they focused on
estimating treatment-specific mean outcomes at fixed times after randomisation
in clinical trial settings, thus their methods are not directly applicable to
marginal regression analyses of irregular longitudinal data. In addition,
Smith et al. (2022) specified their sensitivity parameter in the exponential
tilting function of the observed outcome, instead of through a _selection
function_ incorporated in the model for the visit process, thus it was less
straightforward to anchor the values of their sensitivity parameter to the
observed information (e.g., estimated coefficients of observed covariates in a
Cox model for the visit process). Also using exponential tilting, Wang (2020)
developed sensitivity analyses for clinical trials with informative assessment
times, but a discrete-time framework was used instead and estimation was
through g-computation with fully parametric models. Beyond the clinical trial
settings, there is a lack of research on sensitivity analysis approaches for
IIWEs in marginal regression analyses of irregular longitudinal data. Note
that, apart from IIW, other approaches for handling irregular longitudinal
data, such as joint modelling of the visit and outcome processes Liang et al.
(2009); Sun et al. (2011); Cai et al. (2012), also involve unverifiable
assumptions, and hence cannot avoid the requirement of sensitivity analyses.

In this paper, within the framework of IIWEs, we propose a new sensitivity
analysis approach for marginal regression analyses of irregular longitudinal
data with informative visit times. Specifically, anchored at the visiting at
random assumption and building upon existing IIWEs, we introduce a selection
function into the Cox model for the visit process, similar to the sensitivity
analysis approaches for handling missing longitudinal data Rotnitzky et al.
(1998); Scharfstein et al. (1999); Vansteelandt et al. (2007); Wen and Seaman
(2018). This selection function characterizes the residual dependence of the
visit intensity at time $t$ on the outcome at $t$, given the observed history
up to $t$. A sensitivity parameter (or a vector of sensitivity parameters) is
used to describe the strength of the residual dependence. When the sensitivity
parameter is set at zero, visiting at random is assumed. By varying this
sensitivity parameter over a range of plausible values, the sensitivity of the
substantive conclusions to deviations from the visiting at random assumption
can be assessed. Compared with the methods in Smith et al. (2022), a notable
advantage of the proposed approach is that it can be easily applied to
existing IIWEs using standard software such as the _coxph_ function of the _R_
package _Survival_ (Therneau, 2023).

Furthermore, motivated by improving the robustness and efficiency of IIWEs, we
develop novel balancing weights estimators of regression coefficients by
exactly balancing the covariate distributions that drive the visit and outcome
processes to remove the selection bias after weighting. Balancing weights have
been shown to be more stable than weights obtained by maximum likelihood
estimation (MLE), which dramatically improved the performance of inverse
probability weighted estimators for handling missing data in cross-sectional
settings and dropout in longitudinal studies (e.g., see Graham et al., 2012;
Zubizarreta, 2015; Han, 2016; Tan, 2020; Yiu and Su, 2022). Our simulations
study also show that, under both correct and incorrect model specifications,
the proposed balancing weights estimators perform better than the existing
IIWEs using weights estimated by maximum partial likelihood. The balancing
weights estimators can be easily combined with the proposed sensitivity
analysis approach, which hopefully will enhance the practicality of the IIW
methods for more widespread utilization. To help achieve this, we provide an
_R Markdown_ tutorial to demonstrate the implementation of the proposed
methods (see https://github.com/lisu-stats/IIW_SA).

###  1.2 Motivating data

This research is motivated by data from the University of Toronto Psoriatic
Arthritis (PsA) Clinic cohort. PsA is inflammatory arthritis associated with
the skin disease psoriasis. Manifestations of PsA include the development of
joint activity, which is characterized by the occurrence of pain and/or
swelling in the joints (i.e., active joints). The University of Toronto PsA
clinic is one of the largest cohorts of PsA patients in the world. More than
1000 patients are assessed in the clinic about every 6-12 months, where data
on their demographics, clinical factors, and treatment information are
recorded. Despite the regular scheduling of patients’ visits to the clinic
(6-12 months apart), the actual visit times varied considerably from patient
to patient. Zhu et al. (2017) noted that the gap times between a patient’s
clinic visits appeared to be related to his/her prior disease history, past
visit history and various factors associated with these processes. They
addressed this problem of the PsA cohort data by applying IIW in a parametric
failure time model to examine the association of biologics use with PsA
disease progression. The covariates considered in the visit process model of
Zhu et al. (2017) included demographics, PsA disease duration, erythrocyte
sedimentation rate, histories of various disease activity and progression
variables, treatment use history as well as the median of past visit gap
times. However, as the PsA cohort is clinic-based, patients could visit the
clinic due to a PsA flare at non-scheduled times. In addition, important
socioeconomic factors (e.g., the ability to pay for treatments) that are
indirectly associated with the PsA disease activity and the clinic visits are
not recorded in the PsA cohort. Therefore, there is likely residual dependence
of the visit intensity on the ongoing PsA disease activity. In this paper, we
aim to accommodate the likely informative visit times in the PsA clinic cohort
using the proposed sensitivity analysis approach in a marginal regression
analysis for examining the association between biologics use and PsA active
joint counts over time.

##  2 Methods

###  2.1 Notation, setting and assumptions

We consider a longitudinal study where $n$ patients are enrolled at baseline
(time $0$) and each of the $i=1,\ldots,n$ patients are followed to a maximum
time of $\tau$. Patients’ data were assumed independent and identically
distributed in the study and thus we suppress the patient-specific subscript
$i$ for now and reintroduce it when we describe our sensitivity analysis
approach in Section 2.2. Let $Y(t)$ be a longitudinal outcome measurable at
time $t$ for $0\leq t\leq\tau$. Let $\bm{X}(t)$ denote a vector of exogenous
(possibly time-varying) covariates, which is assumed to be known prior to $t$.
We are interested in a marginal regression model for $Y(t)$,

$g[{\rm
E}\\{Y(t)\mid\bm{X}(t)\\}]=g[\mu\\{t,\bm{X}(t);\bm{\beta}_{0}\\}]=\bm{\beta}_{0}^{\text{T}}\bm{X}(t)$
(1)

for $t\in[0,\tau]$, where $g(\cdot)$ is a known link function and
$\bm{\beta}_{0}$ is a $p$-dimensional vector of regression coefficients.

The marginal regression model in (1) is formulated for the outcome process in
$0\leq t\leq\tau$. However, the outcome measurements are not available
continuously over time but only at visit times. Let $N(t)$ be the counting
process for the number of visits made by time $t$ and
$N(t)=\int_{0}^{t}dN(s)$, where $dN(s)$ is the indicator of a visit made at
time $s$. Let $C$ be the follow-up time measured from the baseline of the
study. Then $\xi(t)=\mbox{I}\\{t<\min(C,\tau)\\}$ is a left-continuous at-risk
process, indicating whether a patient is still under follow-up. We denote the
uncensored visit process by $N^{*}(t)=\xi(t)N(t)$. In addition, let
$\bm{Z}(t)$ be a vector of auxiliary covariates that are associated with the
visit and outcome processes but are excluded from (1) because their
associations with the outcome are not of scientific interest. We assume that
$\bm{Z}(t)$ is known at all times (Bůžková and Lumley, 2009).

Following Bůžková and Lumley (2009) and Pullenayegum and Feldman (2013), we
make the non-informative censoring assumption for the marginal mean of the
outcome

${\rm E}\\{Y(t)\mid\bm{X}(t),t\leq C\leq\tau\\}={\rm
E}\\{\xi(t)Y(t)\mid\bm{X}(t)\\}={\rm E}\\{Y(t)\mid\bm{X}(t)\\},$

such that the marginal mean ${\rm E}\\{Y(t)\mid\bm{X}(t)\\}$ is the same for
those patients whose follow-up is censored at time $t$ and for those who
remain in the study.

Let $H(t)=\\{Y(s):s<t;\bm{X}(s),\bm{Z}(s),N^{*}(s):s\leq t\\}$ include the
outcome process history prior to $t$, the covariate and visit process
histories up to $t$. Let
$O(t)=\\{N^{*}(s)Y(s):s<t;N^{*}(s)\bm{X}(s),\bm{Z}(s),N^{*}(t):s\leq t\\}$
include the _observed_ outcome history prior to $t$, the _observed_ covariate
history and the visit process history up to $t$. We assume that the uncensored
visit process follows

${\rm E}\\{dN^{*}(t)\mid H(t),Y(t),\underline{Y}(t)\\}={\rm E}\\{dN^{*}(t)\mid
H(t),Y(t)\\}={\rm E}\\{dN^{*}(t)\mid O(t),Y(t)\\},$ (2)

where $\underline{Y}(t)=\\{Y(s):s>t\\}$ is the future of the outcome process
beyond $t$. The first equality in (2) is the _non-future dependence_
assumption, which means that conditional on the process history $H(t)$, and
the (possibly unobserved) current value of the outcome process $Y(t)$, the
uncensored visit process is independent of future outcomes beyond time $t$.
The second equality requires that the uncensored visit process be independent
of past unobserved values of the outcome and covariates given observed history
$O(t)$ and the current outcome $Y(t)$. Note that this requirement is _less
restrictive_ compared to the visiting at random assumption which does not
allow the visit process to depend on the current outcome $Y(t)$. By allowing
the visit process to additionally depend on the current outcome $Y(t)$, we are
able to accommodate the common scenario in clinic-based studies or electronic
health record databases where patients’ visits are partly driven by ongoing
disease activities.

Let ${\rm E}\\{dN^{*}(t)\mid O(t),Y(t)\\}=\xi(t)\lambda\\{t,O(t),Y(t)\\}$,
where $\lambda\\{t,O(t),Y(t)\\}$ is the visit intensity at time $t$ that
depends on the observed history and the current outcome. We further make the
positivity assumption for the visit intensity that
$\lambda\\{t,O(t),Y(t)\\}>0$ when $\xi(t)=1$.

###  2.2 Sensitivity analysis for inverse intensity weighted estimators with
informative visit times

Under the assumptions made in Section 2.1 and assuming that the marginal
regression model in (1) is correctly specified, $\bm{\beta}_{0}$ can be
consistently estimated by solving weighted GEE with an independence working
correlation structure of the following form

$\sum_{i=1}^{n}\int_{0}^{\tau}U\\{Y_{i}(t),\bm{X}_{i}(t);\bm{\beta}\\}w_{i}(t)dN^{*}_{i}(t)=0,$
(3)

where
$U\\{Y_{i}(t),\bm{X}_{i}(t);\bm{\beta}\\}=D\\{t,\bm{X}_{i}(t);\bm{\beta}\\}\left[Y_{i}(t)-\mu\\{t,\bm{X}_{i}(t);\bm{\beta}\\}\right]$,
$D\\{t,\bm{X}_{i}(t);\bm{\beta}\\}=\left[\partial\mu\\{t,\bm{X}_{i}(t);\bm{\beta}\\}/\partial\bm{\beta}\right]\mbox{Var}^{-1}\left[t,\mu\\{t,\bm{X}_{i}(t);\bm{\beta}\\}\right]$,
$\mbox{Var}\\{t,\mu(t,\bm{X}_{i}(t);\bm{\beta})\\}$ is a conditional variance
function, $w_{i}(t)=s(t)/\lambda\\{t,O_{i}(t),Y_{i}(t)\\}$ and $s(t)$ is a
user-specified function of $t$ Lin et al. (2004). This is justified because

$\begin{split}&{\rm
E}\left[\int_{0}^{\tau}U\\{Y(t),\bm{X}(t);\bm{\beta}_{0}\\}\frac{s(t)}{\lambda\\{t,O(t),Y(t)\\}}dN^{*}(t)\right]\\\
=&{\rm E}\left(\int_{0}^{\tau}{\rm
E}\left[U\\{Y(t),\bm{X}(t);\bm{\beta}_{0}\\}\frac{s(t)}{\lambda\\{t,O(t),Y(t)\\}}N^{*}(t)\,\middle|\,O(t),Y(t)\right]\right)\\\
=&{\rm
E}\left[\int_{0}^{\tau}U\\{Y(t),\bm{X}(t);\bm{\beta}_{0}\\}\frac{s(t)}{\lambda\\{t,O(t),Y(t)\\}}\,{\rm
E}\left\\{N^{*}(t)\mid O(t),Y(t)\right\\}\right]\\\ =&{\rm
E}\left[\int_{0}^{\tau}U\\{Y(t),\bm{X}(t);\bm{\beta}_{0}\\}s(t)\xi(t)\,dt\right]=0.\end{split}$

Intuitively, the purpose of IIW is to create a pseudo-population that is
representative of the potentially observable population from baseline until
time $\tau$ (i.e., the target population), but with a visit process with an
intensity function $s(t)$ that no longer depends on the observed history and
current outcome. However, the true visit intensity
${\lambda}\\{t,O(t),Y(t)\\}$ is unknown. Even if we make the visiting at
random assumption such that
${\lambda}\\{t,O(t),Y(t)\\}={\lambda}\\{t,O(t)\\}$, the visit intensity given
the observed history ${\lambda}\\{t,O(t)\\}$ still needs to be estimated in
practice.

In the IIW literature, ${\lambda}\\{t,O(t)\\}$ is typically estimated by
specifying a semi-parametric Cox model
${\lambda}\\{t,O(t);\bm{\bm{\gamma}}_{0}\\}={\lambda}_{0}(t)\exp\\{\bm{\bm{\gamma}}^{\text{T}}_{0}\widetilde{\bm{Z}}(t)\\}$,
where ${\lambda}_{0}(t)$ is an unspecified baseline intensity function of $t$,
$\widetilde{\bm{Z}}(t)$ are functions of $O(t)$ that may contain interactions
and transformations of the variables Bůžková and Lumley (2009); Pullenayegum
and Feldman (2013). Then, the regression parameter vector
$\bm{\bm{\gamma}}_{0}$ is estimated by solving the score equations of the Cox
partial likelihood

$\sum_{i=1}^{n}\int_{0}^{\tau}\left\\{\widetilde{\bm{Z}}_{i}(t)-\frac{\sum_{l=1}^{n}\xi_{l}(t)\widetilde{\bm{Z}}_{l}(t)\exp\\{\bm{\bm{\gamma}}^{\text{T}}\widetilde{\bm{Z}}_{l}(t)\\}}{\sum_{l=1}^{n}\xi_{l}(t)\exp\\{\bm{\bm{\gamma}}^{\text{T}}\widetilde{\bm{Z}}_{l}(t)\\}}\right\\}dN^{*}_{i}(t)=0.$
(4)

${{\lambda}}_{0}(t)$ can be estimated by smoothing Breslow’s estimate of the
cumulative baseline intensity ${\Lambda}_{0}(t)$ Lin et al. (2004),

$\hat{{\Lambda}}_{0}(t;\hat{\bm{{\gamma}}})=\int_{0}^{t}\frac{\sum_{i=1}^{n}dN^{*}_{i}(s)}{\sum_{i=1}^{n}\xi_{i}(s)\exp\\{\hat{\bm{\gamma}}^{\text{T}}\widetilde{\bm{Z}}_{i}(s)\\}},$
(5)

where $\hat{\bm{\bm{\gamma}}}$ is obtained from solving (4). Alternatively, if
we set $s(t)={\lambda}_{0}(t)$, then $w_{i}(t)$ can be estimated by
$\exp\\{-\hat{\bm{\bm{\gamma}}}^{\text{T}}\widetilde{\bm{Z}}_{i}(t)\\}$
Bůžková and Lumley (2009). The estimator in (3) is consistent when the
visiting at random assumption is satisfied and there is no model
misspecification in the Cox model for ${\lambda}\\{t,O(t)\\}$.

In our setting with informative visit times, however,
${\lambda}\\{t,O(t),Y(t)\\}$ cannot be identified and estimated from the
observed data because $Y(t)$ is only observed when a visit is made. We instead
propose a sensitivity analysis approach by assuming that

${\lambda}\\{t,O(t),Y(t);\bm{\gamma}_{s0},\bm{\phi}\\}={\lambda}_{0}(t)\exp[\bm{\bm{\gamma}}_{s0}^{\text{T}}\widetilde{\bm{Z}}(t)+q\\{O(t),Y(t);\bm{\phi}\\}],$
(6)

where $q\\{O(t),Y(t);\bm{\phi}\\}$ is a known selection function with a known
sensitivity parameter vector $\bm{\phi}$. For example, in the PsA data
analysis reported in Section 4, we specify
$q\\{O(t),Y(t);\bm{\phi}\\}=\phi\log\\{Y(t)+1\\}$ with a single sensitivity
parameter $\phi$ that characterizes the residual dependence of the visit
intensity on the current (log-transformed) PsA active joint count, after
adjusting for the observed history.

Let $Q(t;\bm{\phi})=\exp[-q\\{O(t),Y(t);\bm{\phi}\\}]$. It is easy to see that
if the model in (6) is correctly specified,

${\rm E}\left\\{Q(t;\bm{\phi})dN^{*}(t)\mid
O(t),Y(t)\right\\}=\xi(t){\lambda}_{0}(t)\exp\\{\bm{\bm{\gamma}}^{\text{T}}_{s0}\widetilde{\bm{Z}}(t)\\}.$

Thus, to estimate $\bm{\gamma}_{s0}$, we can replace $dN^{*}_{i}(t)$ by
$Q_{i}(t;\bm{\phi})dN^{*}_{i}(t)$ in the estimating equation in (4),

$\sum_{i=1}^{n}\int_{0}^{\tau}\left\\{\widetilde{\bm{Z}}_{i}(t)-\frac{\sum_{l=1}^{n}\xi_{l}(t)\widetilde{\bm{Z}}_{l}(t)\exp\\{\bm{\bm{\gamma}}_{s}^{\text{T}}\widetilde{\bm{Z}}_{l}(t)\\}}{\sum_{l=1}^{n}\xi_{l}(t)\exp\\{\bm{\bm{\gamma}}_{s}^{\text{T}}\widetilde{\bm{Z}}_{l}(t)\\}}\right\\}Q_{i}(t;\bm{\phi})dN^{*}_{i}(t)=0.$
(7)

Similarly, the Breslow estimator is modified as

$\hat{{\Lambda}}_{0}(t;\hat{\bm{{\gamma}}}_{s})=\int_{0}^{t}\frac{\sum_{i=1}^{n}Q_{i}(s;\bm{\phi})dN^{*}_{i}(s)}{\sum_{i=1}^{n}\xi_{i}(s)\exp\\{\hat{\bm{\gamma}}_{s}^{\text{T}}\widetilde{\bm{Z}}_{i}(s)\\}},$
(8)

where $\hat{\bm{\gamma}}_{s}$ is the solution to (7). Note that
$Q_{i}(t;\bm{\phi})$, or equivalently the selection function
$q\\{O_{i}(t),Y_{i}(t);\bm{\phi}\\}$, only needs to be evaluated when
$dN^{*}_{i}(t)=1$, i.e., when a patient makes a visit. When
$q\\{O_{i}(t),Y_{i}(t);\bm{\phi}\\}=0$, (7) reduce to the form in (4) under
the visiting at random assumption. We then update (3) by setting
$s(t)=\lambda_{0}(t)$ and replacing $w_{i}(t)$ with
$\hat{w}_{i}(t;\bm{\phi})=\exp[-\hat{\bm{\gamma}}_{s}^{\text{T}}\widetilde{\bm{Z}}_{i}(t)]Q_{i}(t;\bm{\phi})$.
For each set of fixed values of $\bm{\phi}$, we can repeat the estimation of
$\bm{\beta}_{0}$ using different sets of weights $\hat{w}_{i}(t;\bm{\phi})$ in
(3) and evaluate the sensitivity of marginal regression analysis results under
informative visit times.

We provide details about how to estimate $\bm{\gamma}_{s0}$ and obtain
$\hat{w}_{i}(t;\bm{\phi})$ using R (R Development Core Team, 2014) in Section
2.4. To distinguish $\hat{w}_{i}(t;\bm{\phi})$ from the covariate balancing
weights proposed in the next section, we call $\hat{w}_{i}(t;\bm{\phi})$ the
‘MLE weights’ for convenience, even though (7) are not score equations for
$\bm{\gamma}_{s0}$.

###  2.3 Balancing weights estimators under informative visit times

To improve the robustness and efficiency of the IIWEs, we also consider an
alternative weight estimation approach by covariate balancing and develop
novel balancing weights estimators of the regression coefficients.
Specifically, we propose covariate balancing inverse intensity weights
(‘balancing weights’ in short) of the form

$W(t;\bm{\gamma}_{\text{b0}},\bm{\phi})=\exp\left[\bm{\gamma}_{\text{b0}}^{\text{T}}h\\{t,O(t)\\}\right]Q(t;\bm{\phi}),$
(9)

where $\bm{\gamma}_{\text{b0}}$ is a vector of unknown parameters,
$h\\{t,O(t)\\}$ is a vector of functionals of the observed history $O(t)$,
which can contain functions of time $t$ and the interactions between $t$ and
the observed history variables.

We propose the following balancing conditions for estimating
$\bm{\gamma}_{\text{b0}}$ and hence $W(t;\bm{\gamma}_{\text{b0}},\bm{\phi})$,

$\displaystyle\sum_{i=1}^{n}\int_{0}^{\tau}h\\{t,O_{i}(t)\\}\left\\{W_{i}(t;\bm{\gamma}_{\text{b}},\bm{\phi})dN^{*}_{i}(t)-\frac{\xi_{i}(t)\sum_{l=1}^{n}Q_{l}(t;\bm{\phi})dN^{*}_{l}(t)}{\sum_{l=1}^{n}\xi_{l}(t)\exp\\{\hat{\bm{\gamma}}_{s}^{\text{T}}\widetilde{\bm{Z}}_{l}(t)\\}}\right\\}=0,$
(10)

where
$\frac{\sum_{l=1}^{n}Q_{l}(t;\bm{\phi})dN^{*}_{l}(t)}{\sum_{l=1}^{n}\xi_{l}(t)\exp\\{\hat{\bm{\gamma}}_{s}^{\text{T}}\widetilde{\bm{Z}}_{l}(t)\\}}$
is the increment of the modified Breslow estimator in (8) at time $t$, which
is used as an approximation of ${\lambda}_{0}(t)dt$. Basically, (10) ensures
that weighting $dN^{*}_{i}(t)$ with the balancing weights
$W_{i}(t;\bm{\gamma}_{\text{b}},\bm{\phi})$ creates a pseudo-population that
is representative of a target _at-risk_ population at $t$ with a visit
intensity equal to the uncensored baseline intensity $\xi(t){\lambda}_{0}(t)$,
in terms of the functionals of observed history variables $h\\{t,O_{i}(t)\\}$.
In other words, the distributions of $h\\{t,O_{i}(t)\\}$ over time are
_exactly_ balanced between the weighted observed sample and the target at-risk
population with uncensored visit intensity $\xi(t){\lambda}_{0}(t)$. Thus the
selection bias due to $h\\{t,O_{i}(t)\\}$ is removed after weighting.

Assuming that the true visit intensity follows the model in (6), it is easy to
show that the expectations of the left-hand side of (10) are equal to zero if
$W_{i}(t;\bm{\gamma}_{\text{b}},\bm{\phi})$ is replaced by the true inverse
intensity weight
$\exp[-{\bm{\gamma}}_{s0}^{\text{T}}\widetilde{\bm{Z}}_{i}(t)]Q_{i}(t;\bm{\phi})$
and $\exp\\{\hat{\bm{\gamma}}_{s}^{\text{T}}\widetilde{\bm{Z}}_{l}(t)\\}$ is
replaced by $\exp\\{\bm{\gamma}_{s0}^{\text{T}}\widetilde{\bm{Z}}_{l}(t)\\}$.
Thus the estimator of the balancing weights in (10) is consistent for the true
inverse intensity weights if the weight model in (9) includes the correct
covariates $\widetilde{\bm{Z}}_{i}(t)$, the selection function is correctly
specified and the Cox model whose estimates are used in the modified Breslow
estimator is correctly specified. As a result, the IIWE in (3) using the
balancing weights is consistent for $\bm{\beta}_{0}$.

For the choice of $h\\{t,O_{i}(t)\\}$, the same set of variables
$\widetilde{\bm{Z}}_{i}(t)$ as in (6) can be included. In addition, to prevent
extreme weights, we recommend that $h\\{t,O_{i}(t)\\}$ contains $1$ such that
(10) imposes a constraint on the sum of the balancing weights. Moreover, we
can add indicator variables that define a partition of the study follow-up
period to $h\\{t,O_{i}(t)\\}$ such that the sums of the balancing weights are
bounded within the partitioned periods. Ideally, for unbiased estimation of
$\bm{\beta}_{0}$, we would like to impose balancing conditions in (10) to
ensure that the distributions of observed variables are _exactly_ balanced for
all potential visit times. Clearly, this is not possible since $t$ is
continuous and an infinite number of constraints are required. Instead, we
recommend including the interactions between time $t$ (either as a continuous
variable or as partitioned time periods) and $\widetilde{\bm{Z}}_{i}(t)$ for
balancing to provide some parsimony for the conditions in (10).

###  2.4 Implementation

For fixed $\bm{\phi}$, the estimation of $\bm{\gamma}_{s0}$ using (7) can be
implemented by the coxph function in the R package survival. We use a trick by
defining $Q_{i}(t;\bm{\phi})$ as case weights if $dN_{i}^{*}(t)=1$, while the
case weights are equal to $1$ if $dN_{i}^{*}(t)=0$. We add an offset term in
the model formula of the coxph function such that if $dN_{i}^{*}(t)=1$ the
offset is $-\log\\{Q_{i}(t;\bm{\phi})\\}$ and if $dN_{i}^{*}(t)=0$ the offset
is $0$. Using these offsets prevents the coxph function to recalculate the
weighted sum of $\widetilde{\bm{Z}}_{i}(t)$ in score functions of the Cox
partial likelihood using the case weights we just defined. As a results, the
coxph function can implement the estimation using the estimating equation in
(4). The balancing weights in (10) are estimated using the R package nleqslv
(Hasselman, 2023). The estimation of $\bm{\beta}_{0}$ can be achieved using
the glm or lm function in R with the MLE weights or the balancing weights.
Finally, we use nonparametric bootstrap for inference, while for small samples
jackknife can also be applied. An R Markdown tutorial to demonstrate the
implementation of the proposed can be found at https://github.com/lisu-
stats/IIW_SA.

##  3 Simulation

We conduct a simulation study to assess the performance of the proposed
methods in finite samples. Specifically, we compare the IIWEs with the MLE
weights and the balancing weights estimators under the combinations of the
following scenarios: (1) The true visit process depends on baseline and time-
varying covariates as well as the current outcome, but the selection function
is omitted (i.e., assuming visiting at random) or included (i.e.,
accommodating informative visit times) when estimating the weights. (2) The
visit process is highly or moderately dependent on time-varying covariates.
(3) Correct or incorrect (transformed) covariates are included in the visit
process models for estimating the weights.

###  3.1 Data generating mechanism

The data generating mechanism for one patient is summarized in Table 1. We
omit the subscript $i$ for clearer presentation. Both longitudinal count and
continuous outcomes are considered as the analysis of the PsA cohort data in
Section 4 focuses on a longitudinal count outcome.

In this simulation setup, we have one baseline covariate $X$, two time-varying
covariates $Z_{1}(t)$ and $Z_{2}(t)$ as well as the time variable $t$, which
affect visiting at $t$ and also the outcome $Y(t)$. The dependence of the
visit process on the time-varying covariates is characterized by the parameter
$\gamma_{z}$ in the visit process model, where we set $\gamma_{z}=0.5,1.25$ to
allow both moderate and high levels of such dependence. We use the Bernoulli
distribution to approximate simulated data from a Cox model for the visit
process since the visit/event rate is low.

In addition, the current outcome affects the visit process. We let
$Y^{*}(t)=Y(t)$ if the outcome is continuous and $Y^{*}(t)=\log\\{Y(t)+1\\}$
for the case of the longitudinal count outcome. The dependence on $Y^{*}(t)$
is characterized by the sensitivity parameter $\phi$. We set
$\phi=0.15,0.3,0.45,0.6$ such that those patients with larger outcome values
are more likely to make a visit regardless of their observed histories.

For the longitudinal continuous outcome, our aim is to estimate the regression
parameters $\beta_{1}$ and $\beta_{2}$ in the marginal regression model ${\rm
E}\\{Y(t)\mid X\\}=\beta_{0}+\beta_{1}X+\beta_{2}t$. The true values of
$\beta_{1}$ and $\beta_{2}$ are $-4.5$ and $-0.5$, respectively. For the
longitudinal count outcome, we are interested in the parameters
$\tilde{\beta}_{1}$ and $\tilde{\beta}_{2}$ in the log-linear model for the
outcome mean $\log[{\rm E}\\{Y(t)\mid
X\\}]=\tilde{\beta}_{0}+\tilde{\beta}_{1}X+\tilde{\beta}_{2}t$. The true
values of $\tilde{\beta}_{1}$ and $\tilde{\beta}_{2}$ are $-1$ and $-0.5$,
respectively.

Data from each patient are generated independently. We simulate 1000 data sets
with different sample sizes ($n=200,500,1000$) for each scenario with
different values of $\gamma_{z}$ and $\phi$ for both types of outcomes.

###  3.2 Model specification and estimators

We consider the performance of the IIWEs with the MLE weights and the
balancing weights estimators under both correct and incorrect model
specifications for the visit process. There are two types of model
misspecification. The first is to assume visiting at random but in fact, the
visit times are informative with $\phi>0$. The second is functional form
misspecification by including transformed time-varying covariates
$Z^{*}_{1}(t)$, $Z^{*}_{2}(t)$ and their interaction (see Table 1), instead of
$Z_{1}(t)$, $Z_{2}(t)$ and $Z_{1}(t)Z_{2}(t)$, in the models for estimating
the weights.

Specifically, we will evaluate the performance of the IIWEs using both the MLE
weights and the balancing weights estimators, under the following scenarios:

  1. 1)

The selection function $\phi Y^{*}(t)$ is omitted, correct covariates
$Z_{1}(t)$, $Z_{2}(t)$, $Z_{1}(t)Z_{2}(t)$, and $X$ are included.

  2. 2)

The selection function $\phi Y^{*}(t)$ is omitted, transformed covariates
$Z^{*}_{1}(t)$, $Z^{*}_{2}(t)$, $Z^{*}_{1}(t)Z_{2}^{*}(t)$, and $X$ are
included.

  3. 3)

The selection function $\phi Y^{*}(t)$, correct covariates $Z_{1}(t)$,
$Z_{2}(t)$, $Z_{1}(t)Z_{2}(t)$, and $X$ are included.

  4. 4)

The selection function $\phi Y^{*}(t)$, transformed covariates $Z^{*}_{1}(t)$,
$Z^{*}_{2}(t)$, $Z^{*}_{1}(t)Z_{2}^{*}(t)$, and $X$ are included.

The MLE weights and the balancing weights are estimated using the methods
described in Section 2. When the current outcome $Y^{*}(t)$ is accommodated
when estimating the weights, we use the true value of the sensitivity
parameter $\phi$ if the correct covariates are included. When the transformed
covariates are instead used, we fix $\phi$ at the limiting value calculated by
fitting a Cox model for the visit process with transformed covariates and
$Y^{*}(t)$ ($\forall t$) to a large data set with $n=1\times 10^{5}$. For
estimating the balancing weights, we include $1$, the visit times $t$ (treated
as a continuous variable) and their interactions with the baseline and time-
varying covariates in the set of variables for balancing.

Finally, the models for the marginal mean of $Y(t)$ are correctly specified,
where the weighted GEEs in (3) with different sets of weights are used to
estimate the parameters $\beta_{0}$, $\beta_{1}$, $\beta_{2}$,
$\tilde{\beta}_{0}$, $\tilde{\beta}_{1}$ and $\tilde{\beta}_{2}$. For
comparison, we also perform the analysis that uses the complete data at all
possible visit times and the naïve analysis that uses the observed data
without weighting.

###  3.3 Results

Figure 1 presents the empirical bias and root mean squared error (RMSE) of the
IIWEs with the MLE weights and the balancing weights estimators for
$\beta_{1}$ and $\beta_{2}$ in the marginal regression model of a longitudinal
continuous outcome when the sample size $n=200$, the sensitivity parameter
$\phi$ was set at various values and the visit process was _highly_ dependant
on time-varying covariates. It is clear that the estimators for $\beta_{1}$
with the balancing weights had smaller RMSEs than their counterparts with the
MLE weights in all scenarios compared, while they also had smaller biases (in
most of the scenarios) or similar biases. Notably, the estimator of
$\beta_{1}$ with the MLE weights even had larger bias and RMSE than the naïve
estimator, when $\phi=0.15$ and incorrect covariates are included along with
the selection function.

The contrast of the performance of the IIWEs with the MLE weights and the
balancing weights estimators is more striking for $\beta_{2}$, where the
balancing weights estimators performed consistently better in terms of bias
and RMSE than their counterparts with the MLE weights. We also note that the
estimators with the MLE weights had large biases and RMSEs even if the visit
model was correctly specified with correct covariates and true selection
function, and their performance was even worse than the naïve estimator. This
is not surprising because the visit process was highly dependent on time-
varying covariates, thus extreme MLE weights were more likely to occur
(especially under model misspecification), which exacerbated the performance
of the corresponding IIWEs. Covariate balancing weights have been shown that
they can bound the RMSEs of inverse probability weighted estimators under
model misspecification in the settings with point treatments or missing data
in cross-sectional studies (Tan, 2020). Our simulation results confirmed that
the proposed balancing weights also reduced the RMSEs of the weighted
estimators when there was model misspecification and the visit process was
highly dependent on the time-varying covariates.

Interestingly, for both $\beta_{1}$ and $\beta_{2}$, the large RMSEs of the
estimators with the MLE weights were increasing as the value of the
sensitivity parameter $\phi$ decreased. This could be due to that the
variation in the visit process was dominated by the time-varying covariates
when $\phi$ was small, thus leading to more extreme MLE weights. This finding
also suggests that the balancing weights estimators have the most potential
for improved performance when the visit process is highly dependent on time-
varying covariates but weakly dependent on the current outcome (i.e. when the
visiting at random is nearly satisfied).

Figure 2 presents the results when the visit process was _moderately_
dependent on the time-varying covariates with the sample size $n=200$. We also
note that when correct covariates were included and the value of the
sensitivity parameter $\phi$ was not large ($\phi=0.15,0.3$), the biases and
RMSEs of the estimators with the MLE weights were much reduced compared to the
setting with the visit process highly dependent on time-varying covariates.
But when there was model misspecification (omission of the selection function
and inclusion of incorrect covariates), the reduction of the biases and RMSEs
of the estimators with the MLE weights was not obvious. On the other hand,
their counterparts with the balancing weights still performed consistently
better regardless of the associations of the visit process with the time-
varying covariates. This suggests that the performance of the balancing
weights estimators was less sensitive to the true setup of the visit process.

Similar results for the longitudinal continuous outcome with the sample size
$n=500$ and $n=1000$ can be found in Figures 1-4 of the Supplementary
Material.

Results for the longitudinal count outcome can be found in Figures 5-10 of the
Supplementary Material. When the visit process was highly dependent on time-
varying covariates, the balancing weights estimators still had much better
performance than their counterparts with the MLE weights. When the visit
process was moderately dependent on time-varying covariates, the improvements
brought by the balancing weights were less obvious.

###  3.4 Summary

Overall, our simulations suggest that the balancing weights estimators have
better performance than the IIWEs with the MLE weights under both correct and
incorrect model specifications. The balancing weights have the most potential
for improving IIWEs when the true visit process is highly dependent on time-
varying covariates but weakly dependent on the current outcome.

##  4 Analysis of the PsA clinic data

###  4.1 The PsA Clinic cohort

Recall that we are interested in examining the associations between biologics
use and PsA active joint counts in the PsA clinic cohort. Following Zhu et al.
(2017), we created a sub-cohort by including a PsA patient who started his/her
first course of biologics on a specific calendar day and another randomly
selected patient who had never taken biologics as of the same calendar day and
then assembling data from all such patients during the year 1981-2014. This
sub-cohort contained 414 patients, with the baseline defined at the calendar
days when the patients started their first courses of biologics. The cut-off
date of this sub-cohort was December 2014. Overall, the patients with and
without baseline biologics use had an average of 9.75 and 9.31 clinic visits,
with mean inter-visit gap times of 7.29 and 7.44 months (standard deviation
3.07 and 3.45), respectively.

###  4.2 Models and estimators

Since substantial over-dispersion was observed for active joint counts due to
a large proportion (about 50$\%$) of zero counts, we specified a negative
binomial model with the log link for the active joint count $Y_{i}(t)$, where
the scale of time $t$ was in calendar years. The positive active joint counts
ranged from 1-59, with about 10$\%$ of the counts being above 12. For
covariates, we included the indicators of the first year, second year and so
on until after the fifth year since baseline, the biologic use at baseline
(yes/no) and the interactions between baseline biologics use and the yearly
indicators. In addition, we included indicators for calendar year periods
$[1980,2008)$, $[2008,2012)$ and $[2012,2015)$ to account for changing
clinical practice, where roughly similar numbers of visits were made during
these time periods. Finally, the following baseline covariates were considered
in the marginal regression model: gender, PsA disease duration, age, use of
non-steroidal anti-inflammatory drugs (NSAIDs) and disease-modifying
antirheumatic drugs (DMARDs), erythrocyte sedimentation rate (ESR) with square
root transformation, active joint count and damaged joint count (transformed
using $\log(x+1)$).

For the Cox model for the visit process, we included the following variables:
gender, PsA disease duration at baseline, ESR at baseline, active joint count
and damaged joint count at baseline, PsA duration and age at $t$; the
interactions between calendar year periods and the following variables—ESR at
the last visit before $t$ (with square root transformation), active joint
count and damaged joint count at the last visit before $t$ (transformed using
$\log(x+1)$), NSAID use, DMARD use, biologics use, the median of previous
inter-visit times at $t$. We included the same set of variables for estimating
the balancing weights, where additionally the indicators of calendar year
periods were also included for balancing.

Note that for important clinical variables such as the active joint count and
damaged joint count, it may be tempting to incorporate the available
information for the current values of the clinical variables (e.g., the active
joint count at time $t$) into a covariate of the visit intensity model by
applying the last observation-carry-forward method only to the patients who
_did not_ make the clinic visits at $t$ and using the current values of the
clinical variables at $t$ for those who did make the visit at $t$. However,
this approach may provide an unrealistic imputation of the clinical variables
at time $t$, particularly when the clinical variables could dramatically
change over time. Thus this approach is unlikely to adequately address the
selection bias problem from informative visit times. Recall that our proposed
sensitivity analysis approach avoids the need to make assumptions on the
active joint count process between visits since the current active joint count
(when known) is included in the selection function rather than as a covariate
in the visit intensity model.

In our analysis, the main parameters of interest were the conditional
associations between biologics use and active joint counts over time. Let
$\beta_{j}$ ($j=1,\ldots,6$) denote the conditional associations between
biologics use and active joint counts in the $j$th year since baseline. We
applied three estimators to estimate $\beta_{1},\ldots,\beta_{6}$, while
adjusting for other covariates listed above: (1) the GEE estimators without
weighting (the naïve estimator); (2) the IIWEs with the MLE weights; (3) the
balancing weights estimators.

For sensitivity analysis, we set the selection function as
$q\\{O(t),Y(t);\bm{\phi}\\}=\phi\log\\{Y(t)+1\\}$, where $\phi$ was the
sensitivity parameter set at $0$, $0.5$ and $1$. By letting $\phi>0$, we
assumed that the patients currently having a larger number of active joints
were more likely to make a clinic visit.

Standard errors of the parameter estimates were obtained by jackknife, as we
had a large number of variables included in the analysis and this caused
convergence problems when implementing the nonparametric bootstrap.

###  4.3 Results

Table 2 presented the results of the fitted Cox model for the visit process in
the PsA cohort, assuming the visiting at random assumption was satisfied
(i.e., $\phi=0$). It appears that patients with more active joints at the
previous visit were more likely to visit, but those with more damaged joints
at the previous visit were less likely to visit. Patients with longer median
previous inter-visit gap times were less likely to visit. Finally, visits
occurred more frequently as time elapsed since baseline.

Table 3 presented the estimates and standard errors of
$\beta_{1},\ldots,\beta_{6}$ using the naïve estimator, the IIWEs with the MLE
weights and the balancing weights estimators under the different values of the
sensitivity parameter. We note that the results of the IIWEs using the MLE
weights and the balancing weights estimators were almost identical, which
might be because the estimated MLE weights and balancing weights were not
drastically different. For example, when $\phi=0$, the estimated MLE weights
had a minimum of 0.24, a median of 0.63, and a maximum of 4.56, while the
minimum, median and maximum of the balancing weights were 0.24, 0.63 and 3.51,
respectively. Also, the estimated hazard ratios in the fitted Cox model from
Table 2 suggest that the visit process was not highly dependent on the
observed history variables. Therefore, the IIWEs with the MLE weights and the
balancing weights estimators had similar results.

The results based on the naïve estimator suggested that patients who took
biologics at baseline had lower active joint counts around 2-5 years in
follow-up, while the results based on the IIWEs assuming $\phi=0$ were
similar. When the sensitivity parameter was set at $\phi=0.5$, the negative
associations between biologics use and active joint counts were more prominent
throughout the follow-up. This pattern continued when $\phi=1$, where the
magnitude of the negative associations was even larger; and for the first
year, the fourth year and the fifth year, the associations reached the $5\%$
significance level. Since the active joint count is an important clinical
marker for PsA disease progression, our findings could inform further cost-
benefit analyses of PsA treatment options (D’Angiolella et al., 2018).

##  5 Conclusion and discussion

In this paper, we proposed a new sensitivity analysis approach for IIWEs with
informative visit times in marginal regression analyses of irregular
longitudinal data. In particular, we developed novel balancing weights
estimators of regression coefficients that had improved performance
demonstrated in our simulation study. Our proposed sensitivity analysis
approach can be implemented using standard software and we provided an R
Markdown tutorial to facilitate the practical implementation of our methods.

The proposed sensitivity analysis approach can be extended to augmented IIWEs,
where the standard IIWEs with the MLE weights were combined with an outcome
imputation model to form doubly robust estimators (Pullenayegum and Feldman,
2013). Apart from the selection function included in the model for the visit
process, the outcome imputation model for the augmented IIWEs needs to be
fitted by weighting each outcome observation with $Q(t;\bm{\phi})$. This is
similar to the outcome imputation included in the augmented inverse
probability weighted estimators for non-ignorable non-monotone missing data
(Vansteelandt et al., 2007; Wen and Seaman, 2018).

In addition, there are several future research directions for the proposed
methods. First, the specification of the selection function and the choice of
the values of the sensitivity parameter need to be anchored at the observed
information, following substantive input elicited from clinical experts. This
warrants further research. Second, our methods might be computationally
infeasible for large data sets because it requires a Cox model to be fitted.
Moreover, it is time-consuming to use resampling methods such as bootstrap and
jackknife for inference and repeat analyses at different values of the
sensitivity parameter. Therefore, it is important to improve the computational
efficiency of the proposed methods in future work. Third, aiming to reduce
bias, the balancing weights estimators could become inefficient when a large
number of variables are balanced, some of which explain a small amount of
variation of the outcome process. Therefore it would be worth considering a
more parsimonious set of variables to balance, so as to trade off a small
increase in bias for a large increase in efficiency. This is sensible in our
sensitivity analysis setting because the balancing weights are for addressing
the selection bias from the observed variables only and cannot handle the
selection bias from the possibly unobserved current outcome. We could consider
an approach similar to the recently proposed approximate balancing weight
methods (Chattopadhyay et al., 2020) or we could balance a summary score of
the variables that are predictive of the outcome, e.g., the predicted outcome
from an outcome imputation model; see Han (2016) for a similar approach for
longitudinal data with ignorable dropout.

## Acknowledgements

The authors would like to thank Dr Brian Tom for helpful discussions on the
PsA data set and Dr Yayuan Zhu for providing the R code to create the PsA data
set. This research was funded by the Medical Research Council [Unit programme
numbers: MC_UU_00002/8, MC_UU_00002/10, MC_UU_00002/15].

## References

  * Bůžková and Lumley (2009) Bůžková, P. and Lumley, T. (2009) Semiparametric modeling of repeated measurements under outcome-dependent follow-up.  Statistics in Medicine, 28, 987–1003.
  * Cai et al. (2012) Cai, N., Lu, W. and Zhang, H. H. (2012) Time-varying latent effect model for longitudinal data with informative observation times.  Biometrics, 68, 1093–1102.  URL: https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1541-0420.2012.01794.x.
  * Chattopadhyay et al. (2020) Chattopadhyay, A., Hase, C. H. and Zubizarreta, J. R. (2020) Balancing vs modeling approaches to weighting in practice.  Statistics in Medicine, 39, 3227–3254.  URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.8659.
  * D’Angiolella et al. (2018) D’Angiolella, L. S., Cortesi, P. A., Lafranconi, A., Micale, M., Mangano, S., Cesana, G. and Mantovani, L. G. (2018) Cost and cost effectiveness of treatments for psoriatic arthritis: A systematic literature review.  PharmacoEconomics, 36, 567–589.
  * Graham et al. (2012) Graham, B. S., Campos de Xavier Pinto, C. and Egel, D. (2012) Inverse probability tilting for moment condition models with missing data.  Rev. Econ. Stud., 79, 1053–1079.
  * Han (2016) Han, P. (2016) Intrinsic efficiency and multiple robustness in longitudinal studies with drop-out.  Biometrika, 103, 683–700.
  * Hasselman (2023) Hasselman, B. (2023) nleqslv: Solve Systems of Nonlinear Equations.  URL: https://cran.r-project.org/package=nleqslv.  R package version 3.3.4, https://cran.r-project.org/package=nleqslv.
  * Liang et al. (2009) Liang, Y., Lu, W. and Ying, Z. (2009) Joint modeling and analysis of longitudinal data with informative observation times.  Biometrics, 65, 377–384.
  * Lin et al. (2004) Lin, H., Scharfstein, D. O. and Rosenheck, R. A. (2004) Analysis of longitudinal data with irregular, outcome-dependent follow-up.  Journal of the Royal Statistical Society, Series B (Statistical Methodology), 66, 791–813.
  * Pullenayegum and Feldman (2013) Pullenayegum, E. M. and Feldman, B. M. (2013) Doubly robust estimation, optimally truncated inverse-intensity weighting and increment-based methods for the analysis of irregularly observed longitudinal data.  Statistics in Medicine, 15, 1054–1072.
  * Pullenayegum and Lim (2016) Pullenayegum, E. M. and Lim, L. S. (2016) Longitudinal data subject to irregular observation: A review of methods with a focus on visit processes, assumptions, and study design.  Statistical Methods in Medical Research, 25, 2992–3014.
  * R Development Core Team (2014) R Development Core Team (2014) R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Rotnitzky et al. (1998) Rotnitzky, A., Robins, J. and Scharfstein, D. (1998) Semiparametric regression for repeated outcomes with nonignorable nonresponse.  Journal of the American Statistical Association, 93, 1321–1339.  URL: https://www.tandfonline.com/doi/abs/10.1080/01621459.1998.10473795.
  * Scharfstein et al. (1999) Scharfstein, D., Rotnitzky, A. and Robins, J. (1999) Adjusting for nonignorable drop-out using semiparametric nonresponse models.  Journal of the American Statistical Association, 94, 1096–1120.  URL: https://www.tandfonline.com/doi/abs/10.1080/01621459.1999.10473862.
  * Smith et al. (2022) Smith, B., Yang, S., Apter, A. J. and Scharfstein, D. O. (2022) Semi-parametric sensitivity analysis for trials with irregular and informative assessment times.  https://arxiv.org/abs/2204.11979.
  * Sun et al. (2011) Sun, L., Song, X. and Zhou, J. (2011) Regression analysis of longitudinal data with time-dependent covariates in the presence of informative observation and censoring times.  Journal of Statistical Planning and Inference, 141, 2902–2919.
  * Tan (2020) Tan, Z. (2020) Regularized calibrated estimation of propensity scores with model misspecification and high-dimensional data.  Biometrika, 107, 137–158.
  * Therneau (2023) Therneau, T. M. (2023) A Package for Survival Analysis in R.  URL: https://CRAN.R-project.org/package=survival.  R package version 3.5-5.
  * Vansteelandt et al. (2007) Vansteelandt, S., Rotnitzky, A. and Robins, J. (2007) Estimation of regression models for the mean of repeated outcomes under nonignorable nonmonotone nonresponse.  Biometrika, 94, 841–860.  URL: http://www.jstor.org/stable/20441421.
  * Wang (2020) Wang, Z. (2020) Global sensitivity analysis for randomized trials with informative assessment times: A fully parametric approach.  URL: http://jhir.library.jhu.edu/handle/1774.2/62670.  Master’s thesis, Johns Hopkins University.
  * Wen and Seaman (2018) Wen, L. and Seaman, S. (2018) Semi-parametric methods of handling missing data in mortal cohorts under non-ignorable missingness.  Biometrics, 74, 1427–1437.  URL: https://onlinelibrary.wiley.com/doi/abs/10.1111/biom.12891.
  * Yiu and Su (2022) Yiu, S. and Su, L. (2022) Joint calibrated estimation of inverse probability of treatment and censoring weights for marginal structural models.  Biometrics, 78, 115–127.
  * Zhu et al. (2017) Zhu, Y., Lawless, J. F. and Cotton, C. A. (2017) Estimation of parametric failure time distributions based on interval‐censored data with irregular dependent follow‐up.  Statistics in Medicine, 36, 1548–1567.
  * Zubizarreta (2015) Zubizarreta, J. R. (2015) Stable weights that balance covariates for estimation with incomplete outcome data.  Journal of the American Statistical Association, 110, 910–22.

Figure 1: Empirical bias and root mean squared error (RMSE) of the IIWEs with the MLE weights and the balancing weights estimators for $\beta_{1}$ (group effect) and $\beta_{2}$ (time effect) in the marginal regression model of a longitudinal continuous outcome, when the sample size $n=200$, the sensitivity parameter $\phi$ is set at $0.15,0.30,0.45,0.60$ and the visit process is _highly_ dependant on time-varying covariates with $\gamma_{z}=1.25$. The naïve analysis without weighting and the analysis based on complete data are also presented (note that the results from the complete data analyses are not visible due to the small values close to zeros).  Figure 2: Empirical bias and root mean squared error (RMSE) of the IIWEs with the MLE weights and the balancing weights estimators for $\beta_{1}$ (group effect) and $\beta_{2}$ (time effect) in the marginal regression model of a longitudinal continuous outcome, when the sample size $n=200$, the sensitivity parameter $\phi$ is set at $0.15,0.30,0.45,0.60$ and the visit process is _moderately_ dependant on time-varying covariates with $\gamma_{z}=0.5$. The naïve analysis without weighting and the analysis based on complete data are also presented (note that the results from the complete data analyses are not visible due to the small values close to zeros).  Table 1: Data generating mechanism for the simulations Possible visit times: | $t=0.01,0.02,\ldots,5$
---|---
Baseline covariates: | $X\stackrel{{\scriptstyle i.i.d}}{{\sim}}Bernoulli(0.5)$
Time varying covariates: | $Z_{1}(t),Z_{2}(t)\stackrel{{\scriptstyle i.i.d}}{{\sim}}N(-X,1)$
Continuous outcome: |  $Y(t)=5+Z_{1}(t)+Z_{2}(t)-0.5Z_{1}(t)Z_{2}(t)-2X-0.5t+\epsilon(t)$, $\epsilon(t)\stackrel{{\scriptstyle i.i.d}}{{\sim}}N(0,0.25)$
Count outcome: |  $Y(t)\stackrel{{\scriptstyle i.i.d}}{{\sim}}\mbox{NegBin}\\{\mu(t),\theta\\}$ with p.m.f. $\displaystyle\frac{\bm{\Gamma}(y+1/\theta)}{\bm{\Gamma}(1/\theta)y!}\left\\{\frac{\theta\mu(t)}{1+\theta\mu(t)}\right\\}\left\\{\frac{1}{1+\theta\mu(t)}\right\\}^{1/\theta},$
|  mean
$\mu(t)=\exp\\{1.69+Z_{1}(t)+Z_{2}(t)-0.5Z_{1}(t)Z_{2}(t)+0.67X-0.5t\\}$,
|  over-dispersion parameter $\theta=0.5$
Visiting process: |  $dN(t)\sim\text{Bernoulli}\\{\pi(t)\\}$,
|
$\pi(t)=\min\left[1,\exp\left\\{-3.05-2t+\gamma_{z}Z_{1}(t)+\gamma_{z}Z_{2}(t)+0.5Z_{1}(t)Z_{2}(t)+X+\phi
Y(t)^{*}\right\\}\right]$,
|  $Y(t)^{*}=Y(t)$ if $Y(t)$ is continuous,
|  $Y(t)^{*}=\log\\{Y(t)+1\\}$ if $Y(t)$ is the count outcome.
|  Sensitivity parameter: $\phi=0.15,0.3,0.45,0.6$
|  Dependence on covariates: $\gamma_{z}=0.5,1.25$
Transformed covariates: | $Z^{*}_{1}(t)=Z_{1}(t)-Z_{2}(t)$
|  $Z^{*}_{2}(t)=Z_{2}(t)+e(t)$, $e(t)$ $\stackrel{{\scriptstyle
i.i.d}}{{\sim}}$ $N(0,0.01)$
Table 2: Estimated hazard ratios with 95% confidence intervals and p-values from the fitted Cox model for the visit process in the PsA clinic cohort, assuming the visiting at random assumption was satisfied (i.e., $\phi=0$).  | Hazard ratio | 95 % CI | p-value
---|---|---|---
Gender (male) | 1.02 | 0.96 | 1.09 | 0.520
PsA disease duration at baseline | 1.07 | 1.05 | 1.09 |  $<$0.001
Damaged joints at baseline | 1.06 | 0.99 | 1.13 | 0.104
Active joints at baseline | 0.98 | 0.95 | 1.02 | 0.422
ESR at baseline | 1.00 | 0.98 | 1.02 | 0.870
1-2 years (vs. 0-1 year) since baseline | 2.62 | 2.19 | 3.13 |  $<$0.001
2-3 years (vs. 0-1 year) since baseline | 2.68 | 2.22 | 3.24 |  $<$0.001
3-4 years (vs. 0-1 year) since baseline | 3.14 | 2.57 | 3.83 |  $<$0.001
4-5 years (vs. 0-1 year) since baseline | 3.36 | 2.71 | 4.17 |  $<$0.001
$\geq 5$ years (vs. 0-1 year) since baseline | 4.63 | 3.76 | 5.70 |  $<$0.001
0-1 years since baseline $*$ biologics use | 1.09 | 0.90 | 1.33 | 0.355
1-2 years since baseline $*$ biologics use | 1.21 | 1.02 | 1.45 | 0.033
2-3 years since baseline $*$ biologics use | 1.24 | 1.02 | 1.49 | 0.028
3-4 years since baseline $*$ biologics use | 1.23 | 1.01 | 1.50 | 0.038
4-5 years since baseline $*$ biologics use | 1.35 | 1.09 | 1.67 | 0.007
$\geq 5$ years since baseline $*$ biologics use | 1.13 | 1.00 | 1.27 | 0.046
PsA disease duration at $t$ | 0.94 | 0.92 | 0.95 |  $<$0.001
Age at $t$ | 1.00 | 1.00 | 1.00 | 0.619
Year $[1980,2008]$ $*$ latest ESR prior to $t$ | 1.00 | 0.97 | 1.04 | 0.940
Year $[2008,2012)$ $*$ latest ESR prior to $t$ | 1.00 | 0.97 | 1.03 | 0.868
Year $[2012,2015)$ $*$ latest ESR prior to $t$ | 0.99 | 0.95 | 1.03 | 0.700
Year $[1980,2008]$ $*$ latest active joint count prior to $t$ | 1.05 | 1.00 | 1.11 | 0.054
Year $[2008,2012)$ $*$ latest active joint count prior to $t$ | 1.06 | 1.01 | 1.12 | 0.023
Year $[2012,2015)$ $*$ latest active joint count prior to $t$ | 1.07 | 0.99 | 1.14 | 0.075
Year $[1980,2008]$ $*$ latest damaged joint count prior to $t$ | 0.96 | 0.89 | 1.03 | 0.281
Year $[2008,2012)$ $*$ latest damaged joint count prior to $t$ | 0.94 | 0.87 | 1.01 | 0.077
Year $[2012,2015)$ $*$ latest damaged joint count prior to $t$ | 0.96 | 0.89 | 1.03 | 0.271
Year $[1980,2008]$ $*$ latest NSAID use prior to $t$ | 1.11 | 0.97 | 1.27 | 0.114
Year $[2008,2012)$ $*$ latest NSAID use prior to $t$ | 0.99 | 0.89 | 1.10 | 0.874
Year $[2012,2015)$ $*$ latest NSAID use prior to $t$ | 0.99 | 0.87 | 1.12 | 0.839
Year $[1980,2008]$ $*$ latest DMARD use prior to $t$ | 1.05 | 0.92 | 1.20 | 0.455
Year $[2008,2012)$ $*$ latest DMARD use prior to $t$ | 1.09 | 0.98 | 1.22 | 0.124
Year $[2012,2015)$ $*$ latest DMARD use prior to $t$ | 1.07 | 0.94 | 1.21 | 0.326
Year $[1980,2008]$ $*$ biologics use | 0.73 | 0.64 | 0.83 |  $<$0.001
Year $[2008,2012)$ $*$ biologics use | 0.80 | 0.71 | 0.91 |  $<$0.001
Year $[2012,2015)$ $*$ biologics use | 0.82 | 0.71 | 0.94 | 0.004
Year $[1980,2008]$ $*$ median inter-visit time by $t$ | 0.60 | 0.46 | 0.79 |  $<$0.001
Year $[2008,2012)$ $*$ median inter-visit time by $t$ | 0.40 | 0.26 | 0.61 |  $<$0.001
Year $[2012,2015)$ $*$ median inter-visit time by $t$ | 0.42 | 0.23 | 0.75 | 0.004
Table 3: Estimates (EST) and jackknife standard errors (SE) of the conditional associations between biologics use and active joint counts by years, $\beta_{1},\ldots,\beta_{6}$, based on the naïve estimator and the IIWEs with the MLE weights and the balancing weights estimators under different values of the sensitivity parameter ($\phi=0,0.5,1$) in the PsA clinic data analysis.  |  | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$ | $\beta_{4}$ | $\beta_{5}$ | $\beta_{6}$
---|---|---|---|---|---|---|---
|  | EST | SE | EST | SE | EST | SE | EST | SE | EST | SE | EST | SE
Naïve estimator | 0.12 | 0.19 | 0.08 | 0.22 | -0.11 | 0.28 | -0.27 | 0.25 | -0.38 | 0.29 | 0.02 | 0.27
IIWE, $\phi=0$ |  |  |  |  |  |  |  |  |  |  |  |
MLE wts. | 0.13 | 0.22 | 0.14 | 0.25 | -0.09 | 0.31 | -0.23 | 0.26 | -0.48 | 0.29 | 0.01 | 0.30
Balancing wts. | 0.13 | 0.21 | 0.12 | 0.24 | -0.12 | 0.31 | -0.23 | 0.26 | -0.48 | 0.29 | 0.02 | 0.29
IIWE, $\phi=0.5$ |  |  |  |  |  |  |  |  |  |  |  |
MLE wts. | -0.10 | 0.18 | -0.02 | 0.22 | -0.18 | 0.26 | -0.41 | 0.26 | -0.54 | 0.28 | -0.10 | 0.28
Balancing wts. | -0.10 | 0.18 | -0.02 | 0.22 | -0.18 | 0.26 | -0.42 | 0.26 | -0.53 | 0.28 | -0.08 | 0.29
IIWE, $\phi=1$ |  |  |  |  |  |  |  |  |  |  |  |
MLE wts. | -0.36 | 0.18 | -0.18 | 0.20 | -0.30 | 0.22 | -0.48 | 0.25 | -0.58 | 0.27 | -0.22 | 0.26
Balancing wts. | -0.35 | 0.18 | -0.16 | 0.21 | -0.31 | 0.23 | -0.51 | 0.26 | -0.56 | 0.28 | -0.18 | 0.27