Ordinal logistic with Panel data and violation of the proportional odds assumption

[Nikolaos Pandis]

Hi,
I have ordinal data with repeated measurements and I would like to fit a random effects model but the PO assumption does not hold.
Is there a way to fit such a model in Stata?
Thank you,
Nikos

[Joseph Coveney]

Yes, but you'll have to program the likelihood yourself for the particular type of ordered-categorical model that you want to fit.

help menl

A search of the Internet might turn up SAS code for the task (using PROC NLMIXED) that you could adapt, if that's any consolation.

How do you know that the proportional odds assumption is so severely violated?

[Nikolaos Pandis]

Dear Joseph,
Thank you for the response.
We have verified PO violation in SAS. We can also see clearly from the evolutions in the plots the big deviations. We are using proc GLIMMIX with various options such as using a baseline model and we are trying also proc NLMIXED wihout a real solution. SAS has the option of unequal slopes but the models do not converge.
I thought I give Stata a try.
We have fitted also 2 separate logistic models like 2 vs 1 and 3 vs 1. The sum of the -2loglikelihoods is smaller (more negative-worse fit) than the one after fitting a single PO model assuming no PO violation. Is it proper to add the log-lik from the 2 models and compare it with the single model, I wonder?
Best ,
Nikos

[Joseph Coveney]

. . . we are trying also proc NLMIXED wihout a real solution.

A quick search turns up some code here (Slide 60) and here (Slide 51) that you could look into if you're planning to continue with SAS's PROC NLMIXED (change the initial values etc.), or adapt for Stata's menl. (You could also just get his SuperMix software package.)

SAS has the option of unequal slopes but the models do not converge. I thought I give Stata a try.

You could use the user-written gologit2

search gologit2

and try its cluster() option to fit a marginal model with adjusted coefficient standard errors.

We have fitted also 2 separate logistic models like 2 vs 1 and 3 vs 1. The sum of the -2loglikelihoods is smaller (more negative-worse fit) than the one after fitting a single PO model assuming no PO violation. Is it proper to add the log-lik from the 2 models and compare it with the single model, I wonder?

How about just comparing the magnitudes of the two models' corresponding regression coefficients? If they are not too different given the scientific context, then I would just go with the proportional odds model. Otherwise, I'd interpret and discuss using both models' coefficients.

[Nikolaos Pandis]

Thank you Joseph,

We need to fit a random effects (intercept and slope model).
We have managed to fit a marginal model like GEE.
I have been able to fit the model with xtlogit but assumes PO and not sure if possible to add RE intercept and slope?
I need to check how with xtologit I could get as close as possible to the desired RE model.
We have managed also to fit an adjacent type model in SAS 0 vs 1+2 and 0+1 vs 2 ( as 2 separate models) and looks ok but we get some negative probabilitis which however are for about 20 obs and the largest was - 0.0325 with several much smaller. Then we get zero variance for the intercept for one of the 2 models and then choosing the unbound option (SAS) with Laplace and we get convergence.

xtologit weightcat i.treatmen##c.month i.treatmen##c.month2

Fitting comparison model:

Iteration 0: log likelihood = -1196.5581
Iteration 1: log likelihood = -1077.2717
Iteration 2: log likelihood = -1075.0014
Iteration 3: log likelihood = -1074.9959
Iteration 4: log likelihood = -1074.9959

Refining starting values:

Grid node 0: log likelihood = -951.61366

Fitting full model:

Iteration 0: log likelihood = -951.61366
Iteration 1: log likelihood = -899.82211
Iteration 2: log likelihood = -893.82354
Iteration 3: log likelihood = -893.25404
Iteration 4: log likelihood = -893.24881
Iteration 5: log likelihood = -893.2488

Random-effects ordered logistic regression Number of obs = 1,127
Group variable: patid Number of groups = 161

Random effects u_i ~ Gaussian Obs per group:
min = 7
avg = 7.0
max = 7

Integration method: mvaghermite Integration pts. = 12

Wald chi2(5) = 228.48
Log likelihood = -893.2488 Prob > chi2 = 0.0000

-----------------------------------------------------------------------------------
weightcat | Coef. Std. Err. z P>|z| [95% Conf. Interval]
------------------+----------------------------------------------------------------
|
treatmen |
New Treatment | -0.053 0.458 -0.12 0.908 -0.951 0.846
month | -0.052 0.164 -0.32 0.750 -0.373 0.269
|
treatmen#c.month |
New Treatment | -2.898 0.276 -10.51 0.000 -3.438 -2.357
|
month2 | -0.028 0.026 -1.08 0.281 -0.079 0.023
|
treatmen#c.month2 |
New Treatment | 0.571 0.047 12.24 0.000 0.480 0.662
------------------+----------------------------------------------------------------
/cut1 | -3.616 0.342 -4.285 -2.946
/cut2 | -0.219 0.309 -0.824 0.387
------------------+----------------------------------------------------------------
/sigma2_u | 4.172 0.672 3.042 5.722
-----------------------------------------------------------------------------------
LR test vs. ologit model: chibar2(01) = 363.49 Prob >= chibar2 = 0.0000

.

Fitting the model using meologit does not converge. Also using 2 separate models one would converge (1 vs 0) and the other would not (2 vs 1).
"cannot compute an improvement -- discontinuous region encountered"

My data looks like this:
tab month weightcat,all

| weightcat
month | 0 1 2 | Total

0 | 0 76 85 | 161
1 | 25 91 45 | 161
2 | 50 65 46 | 161
3 | 59 53 49 | 161
4 | 38 64 59 | 161
5 | 35 48 78 | 161
6 | 30 46 85 | 161

Total | 237 443 447 | 1,127


I have fitted the model below and says no PO violation:

. xi:gologit2 weightcat treatmen month month2 i.treatmen*month i.treatmen*month2,autofit cluster() or
i.treatmen _Itreatmen_0-1 (naturally coded; _Itreatmen_0 omitted)
i.treat~n*month _ItreXmonth_# (coded as above)
i.trea~n*month2 _ItreXmontha# (coded as above)
note: _Itreatmen_1 dropped because of collinearity
note: month dropped because of collinearity
note: _Itreatmen_1 dropped because of collinearity
note: month2 dropped because of collinearity

------------------------------------------------------------------------------
Testing parallel lines assumption using the .05 level of significance...

Step 1: Constraints for parallel lines imposed for _ItreXmontha1 (P Value = 0.0717)
note: _Itreatmen_1 dropped because of collinearity
note: month dropped because of collinearity
note: _Itreatmen_1 dropped because of collinearity
note: month2 dropped because of collinearity
Step 2: Constraints for parallel lines imposed for treatmen (P Value = 0.1780)
note: _Itreatmen_1 dropped because of collinearity
note: month dropped because of collinearity
note: _Itreatmen_1 dropped because of collinearity
note: month2 dropped because of collinearity
Step 3: Constraints for parallel lines imposed for _ItreXmonth_1 (P Value = 0.1249)
note: _Itreatmen_1 dropped because of collinearity
note: month dropped because of collinearity
note: _Itreatmen_1 dropped because of collinearity
note: month2 dropped because of collinearity
Step 4: Constraints for parallel lines are not imposed for
month (P Value = 0.00015)
month2 (P Value = 0.01909)

Wald test of parallel lines assumption for the final model:

( 1) [0]_ItreXmontha1 - [1]_ItreXmontha1 = 0
( 2) [0]treatmen - [1]treatmen = 0
( 3) [0]_ItreXmonth_1 - [1]_ItreXmonth_1 = 0

chi2( 3) = 7.08
Prob > chi2 = 0.0694

An insignificant test statistic indicates that the final model
does not violate the proportional odds/ parallel lines assumption

If you re-estimate this exact same model with gologit2, instead
of autofit you can save time by using the parameter

pl(_ItreXmontha1 treatmen _ItreXmonth_1)

------------------------------------------------------------------------------

Generalized Ordered Logit Estimates Number of obs = 1,127
Wald chi2(7) = 222.44
Prob > chi2 = 0.0000
Log likelihood = -1049.2255 Pseudo R2 = 0.1231

( 1) [0]_ItreXmontha1 - [1]_ItreXmontha1 = 0
( 2) [0]treatmen - [1]treatmen = 0
( 3) [0]_ItreXmonth_1 - [1]_ItreXmonth_1 = 0
-------------------------------------------------------------------------------
weightcat | Odds Ratio Std. Err. z P>|z| [95% Conf. Interval]
--------------+----------------------------------------------------------------
0 |
treatmen | 0.958 0.261 -0.16 0.874 0.561 1.635
month | 0.518 0.115 -2.96 0.003 0.335 0.800
month2 | 1.044 0.034 1.34 0.181 0.980 1.113
_ItreXmonth_1 | 0.161 0.035 -8.30 0.000 0.104 0.247
_ItreXmontha1 | 1.435 0.052 9.89 0.000 1.336 1.541
_cons | 32.536 11.040 10.26 0.000 16.732 63.268
--------------+----------------------------------------------------------------
1 |
treatmen | 0.958 0.261 -0.16 0.874 0.561 1.635
month | 1.103 0.157 0.69 0.491 0.835 1.458
month2 | 0.976 0.022 -1.06 0.289 0.933 1.021
_ItreXmonth_1 | 0.161 0.035 -8.30 0.000 0.104 0.247
_ItreXmontha1 | 1.435 0.052 9.89 0.000 1.336 1.541
_cons | 0.884 0.160 -0.68 0.498 0.620 1.262
-------------------------------------------------------------------------------

However, fitting omodel shows PO violation:
xi:omodel logit weightcat i.treatmen month i.treatmen*month i.treatmen*month2
i.treatmen _Itreatmen_0-1 (naturally coded; _Itreatmen_0 omitted)
i.treat~n*month _ItreXmonth_# (coded as above)
i.trea~n*month2 _ItreXmontha# (coded as above)

note: _Itreatmen_1 dropped because of collinearity
note: month dropped because of collinearity
note: _Itreatmen_1 dropped because of collinearity
Iteration 0: log likelihood = -1196.5581
Iteration 1: log likelihood = -1077.2717
Iteration 2: log likelihood = -1075.0082
Iteration 3: log likelihood = -1074.9959
Iteration 4: log likelihood = -1074.9959

Ordered logit estimates Number of obs = 1127
LR chi2(5) = 243.12
Prob > chi2 = 0.0000
Log likelihood = -1074.9959 Pseudo R2 = 0.1016

------------------------------------------------------------------------------
weightcat | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_Itreatmen_1 | .135386 .2573981 0.53 0.599 -.369105 .6398769
month | .0299755 .1365244 0.22 0.826 -.2376074 .2975583
_ItreXmon~_1 | -1.900323 .2146415 -8.85 0.000 -2.321013 -1.479634
month2 | -.0239881 .0221225 -1.08 0.278 -.0673474 .0193713
_ItreXmon~a1 | .3680163 .0357958 10.28 0.000 .2978579 .4381748
-------------+----------------------------------------------------------------
_cut1 | -2.136234 .18599 (Ancillary parameters)
_cut2 | -.0495934 .1709424
------------------------------------------------------------------------------

Approximate likelihood-ratio test of proportionality of odds
across response categories:
chi2(5) = 59.33
Prob > chi2 = 0.0000

Then I get collinearity warnings from gologit2 and omodel but not from xtologit.
We do need the qaudratic term though as the evolution follows at first descending and then ascending curve.
I do not know if data is not sufficient or the pattern of the data.
Perhaps, not possible.
Thank you,
Nikos