Multivariate logistic regression for multiple events data

[Mengmeng Li]

Hi,

I am writing to know if using stata commend -logistic depvar covariate1 covariate2 .....- is still appropriate for my case? Or should I revise the code as necessary? Please refer to the sample data below.

The dependent variable is the event (e1-e25) outcome which is code as 0 and 1. From the sample data we can see that not all of the observation have a complete values for all 25 events. Pretty much every one of them has missing values for some events out of the 25. It is due to the structure of the database that was designed for data collection. It is not indicating that data is actually missing. It happens only because not each observation has 25 events, some of them has only 6 events, and some of them has 20, for example.
The independent variable in the sample data is only age and gonad_total. There are others but were not included here just wanted to keep the example as simple as possible. Both age and gonad_total are continuous. But later, I am planning to turn gonad_total into a categorical variable (0, 1, 2).

I did have given the dataset a try to perform the most simple multivariate logistic regression (2 covariates) and only included these three variables. What I did first was to -reshape- the dataset into a long one (number of observation rocked from 320 to over 8,000 due to those events). Coded e1-e25 to dichotomous variable (turned 1 and 2 to 0, and turned 3 to 1), coded gonad_total into categorical variable (0,1,2). Left age as continuous as it was. Then based on the long format, I ran the logistic regression code (-logistic e age gonad_total-) and got the output. I would like to verify if in this case Stata knows to count all the "0" results and "1" results of the depvar from all the events (8,000 in this case) among all 320 observations and fit it into logistic regression model as it is usually performed on a wide shape dataset (no multiple events per observation)? How Stata handle the "missing value" (not the real missing ones in my case) for the logistic model?

 
* Example generated by -dataex-. To install: ssc install dataex
clear
input float(age e1 e2 e3 e4 e5 e6 e7 e8 e9 e10 e11 e12 e13 e14 e15 e16 e17 e18 e19 e20 e21 e22 e23 e24 e25 gonad_total)
40 . 3 3 . . . . . . . . . . . . . . . . . . . . . .   3600
36 2 2 3 3 3 3 3 3 . . . . . . . . . . . . . . . . .   4950
37 1 3 2 3 3 3 3 3 3 3 . . . . . . . . . . . . . . .   3375
34 3 3 3 3 3 . . . . . . . . . . . . . . . . . . . .   4050
38 3 2 2 3 3 3 1 1 2 3 3 3 3 3 . . . . . . . . . . .   1950
38 2 3 3 3 3 3 . . . . . . . . . . . . . . . . . . .   4230
41 3 3 3 3 3 3 3 3 . . . . . . . . . . . . . . . . .   2700
37 1 3 . . . . . . . . . . . . . . . . . . . . . . .   5400
31 2 3 3 3 3 3 3 3 3 2 . . . . . . . . . . . . . . .   1350
43 3 . . . . . . . . . . . . . . . . . . . . . . . .      0
35 1 1 2 2 3 3 3 2 . . . . . . . . . . . . . . . . .   3600
34 3 3 3 3 2 1 2 3 2 3 . . . . . . . . . . . . . . .   1425
37 1 3 3 3 . . . . . . . . . . . . . . . . . . . . .   1400
37 1 2 2 2 . . . . . . . . . . . . . . . . . . . . .   3600
40 3 3 3 3 3 3 . . . . . . . . . . . . . . . . . . .   2625
36 2 1 1 2 2 1 3 3 3 3 3 . . . . . . . . . . . . . .   2150
37 2 3 3 3 3 3 3 3 3 3 3 3 3 2 2 1 3 . . . . . . . .   2175
36 1 2 2 2 1 1 3 2 3 2 2 3 3 3 3 . . . . . . . . . .   6600
34 2 3 3 3 2 3 . . . . . . . . . . . . . . . . . . .   4500
36 3 . . . . . . . . . . . . . . . . . . . . . . . .   4950
41 2 1 3 3 3 . . . . . . . . . . . . . . . . . . . .   4950
36 3 3 3 3 3 3 3 3 3 3 3 3 . . . . . . . . . . . . .   2400
34 2 3 3 2 1 3 3 3 3 3 3 3 3 3 . . . . . . . . . . . 1837.5
38 2 2 3 3 . . . . . . . . . . . . . . . . . . . . .   2775
37 2 2 1 1 1 2 1 2 1 1 3 3 3 3 3 2 2 2 1 3 3 3 2 2 .   2925
36 2 2 3 3 3 3 2 2 2 . . . . . . . . . . . . . . . .   1600
40 2 2 2 . . . . . . . . . . . . . . . . . . . . . .   4500
32 2 2 2 1 1 3 3 2 3 3 3 2 . . . . . . . . . . . . .   2925
37 3 2 3 3 . . . . . . . . . . . . . . . . . . . . .   2925
38 1 2 3 3 . . . . . . . . . . . . . . . . . . . . .   3375
39 2 1 3 2 1 1 3 3 3 3 3 2 . . . . . . . . . . . . .   4500
35 2 2 1 3 3 3 . . . . . . . . . . . . . . . . . . .   4050
34 1 1 1 . . . . . . . . . . . . . . . . . . . . . .   4500
34 3 3 . . . . . . . . . . . . . . . . . . . . . . .   4500
33 2 1 1 3 2 . . . . . . . . . . . . . . . . . . . .   3600
33 3 3 3 3 3 3 3 3 2 3 . . . . . . . . . . . . . . .   2325
30 2 2 3 3 3 3 3 3 3 1 . . . . . . . . . . . . . . .   2700
35 2 3 2 3 3 . . . . . . . . . . . . . . . . . . . .   3000
39 3 3 . . . . . . . . . . . . . . . . . . . . . . .   4050
39 . 2 1 3 2 2 1 1 . . . . . . . . . . . . . . . . .   4500
35 3 3 3 3 3 . . . . . . . . . . . . . . . . . . . .   5400
37 1 1 1 1 1 2 2 2 2 3 3 3 . . . . . . . . . . . . .   2775
33 2 3 3 3 3 3 3 3 3 3 . . . . . . . . . . . . . . .   2475
40 2 3 . . . . . . . . . . . . . . . . . . . . . . .   4560
40 2 3 3 3 3 3 3 . . . . . . . . . . . . . . . . . .   3075
35 1 1 2 1 3 3 . . . . . . . . . . . . . . . . . . .   1800
38 1 1 3 2 1 3 3 3 3 3 3 3 3 3 3 2 3 . . . . . . . .   3000
37 3 3 3 3 3 . . . . . . . . . . . . . . . . . . . .   1575
40 2 3 2 2 3 3 2 . . . . . . . . . . . . . . . . . .   3600
36 3 3 . . . . . . . . . . . . . . . . . . . . . . .   5400
40 1 2 2 1 1 3 1 1 . . . . . . . . . . . . . . . . .   3750
32 2 3 3 3 3 3 3 3 3 3 3 . . . . . . . . . . . . . .   3375
42 2 2 2 3 3 . . . . . . . . . . . . . . . . . . . .   4500
36 2 3 3 3 3 3 3 3 3 . . . . . . . . . . . . . . . .   4050
30 2 2 2 2 2 2 2 3 3 3 3 2 3 3 . . . . . . . . . . .   3000
36 1 1 3 2 2 2 3 2 3 3 3 . . . . . . . . . . . . . .   3750
31 2 3 2 3 3 3 3 . . . . . . . . . . . . . . . . . .   3675
37 1 1 1 1 1 3 3 3 3 3 3 3 3 3 2 2 2 2 . . . . . . .   3000
31 1 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 . . . . . . . .   1800
37 3 3 2 3 3 3 2 3 3 3 2 3 3 3 3 . . . . . . . . . .   3900
33 3 3 3 3 3 . . . . . . . . . . . . . . . . . . . .   2475
41 3 2 3 3 3 3 . . . . . . . . . . . . . . . . . . .   4950
37 2 2 2 3 3 . . . . . . . . . . . . . . . . . . . .   4050
41 2 2 2 2 1 3 3 3 2 3 3 2 3 3 3 3 3 3 . . . . . . .   1661
42 1 3 3 3 3 . . . . . . . . . . . . . . . . . . . .   4500
40 2 3 3 3 3 3 3 3 3 3 3 3 3 . . . . . . . . . . . .   2700
37 2 2 2 3 3 1 3 3 3 3 3 3 3 . . . . . . . . . . . .   4500
32 3 3 . . . . . . . . . . . . . . . . . . . . . . .   3600
40 . 1 1 3 . . . . . . . . . . . . . . . . . . . . .   4500
37 3 2 3 2 3 3 3 3 3 3 3 2 3 2 . . . . . . . . . . .   1800
36 2 2 2 2 2 1 3 . . . . . . . . . . . . . . . . . .   4050
35 3 3 3 3 2 2 2 . . . . . . . . . . . . . . . . . .   2900
41 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 . . . . . . . . . .   5400
36 2 2 1 . 3 3 1 1 1 2 . . . . . . . . . . . . . . .   4950
40 1 3 3 . . . . . . . . . . . . . . . . . . . . . .   4950
39 3 3 3 3 . . . . . . . . . . . . . . . . . . . . .   2100
35 2 1 3 3 3 3 3 3 2 3 2 2 2 3 3 3 3 3 3 3 . . . . .   3600
30 1 3 3 3 3 3 3 3 3 3 3 . . . . . . . . . . . . . .   1000
33 1 1 1 2 2 2 3 3 3 3 2 . . . . . . . . . . . . . .   2700
32 1 2 2 1 2 2 2 1 2 1 3 3 3 3 2 2 1 3 2 . . . . . .   2625
35 2 3 3 3 3 . . . . . . . . . . . . . . . . . . . .   4500
34 . . . . . . . . . . . . . . . . . . . . . . . . .      .
34 1 3 3 3 . . . . . . . . . . . . . . . . . . . . .   3150
33 2 2 . . . . . . . . . . . . . . . . . . . . . . .   3075
33 2 2 1 1 1 1 2 3 3 3 3 3 3 3 3 3 3 3 2 2 2 . . . .   1950
36 1 3 . . . . . . . . . . . . . . . . . . . . . . .   4500
36 3 3 3 3 3 1 3 2 1 . . . . . . . . . . . . . . . .   2700
43 1 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 . . . . . .   4950
41 3 2 3 3 . . . . . . . . . . . . . . . . . . . . .   3375
35 1 1 3 1 1 2 3 3 3 3 2 3 3 2 2 . . . . . . . . . .   4050
29 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 2 2 3 3   1700
35 1 1 1 . . . . . . . . . . . . . . . . . . . . . .   5400
36 1 2 1 3 2 . . . . . . . . . . . . . . . . . . . .   2550
39 3 3 1 3 3 3 3 . . . . . . . . . . . . . . . . . .   3075
40 2 1 1 1 3 3 2 . . . . . . . . . . . . . . . . . .   2325
38 3 3 3 2 3 3 3 3 3 3 3 . . . . . . . . . . . . . .   4950
34 1 3 3 3 3 3 3 3 . . . . . . . . . . . . . . . . .   4500
43 3 3 . . . . . . . . . . . . . . . . . . . . . . .   4050
41 1 1 2 2 2 2 2 2 3 3 . . . . . . . . . . . . . . .   4050
34 1 2 3 2 2 2 2 3 3 2 2 3 2 . . . . . . . . . . . .   3600
end

Thanks so much for helping clarify my mind.

Regards,
Mengmeng

[Joseph Coveney]

The dependent variable is the event (e1-e25) outcome which is code as 0 and 1.

From your listing, it looks as if the data are 1, 2, 3 or missing, and not 0 or 1. Are these the number of successes out of a fixed number (e.g., 3) Bernoulli trials? If not, then what are they? If so, then why aren't there at least an occasional zero successes out of three trials?

From the sample data we can see that not all of the observation have a complete values for all 25 events. Pretty much every one of them has missing values for some events out of the 25. It is due to the structure of the database that was designed for data collection. It is not indicating that data is actually missing. It happens only because not each observation has 25 events, some of them has only 6 events, and some of them has 20, for example.

You'll need to explain just what these "events" are. With a variable number of opportunities for an "event" to occur, logistic regression might not be a suitable method, or at least you might need an offset variable or something.

How Stata handle the "missing value" (not the real missing ones in my case) for the logistic model?

logistic ignores observations where the dependent variable is missing-valued (>= .).

[Mengmeng Li]

From your listing, it looks as if the data are 1, 2, 3 or missing, and not 0 or 1. Are these the number of successes out of a fixed number (e.g., 3) Bernoulli trials? If not, then what are they? If so, then why aren't there at least an occasional zero successes out of three trials?

You'll need to explain just what these "events" are. With a variable number of opportunities for an "event" to occur, logistic regression might not be a suitable method, or at least you might need an offset variable or something.

logistic ignores observations where the dependent variable is missing-valued (>= .).

Hi Joseph,

Thanks for your kindly reply.

For your first question, event variables were coded as 1, 2, 3 as you saw in the sample. However, before performing -logistic-, I recoded them to 0 or 1 (1 and 2 to 0, 3 to 1). They are actually the by-products produced by each observation at one time under the same condition per observation (however conditions are different among different observations, indicated by multiple covariates). 1,2,3 does not mean a "yes" or "no". It actually indicates the different degrees of the same one aspect regarding the by-product. Different observation has different number of by-products. Some have 5 by-products, some have less or more.

For you second question, please see above. What do you mean by "offset variable"? How to incorporate it into the logistic regression model? Or should I use other models instead?

Thanks so much!!!

Regards,
Mengmeng

[Joseph Coveney]

I don't know what you're trying to do, what research question you're trying to answer, but I recommend against coarsening the outcome variable if it's avoidable. There are Stata estimation commands that will accommodate ordered-categorical outcomes. I show an example below that uses one and mentions another that you could consider. (After fitting a model using an estimation command for an ordered-categorical outcome, you can get your 0/1 prediction—if that's helpful for interpretation—by predicting the probability of outcome == 3. See the help file for postestimation commands associated with the estimation command that I use in the illustration below.)

In example code below, I assume that the values of the by-products are exchangeable, that there's no particular information in the indexes 1 through 25. Given how closely you're holding the cards to your chest, the models below are perforce entirely exploratory. But you can run the code and see whether it helps clarify whatever relationship it is between whatever phenomena they are that you're looking into.

version 14.2

clear *
set more off

input double(age e1 e2 e3 e4 e5 e6 e7 e8 e9 e10 e11 e12 e13 ///
    e14 e15 e16 e17 e18 e19 e20 e21 e22 e23 e24 e25 gonad_total)
40 . 3 3 . . . . . . . . . . . . . . . . . . . . . .   3600
36 2 2 3 3 3 3 3 3 . . . . . . . . . . . . . . . . .   4950
37 1 3 2 3 3 3 3 3 3 3 . . . . . . . . . . . . . . .   3375
34 3 3 3 3 3 . . . . . . . . . . . . . . . . . . . .   4050
38 3 2 2 3 3 3 1 1 2 3 3 3 3 3 . . . . . . . . . . .   1950
38 2 3 3 3 3 3 . . . . . . . . . . . . . . . . . . .   4230
41 3 3 3 3 3 3 3 3 . . . . . . . . . . . . . . . . .   2700
37 1 3 . . . . . . . . . . . . . . . . . . . . . . .   5400
31 2 3 3 3 3 3 3 3 3 2 . . . . . . . . . . . . . . .   1350
43 3 . . . . . . . . . . . . . . . . . . . . . . . .      0
35 1 1 2 2 3 3 3 2 . . . . . . . . . . . . . . . . .   3600
34 3 3 3 3 2 1 2 3 2 3 . . . . . . . . . . . . . . .   1425
37 1 3 3 3 . . . . . . . . . . . . . . . . . . . . .   1400
37 1 2 2 2 . . . . . . . . . . . . . . . . . . . . .   3600
40 3 3 3 3 3 3 . . . . . . . . . . . . . . . . . . .   2625
36 2 1 1 2 2 1 3 3 3 3 3 . . . . . . . . . . . . . .   2150
37 2 3 3 3 3 3 3 3 3 3 3 3 3 2 2 1 3 . . . . . . . .   2175
36 1 2 2 2 1 1 3 2 3 2 2 3 3 3 3 . . . . . . . . . .   6600
34 2 3 3 3 2 3 . . . . . . . . . . . . . . . . . . .   4500
36 3 . . . . . . . . . . . . . . . . . . . . . . . .   4950
41 2 1 3 3 3 . . . . . . . . . . . . . . . . . . . .   4950
36 3 3 3 3 3 3 3 3 3 3 3 3 . . . . . . . . . . . . .   2400
34 2 3 3 2 1 3 3 3 3 3 3 3 3 3 . . . . . . . . . . . 1837.5
38 2 2 3 3 . . . . . . . . . . . . . . . . . . . . .   2775
37 2 2 1 1 1 2 1 2 1 1 3 3 3 3 3 2 2 2 1 3 3 3 2 2 .   2925
36 2 2 3 3 3 3 2 2 2 . . . . . . . . . . . . . . . .   1600
40 2 2 2 . . . . . . . . . . . . . . . . . . . . . .   4500
32 2 2 2 1 1 3 3 2 3 3 3 2 . . . . . . . . . . . . .   2925
37 3 2 3 3 . . . . . . . . . . . . . . . . . . . . .   2925
38 1 2 3 3 . . . . . . . . . . . . . . . . . . . . .   3375
39 2 1 3 2 1 1 3 3 3 3 3 2 . . . . . . . . . . . . .   4500
35 2 2 1 3 3 3 . . . . . . . . . . . . . . . . . . .   4050
34 1 1 1 . . . . . . . . . . . . . . . . . . . . . .   4500
34 3 3 . . . . . . . . . . . . . . . . . . . . . . .   4500
33 2 1 1 3 2 . . . . . . . . . . . . . . . . . . . .   3600
33 3 3 3 3 3 3 3 3 2 3 . . . . . . . . . . . . . . .   2325
30 2 2 3 3 3 3 3 3 3 1 . . . . . . . . . . . . . . .   2700
35 2 3 2 3 3 . . . . . . . . . . . . . . . . . . . .   3000
39 3 3 . . . . . . . . . . . . . . . . . . . . . . .   4050
39 . 2 1 3 2 2 1 1 . . . . . . . . . . . . . . . . .   4500
35 3 3 3 3 3 . . . . . . . . . . . . . . . . . . . .   5400
37 1 1 1 1 1 2 2 2 2 3 3 3 . . . . . . . . . . . . .   2775
33 2 3 3 3 3 3 3 3 3 3 . . . . . . . . . . . . . . .   2475
40 2 3 . . . . . . . . . . . . . . . . . . . . . . .   4560
40 2 3 3 3 3 3 3 . . . . . . . . . . . . . . . . . .   3075
35 1 1 2 1 3 3 . . . . . . . . . . . . . . . . . . .   1800
38 1 1 3 2 1 3 3 3 3 3 3 3 3 3 3 2 3 . . . . . . . .   3000
37 3 3 3 3 3 . . . . . . . . . . . . . . . . . . . .   1575
40 2 3 2 2 3 3 2 . . . . . . . . . . . . . . . . . .   3600
36 3 3 . . . . . . . . . . . . . . . . . . . . . . .   5400
40 1 2 2 1 1 3 1 1 . . . . . . . . . . . . . . . . .   3750
32 2 3 3 3 3 3 3 3 3 3 3 . . . . . . . . . . . . . .   3375
42 2 2 2 3 3 . . . . . . . . . . . . . . . . . . . .   4500
36 2 3 3 3 3 3 3 3 3 . . . . . . . . . . . . . . . .   4050
30 2 2 2 2 2 2 2 3 3 3 3 2 3 3 . . . . . . . . . . .   3000
36 1 1 3 2 2 2 3 2 3 3 3 . . . . . . . . . . . . . .   3750
31 2 3 2 3 3 3 3 . . . . . . . . . . . . . . . . . .   3675
37 1 1 1 1 1 3 3 3 3 3 3 3 3 3 2 2 2 2 . . . . . . .   3000
31 1 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 . . . . . . . .   1800
37 3 3 2 3 3 3 2 3 3 3 2 3 3 3 3 . . . . . . . . . .   3900
33 3 3 3 3 3 . . . . . . . . . . . . . . . . . . . .   2475
41 3 2 3 3 3 3 . . . . . . . . . . . . . . . . . . .   4950
37 2 2 2 3 3 . . . . . . . . . . . . . . . . . . . .   4050
41 2 2 2 2 1 3 3 3 2 3 3 2 3 3 3 3 3 3 . . . . . . .   1661
42 1 3 3 3 3 . . . . . . . . . . . . . . . . . . . .   4500
40 2 3 3 3 3 3 3 3 3 3 3 3 3 . . . . . . . . . . . .   2700
37 2 2 2 3 3 1 3 3 3 3 3 3 3 . . . . . . . . . . . .   4500
32 3 3 . . . . . . . . . . . . . . . . . . . . . . .   3600
40 . 1 1 3 . . . . . . . . . . . . . . . . . . . . .   4500
37 3 2 3 2 3 3 3 3 3 3 3 2 3 2 . . . . . . . . . . .   1800
36 2 2 2 2 2 1 3 . . . . . . . . . . . . . . . . . .   4050
35 3 3 3 3 2 2 2 . . . . . . . . . . . . . . . . . .   2900
41 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 . . . . . . . . . .   5400
36 2 2 1 . 3 3 1 1 1 2 . . . . . . . . . . . . . . .   4950
40 1 3 3 . . . . . . . . . . . . . . . . . . . . . .   4950
39 3 3 3 3 . . . . . . . . . . . . . . . . . . . . .   2100
35 2 1 3 3 3 3 3 3 2 3 2 2 2 3 3 3 3 3 3 3 . . . . .   3600
30 1 3 3 3 3 3 3 3 3 3 3 . . . . . . . . . . . . . .   1000
33 1 1 1 2 2 2 3 3 3 3 2 . . . . . . . . . . . . . .   2700
32 1 2 2 1 2 2 2 1 2 1 3 3 3 3 2 2 1 3 2 . . . . . .   2625
35 2 3 3 3 3 . . . . . . . . . . . . . . . . . . . .   4500
34 . . . . . . . . . . . . . . . . . . . . . . . . .      .
34 1 3 3 3 . . . . . . . . . . . . . . . . . . . . .   3150
33 2 2 . . . . . . . . . . . . . . . . . . . . . . .   3075
33 2 2 1 1 1 1 2 3 3 3 3 3 3 3 3 3 3 3 2 2 2 . . . .   1950
36 1 3 . . . . . . . . . . . . . . . . . . . . . . .   4500
36 3 3 3 3 3 1 3 2 1 . . . . . . . . . . . . . . . .   2700
43 1 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 . . . . . .   4950
41 3 2 3 3 . . . . . . . . . . . . . . . . . . . . .   3375
35 1 1 3 1 1 2 3 3 3 3 2 3 3 2 2 . . . . . . . . . .   4050
29 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 2 2 3 3   1700
35 1 1 1 . . . . . . . . . . . . . . . . . . . . . .   5400
36 1 2 1 3 2 . . . . . . . . . . . . . . . . . . . .   2550
39 3 3 1 3 3 3 3 . . . . . . . . . . . . . . . . . .   3075
40 2 1 1 1 3 3 2 . . . . . . . . . . . . . . . . . .   2325
38 3 3 3 2 3 3 3 3 3 3 3 . . . . . . . . . . . . . .   4950
34 1 3 3 3 3 3 3 3 . . . . . . . . . . . . . . . . .   4500
43 3 3 . . . . . . . . . . . . . . . . . . . . . . .   4050
41 1 1 2 2 2 2 2 2 3 3 . . . . . . . . . . . . . . .   4050
34 1 2 3 2 2 2 2 3 3 2 2 3 2 . . . . . . . . . . . .   3600
end

generate int pid = _n
quietly reshape long e, i(pid) j(seq)
quietly drop if mi(e)
quietly bysort pid (seq): replace seq = _n

foreach var of varlist age gonad_total {
    summarize `var', meanonly
    generate double `=substr("c`var'", 1, 3)' = `var' - r(mean)
}

ologit e c.cgo c.cag c.cgo#c.cgo c.cag#c.cag c.cgo#c.cag, vce(cluster pid) nolog
/* alternative: meologit e c.cgo c.cag c.cgo#c.cgo c.cag#c.cag c.cgo#c.cag || pid: , nolog */
predict double xb, xb

foreach var of varlist cag cgo {
    quietly centile `var', centile(10(10)90)
    forvalues i = 1/`=r(n_cent)' {
        local `var'_list ``var'_list' `=round(r(c_`i'), 1)'
    }
}
di `cag_list'
margins , predict(xb) atmeans at(cag = (`cag_list'))
marginsplot ,  plotopts(lcolor(black) mcolor(black)) ciopts(lcolor(black)) ///
    level(50) xtitle(Mean-centered Age (d)) ylabel( , angle(horizontal) nogrid)
pause on
pause

margins , predict(xb) atmeans at(cgo = (`cgo_list'))
marginsplot ,  plotopts(lcolor(black) mcolor(black)) ciopts(lcolor(black)) ///
    level(50) xtitle(Mean-centered Gonad Tally) xlabel( , angle(45)) ///
    ylabel( , angle(horizontal) nogrid)

quietly bysort pid (seq): keep if _n == _N
nbreg seq c.cag c.cgo c.cag#c.cag c.cgo#c.cgo c.cag#c.cgo, nolog

exit