Analyzing zero-inflated clustered longitudinal ordinal outcomes using GEE-type models with an application to dental fluorosis studies

Generalized Estimating Equations (GEEs) proposed by Liang and Zeger [1986] extend the Generalized Linear Model (GLM) framework to analyze longitudinal and clustered data. GEEs are designed to handle correlated observations; this makes them particularly useful in situations where repeated measures are collected on the same subject. Observations in these subject-specific “clusters” are not necessarily independent. GEEs use quasi-likelihood methods, focusing on the mean structure and correlation rather than fully specifying the joint distribution of the data. Parameter estimates in are obtained iteratively until convergence. In a setting with repeated measures on $N$ subjects, a GEE is of the following form:

where $\bm{\beta}$ are parameters, $\mathbf{Y_{i}}$ is the vector of responses for subject $i$, $\bm{\mu_{i}}=E(\mathbf{Y_{i}})=\mathbf{X_{i}^{\prime}\beta}$ is the mean response, $\mathbf{X_{i}}$ denotes covariates for subject $i$, $\mathbf{D_{i}}=\frac{\partial\bm{\mu_{i}}}{\partial\bm{\beta}}$ is the matrix of derivatives, $\mathbf{V_{i}}$ is the working covariance matrix. A working correlation structure $\mathbf{V_{i}}$ is specified for repeated measurements within cluster $i$. GEEs provide consistent parameter estimates despite the misspecification of $\mathbf{V_{i}}$. This robustness is valuable when the true correlation structure is unknown or complex.

A hurdle model is a two-part model used to handle datasets with a substantial number of zero observations. It consists of a zero model and a conditional non-zero outcome model. Our modeling approach in the subsequent sections is inspired by the hurdle model, where all “zeroes” are viewed as originating from a single source.

In the context of our study, the zero model, which we will refer to as our “presence” model, determines whether fluorosis is present or absent. The “non-zero” model conditioning on a non-zero FRI score outcome (indicating the presence of some level of observable fluorosis) is known as the “severity” model. This model specifically quantifies the degree of fluorosis in affected teeth or zones.

In Section (2.3), we will introduce presence models that address all zero/non-zero FRI scores, and severity models that handles the non-zero outcomes. A combined modeling scheme that connects both of these components is then proposed.

The presence model utilizes $W_{P,itjk}=I(Y_{itjk}\neq 0)$, the indicator variable for the presence of fluorosis. In general, it can be expressed through McCullagh’s model [McCullagh, 1980]:

We will illustrate the case where $F(.)$ is the logistic CDF, that is, $F(t)=\frac{\exp(t)}{1+\exp(t)}$. This specification corresponds to the case of logistic regression. Accordingly, we have

Corresponding to each time point $t=1,\cdots,T$, parameter estimates for $\alpha_{P,t}$ and $\bm{\beta_{P,t}}$ are obtained from the following GEE; by specifying cluster correlation matrices, this approach accounts for correlations in the observations, an improvement over ordinary logistic regression.

Here, $\bm{W_{P,i_{p},t}}$ denotes the vector of $W_{P,tijk}$ for fixed $t$ and $i$, stacked over the levels $j=1,\cdots,J$ and $k=1,\cdots,K$. $\bm{\mu_{P,i_{p},t}=E(W_{P,i_{p},t}})$, which is also a vector of $\mu_{P,tijk}$ for fixed $t$ and $i$, stacked over the levels of $j=1,\cdots,J$ and $k=1,\cdots,K$. $\bm{V_{p,i_{p},t}=A^{1/2}_{p,i_{p},t}R_{p,i_{p},t}A^{1/2}_{p,i_{p},t}}$, where $\bm{R_{p,i_{p},t}}$ denotes the working cluster correlation matrix of cluster $i$ at time of dental visit $t$ and $\bm{A_{p,i_{p},t}}$ denotes $Var(\bm{W_{P,i_{p},t}})$, the variance-covariance matrix of $\bm{W_{P,i_{p},t}}$. The entries of $\bm{R_{p,i_{p},t}}$ are specific to the chosen working correlation structure, and are detailed in Section 2.4. Additional details on the GEE components can be found in the Appendix (Section 6).

The severity model, modelling $W_{S,itjk}$, uses a proportional odds framework to model the conditional probability of different severity levels. The model is specified as

so

Estimates for $\alpha_{S,t,1|2},\cdots,\alpha_{S,t,L-1|L},$ and $\bm{\beta_{S,t}}$ are obtained from the following GEEs:

where $\bm{Z_{S,it}=(Z_{S,t,111},\cdots,Z_{S,i,tjk},\cdots,Z_{S,t,NJK}})^{T}$, with $\bm{Z_{S,i,tjk}}=(Z_{S,i,tjk,1},Z_{S,i,tjk,2},Z_{S,i,tjk,3})^{T}$, where $Z_{S,i,tjk,l}=I[W_{S,itjk}=l]$, $l=1,\cdots,L$. $\bm{\pi_{S,it}}=E[\bm{Z_{S,it}}]=(\pi_{S,t,111},\cdots,\pi_{S,i,tjk},\cdots,\pi_{S,t,NJK})^{T}$ with $\bm{\pi_{S,i,tjk}}=(\pi_{S,i,tjk,1},\pi_{S,i,tjk,2},\pi_{S,i,tjk,3})^{T}$, where $\pi_{S,i,tjk,l}=P[W_{S,itjk}=l]=\mu_{S,i,tjk,l+1}-\mu_{S,i,tjk,l}$, $l=1,\cdots,L-1$.

$\bm{V_{S,it}=A^{1/2}_{S,it}R_{S,it}A^{1/2}_{S,it}}$, where $\bm{R_{S,it}}$ denotes the working cluster correlation matrix of cluster $i$ at time of dental visit $t$ and $\bm{A_{S,i_{p},t}}$ denotes $Var(\bm{Z_{S,it}})$, the variance-covariance matrix of $\bm{Z_{S,it}}$. The entries of $\bm{R_{S,it}}$, specific to the chosen working correlation structure, are detailed in Section 2.4. Additional details on the GEE components can be found in the Appendix (Section 6).

We outline a combined modelling approach that integrates the presence and severity components by assuming a relationship $\mathbf{\beta_{S,t}}=\gamma_{t}\mathbf{\beta_{P,t}}$ between their regression coefficients. The parameter $\mathbf{\beta_{P,t}}$ signifies the effect that predictors exert on dental fluorosis in the presence piece. The scalar parameter $\gamma_{t}$ serves as an amplifying factor in the severity piece, specifically augmenting the effect of predictors within the model’s severity dimension, as $\mathbf{\beta_{S,t}}=\gamma_{t}\mathbf{\beta_{P,t}}$. Specifically, we model

In a two-step approach, the parameter $\gamma_{t}$ is first estimated from separate estimates of presence and severity parameters. As the intercepts $\alpha_{P,t},\alpha_{S,t,1|2},\cdots,\alpha_{S,t,L|L-1}$ are ancillary for the slope parameters, they are estimated from separate model pieces (presence and severity). Next, the shared parameters $\beta_{P,t,1},\cdots,\beta_{P,t,q}$ are updated from combined estimating equations derived from averaging the GEEs of both models.

The intercepts $\alpha_{P,t},\alpha_{S,t,1|2},\cdots,\alpha_{S,t,L|L-1}$ and $\gamma_{t}$ estimated from separate model pieces. Specifically, ${\gamma}_{t}$ is estimated as the ratio of the sum $\sum_{g=1}^{q}\hat{\beta}_{S,t,g}$ of the coefficients of the severity model to sum $\sum_{g=1}^{q}\hat{\beta}_{P,t,g}$ of the coefficients of the presence model, that is, ${\hat{\gamma}}_{t}=\frac{\sum_{g=1}^{q}\hat{\beta}_{S,t,g}}{\sum_{g=1}^{q}\hat{\beta}_{P,t,g}}$, $g=1,\cdots,q$. Then, only the parameters in $\bm{\beta_{P,t}}$ remain to be estimated. As $\mathbf{\beta_{S,t}}$ depends on $\gamma_{t}$ and $\mathbf{\beta_{P,t}}$ and is estimated as $\mathbf{\hat{\beta}_{S,t}}=\hat{\gamma}_{t}\mathbf{\hat{\beta}_{P,t}}$ once $\gamma_{t}$ and $\mathbf{\beta_{P,t}}$ are estimated.

The parameters in $\bm{\beta_{P,t}}$ is derived by estimating the following combined estimating equation:

where $\bm{\beta_{S,t}}=\gamma_{t}\bm{\beta_{P,t}}$. Applying the modified Newton-Raphson method, we update $\bm{\beta_{P,t}}$ for $t=1,\cdots,T$ iteratively until convergence:

The correlation parameters are estimated according to the following schemes proposed by Liang and Zeger [1986].

Let $r_{P,tijk}=\frac{w_{P,tijk}-\mu_{P,tijk}}{\sqrt{Var(W_{P,tijk})}}$ and $r_{S,tijk,l}=\frac{Z_{S,tijk,l}-\pi_{S,tijk,l}}{\sqrt{Var(Z_{S,tijk,l})}}$.

We define

and

. Then, for the exchangeable cluster correlation structure,

$\hat{\rho}_{P,t,exch}=\frac{1}{\sum_{i=1}^{N}\sum_{j,j^{\prime}=1}^{J_{P,ti}}\sum_{\begin{subarray}{c}k,k^{\prime}=1\\ j\neq j,k\neq k\end{subarray}}^{K_{P,ti}}1}{\phi}_{P,t}\sum_{i=1}^{N}\sum_{\begin{subarray}{c}j,j^{\prime}=1\\ j\neq j^{\prime}\end{subarray}}^{J_{P,ti}}\sum_{\begin{subarray}{c}k,k^{\prime}=1\\ k\neq k\end{subarray}}^{K_{P,ti}}\hat{r}_{P,tijk}\hat{r}_{P,tij^{\prime}k^{\prime}}$

$\hat{\rho}_{S,t,exch}=\frac{1}{\sum_{i=1}^{N}\sum_{l,l^{\prime}=1}^{L}\sum_{j,j^{\prime}=1}^{J_{S,ti}}\sum_{\begin{subarray}{c}k,k^{\prime}=1\\ l\neq l^{\prime},j\neq j^{\prime},k\neq k\end{subarray}}^{K_{S,ti}}1}{\phi}_{S,t}\sum_{i=1}^{N}\sum_{l,l^{\prime}=1}^{L}\sum_{j,j^{\prime}=1}^{J_{S,ti}}\sum_{\begin{subarray}{c}k,k^{\prime}=1\\ l\neq l^{\prime},j\neq j^{\prime},k\neq k\end{subarray}}^{K_{S,ti}}\hat{r}_{S,tijk,l}\hat{r}_{S,tij^{\prime}k^{\prime},l^{\prime}}$

For the AR(1) cluster correlation structure, $\rho_{P,t,AR1}$ is estimated from the slope of the regression of $\log(\hat{r}_{P,tijk}\hat{r}_{P,tij^{\prime}k^{\prime}})$ on $|j-j^{\prime}|$. Likewise, $\rho_{S,t,AR1}$ is estimated from the slope of the regression of $\log(\hat{r}_{S,tijk,l}\hat{r}_{S,tij^{\prime}k^{\prime},l})$ on $|j-j^{\prime}|$.

For the jackknifed variance-covariance structure, we define
$\hat{\rho}_{P,t,jk,j^{\prime}k^{\prime},jackknife(-i)}=\frac{1}{\sum_{\begin{subarray}{c}i^{*},i^{**}=1\\ i^{*},i^{**}\neq i\end{subarray}}^{N}\sum_{j,j^{\prime}=1}^{J_{P,ti}}\sum_{k,k^{\prime}=1}^{K_{P,ti}}1}{\phi}_{P,t,-i}\sum_{\begin{subarray}{c}i^{*},i^{**}=1\\ i^{*},i^{**}\neq i\end{subarray}}^{N}\sum_{j,j^{\prime}=1}^{J_{P,ti}}\sum_{k,k^{\prime}=1}^{K_{P,ti}}\hat{r}_{P,t{i^{*}}jk}\hat{r}_{P,t{i^{**}}j^{\prime}k^{\prime}}$, where $\phi_{P,t,-i}$ corresponds to the quantity in 13 computed leaving out cluster $i$.

Then $\hat{\rho}_{P,t,jk,j^{\prime}k^{\prime},jackknife}=\frac{1}{N}\sum_{i=1}^{N}\hat{\rho}_{P,t,jk,j^{\prime}k^{\prime},jackknife(-i)}$. Likewise, $\hat{\rho}_{S,t,jk,j^{\prime}k^{\prime},jackknife}=\frac{1}{N}\sum_{i=1}^{N}\hat{\rho}_{S,t,jk,j^{\prime}k^{\prime},jackknife(-i)}$, where

with $\phi_{S,t,-i}$ corresponding to the quantity in 14 computed leaving out cluster $i$.