For my advanced research design course this semester I have been providing code snippets in Stata and R. This is the first time I’ve really sat down and programmed extensively in Stata, and this is a followup to produce some of the same plots and model fit statistics for group based trajectory statistics as this post in R. The code and the simulated data I made to reproduce this analysis can be downloaded here.

First, for my own notes my version of Stata is on a server here at UT Dallas. So I cannot simply go

```
net from http://www.andrew.cmu.edu/user/bjones/traj
net install traj, force
```

to install the group based trajectory code. First I have this set of code in the header of all my do files now

```
*let stata know to search for a new location for stata plug ins
adopath + "C:\Users\axw161530\Documents\Stata_PlugIns"
*to install on your own in the lab it would be
net set ado "C:\Users\axw161530\Documents\Stata_PlugIns"
```

So after running that then I can install the `traj`

command, and Stata will know where to look for it!

Once that is taken care of after setting the working directory, I can simply load in the csv file. Here my first variable was read in as `Ã¯id`

instead of `id`

(I’m thinking because of the encoding in the csv file). So I rename that variable to `id`

.

```
*Load in csv file
import delimited GroupTraj_Sim.csv
*BOM mark makes the id variable weird
rename Ã¯id id
```

Second, the `traj`

model expects the data in wide format (which this data set already is), and has counts in `count_1`

, `count_2`

… `count_10`

. The `traj`

command also wants you to input a time variable though to, which I do not have in this file. So I create a set of `t_`

variables to mimic the counts, going from 1 to 10.

```
*Need to generate a set of time variables to pass to traj, just label 1 to 10
forval i = 1/10 {
generate t_`i' = `i'
}
```

Now we can estimate our group based models, and we get a pretty nice default plot.

```
traj, var(count_*) indep(t_*) model(zip) order(2 2 2) iorder(0)
trajplot
```

Now for absolute model fit statistics, there are the average posterior probabilities, the odds of correct classification, and the observed classification proportion versus the expected classification proportion. Here I made a program that I will surely be ashamed of later (I should not brutalize the data and do all the calculations in matrix), but it works and produces an ugly table to get us these stats.

```
*I made a function to print out summary stats
program summary_table_procTraj
preserve
*now lets look at the average posterior probability
gen Mp = 0
foreach i of varlist _traj_ProbG* {
replace Mp = `i' if `i' > Mp
}
sort _traj_Group
*and the odds of correct classification
by _traj_Group: gen countG = _N
by _traj_Group: egen groupAPP = mean(Mp)
by _traj_Group: gen counter = _n
gen n = groupAPP/(1 - groupAPP)
gen p = countG/ _N
gen d = p/(1-p)
gen occ = n/d
*Estimated proportion for each group
scalar c = 0
gen TotProb = 0
foreach i of varlist _traj_ProbG* {
scalar c = c + 1
quietly summarize `i'
replace TotProb = r(sum)/ _N if _traj_Group == c
}
*This displays the group number, the count per group, the average posterior probability for each group,
*the odds of correct classification, and the observed probability of groups versus the probability
*based on the posterior probabilities
list _traj_Group countG groupAPP occ p TotProb if counter == 1
restore
end
```

This should work after any model as long as the naming conventions for the assigned groups are `_traj_Group`

and the posterior probabilities are in the variables `_traj_ProbG*`

. So when you run

`summary_table_procTraj`

You get this ugly table:

```
| _traj_~p countG groupAPP occ p TotProb |
|------------------------------------------------------------|
1. | 1 103 .937933 43.57431 .2575 .2573651 |
104. | 2 161 .9935606 229.0434 .4025 .4064893 |
265. | 3 136 .9607248 47.48378 .34 .3361456 |
```

The `groupAPP`

are the average posterior probabilities – here you can see they are all quite high. `occ`

is the odds of correct classification, and again they are all quite high. `p`

is the proportion in each group based on the assignments for the maximum posterior probability, and the `TotProb`

are the expected number based on the sums of the posterior probabilities. `TotProb`

should be the same as in the Group Membership part at the bottom of the `traj`

model. They are close (to 5 decimals), but not exactly the same (and I do not know why that is the case).

Next, to generate a plot of the individual trajectories, I want to reshape the data to long format. I use `preserve`

in case I want to go back to wide format later and estimate more models. If you need to look to see how the `reshape`

command works, type `help reshape`

at the Stata prompt. (Ditto for help on all Stata commands.)

```
preserve
reshape long count_ t_, i(id)
```

To get the behavior I want in the plot, I use a scatter plot but have them connected via `c(L)`

. Then I create small multiples for each trajectory group using the `by()`

option. Before that I slightly jitter the count data, so the lines are not always perfectly overlapped. I make the lines thin and grey — I would also use transparency but Stata graphs do not support this.

```
gen count_jit = count_ + ( 0.2*runiform()-0.1 )
graph twoway scatter count_jit t_, c(L) by(_traj_Group) msize(tiny) mcolor(gray) lwidth(vthin) lcolor(gray)
```

I’m too lazy at the moment to clean up the axis titles and such, but I think this plot of the individual trajectories should always be done. See *Breaking Bad: Two Decades of Life-Course Data Analysis in Criminology, Developmental Psychology, and Beyond* (Erosheva et al., 2014).

While this fit looks good, this is not the correct number of groups given how I simulated the data. I will give those trying to find the right answer a few hints; none of the groups have a higher polynomial than 2, and there is a constant zero inflation for the entire sample, so `iorder(0)`

will be the correct specification for the zero inflation part. If you take a stab at it let me know, I will fill you in on how I generated the simulation.

## mario spiezio

/ April 19, 2018Dear Prof Wheeler, thanks a lot for the code. It is really helpful. I was wondering whether it is possible to use margins after after traj to compute average marginal effects of the covariates included in risk(). Thank you.

## apwheele

/ April 19, 2018I don’t think so, it appears you will need to calculate that effect yourself. I’m not sure exactly how you would do it with the mixture model.

## Dimi

/ May 7, 2018Dear Prof Wheeler,

Thank you for the code.

For my research, I am trying to use the traj plugin for a survival analysis. As I am not a statistician (I am a physician) you can imagine I need a little help.

My initial issue is that, although I do not have missing values in the long format, in the wide format (due to the fact that some patients died during the follow-up) I do. When I run the code, these patients with “missing values” are dropped from the analysis. Do you know how to overcome this problem?

All my best,

Dimi

## apwheele

/ August 13, 2018I am not sure the best way to deal with missing data in such models.

## lori

/ August 12, 2018can you covary outcome group when you have a risk variable?

## apwheele

/ August 13, 2018Not exactly sure what you mean — are you asking about a joint trajectory model with multiple outcomes and predicting probability of being in a particular trajectory?

Also to see various examples you can type

help traj

into Stata and it provides various examples.

## lori

/ August 13, 2018I have an assessment that is completed across 4 different time points – 12, 18, 24, and 36 months. I am putting scores on a different assessment completed at 6 months as a risk variable. Can I also covary outcome group along with this risk marker to control for diagnosis at 36 months – the outcome group)?

## apwheele

/ August 13, 2018I think that would just be something like

traj, var(y1-y4) model(logit) order(2 2 2 2) risk(pre)

if I am understanding correctly. [Sorry if I am not!]

## lori

/ August 13, 2018My formulas is traj, var(abc*) indep(time*) model(cnorm) min(53) max(142) order(3) risk(mu6cmss)

When i try to covary for outcome group at the last time point, I get:

. traj, var(abc*) indep(time*) model(cnorm) min(53) max(142) order(3) risk(mu6cmss) cov(dxgroup)

option cov() not allowed

r(198);

. traj, var(abc*) indep(time*) model(cnorm) min(53) max(142) order(3) risk(mu6cmss) tcov(dxgroup)

The number of variables in tcov1 and var1 must match or be a multiple.

r(198);

var has 4 time points, dxgroup has 1 time point

## apwheele

/ August 13, 2018tcov is for time-varying covariates. If it is not time varying, it can be used to predict what trajectory group a person in likely to be in, but not the shape of the trajectory over time. Not time-varying covariates should probably go in the “risk()” option in most cases then.

What exactly is “dxgroup”?

## lori

/ August 13, 2018dxgroup is the outcome group. We collect data at 12, 18, 24, and 36 months. Then at 36 months, we do a diagnostic check to see if any of our kids end up with a diagnosis of autism. I want to be able to add this into my formula to covariate out any effects of outcome group to check if my predictor variable (an assessment completed at 6 months of age) can predict trajectory membership above and beyond outcome group.

## apwheele

/ August 13, 2018Given that description your formula should be:

traj, var(abc*) indep(time*) model(cnorm) min(53) max(142) order(3) risk(mu6cmss dxgroup)

with only one group though this is probably not identified. So you will need to change it to something like:

“order(2 2)” for two trajectory groups or

“order(2 2 2)” for three trajectory groups etc.

(With only four time points I would not do a cubic, about the best you can do is a quadratic).

## lori

/ August 13, 2018If you have two risk variables, does it matter what order you put them in, because if i put the 6 month assessment first followed by dxgroup, i get different information.

## apwheele

/ August 13, 2018If you mean should

“risk(mu6cmss dxgroup)”

give a different result than

“risk(dxgroup mu6cmss)”

it should not. One caveat is that the labelling of the trajectory groups is arbitrary, so they could in theory converge to the same estimates, but different groups get different labels (and subsequently the coefficients predicting trajectory membership may then be altered).

Given that the mixture models have difficulty converging that could also be a culprit as well.

This is a bit awkward to give so many comments, just send me an email and it will be easier for me to respond.

## Yong Kyu Lee

/ August 16, 2018Thank you for the thorough review.

I had a great help from it.

I usually use SAS, so I am not familiar with STATA, but by some political reason I have to do GBTM in STATA.

And I am wondering a way to get a file, as in ‘SAS ods output’, of the result of this GBTM which is ID and designated group by posterior probability.

## apwheele

/ August 17, 2018You can just save the results to a csv file. After you run the traj command, something like:

********************

preserve

keep ID _traj*

outsheet using “TrajResults.csv:, replace comma

restore

********************