 Research
 Open access
 Published:
A biasreduced generalized estimating equation approach for proportional odds models with smallsample longitudinal ordinal data
BMC Medical Research Methodology volume 24, Article number: 140 (2024)
Abstract
Background
Longitudinal ordinal data are commonly analyzed using a marginal proportional odds model for relating ordinal outcomes to covariates in the biomedical and health sciences. The generalized estimating equation (GEE) consistently estimates the regression parameters of marginal models even if the working covariance structure is misspecified. For smallsample longitudinal binary data, recent studies have shown that the bias of regression parameters may result from the GEE and have addressed the issue by applying Firth’s adjustment for the likelihood score equation to the GEE as if generalized estimating functions were likelihood score functions. In this manuscript, for the proportional odds model for longitudinal ordinal data, the smallsample properties of the GEE were investigated, and a biasreduced GEE (BRGEE) was derived.
Methods
By applying the adjusted function originally derived for the likelihood score function of the proportional odds model to the GEE, we produced the BRGEE. We investigated the smallsample properties of both GEE and BRGEE through simulation and applied them to a clinical study dataset.
Results
In simulation studies, the BRGEE had a bias closer to zero, smaller root mean square error than the GEE with coverage probability of confidence interval near or above the nominal level. The simulation also showed that BRGEE maintained a type I error rate near or below the nominal level.
Conclusions
For the analysis of longitudinal ordinal data involving a small number of subjects, the BRGEE is advantageous for obtaining estimates of the regression parameters of marginal proportional odds models.
Introduction
Longitudinal ordinal data are frequently collected in biomedical and health science studies. In such a study, each subject is observed repeatedly over a period, and ordinal responses of interest at each observation are recorded. Because repeated responses from the same subject are usually correlated, a straightforward application of generalized linear models for a single response variable to longitudinal data is not appropriate. There are numerous approaches for extending generalized linear models to longitudinal data. One approach for extending generalized linear models to longitudinal data is a class of regression models where the model for the mean response at each observation does not incorporate dependence on any random effects or previous responses. The model is known as a marginal model. An alternative approach for accounting for the withinsubject association is via the introduction of random effects. The model is known as the generalized linear mixed effects model. Due to the interpretation of their regression coefficients, the former are often referred to as “populationaverage models,” and the latter are referred to as “subjectspecific models” [1]. Assuming a situation where the target of inference is the populationlevel summary described as a basis for comparison between treatment conditions in ICH E9 (R1), “Addendum on Estimands and Sensitivity Analysis in Clinical Trials” to the guideline on Statistical Principles for Clinical Trials [2], we focus on marginal models here. For estimation of the marginal model parameters, assumptions about the joint distribution of the responses are not necessary. The avoidance of distributional assumptions leads to a method of estimation known as the generalized estimating equation (GEE) [3]. An appealing property of the GEE estimator is that it is a consistent estimator even if the assumed model for the covariance among the repeated measures is not correctly specified. In addition, we can obtain valid variance estimate for regression coefficients by utilizing the empirical or sandwich estimator. The sandwich estimator possesses a notable property as it demonstrates robustness by providing valid variance estimate even when the working covariance among the repeated measures is not correct. The GEE requires only that the model for the mean response, the link function and the variance of each response be correct.
For binary response data, recent studies have investigated the smallsample property of GEE estimates of the regression parameter and the biasreduced estimates obtained by applying Firth’s adjustment for the likelihood score equation [4] to the GEE as if generalized estimating functions were likelihood score functions [5,6,7,8]. Paul and Zhang [5] found that GEE estimates yielded biased estimates of the regression parameters when the sample size was small and that biasreduced estimates improved bias and mean squared error (MSE) by simulation studies. Mondol and Rahman [6] showed that biasreduced GEE (also referred to as penalized GEE) achieved convergence and provided finite estimates in the presence of separation. Geroldinger et al. [7] reported that the biasreduced estimation improved convergence compared to the standard GEE with a similar or even better performance in terms of the accuracy of estimates. Gosho et al. [8] reviewed biasreduced GEEs and modified covariance estimators and evaluated their performance in sparse binary data from smallsample longitudinal studies. In a study investigating marginal models in small samples, Bie et al. [9] reported that both GEE and marginalized multilevel models exhibit smallsample bias when correct correlation structure is adopted. Greenland [10] states that “the potential for smallsample bias in results from asymptotic procedures needs to be checked more routinely than is current practice.”
The proportional odds model is the commonly used model for relating ordinal outcome to covariates. The model is also referred to as the proportional odds version of the cumulative logit model [11]. The proportional odds model assumes that the association of each covariate with the outcome is represented by a single odds ratio. The purpose of this paper is to obtain the biasreduced GEE (BRGEE) estimates of the proportional odds model for longitudinal ordinal data and to investigate the smallsample properties of both the GEE and the BRGEE through simulation. Kosmidis and Firth [12] derived general expressions for the adjusted score equations for the general class of multivariate models. Using the adjusted score equations exploited in Kosmidis and Firth [12] for the multivariate generalized linear model, Kosmidis [13] derived the adjusted score equation for the proportional odds model. We apply the adjusted likelihood score equation for the proportional odds model to the GEE to obtain the BRGEE estimate as a solution to an adjusted estimating equation, which is induced by adding the adjustment function.
“Methods” section provides a brief summary of the BRGEE estimates for longitudinal ordinal data. In “Simulation study” section, the results of the simulation study are presented and the smallsample properties of the GEE and the BRGEE are investigated. In “Example” section, we apply the BRGEE to the data of postoperative pain after laparoscopic cholecystectomy illustrated in Lumley [14]. A discussion of our findings follows in “Discussion” section.
Methods
Notation and GEE for the proportional odds model
Suppose a study includes \(N\) subjects with \({n}_{i}\) repeated observations of the \(K\) multinomial ordered categories for the ith subject, \(i=1,\dots ,N\). For simplicity, we assume equal repetition, \({n}_{i}=n\). Let \({{\varvec{Y}}}_{it}\) represent responses as a vector \({{\varvec{Y}}}_{it}={\left({Y}_{it1},\dots ,{Y}_{itq}\right)}^{T}\) of \(q=K1\) dummy variables, where \({Y}_{itk}=1\) for the observation at occasion \(t (t=1,\dots ,n)\) of subject \(i (i=1,\dots ,N)\) falling in category \(k (k=1,\dots ,q)\) and \({Y}_{itk}=0\) otherwise. The vector of marginal means or response probabilities of \({{\varvec{Y}}}_{it}\) is \({{\varvec{\pi}}}_{it}={\left({\pi }_{it1},\dots ,{\pi }_{itq}\right)}^{T}\). The proportional odds model links the cumulative probability \({\gamma }_{itk}={\pi }_{it1}+\dots +{\pi }_{itk}\) to a \(p\)vector of covariates \({\mathbf{x}}_{it}\) via the following relationship:
where \({\varvec{\delta}}={\left({\beta }_{01},\dots ,{\beta }_{0q},{{\varvec{\beta}}}^{T}\right)}^{T}\) is a \(\left(q+p\right)\)vector of model parameters with \({\beta }_{01}<\dots <{\beta }_{0q}\) and \({\varvec{\beta}}={\left({\beta }_{1},\dots ,{\beta }_{p}\right)}^{T}\). The linear term \({\eta }_{itk}={\beta }_{0k}+{\text{x}}_{it}^{T}{\varvec{\beta}}\) is determined by \({\eta }_{itk}={\sum }_{r=1}^{p+q}{\delta }_{r}{z}_{itkr}\), where \({z}_{itkr}\) is the \(\left(k,r\right)\) th element of the \(q\times \left(q+p\right)\) matrix,
Let \({{\varvec{Y}}}_{i}={\left({{\varvec{Y}}}_{i1}^{T},\dots ,{{\varvec{Y}}}_{in}^{T}\right)}^{T}\) and \({{\varvec{\pi}}}_{i}={\left({{\varvec{\pi}}}_{i1}^{T},\dots ,{{\varvec{\pi}}}_{in}^{T}\right)}^{T}\) denote responses and response probabilities, respectively, for subject \(i\). Then, a GEE for consistent estimation of \({\varvec{\delta}}\) is given in the form,
where \({{\varvec{D}}}_{i}=D\left({{\varvec{\pi}}}_{i};{\varvec{\delta}}\right)\) is the \(qn\times p\) Jacobian of \({{\varvec{\pi}}}_{i}\) with respect to \({\varvec{\updelta}}\) and \({{\varvec{V}}}_{i}={{\varvec{V}}}_{i}({\varvec{\delta}},\alpha )\) is a ‘working’ covariance matrix indexed by \({\varvec{\updelta}}\) and an association parameter \(\alpha\). Further details on the GEE for ordinal data can be found in Agresti [11], Heagerty and Zeger [15], Fahrmeir and Pritscher [16], Lipsitz et al. [17], Liang et al. [18], and Williamson et al. [19].
Let \({{\varvec{V}}}_{itt}={\text{Var}}\left({{\varvec{Y}}}_{it}\right)=\text{diag}\left({{\varvec{\pi}}}_{it}\right){{\varvec{\pi}}}_{it}{{\varvec{\pi}}}_{it}^{T}\) denote the \(q\times q\) multinomial covariance matrix for \({{\varvec{Y}}}_{it}\). For \(t\ne {t}{\prime}\), \({{\varvec{V}}}_{it{t}{\prime}}={\text{Cov}}\left({{\varvec{Y}}}_{it},{{\varvec{Y}}}_{i{t}{\prime}}\right)\) is the covariance matrix between two different occasions in the same subject. The simplest model for the working covariance matrix is the independent model, \({{\varvec{V}}}_{it{t}{\prime}}=\bf{0}\). Lumley [14] described working covariance structures \({{\varvec{V}}}_{it{t}{\prime}}\) and computational methods for ordinal data based on the global odds ratio that allow the GEE to be used for smaller sets of ordinal data and with less effort expended on modeling associations. For each pair of categories \(c\) and \({c}{\prime}\) of two ordinal responses \({R}_{it}\) and \({R}_{i{t}{\prime}}\) in the same subject, the global odds ratio is defined as
The covariance of \({Y}_{itc}\) and \({Y}_{i{t}{\prime}{c}{\prime}}\) is given by \(\text{cov}\left({Y}_{itc}, {Y}_{i{t}{\prime}{c}{\prime}}\right)=E\left({Y}_{itc}{Y}_{i{t}{\prime}{c}{\prime}}\right)E\left({Y}_{itc}\right)E\left({Y}_{i{t}{\prime}{c}{\prime}}\right)\), where the second term is the product of \({\pi }_{itc}\) and \({\pi }_{i{t}{\prime}{c}{\prime}}\). The first term \(E\left({Y}_{itc}{Y}_{i{t}{\prime}{c}{\prime}}\right)\) can be expressed by the joint cumulative probabilities \({\gamma }_{itc{t}{\prime}{c}{\prime}}=\text{Pr}\left({R}_{it}\le c, {R}_{i{t}{\prime}}\le {c}{\prime}\right)\):
The joint cumulative probabilities \({\gamma }_{itc{t}{\prime}{c}{\prime}}\) can be expressed in terms of global odds ratio and marginal cumulative probabilities:
where \(\kappa =1+({\gamma }_{itc}+{\gamma }_{i{t}{\prime}{c}{\prime}})({\Psi }_{itc{t}{\prime}{c}{\prime}}1)\). A simple estimate of the crude global odds ratio is calculated by the weighted mean of the log odds ratios with weights inversely proportional to the variances for every possible pair of the two time points \(\left(t,{t}{\prime}\right)\), where the odds ratios are computed based on the \(2\times 2\) tables obtained from the \({q}^{2}\) ways of collapsing the row and column classifications into dichotomies for a \(K\times K\) table of the variables \({R}_{it}\) and \({R}_{i{t}{\prime}}\) over all subjects. Hines [20] compared the various models for \({{\varvec{V}}}_{it{t}{\prime}}\) both in terms of efficiency and computational stability. She showed that the performance of Lumley’s crude global odds ratio method is satisfactory. Therefore, our approach is based on the working covariance structure of the independence or exchangeable model, where the crude global odds ratio is constant for all pairs of times.
Biasreduced GEE for the proportional odds model
Firth [4] showed that the solution of the following adjusted score equation results in an estimator that is free from the firstorder term in the asymptotic expansion of the bias in the MLE of the regression parameter:
with \({\varvec{A}}\left({\varvec{\theta}}\right)={\varvec{I}}({\varvec{\theta}}){\varvec{b}}({\varvec{\theta}})\), where \({\varvec{\theta}}={({\theta }_{1},\dots ,{\theta }_{p})}^{T}\) is a \(p\)vector of regression parameters from the generalized linear model, \({\varvec{U}}\left({\varvec{\theta}}\right)=\frac{\partial }{\partial{\varvec{\theta}}}l({\varvec{\theta}})=\bf{0}\) is the standard score equation based on the loglikelihood function \(l({\varvec{\theta}})\), \({\varvec{I}}({\varvec{\theta}})\) is the expected information matrix for \({\varvec{\theta}}\), and \({\varvec{b}}({\varvec{\theta}})\) is the first term in the asymptotic expansion of the bias of the MLE.
Kosmidis and Firth [12] derived general expressions for the adjusted score functions for the general class of multivariate models, which include multivariate generalized linear models. Using the adjustment functions exploited in Kosmidis and Firth [12] for the score functions of the multivariate generalized linear model, Kosmidis [13] formulated the adjustment function of the score function of the proportional odds model. Considering that the GEE is equivalent to the generalized linear model score function under the identity working covariance structure, we treat \({\varvec{U}}\left({\varvec{\delta}},\alpha \right)\) as if it were a likelihood score function and apply the adjustment function in the score function to the GEE. Then, the adjustment function for the generalized estimating function for the proportional odds model is given by
where \({{\varvec{Z}}}_{i}={\left({{\varvec{Z}}}_{i1}^{T},\dots ,{{\varvec{Z}}}_{in}^{T}\right)}^{T}\), \({{\varvec{\eta}}}_{i}={({{\varvec{\eta}}}_{i1}^{T},\dots ,{{\varvec{\eta}}}_{in}^{T})}^{T}\) with elements of \({{\varvec{\eta}}}_{it}={({\eta }_{it1},\dots ,{\eta }_{itq})}^{T}\), \(D\left({{\varvec{\pi}}}_{i};{{\varvec{\eta}}}_{i}\right)\) is the \(qn\times qn\) Jacobian of \({{\varvec{\pi}}}_{i}\) with respect to \({{\varvec{\eta}}}_{i}\), \({\varvec{F}}={\sum }_{i=1}^{N}{(D\left({{\varvec{\pi}}}_{i};{{\varvec{\eta}}}_{i}\right){{\varvec{Z}}}_{i})}^{T}{{\varvec{V}}}_{i}^{1}D\left({{\varvec{\pi}}}_{i};{{\varvec{\eta}}}_{i}\right){{\varvec{Z}}}_{i}\), \({D}^{2}\left({{\varvec{\pi}}}_{i};{{\varvec{\eta}}}_{i}\right)\) is the \({(qn)}^{2}\times qn\) matrix with the \(s\)th block as the Hessian of \({\pi }_{is}\) with respect to \({{\varvec{\eta}}}_{i} (s=1,\dots ,qn)\), \({{\varvec{I}}}_{qn}\) is the \(qn\times qn\) identity matrix, \({z}_{isr}\) is the \(\left(s,r\right)\) th element of \({{\varvec{Z}}}_{i}\), \({{\varvec{V}}}_{i}\) is the covariance matrix of the vector \({{\varvec{Y}}}_{i}\), and \({\left(D\left({{\varvec{\pi}}}_{i};{{\varvec{\eta}}}_{i}\right){{\varvec{V}}}_{i}^{1}\right)}_{s}\) is the \(s\)th row of \(D\left({{\varvec{\pi}}}_{i};{{\varvec{\eta}}}_{i}\right){{\varvec{V}}}_{i}^{1}\). The matrix \(D\left({{\varvec{\pi}}}_{i};{{\varvec{\eta}}}_{i}\right)\) is the block diagonal matrix whose \(t\) th diagonal element is
where \({g}_{itk}=g\left({\eta }_{itk}\right)\) with \(g\left(\eta \right)=dG(\eta )/d\eta\) and \(G\left(\eta \right)=\text{exp}(\eta )/(1+\text{exp}(\eta ))\); hence, \({g}_{itk}=G\left({\eta }_{itk}\right)(1G\left({\eta }_{itk}\right))={\gamma }_{ikt}\left(1{\gamma }_{ikt}\right)\).
The matrix \({D}^{2}\left({{\varvec{\pi}}}_{i};{{\varvec{\eta}}}_{i}\right)\) is the \({(qn)}^{2}\times qn\) matrix with the \(\left(s,u\right)\) th element given by:
where \(v=\left(t1\right)q+k \left(t=1,\dots ,n;k=1,\dots ,q\right)\) and \(w=\left(t1\right)q+k \left(t=1,\dots ,n;k=1,\dots ,q1\right)\).
The \(r\) th element of the biasreduced generalized estimating function is \({U}_{r}^{*}\left({\varvec{\delta}},\alpha \right)={U}_{r}\left({\varvec{\delta}},\alpha \right)+{A}_{r}\left({\varvec{\delta}},\alpha \right)\). The BRGEE estimate \(\widetilde{{\varvec{\delta}}}\) is such that \({U}_{r}^{*}\left(\widetilde{{\varvec{\delta}}},\widetilde{\alpha }\right)=0\) for every \(r=1,\dots ,q+p\), where \(\widetilde{\alpha }\) is the estimate of global odds ratio. \(\widetilde{{\varvec{\delta}}}\) can be estimated using an iterative method described below:
Step 1: Choose an initial value \({\widetilde{{\varvec{\delta}}}}^{0}\) of \({\varvec{\delta}}\).
Step 2: If the working covariance structure is exchangeable, we calculate the log (crude) global odds ratio \(\widetilde{\alpha }\) [14].
Step 3: Given \({\widetilde{{\varvec{\delta}}}}^{(t)}\) at the tth iteration and \(\widetilde{\alpha }\) (unnecessary if the working covariance structure is independent), the covariance matrix is estimated from the global odds ratio and the predicted cumulative probabilities [14].
Step 4: Given the working covariance matrix \({{\varvec{V}}}_{i}^{1}\left({\widetilde{{\varvec{\delta}}}}^{(t)},\widetilde{\alpha }\right)\), the current estimate \({\widetilde{{\varvec{\delta}}}}^{(t)}\) is updated according to the adjusted GEE using the Newton‒Raphson method given by
with
Step 5: Iterate steps 3 and 4 until a desired convergence criterion is satisfied (for example, \(\text{max}\left{\widetilde{{\varvec{\delta}}}}^{(t+1)}{\widetilde{{\varvec{\delta}}}}^{(t)}\right\le 0.0001\)). At convergence, the estimate of \({\varvec{\delta}}\) is denoted by \(\widetilde{{\varvec{\delta}}}\) and referred to as the BRGEE estimate.
The consistent estimate of the covariance matrix of the GEE regression parameter estimates is given by the sandwich estimator of Liang and Zeger [3],
where \({\widehat{{\varvec{r}}}}_{i}{\widehat{{\varvec{r}}}}_{i}^{T}=({{\varvec{Y}}}_{i}{\widehat{{\varvec{\pi}}}}_{i}){({{\varvec{Y}}}_{i}{\widehat{{\varvec{\pi}}}}_{i})}^{T}\) is used to estimate \({\text{Cov}}({{\varvec{Y}}}_{i})\).
In this manuscript, we calculate the covariance matrix for the BRGEE derived from the estimating function \({{\varvec{U}}}^{\boldsymbol{*}}\left({\varvec{\delta}},\alpha \right)={\varvec{U}}\left({\varvec{\delta}},\alpha \right)+{\varvec{A}}\left({\varvec{\delta}},\alpha \right)\),
with
Mancl and DeRouen [21] described that the biascorrected covariance estimator used with the critical value based on the Fdistribution instead of chisquare produced proper test sizes. Therefore, in the following, we used the critical value of the tdistribution with the number of the sample size minus the number of coefficients in the regression model as the degree of freedom to compute the confidence interval and the statistical significance.
When the sample size is small, the sandwich covariance matrix is expected to underestimate the covariance matrix. To reduce the smallsample bias, several modified variance estimators have been developed [22]. In this manuscript, in addition to \({{\varvec{\Sigma}}}_{{\varvec{G}}{\varvec{E}}{\varvec{E}}}\) and \({{\varvec{\Sigma}}}_{{\varvec{B}}{\varvec{R}}{\varvec{G}}{\varvec{E}}{\varvec{E}}}\), we used the biascorrected covariance estimates of Mancl and DeRouen [21] obtained by substituting the GEE and BRGEE estimates into the following equation in order to calculate the 95% confidence interval,
with
where \({{\varvec{I}}}_{qn}\) is the \(qn\times qn\) identity matrix and \({{\varvec{H}}}_{ii}\) is the matrix calculated by
Simulation study
In this section, we conducted a limited simulation study to investigate the smallsample properties of both the standard GEE (GEE) and the proposed biasreduced GEE (BRGEE) estimation of the regression parameters in a marginal model for correlated ordinal data. We assumed a randomized clinical trial with a parallel group design, where each subject was assigned to either the treatment group or the control group. Balanced correlated ordinal data with \(K=3\) of the total number of response categories and with observations over \(n=4\) time occasions were generated using the algorithm given in Ibrahim et al. [23]. They used the Goodman and Kruskal \(\Gamma\) coefficient as a measurement of the association for the responses of each subject. We modified the algorithm to specify the association of the withinsubject observation by using the global odds ratio. Specifically, we modified the macro developed by Ibrahim et al. [23] to calculate the joint probabilities of responses between two time points, which are used to generate correlated ordinal data, by substituting the global odds ratio and marginal cumulative probability based on the marginal model into expression (1). For the exchangeable correlation, we used a common \(\text{exp}\left(\alpha \right)\) as global odds ratio for all pairs of time points, while for ARtype correlation, we used \(\text{exp}\left(\alpha /\tau \right)\) as global odds ratio based on the time interval \(\tau\) for each pair of time points. We used a proportional odds marginal model both to simulate and to fit the data. We considered the following proportional odds model for the \(i\) th subject measured at the \(t\) th occasion (\(i=1,\dots ,N;t=1,\dots ,4\)):
where \({\gamma }_{itk}\) is the cumulative probability, \({\text{Trt}}_{i}\) is a treatment assignment variable coded as 1 for the first half of the subjects or 0 for the second half of the subjects, and \({\text{Time}}_{it}\) is an occasion coded as 1, 2, 3 and 4. We conducted simulations for four different sample sizes: \(N=20, 30, 40, 50\). We let the \((q+p)\)vector parameter set used with regression model (2) be \({\varvec{\delta}}={\left({\beta }_{01},{\beta }_{02},{\beta }_{1},{\beta }_{2},{\beta }_{3},{\beta }_{4},{\beta }_{5},{\beta }_{6},{\beta }_{7}\right)}^{T}={\left(0.1, 1, 1.2, 0.9, 0.6, 0.3, 0.3, 0.2, 0.1\right)}^{T}\) to evaluate the performance under the alternative hypothesis (Scenario 1). In addition, we let \({\varvec{\delta}}={\left(1.1, 2.2, 0, 2.1, 1.4, 0.7, 0, 0, 0\right)}^{T}\) to evaluate the empirical type I error rate (Scenario 2). For this set of simulations, the marginal response probabilities of the treatment group at the last time point are \(\left(\text{0.75,0.15,0.10}\right)\). The covariance structure used to generate data was an exchangeable structure with the log global odds ratio of \(\alpha\) as \(1, 1.5, 2\) and an autoregressivetype (ARtype) structure with the log global odds ratio of \(\alpha /\leftt{t}{\prime}\right (t\ne {t}{\prime})\) for each \(\alpha\) as \(1, 1.5, 2\). For each parameter configuration, 5,000 simulation replications were performed.
For each scenario, we fit the correct marginal mean model, and the covariate coefficients were estimated by the GEE and the BRGEE with the independent and exchangeable in terms of the crude global odds ratio suggested in Lumley [14] as the working covariance structure. The algorithm convergence criterion was less than or equal to 0.0001, and the maximum iterations allowed were set to 50. The standard errors of the GEE and the BRGEE estimated regression coefficients were estimated by substituting each parameter estimate into \({{\varvec{\Sigma}}}_{{\varvec{G}}{\varvec{E}}{\varvec{E}}}\) and \({{\varvec{\Sigma}}}_{{\varvec{B}}{\varvec{R}}{\varvec{G}}{\varvec{E}}{\varvec{E}}}\), respectively. We defined a model fit as having a convergence problem if the algorithm convergence criterion was not met or if the maximum of the absolute value of the estimated regression parameters was above 10. There were a small number of datasets where the diagonal elements of the sandwich variance estimate were negative when there was a convergence problem. At that time, we replaced it with the modelbased variance to evaluate the performance measures. We focused on the results of the coefficient \({\beta }_{1}\), which represents the treatment effect at the last time point. In this simulation study, we defined the bias as \(\text{bias}\left({\widehat{\beta }}_{1}\right)={\widehat{\beta }}_{1}{\beta }_{1}\), where \({\widehat{\beta }}_{1}={\sum }_{l=1}^{{L}_{s}}{\widehat{\beta }}_{1l}/{L}_{s}\), \({\widehat{\beta }}_{1l}\) is the estimate of \({\beta }_{1}\) for the \(l\) th replication, and \({L}_{s}\) is the number of runs in which the regression parameter was estimated. The root mean square error (RMSE) was defined by \(\text{RMSE}\left({\widehat{\beta }}_{1}\right)=\sqrt{{\sum }_{l=1}^{{L}_{s}}{({\widehat{\beta }}_{1l}{\beta }_{1})}^{2}/{L}_{s}}\). Using the standard error estimate, we constructed 95% confidence intervals by multiplying the standard error with the 2.5th percentile and 97.5th percentile of the tdistribution. Then, we defined the coverage probability of 95% confidence intervals for the parameter \({\beta }_{1}\) to be the proportion of the number of 95% confidence intervals that contain the true parameter value to the total number of runs where the regression parameter was estimated. Furthermore, the empirical type I error rate was defined as the proportion of times \(\left{\widehat{\beta }}_{1l}/{\text{SE}}({\widehat{\beta }}_{1l})\right\ge {t}_{0.975}\) for the null hypothesis, where \({t}_{0.975}\) is the 97.5th percentile of the tdistribution.
In addition, the bias, RMSE and coverage for the specified scenarios with \(\left(K,n\right)=\left(\text{3,6}\right), \left(\text{4,4}\right), \left(\text{4,6}\right)\) were investigated. And, for each scenario, coverage probability of the 95% confidence interval based on the Mancl and DeRouen’s biascorrected covariance estimates were calculated in addition to that based on \({{\varvec{\Sigma}}}_{{\varvec{G}}{\varvec{E}}{\varvec{E}}}\) and \({{\varvec{\Sigma}}}_{{\varvec{B}}{\varvec{R}}{\varvec{G}}{\varvec{E}}{\varvec{E}}}\).
Results for the convergence
First, we report the percentage of simulation replications in which there was evidence of quasicomplete separation as well as the percentage of simulation sets where the model had convergence problems for Scenario 1 and Scenario 2 (see Table S1 for the results with true exchangeable covariance structure (Scenario 1), Table S2 for the results with true exchangeable covariance structure (Scenario 2), Table S3 for the results with true ARtype covariance structure (Scenario 1), and Table S4 for the results with true ARtype covariance structure (Scenario 2) in Additional file). We defined that the simulation datasets had evidence of quasicomplete separation when all categories but one level of 1 or \(K\) of the outcome for either treatment group were zero at a time point, such as Table 1. The model fit by the GEE had the convergence problem in the presence of quasicomplete separation. For instance, the percentage of datasets where GEE estimation had a convergence problem fell within the range 7.30–8.06 when \(N=20\), 1.56–1.66 when \(N=30\), 0.28–0.40 when \(N=40\), and 0.06–0.08 when \(N=50\) for Scenario 1 with true exchangeable covariance structure. In contrast, it was 0.00% for all settings for the BRGEE. The performance of the estimates of \({\beta }_{1}\) was similar for both true covariance structures. Therefore, here, we present the results of the estimates of \({\beta }_{1}\) in the true exchangeable covariance structure only. See Additional file for the results in the ARtype covariance structure.
Results for the performance under the alternative hypothesis (Scenario 1)
In this section, we first report the distribution of the estimates. Then, we summarize the results for the simulation in terms of bias, RMSE, and coverage probability of the 95% confidence interval.
Figure 1 shows the box and whisker plot for estimates of \({\beta }_{1}\) with an exchangeable (correctly specified) and an independent (misspecified) working covariance for the true exchangeable covariance structure. The distribution of the GEE estimates had substantially large outliers for datasets with evidence of quasicomplete separation. On the other hand, the BRGEE did not have such outliers, regardless of the value of the global odds ratio and the working covariance structure for each sample size (\(N\le 50\)). See Fig. S1 for the results in the ARtype covariance structure.
Figure 2 shows the bias of the estimates of \({\beta }_{1}\). Although it decreased as the sample size increased, bias of the GEE estimates still existed when \(N=50\). On the other hand, the bias of the BRGEE estimates was close to zero regardless of the value of the global odds ratio and the working covariance structure for each sample size (\(N\le 50\)). See Fig. S2 for the results in the ARtype covariance structure.
Figure 3 shows the RMSE of the estimates of \({\beta }_{1}\). The RMSE of the BRGEE estimates was smaller than that of the GEE estimates regardless of the value of the global odds ratio and the working covariance structure for each sample size (\(N\le 50\)). See Fig. S3 for the results in the ARtype covariance structure.
Figure 4 shows the coverage probabilities of the 95% confidence interval of \({\beta }_{1}\). For \(N\)=20, the coverage probability of the 95% confidence interval based on the GEE was below the nominal level of 95%; on the other hand, that based on the BRGEE was close to 95%. In general, as the number of subjects increased (\(N\ge\) 30), the coverage probabilities of both methods became similar and slightly above 95%. See Fig. S4 for the results in the ARtype covariance structure. When the 95% confidence interval based on the BRGEE was estimated by substituting the parameter estimate into \({{\varvec{\Sigma}}}_{{\varvec{G}}{\varvec{E}}{\varvec{E}}}\) instead of \({{\varvec{\Sigma}}}_{{\varvec{B}}{\varvec{R}}{\varvec{G}}{\varvec{E}}{\varvec{E}}}\), undercoverage of the 95% confidence interval based on the BRGEE was observed for \(N\)=20 but not for \(N\ge\) 30. See Fig. S6 for the results in the exchangeable covariance structure, and Fig. S7 for those in the ARtype covariance structure.
The bias and RMSE for the scenarios with \(\left(K,n\right)=\left(\text{3,6}\right), \left(\text{4,4}\right), \left(\text{4,6}\right)\) were similar to the results when \(K=3\) and \(n=4\). The differences in coverage probabilities between GEE and BRGEE exhibited similar characteristics for \(\left(K,n\right)=\left(\text{3,4}\right)\text{ and }\left(3,6\right)\). Specifically, when \(N=20\), the coverage probabilities based on the GEE fell below the nominal level while those based on the BRGEE exceeded it. However, when \(K=4\), such differences in coverage probabilities between GEE and BRGEE were not observed. See Fig. S10 to Fig. S27 in Additional file. The coverage probability of the 95% confidence interval based on the Mancl and DeRouen’s biascorrected covariance estimates calculated by substituting the BRGEE estimates generally remained above the nominal confidence level, similar to that of the 95% confidence interval based on \({{\varvec{\Sigma}}}_{{\varvec{B}}{\varvec{R}}{\varvec{G}}{\varvec{E}}{\varvec{E}}}\). See Table. S5 to Table. S12 in Additional file.
Results for the type I error rate under the null hypothesis (Scenario 2)
Figure 5 presents the empirical type I error rate of the ttest under the null hypothesis \({H}_{0}:{\beta }_{1}=0\) resulting from the use of the GEE or the BRGEE. For \(N\)=20, although the type I error rate for the GEE was somewhat inflated, that of the BRGEE was approximately 0.05. In general, as the number of subjects increased (\(N\ge\) 30), the type I error rates of both methods became similar and slightly conservative. See Fig. S5 for the results in the ARtype covariance structure. For the empirical type I error rate of the ttest under the null hypothesis, \({H}_{0}:{\beta }_{1}=0\), when the standard error for the BRGEE as well as the GEE was estimated using the standard sandwich variance estimator, inflation of the type I error rate was observed for \(N\)=20 but not for \(N\ge\) 30. See Fig. S8 for the results in the exchangeable covariance structure, and Fig. S9 for those in the ARtype covariance structure.
Example
In this section, we applied the proposed method to data from a randomized clinical trial of patients with postoperative pain after laparoscopic cholecystectomy. One of the aims of the study was to determine whether the use of an abdominal suction drain reduces shoulder tip pain after laparoscopic surgery. Patients rated their shoulder pain levels on a visual analog scale in the morning and afternoon of the first 3 days after the operation. The dataset from Lumley [14] contained 41 patients (22 patients in the active treatment group and 19 patients in the control group) with 6 time points and 5 ordered pain score categories per patient. The pain in the control group peaked on the afternoon of day 2 [24]. Therefore, in this manuscript, we analyzed the data of 4 time points (Day 1 (am), Day 1 (pm), Day 2 (am), and Day 2 (pm)) of all 6 time points and estimated the treatment effect at Day 2 (pm). The time points for each morning and afternoon from Day 1 to Day 2 are designated as time = 1 to time = 4 in sequential order. To estimate the treatment effect at the last time point, we fit the marginal model of the proportional odds model with 9 covariates: 8 binary indicators (treatment, sex, time point, interaction of treatment and time point) and one continuous variable (age),
where \({\text{Sex}}_{i}=1\) denotes that the subject is male, \({\text{Trt}}_{i}=1\) if the subject receives the active treatment, and 0 otherwise. For this model, we estimated the regression parameter by using the GEE and the BRGEE, assuming an independent and exchangeable working covariance structure. In addition, we calculated the 95% confidence interval by using the tdistribution. In this model, the regression parameter \({\beta }_{1}\) represents the treatment effect at the last time point \(\left(\text{Time}=4\right)\).
The odds ratio of favorable response comparing the active treatment group with the control group and the 95% confidence intervals by using the GEE with independent working covariance structure, the GEE with exchangeable working covariance structure, the BRGEE with independent working covariance structure, and the BRGEE with exchangeable working covariance structure were 14.0 (3.6, 54.5), 12.8 (3.5, 47.0), 12.1 (3.6, 41.1) and 10.6 (3.2, 34.8), respectively. Each of the odds ratios can be interpreted as a common odds ratio derived from 4 logistic regression models for each possible binary dichotomization (1 vs. 2,3,4,5; 1,2 vs. 3,4,5; 1,2,3 vs. 4,5; 1,2,3,4 vs. 5) of the outcome. Overall, the estimate was smaller for the BRGEE than for the GEE. And the results from any methods suggested that subjects with the active treatment were likely to have a more favorable response when compared to subjects with the control treatment.
Discussion
In this paper, we have investigated a BRGEE to obtain GEE estimates of the proportional odds model for longitudinal ordinal data. The BRGEE was derived by applying adjusted score equations for the proportional odds model in Kosmidis [13] to the GEE as if it were equivalent to the likelihood score equation. A similar approach for smallsample bias reduction in the GEE estimate for clustered binary data was investigated in recent studies (Paul and Zhang [5], Mondol and Rahman [6], Geroldinger [7], Gosho et al. [8]).
In the simulations that are reported here, the BRGEE performs better than the standard GEE for smallsample by providing finite estimates in the presence of quasicomplete separation. This finding is consistent with that of Mondol and Rahman [6].
With the above consideration, in the estimation of the regression parameter when the proportional odds model is applied to longitudinal ordinal data in a smallsample size of 50 or less, it is advantageous to apply the BRGEE instead of the standard GEE. However, these findings are subject to the following limitations. First, our simulation setting is limited. In fact, we did not consider simulations with extreme parameter settings where datasets with no observations for a specific category throughout all time points are generated. Such settings would lead to conditional performance evaluation, and therefore we decided to avoid such configurations in this study. Second, we assumed we had complete data, so further study may be required to evaluate the properties of the BRGEE when there are missing data. In addition, the biascorrected covariance estimator for GEE with longitudinal ordinal data may also be in need of further investigation. The bias correction methods for the variance–covariance matrix in smallsample GEE estimation have primarily been evaluated assuming binary or count response data. Considering ordinal response data, we believe that new insights can be obtained by applying the previously proposed methods as well as the method proposed in this study and conducting a detailed evaluation of operational characteristics. As alternative models when proportional odds model is not appropriate, we have other types of logit including the adjacentcategories or continuationratio. Furthermore, for applying the stereotype model to longitudinal data, estimation methods using GEE was considered [25]. It is worth considering the application of biasreduction methods to these models.
Availability of data and materials
SAS code for the BRGEE can be found at https://github.com/yukiotada0/brgeegor.git. Data on postoperative pain after laparoscopic cholecystectomy can be found in Lumley [14].
Abbreviations
 GEE:

Generalized estimating equations
 BRGEE:

Biasreduced generalized estimating equations
 MSE:

Mean squared error
 RMSE:

Root mean squared error
References
Fitzmaurice GM, Laird NM, Ware JH. Applied longitudinal analysis. 2nd ed. New Jersey: John Wiley & Sons; 2011.
International Conference on Harmonization. ICH E9(R1)  Addendum on estimands and sensitivity analysis in clinical trials to the guideline on statistical principles for clinical trials. 2019. https://database.ich.org/sites/default/files/E9R1_Step4_Guideline_2019_1203.pdf. Accessed 25 Mar 2023.
Liang KY, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika. 1986;73(1):13–22. https://doi.org/10.1093/biomet/73.1.13.
Firth D. Bias reduction of maximum likelihood estimates. Biometrika. 1993;80(1):27–38. https://doi.org/10.1093/biomet/80.1.27.
Paul S, Zhang X. Small sample GEE estimation of regression parameters for longitudinal data. Statist Med. 2014;33(22):3869–81. https://doi.org/10.1002/sim.6198.
Mondol MH, Rahman MS. Biasreduced and separationproof GEE with small or sparse longitudinal binary data. Statist Med. 2019;38(14):2544–60. https://doi.org/10.1002/sim.8126.
Geroldinger A, Blagus R, Ogden H, et al. An investigation of penalization and data augmentation to improve convergence of generalized estimating equations for clustered binary outcomes. BMC Med Res Methodol. 2022;22(168). https://doi.org/10.1186/s12874022016416.
Gosho M, Ishii R, Noma H, Maruo K. A comparison of biasadjusted generalized estimating equations for sparse binary data in smallsample longitudinal studies. Statist Med. 2023;42(15):2711–27. https://doi.org/10.1002/sim.9744.
Bie R, Haneuse S, Huey N, Schildcrout J, McGee G. Fitting marginal models in small samples: a simulation study of marginalized multilevel models and generalized estimating equations. Statist Med. 2021;40:5298–312. https://doi.org/10.1002/sim.9126.
Greenland S. Smallsample bias and corrections for conditional maximumlikelihood oddsratio estimators. Biostatistics. 2000;1(1):113–22. https://doi.org/10.1093/biostatistics/1.1.113.
Agresti A. Analysis of Ordinal Categorical Data. 2nd ed. New Jersey: John Wiley & Sons; 2010.
Kosmidis I, Firth D. Bias reduction in exponential family nonlinear models. Biometrika. 2009;96(4):793–804. https://doi.org/10.1093/biomet/asp055.
Kosmidis I. Improved estimation in cumulative link models. J Royal Stat Soc Ser B Methodol. 2014;76(1):169–96. https://doi.org/10.1111/rssb.12025.
Lumley T. Generalized estimating equations for ordinal data: a note on working correlation structures. Biometrics. 1996;52(1):354–61. https://doi.org/10.2307/2533173.
Heagerty PJ, Zeger SL. ginal regression models for clustered ordinal measurements. J Am Stat Assoc. 1996;91(435):1024–36. https://doi.org/10.1080/01621459.1996.10476973.
Fahrmeir L, Pritscher L. Regression analysis of forest damage by marginal models for correlated ordinal responses. Environ Ecol Stat. 1996;3(3):257–68. https://doi.org/10.1007/BF00453014.
Lipsitz SR, Kim K, Zhao L. Analysis of repeated categorical data using generalized estimating equations. Statist Med. 1994;13(11):1149–63. https://doi.org/10.1002/sim.4780131106.
Liang KY, Zeger SL, Qaqish B. Multivariate regression analyses for categorical data. J Royal Stat Soc Ser B Methodol. 1992;54(1):3–40. https://doi.org/10.1111/j.25176161.1992.tb01862.x.
Williamson JM, Kim K, Lipsitz SR. Analyzing bivariate ordinal data using a global odds ratio. J Am Stat Assoc. 1995;90(432):1432–7. https://doi.org/10.2307/2291535.
Hines RJO. Anaysis of clustered polytomous data using generalized estimating equations and working covariance structures. Biometrics. 1997;53(4):1552–6. https://doi.org/10.2307/2533523.
Mancl LA, DeRouen TA. A covariance estimator for GEE with improved smallsample properties. Biometrics. 2001;57(1):126–34. https://doi.org/10.1111/j.0006341x.2001.00126.x.
Wang M, Kong L, Li Z, Zhang L. Covariance estimators for generalized estimating equations (GEE) in longitudinal analysis with small samples. Statist Med. 2016;35:1706–21. https://doi.org/10.1002/sim.6817.
Ibrahim NA, Suliadi S. Generating correlated discrete ordinal data using R and SAS IML. Comp Method Prog Biomed. 2011;104(3):e122–32. https://doi.org/10.1016/j.cmpb.2011.06.003.
Jorgensen JO, Gillies RB, Hunt DR, Caplehorn JR, Lumley T. A simple and effective way to reduce postoperative pain after laparoscopic cholecystectomy. Aust N Z J Surg. 1995;65(7):466–9. https://doi.org/10.1111/j.14452197.1995.tb01787.x.
Spiess M, Fernández D, Nguyen T, Liu I. Generalized estimating equations to estimate the ordered stereotype logit model for panel data. Statist Med. 2020;39:1919–40. https://doi.org/10.1002/sim.8520.
Acknowledgements
Not applicable.
Funding
This research was supported by AMED under Grant/Award Number JP21lk0201702. The funding bodies played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Author information
Authors and Affiliations
Contributions
Y.T. and T.S. contributed to the study conception and design; Y.T. conducted the simulation study and the statistical analysis; Y.T. wrote the draft of the manuscript; and T.S. provided input. All authors have read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Tada, Y., Sato, T. A biasreduced generalized estimating equation approach for proportional odds models with smallsample longitudinal ordinal data. BMC Med Res Methodol 24, 140 (2024). https://doi.org/10.1186/s12874024022596
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12874024022596