Dear All,
I want to estimate the impact of COVID-19 infections on onset of new chronic pain diagnosis. I would like to jointly model having a COVID-19 infection (logistic regression) with time to new chronic pain diagnosis, due to shared risks of both. Here is the dataex of my longitudinal data for 3 patients :
User written command -merlin- by Michael Crowther seems appropriate for the analyses, and I have followed the advice of the author to "keep the data in wide format but within each outcome the observations are long". But now I am running into "cannot compute an improvement -- discontinuous region encountered" errors:
merlin—A unified modeling framework for data analysis and methods development in Stata - Michael J. Crowther, 2020 (sagepub.com)
Any guidance on what the trouble is and ideas for troubleshooting would be much appreciated.
Hoping for help.
Gratefully,
Sumedha
I want to estimate the impact of COVID-19 infections on onset of new chronic pain diagnosis. I would like to jointly model having a COVID-19 infection (logistic regression) with time to new chronic pain diagnosis, due to shared risks of both. Here is the dataex of my longitudinal data for 3 patients :
Code:
. stset stop, id(grpatidtreat) enter(start) failure(d=1) time0(start)
Survival-time data settings
ID variable: grpatidtreat
Failure event: d==1
Observed time interval: (start, stop]
Enter on or after: time start
Exit on or before: failure
--------------------------------------------------------------------------
700,555 total observations
0 exclusions
--------------------------------------------------------------------------
700,555 observations remaining, representing
26,158 subjects
4,039 failures in single-failure-per-subject data
700,555 total analysis time at risk and under observation
At risk from t = 0
Earliest observed entry t = 0
Last observed exit t = 31
. dataex monthlydate grpatidtreat Female age80plus start pain dead COVID d t0 failtime stop _st _d _t _t0 if (grp
> atidtreat==2346371|grpatidtreat==2345602|grpatidtreat==2346479)
----------------------- copy starting from the next line -----------------------
Code:
* Example generated by -dataex-. For more info, type help dataex
clear
input float(monthlydate grpatidtreat) byte(Female age80plus) float(start pain dead COVID d t0 failtime stop) byte(_st _d _t _t0)
714 2345602 1 0 0 1 0 0 0 0 30 1 1 0 1 0
715 2345602 1 0 1 1 0 0 0 0 30 2 1 0 2 1
716 2345602 1 0 2 1 0 0 0 0 30 3 1 0 3 2
717 2345602 1 0 3 1 0 0 0 0 30 4 1 0 4 3
718 2345602 1 0 4 1 0 0 0 0 30 5 1 0 5 4
719 2345602 1 0 5 1 0 0 0 0 30 6 1 0 6 5
720 2345602 1 0 6 1 0 0 0 0 30 7 1 0 7 6
721 2345602 1 0 7 1 0 0 0 0 30 8 1 0 8 7
722 2345602 1 0 8 1 0 0 0 0 30 9 1 0 9 8
723 2345602 1 0 9 1 0 0 0 0 30 10 1 0 10 9
724 2345602 1 0 10 1 0 0 0 0 30 11 1 0 11 10
725 2345602 1 0 11 1 0 0 0 0 30 12 1 0 12 11
726 2345602 1 0 12 1 0 0 0 0 30 13 1 0 13 12
727 2345602 1 0 13 1 0 0 0 0 30 14 1 0 14 13
728 2345602 1 0 14 1 0 0 0 0 30 15 1 0 15 14
729 2345602 1 0 15 1 0 0 0 0 30 16 1 0 16 15
730 2345602 1 0 16 1 0 0 0 0 30 17 1 0 17 16
731 2345602 1 0 17 1 0 0 0 0 30 18 1 0 18 17
732 2345602 1 0 18 1 0 0 0 0 30 19 1 0 19 18
733 2345602 1 0 19 1 0 0 0 0 30 20 1 0 20 19
734 2345602 1 0 20 1 0 0 0 0 30 21 1 0 21 20
735 2345602 1 0 21 1 0 0 0 0 30 22 1 0 22 21
736 2345602 1 0 22 1 0 0 0 0 30 23 1 0 23 22
737 2345602 1 0 23 1 0 0 0 0 30 24 1 0 24 23
738 2345602 1 0 24 1 0 0 0 0 30 25 1 0 25 24
739 2345602 1 0 25 1 0 0 0 0 30 26 1 0 26 25
740 2345602 1 0 26 1 0 0 0 0 30 27 1 0 27 26
741 2345602 1 0 27 1 0 0 0 0 30 28 1 0 28 27
742 2345602 1 0 28 1 0 0 0 0 30 29 1 0 29 28
743 2345602 1 0 29 1 0 0 1 0 30 30 1 1 30 29
714 2346371 1 0 0 0 0 0 0 0 31 1 1 0 1 0
715 2346371 1 0 1 0 0 0 0 0 31 2 1 0 2 1
716 2346371 1 0 2 0 0 0 0 0 31 3 1 0 3 2
717 2346371 1 0 3 0 0 0 0 0 31 4 1 0 4 3
718 2346371 1 0 4 0 0 0 0 0 31 5 1 0 5 4
719 2346371 1 0 5 0 0 0 0 0 31 6 1 0 6 5
720 2346371 1 0 6 0 0 0 0 0 31 7 1 0 7 6
721 2346371 1 0 7 0 0 0 0 0 31 8 1 0 8 7
722 2346371 1 0 8 0 0 0 0 0 31 9 1 0 9 8
723 2346371 1 0 9 0 0 0 0 0 31 10 1 0 10 9
724 2346371 1 0 10 0 0 0 0 0 31 11 1 0 11 10
725 2346371 1 0 11 0 0 0 0 0 31 12 1 0 12 11
726 2346371 1 0 12 0 0 0 0 0 31 13 1 0 13 12
727 2346371 1 0 13 0 0 0 0 0 31 14 1 0 14 13
728 2346371 1 0 14 0 0 0 0 0 31 15 1 0 15 14
729 2346371 1 0 15 0 0 0 0 0 31 16 1 0 16 15
730 2346371 1 0 16 0 0 0 0 0 31 17 1 0 17 16
731 2346371 1 0 17 0 0 0 0 0 31 18 1 0 18 17
732 2346371 1 0 18 0 0 0 0 0 31 19 1 0 19 18
733 2346371 1 0 19 0 0 0 0 0 31 20 1 0 20 19
734 2346371 1 0 20 0 0 0 0 0 31 21 1 0 21 20
735 2346371 1 0 21 0 0 0 0 0 31 22 1 0 22 21
736 2346371 1 0 22 0 0 0 0 0 31 23 1 0 23 22
737 2346371 1 0 23 0 0 0 0 0 31 24 1 0 24 23
738 2346371 1 0 24 0 0 0 0 0 31 25 1 0 25 24
739 2346371 1 0 25 0 0 0 0 0 31 26 1 0 26 25
740 2346371 1 0 26 0 0 0 0 0 31 27 1 0 27 26
741 2346371 1 0 27 0 0 0 0 0 31 28 1 0 28 27
742 2346371 1 0 28 0 0 0 0 0 31 29 1 0 29 28
743 2346371 1 0 29 0 0 1 0 0 31 30 1 0 30 29
744 2346371 1 0 30 0 0 1 0 0 31 31 1 0 31 30
714 2346479 0 1 0 1 0 0 0 0 14 1 1 0 1 0
715 2346479 0 1 1 1 0 0 0 0 14 2 1 0 2 1
716 2346479 0 1 2 1 0 0 0 0 14 3 1 0 3 2
717 2346479 0 1 3 1 0 0 0 0 14 4 1 0 4 3
718 2346479 0 1 4 1 0 0 0 0 14 5 1 0 5 4
719 2346479 0 1 5 1 0 0 0 0 14 6 1 0 6 5
720 2346479 0 1 6 1 0 0 0 0 14 7 1 0 7 6
721 2346479 0 1 7 1 0 0 0 0 14 8 1 0 8 7
722 2346479 0 1 8 1 0 0 0 0 14 9 1 0 9 8
723 2346479 0 1 9 1 0 0 0 0 14 10 1 0 10 9
724 2346479 0 1 10 1 0 0 0 0 14 11 1 0 11 10
725 2346479 0 1 11 1 0 0 0 0 14 12 1 0 12 11
726 2346479 0 1 12 1 0 0 0 0 14 13 1 0 13 12
727 2346479 0 1 13 1 0 0 1 0 14 14 1 1 14 13
end
format %tm monthlydate
------------------ copy up to and including the previous line ------------------
Listed 75 out of 700555 observations
. bys grpatidtreat (start) : replace failtime = . if _n>1
(674397 real changes made, 674397 to missing)
. bys grpatidtreat (start) : replace pain = . if _n>1
(674397 real changes made, 674397 to missing)
. bys grpatidtreat (start) : replace dead = . if _n>1
(674397 real changes made, 674397 to missing)
. dataex monthlydate grpatidtreat Female age80plus start pain dead COVID d t0 failtime stop _st _d _t _t0 if (grp
> atidtreat==2346371|grpatidtreat==2345602|grpatidtreat==2346479)
----------------------- copy starting from the next line -----------------------
Code:
* Example generated by -dataex-. For more info, type help dataex
clear
input float(monthlydate grpatidtreat) byte(Female age80plus) float(start pain dead COVID d t0 failtime stop) byte(_st _d _t _t0)
714 2345602 1 0 0 1 0 0 0 0 30 1 1 0 1 0
715 2345602 1 0 1 . . 0 0 0 . 2 1 0 2 1
716 2345602 1 0 2 . . 0 0 0 . 3 1 0 3 2
717 2345602 1 0 3 . . 0 0 0 . 4 1 0 4 3
718 2345602 1 0 4 . . 0 0 0 . 5 1 0 5 4
719 2345602 1 0 5 . . 0 0 0 . 6 1 0 6 5
720 2345602 1 0 6 . . 0 0 0 . 7 1 0 7 6
721 2345602 1 0 7 . . 0 0 0 . 8 1 0 8 7
722 2345602 1 0 8 . . 0 0 0 . 9 1 0 9 8
723 2345602 1 0 9 . . 0 0 0 . 10 1 0 10 9
724 2345602 1 0 10 . . 0 0 0 . 11 1 0 11 10
725 2345602 1 0 11 . . 0 0 0 . 12 1 0 12 11
726 2345602 1 0 12 . . 0 0 0 . 13 1 0 13 12
727 2345602 1 0 13 . . 0 0 0 . 14 1 0 14 13
728 2345602 1 0 14 . . 0 0 0 . 15 1 0 15 14
729 2345602 1 0 15 . . 0 0 0 . 16 1 0 16 15
730 2345602 1 0 16 . . 0 0 0 . 17 1 0 17 16
731 2345602 1 0 17 . . 0 0 0 . 18 1 0 18 17
732 2345602 1 0 18 . . 0 0 0 . 19 1 0 19 18
733 2345602 1 0 19 . . 0 0 0 . 20 1 0 20 19
734 2345602 1 0 20 . . 0 0 0 . 21 1 0 21 20
735 2345602 1 0 21 . . 0 0 0 . 22 1 0 22 21
736 2345602 1 0 22 . . 0 0 0 . 23 1 0 23 22
737 2345602 1 0 23 . . 0 0 0 . 24 1 0 24 23
738 2345602 1 0 24 . . 0 0 0 . 25 1 0 25 24
739 2345602 1 0 25 . . 0 0 0 . 26 1 0 26 25
740 2345602 1 0 26 . . 0 0 0 . 27 1 0 27 26
741 2345602 1 0 27 . . 0 0 0 . 28 1 0 28 27
742 2345602 1 0 28 . . 0 0 0 . 29 1 0 29 28
743 2345602 1 0 29 . . 0 1 0 . 30 1 1 30 29
714 2346371 1 0 0 0 0 0 0 0 31 1 1 0 1 0
715 2346371 1 0 1 . . 0 0 0 . 2 1 0 2 1
716 2346371 1 0 2 . . 0 0 0 . 3 1 0 3 2
717 2346371 1 0 3 . . 0 0 0 . 4 1 0 4 3
718 2346371 1 0 4 . . 0 0 0 . 5 1 0 5 4
719 2346371 1 0 5 . . 0 0 0 . 6 1 0 6 5
720 2346371 1 0 6 . . 0 0 0 . 7 1 0 7 6
721 2346371 1 0 7 . . 0 0 0 . 8 1 0 8 7
722 2346371 1 0 8 . . 0 0 0 . 9 1 0 9 8
723 2346371 1 0 9 . . 0 0 0 . 10 1 0 10 9
724 2346371 1 0 10 . . 0 0 0 . 11 1 0 11 10
725 2346371 1 0 11 . . 0 0 0 . 12 1 0 12 11
726 2346371 1 0 12 . . 0 0 0 . 13 1 0 13 12
727 2346371 1 0 13 . . 0 0 0 . 14 1 0 14 13
728 2346371 1 0 14 . . 0 0 0 . 15 1 0 15 14
729 2346371 1 0 15 . . 0 0 0 . 16 1 0 16 15
730 2346371 1 0 16 . . 0 0 0 . 17 1 0 17 16
731 2346371 1 0 17 . . 0 0 0 . 18 1 0 18 17
732 2346371 1 0 18 . . 0 0 0 . 19 1 0 19 18
733 2346371 1 0 19 . . 0 0 0 . 20 1 0 20 19
734 2346371 1 0 20 . . 0 0 0 . 21 1 0 21 20
735 2346371 1 0 21 . . 0 0 0 . 22 1 0 22 21
736 2346371 1 0 22 . . 0 0 0 . 23 1 0 23 22
737 2346371 1 0 23 . . 0 0 0 . 24 1 0 24 23
738 2346371 1 0 24 . . 0 0 0 . 25 1 0 25 24
739 2346371 1 0 25 . . 0 0 0 . 26 1 0 26 25
740 2346371 1 0 26 . . 0 0 0 . 27 1 0 27 26
741 2346371 1 0 27 . . 0 0 0 . 28 1 0 28 27
742 2346371 1 0 28 . . 0 0 0 . 29 1 0 29 28
743 2346371 1 0 29 . . 1 0 0 . 30 1 0 30 29
744 2346371 1 0 30 . . 1 0 0 . 31 1 0 31 30
714 2346479 0 1 0 1 0 0 0 0 14 1 1 0 1 0
715 2346479 0 1 1 . . 0 0 0 . 2 1 0 2 1
716 2346479 0 1 2 . . 0 0 0 . 3 1 0 3 2
717 2346479 0 1 3 . . 0 0 0 . 4 1 0 4 3
718 2346479 0 1 4 . . 0 0 0 . 5 1 0 5 4
719 2346479 0 1 5 . . 0 0 0 . 6 1 0 6 5
720 2346479 0 1 6 . . 0 0 0 . 7 1 0 7 6
721 2346479 0 1 7 . . 0 0 0 . 8 1 0 8 7
722 2346479 0 1 8 . . 0 0 0 . 9 1 0 9 8
723 2346479 0 1 9 . . 0 0 0 . 10 1 0 10 9
724 2346479 0 1 10 . . 0 0 0 . 11 1 0 11 10
725 2346479 0 1 11 . . 0 0 0 . 12 1 0 12 11
726 2346479 0 1 12 . . 0 0 0 . 13 1 0 13 12
727 2346479 0 1 13 . . 0 1 0 . 14 1 1 14 13
end
format %tm monthlydate
------------------ copy up to and including the previous line ------------------
Listed 75 out of 700555 observations
User written command -merlin- by Michael Crowther seems appropriate for the analyses, and I have followed the advice of the author to "keep the data in wide format but within each outcome the observations are long". But now I am running into "cannot compute an improvement -- discontinuous region encountered" errors:
merlin—A unified modeling framework for data analysis and methods development in Stata - Michael J. Crowther, 2020 (sagepub.com)
Code:
. merlin (COVID _t Female age80plus , family(bernoulli) timevar(_t)) ///
> (failtime EV[COVID] Female age80plus, family(weibull, failure(pain)) timevar(failtime))
Fitting full model:
Iteration 0: Log likelihood = -1186142.7 (not concave)
Iteration 1: Log likelihood = -187814.37 (not concave)
Iteration 2: Log likelihood = -133445.55 (not concave)
Iteration 3: Log likelihood = -115807.53 (not concave)
Iteration 4: Log likelihood = -112170.1 (not concave)
Iteration 5: Log likelihood = -107129.28 (not concave)
Iteration 6: Log likelihood = -104656.06 (not concave)
Iteration 7: Log likelihood = -103565.47 (not concave)
Iteration 8: Log likelihood = -102470.69
Iteration 9: Log likelihood = -100756.81 (backed up)
Iteration 10: Log likelihood = -98454.518
Iteration 11: Log likelihood = -98272.59
Iteration 12: Log likelihood = -98248.121
Iteration 13: Log likelihood = -98247.977
Iteration 14: Log likelihood = -98247.977
Fixed effects regression model Number of obs = 700,555
Log likelihood = -98247.977
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
COVID: |
_t | .1494616 .0012058 123.95 0.000 .1470983 .151825
Female | -.0788066 .0152507 -5.17 0.000 -.1086975 -.0489157
age80plus | -.1059348 .0186377 -5.68 0.000 -.1424639 -.0694056
_cons | -6.587389 .0317108 -207.73 0.000 -6.649541 -6.525237
-------------+----------------------------------------------------------------
failtime: |
EV[] | -1.575292 .7529404 -2.09 0.036 -3.051028 -.0995561
Female | .1123753 .0321167 3.50 0.000 .0494276 .1753229
age80plus | 1.137956 .0316068 36.00 0.000 1.076008 1.199905
_cons | -5.170396 .0551396 -93.77 0.000 -5.278467 -5.062324
log(gamma) | -.1407635 .0188433 -7.47 0.000 -.1776956 -.1038314
------------------------------------------------------------------------------
.
. merlin (COVID _t Female age80plus M1[grpatidtreat]@1, family(bernoulli) timevar(_t)) ///
> (failtime EV[COVID] Female age80plus, family(weibull, failure(pain)) timevar(failtime))
Fitting fixed effects model:
Fitting full model:
Iteration 0: Log likelihood = -66848.647
Iteration 1: Log likelihood = -52114.54
Iteration 2: Log likelihood = -43047.014
Iteration 3: Log likelihood = -41463.026
Iteration 4: Log likelihood = -38762.77 (not concave)
Iteration 5: Log likelihood = -38016.249 (not concave)
Iteration 6: Log likelihood = -37524.339
Iteration 7: Log likelihood = -37220.634 (not concave)
Iteration 8: Log likelihood = -37018.575
Iteration 9: Log likelihood = -36840.872
Iteration 10: Log likelihood = -36742.457
Iteration 11: Log likelihood = -36659.819
Iteration 12: Log likelihood = -36655.375
Iteration 13: Log likelihood = -36620.241 (backed up)
cannot compute an improvement -- discontinuous region encountered
r(430);
end of do-file
r(430);
Code:
Any guidance on what the trouble is and ideas for troubleshooting would be much appreciated.
Hoping for help.
Gratefully,
Sumedha
