I use your stata command traj to find the group-based trajectories. Usually, the comparison of BIC values from 1 to x groups leads to the decision to select the model with optimal groups.This command gives the optimal number of groups. I have read the help file for many times but i cannot find the way to get BIC values for different numbers of groups.
Does anyone use stata command traj to conduct group-based trajectory modeling? And how can i get the BIC values from models with different number of groups?
I have check the SAS codes. They include an option that I can set the number of groups. However, I have not installed SAS software in my computer.
Thank you!
Yimin
I am attaching two studies that are somewhat like what you want to do. Read the papers and compare with what you have. Look at the references. This should help you get back.to your problem. Please feel free to ask any questions you want to. Best wishes, David Booth
You can access the relevant parameters with the use of ereturn list after each trajectory analysis, i.e. the number of groups is stored in `e(numGroups1)', and the corresponding BIC values are stores in `e(BIC_N_data)' and `e(BIC_n_subjects)'.
If you write a loop using the forvalues command, you can easily set the range you want to test, store the values in a local macro and ultimately display them at the end.
I have written a concise example that is built upon the first example in the help file:
use "http://www.andrew.cmu.edu/user/bjones/traj/data/montreal_sim.dta", clear
*SET PARAMETERS HERE (i.e. number of polynoms, max. number of groups)
local polynom 1
local maxnum 5
matrix a = J(`maxnum',3,.)
matrix colnames a = "Groups" "BIC(N)" "BIC(panels)"
local order `polynom'
forvalues i = 1/`maxnum' {
local groupsnum `:list sizeof order'
quietly traj, var(qcp*op) indep(age*) model(cnorm) min(0) max(10) order(`order')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= e(BIC_N_data)
matrix a[`i',3]= e(BIC_n_subjects)
local order `order' `polynom'
}
matrix list a
Dear David and Jan, thank you for your answers.
Jan Helmdag I have run the codes you provided. It is useful. However, It can only compare models with same polynomials and find the optimal number of groups. For example, if i set the number of polynomial as 1, these codes will give a comparison of models with linear terms.
I wonder if there is a way to write a loop that can list the BIC values of all models with the following set:
the number of groups: from 2 to 5; the number of polynomials: from 1 to 4.
It costs me a lot of time to find the reasonable combinations of the numbers of polynomials. For example, when I set the number of polynomials as 1 and find the optimal number of groups is 3, then I have to compare 3-group models with different combinations of polynomials, like (1 1 2), ( 1 2 2 ), ( 1 3 2) ....Although the shapes of trajectories in the graph and the p values can help me to exclude some combinations, it is still a huge task.
Thank you very much!
Yimin Mao: There is a solution to that particular problem and requires a bit of programming and applying loops. It may not be the most elegant of all solutions, but it does trick. Be aware that applying permutations of polynomials in the trajectory framework consumes lots of time – depending on the size of the dataset you're using. To circumvent this problem, I always use test datasets (approx. 10 percent of original size) for bigger datasets. Also, feel free to recommend my answers if they helped you.
*Permutations of polynomials and group numbers
use "http://www.andrew.cmu.edu/user/bjones/traj/data/montreal_sim.dta", clear
local polynom 4
local maxnum: display `polynom'^1+`polynom'^2+`polynom'^3+`polynom'^4
matrix a = J(`maxnum',4,.)
matrix colnames a = "Groups" "Polynomials" "BIC(N)" "BIC(panels)"
local i 1
local bicn: display -10^99
local bicp: display -10^99
local e(BIC_N_data): display -10^99
local e(BIC_N_subjects): display -10^99
forvalues a = 1/`polynom' {
quietly traj, var(qcp*op) indep(age*) model(cnorm) min(0) max(10) order(`a')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'
}
local i `++i'
forvalues b = 1/`polynom' {
quietly traj, var(qcp*op) indep(age*) model(cnorm) min(0) max(10) order(`a' `b')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'`b'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'`b'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'`b'
}
local i `++i'
forvalues c = 1/`polynom' {
quietly traj, var(qcp*op) indep(age*) model(cnorm) min(0) max(10) order(`a' `b' `c')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'`b'`c'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'`b'`c'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'`b'`c'
}
local i `++i'
forvalues d = 1/`polynom' {
quietly traj, var(qcp*op) indep(age*) model(cnorm) min(0) max(10) order(`a' `b' `c' `d')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'`b'`c'`d'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'`b'`c'`d'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'`b'`c'`d'
}
local i `++i'
}
}
}
}
matrix list a
display "Best solution BIC(n): " `solutionn'
display "Best solution BIC(p): " `solutionp'
Hi, thank you Jan Helmdag for your great reply. I have used your code, in Stata, to run iterations of traj on my dataset. I also logged all output to later extract all error messages (such as "variance matrix is nonsymmetric or highly singular" and so on) and save everything in a matrix. Your code gave me a really good start. Thanks.
I have another question regarding traj. I have fitted a model with 5 trajectories. I now want to supply different start values to my model in order to be sure I have found the global, and not the local, maxima in my model. Then I want to run the model with different start values, for let say 200 times, and then plot and compare the BICs and AICs. Depending on how they differ I can then draw the conclusion if I have the local or global maxima.
I've read (https://www.andrew.cmu.edu/user/bjones/example.htm) that proc traj in Stata and SPSS uses the following to provide start values: "TRAJ procedure provides default start values by calculating group intercepts spaced based on the range or standard deviation of the dependent variables".
Anyone knows how to randomly assign different start values to my model? Using "seed" does not work since proc TRAJ calculates default start values every time so it ends up giving the same start values regardless of seed settings.
All the best,
Jesper Eriksson
Jesper Eriksson: I think the trajstart command does exactly what you need. Maybe check if you have the most recent version of traj installed.
Hello everyone,
Thanks for the great thread Jan Helmdag . I have been playing around with the codes provided here but I am consistently receiving error messages.
"Warning: variance matrix is nonsymmetric or highly singular"
+ the code takes an extremely long time to run. Im not sure if this is some error from my part or the trajectory piece takes up a lot of space.
Thoughts?
Thank you in advance!
Hi Alexandra,
yes, unfortunately the trajectory analysis can take up huge amounts of computing time. I use a dedicated work station for my analyses and always make sure that the the entire code runs beforehand. I am always using sub-samples of my data to check if a code runs through, since this saves me lots of time. If you do not have access to an additional computer (or a computing unit/cluster/server provided by your institution) , I recommend running the analyses over night and use a log file to review the entire output once it is finished.
When you get the error message that the variance matrix is non-symmetric or highly singular, the traj procedure fails to compute point estimates and standard errors. In my experience, this typically happens when you try to estimate too many groups or a specific combination of polynomials that just does not fit the data.
Jan Helmdag thank you so much for your quick and detailed response. I really appreciate it.
When trying to understand the combination of polynomials to fit the data, is there anything you recommend I try?
Does this also include possibly testing different types of models ( poisson vs normal?).
The answers to those issues heavily depend on your data and its distribution. A very comprehensive, yet easy to understand introduction on trajectory analyses that covers all those issues can be found in Nagin (2005). Chapter 4 of the book is dedicated to model selection and the entire second part is concerned with its subsequent generalization. I highly recommend giving the book a try.
Nagin, D. S. (2005). Group-Based Modeling of Development. Harvard University Press, Cambridge.
Jan Helmdag thank you! I will try that. you have been super helpful
I am running GBTM in Stata software for my count data using the zip model. Upon running the traj syntax i.e., “traj, var(y*) indep(t*) model(zip) order(x x x), a particular error pops up i.e., ‘Likelihood could not be computed at start values’. i was advised that there is an issue with the start values in your data and need to generate alternative start values to over-ride the default start values.
After reading the above messages and reviewing the notes on https://www.andrew.cmu.edu/user/bjones/strtxmpl1.htm; I tried to do two things:
1. I ran the syntax “traj, var(y*) indep(t*) model(zip) order(x x x), detail
2. Instead of considering the zip model for count data, I used the cnorm model by putting min(0) and max(999) values and ran the syntax "traj, var(y*) indep(t*) model(cnorm) min(0) max(999) order(x x x), detail.
Following these syntax's, the estimate parameters were surely generated but i got an error message i.e., "Warning: variance matrix is nonsymmetric or highly singular".
Can anyone please provide me some useful advice on how to generate new starting values so that i can produce the correct trajectories.
Thanks for your kind help.
Thank you Jan Helmdag for the code. Jesper Eriksson I would be very interested in seeing your developments, that you describe in your first paragraph, on the previously posted code. Would you mind sharing the code?
Jan Helmdag: Are you familiar with the exact meaning of the sigma(#) option in the trajstart command? And do you recommend a particular value for this option? I would like to generate new random starting values for my estimation with traj to prevent convergence on a local maximum.
Antonino Polizzi: Sigma defines the standard deviation of your starting values. By default, the value is arbitrarily set at .1. You can trace such information using the command viewsource trajstart.ado and there you can see how it is implemented when looking up the local called `sigma'.
Regarding a recommendation on how to set sigma: I don't actually know a definite answer to that. I guess it's what you assume the true population value for the standard deviation actually is.
Jan Helmdag and Yimin Mao Please could you tell an quick approach to test the best fitting in Multigroup Based trajectory modeling ?
Jan Helmdag - am I correct in observing this code works for potential solutions of up to 4 groups? Is there an easy way to expand the range of possible groups to consider here? I tried making the following changes for up to 8 groups, but I cannot get a large enough matsize to make this work...so I'm wondering if I did something wrong along the way. Thanks in advance.
local polynom 8
local maxnum: display `polynom'^1+`polynom'^2+`polynom'^3+`polynom'^4 +`polynom'^5 +`polynom'^6 +`polynom'^7 +`polynom'^8
matrix a = J(`maxnum',8,.)
matrix colnames a = "Groups" "Polynomials" "BIC(N)" "BIC(panels)"
local i 1
local bicn: display -10^99
local bicp: display -10^99
local e(BIC_N_data): display -10^99
local e(BIC_N_subjects): display -10^99
forvalues a = 1/`polynom' {
quietly traj, var(relhw*) indep(tyr*) model(cnorm) min(0) max(2) order(`a')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'
}
local i `++i'
forvalues b = 1/`polynom' {
quietly traj, var(relhw*) indep(tyr*) model(cnorm) min(0) max(2) order(`a' `b')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'`b'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'`b'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'`b'
}
local i `++i'
forvalues c = 1/`polynom' {
quietly traj, var(relhw*) indep(tyr*) model(cnorm) min(0) max(2) order(`a' `b' `c')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'`b'`c'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'`b'`c'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'`b'`c'
}
local i `++i'
forvalues d = 1/`polynom' {
quietly traj, var(relhw*) indep(tyr*) model(cnorm) min(0) max(2) order(`a' `b' `c' `d')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'`b'`c'`d'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'`b'`c'`d'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'`b'`c'`d'
}
local i `++i'
forvalues e = 1/`polynom' {
quietly traj, var(relhw*) indep(tyr*) model(cnorm) min(0) max(2) order(`a' `b' `c' `d' `e')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'`b'`c'`d'`e'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'`b'`c'`d'`e'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'`b'`c'`d'`e'
}
local i `++i'
forvalues f = 1/`polynom' {
quietly traj, var(relhw*) indep(tyr*) model(cnorm) min(0) max(2) order(`a' `b' `c' `d' `e' `f')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'`b'`c'`d'`e'`f'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'`b'`c'`d'`e'`f'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'`b'`c'`d'`e'`f'
}
local i `++i'
forvalues g = 1/`polynom' {
quietly traj, var(relhw*) indep(tyr*) model(cnorm) min(0) max(2) order(`a' `b' `c' `d' `e' `f' `g')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'`b'`c'`d'`e'`f'`g'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'`b'`c'`d'`e'`f'`g'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'`b'`c'`d'`e'`f'`g'
}
local i `++i'
forvalues h = 1/`polynom' {
quietly traj, var(relhw*) indep(tyr*) model(cnorm) min(0) max(2) order(`a' `b' `c' `d' `e' `f' `g' `h')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'`b'`c'`d'`e'`f'`g'`h'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'`b'`c'`d'`e'`f'`g'`h'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'`b'`c'`d'`e'`f'`g'`h'
}
local i `++i'
}
}
}
}
matrix list a
display "Best solution BIC(n): " `solutionn'
display "Best solution BIC(p): " `solutionp'
log close
Tom Kemeny: I think I know what you're aiming at and what might causes the problem. However, the problem you're describing is multi-layered and I have to break it down into four different issues:
Firstly, depending on the version and flavor (i.e. MP/SE/IC) you're using, Stata allows for only a limited matrix size. The standard value is set to 400 to allow for less memory consumption. If you want to change that value, just type in set matsize [value] to get the maximum matrix size. Also look at the help file by typing help set matsize.
Secondly, there is a line in my code that I should have commented on for better transparency. The line where it states “local polynom 4” (in your example local polynom 8) sets the maximum degree of polynomials that is tested to four (i.e. quartic); not the number of groups. The latter is defined by the loops themselves that follow – as you have already figured out.
Thirdly, the traj command allows for a maximum polynomial degree of five (i.e. quintic), at the moment. So the maximum number that can be defined for the local mentioned above can only be set to a maximum of five.
Fourthly, taking all things together, you may want to check how big the matrix will be that you are going to produce. For example, when you want to test for third-degree polynomials and eight groups, you'll end up with 9,840 models. By just changing the degree of polynomials to four (i.e. quartic), you'll end up with 87,380 models. This not just clearly demands more computational power and patience on your side, but also goes beyond the matrix size limit I mentioned earlier.
Ultimately, you may have to think of a different solution if you want the eight-group solution and a higher degree of polynomials. The first thing that comes to my mind is to simply get rid of the idea of displaying all possible solutions at once and may calculate them separately.
Radhouene Doggui: On model selection, there is a helpful chapter in the book Group-Based Modeling of Development by Daniel S. Nagin (2005), Chapter 4 Model selection, pp. 61–77.
I have also written a code on how to calculate the Bayes factor for the example that is provided in the book on pages 68–9.
input groups bic bicb
2 -7460.99 -7466.56
3 -7276.69 -7285.97
4 -7269.63 -7282.63
5 -7274.83 -7291.53
6 -7271.98 -7292.40
7 -7278.08 -7302.21
end
quietly gen double nominator = exp(bicb-bicb[_N])
quietly sum nominator
quietly gen double denominator = `r(sum)'
quietly gen double prob_correct = nominator/denominator
format %3.2f prob_correct
list, clean
Jan Helmdag - thanks for your earlier response. I tried adding an additional potential group, using the code below. I would appreciate any thoughts you might have - I get the following errors:
Warning: variance matrix is nonsymmetric or highly singular
conformability error
My syntax is as follows:
// permutations
*log using traj_permutations, replace
local polynom 4
local maxnum: display `polynom'^1+`polynom'^2+`polynom'^3+`polynom'^4
matrix a = J(`maxnum',4,.)
matrix colnames a = "Groups" "Polynomials" "BIC(N)" "BIC(panels)"
local i 1
local bicn: display -10^99
local bicp: display -10^99
local e(BIC_N_data): display -10^99
local e(BIC_N_subjects): display -10^99
forvalues a = 1/`polynom' {
quietly traj, var(relhw*) indep(tyr*) model(cnorm) min(0) max(2) order(`a')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'
}
local i `++i'
forvalues b = 1/`polynom' {
quietly traj, var(relhw*) indep(tyr*) model(cnorm) min(0) max(2) order(`a' `b')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'`b'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'`b'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'`b'
}
local i `++i'
forvalues c = 1/`polynom' {
quietly traj, var(relhw*) indep(tyr*) model(cnorm) min(0) max(2) order(`a' `b' `c')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'`b'`c'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'`b'`c'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'`b'`c'
}
local i `++i'
forvalues d = 1/`polynom' {
quietly traj, var(relhw*) indep(tyr*) model(cnorm) min(0) max(2) order(`a' `b' `c' `d')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'`b'`c'`d'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'`b'`c'`d'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'`b'`c'`d'
}
local i `++i'
forvalues e = 1/`polynom' {
quietly traj, var(relhw*) indep(tyr*) model(cnorm) min(0) max(2) order(`a' `b' `c' `d' `e')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'`b'`c'`d'`e'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'`b'`c'`d'`e'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'`b'`c'`d'`e'
}
local i `++i'
}
}
}
}
}
matrix list a
display "Best solution BIC(n): " `solutionn'
display "Best solution BIC(p): " `solutionp'
Tom Kemeny: I've rewritten the code and it works with the example dataset I'm using. Now the loops check for 1,364 combinations (5 groups, quartic polynomials). The first lines are now easier to modify and hopefully a bit more straightforward: Now the number of groups and the degree of polynomials are defined as locals and there's also a line where the trajectory model itself is defined (see code in bold). The advantage is that you don't have to modify multiple lines when changing the model this way.
use "http://www.andrew.cmu.edu/user/bjones/traj/data/montreal_sim.dta", clear
*Define groups and degree of polynomials
local polynom 4
local groups 5
forvalues i=1/`groups' {
local maxnum `maxnum'+`polynom'^`i'
}
*Define traj model (except number of groups and order of polynomials
local traj "var(qcp*op) indep(age*) model(cnorm) min(0) max(10)"
*Define matrix and BIC values
matrix a = J(`maxnum',4,.)
matrix colnames a = "Groups" "Polynomials" "BIC(N)" "BIC(panels)"
local i 1
local bicn: display -10^99
local bicp: display -10^99
local e(BIC_N_data): display -10^99
local e(BIC_N_subjects): display -10^99
*Add additional loop here, if necessary
forvalues a = 1/`polynom' {
quietly traj, `traj' order(`a')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'
}
local i `++i'
forvalues b = 1/`polynom' {
quietly traj, `traj' order(`a' `b')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'`b'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'`b'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'`b'
}
local i `++i'
forvalues c = 1/`polynom' {
quietly traj, `traj' order(`a' `b' `c')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'`b'`c'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'`b'`c'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'`b'`c'
}
local i `++i'
forvalues d = 1/`polynom' {
quietly traj, `traj' order(`a' `b' `c' `d')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'`b'`c'`d'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'`b'`c'`d'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'`b'`c'`d'
}
local i `++i'
forvalues e = 1/`polynom' {
quietly traj, `traj' order(`a' `b' `c' `d' `e')
matrix a[`i',1]= e(numGroups1)
matrix a[`i',2]= `a'`b'`c'`d'`e'
matrix a[`i',3]= e(BIC_N_data)
matrix a[`i',4]= e(BIC_n_subjects)
if `e(BIC_N_data)' > `bicn' {
local bicn `e(BIC_N_data)'
local solutionn `a'`b'`c'`d'`e'
}
if `e(BIC_n_subjects)' > `bicp' {
local bicp `e(BIC_n_subjects)'
local solutionp `a'`b'`c'`d'`e'
}
local i `++i'
}
}
}
}
}
matrix list a
display "Best solution BIC(n): " `solutionn'
display "Best solution BIC(p): " `solutionp'
Thanks Jan Helmdag - that works for me quite nicely., I've been going through Nagin's book, which is surprisingly clear and readable. One question I had I'm hoping you might answer: it seems as if everyone's examples feature BIC values that are negative, and fairly close together. But for my current project they are always positive, and there can be some quite large differences between BICs for an additional potential group. So interpreting BIC gets difficult - I cannot, for instance, exponentiate the difference between two potential BICs (as in Nagin's equation 4.2, p.70), as the difference is undefined. Any thoughts?
Tom Kemeny: I am not sure whether I understand the issue you're having correctly, but if it is just the amount of decimal places that's causing you trouble, you have to modify the third-last line of code to display more decimal places in the matrix, e.g.:
. matrix list a, format(%20.0g)
Jan Helmdag - sorry if I was not clear enough. Your code works fine for me. My question is more of a general GBTM question. Across traj models using my data, with any number of groups and polynomial orders, I consistently get (a) BIC values with positive signs; and (b) differences between BIC values (say comparing a 4-group to a 5-group model) that are huge - routinely in the thousands. Below is a little snippet of output from your syntax, to give you an idea:
Groups Polynomials BIC(N) BIC(panels)
r1024 1 4 17771.853 17778.441
r1025 2 41 34507.903 34517.786
r1026 3 411 41004.974 41018.151
r1027 4 4111 40786.743 40803.214
r1028 5 41111 46713.198 46732.964
r1029 5 41112 44219.926 44240.79
r1030 5 41113 46201.607 46223.569
r1031 5 41114 47269.597 47292.657
I'm not sure if the sign on BIC matters - it may not, though I am not certain why my BIC values are positive while essentially all examples in Nagin and elsewhere are all negative.
But the large differences between BICs across models seem to matter in the sense that it makes it impossible to compare model fit using the method Nagin recommends (p.70), as exponentiating a difference between BICs of say 6,000 (the approximate difference between BIC on r1028 and the maximum (r1031)) does not yield an actual number. I could simply assume the largest BIC is the correct one. But if the largest BIC was automatically the right choice, why bother with model fit comparisons like Nagin's equation 4.2:
p_j = e^(bic_j - bic_max)/sigma_j(e^(bic_j - bic_max))
Jan Helmdag Thanks Jan for very helpful answers. I have followed your code, but some models in my data can not calculated the Entropy (Entropy = unknown function /log() ). And when the loop meet that model, it will be stop immediately, so I could not get the results. So Do you have any suggestions for this problem. Thank you!
Jan Helmdag Thank you Jan for your wonderful code. I have to admit I am a fan!
Do you by any chance have a similar code to get the BIC values for a multi trajectory model?
e.g.
traj, multgroups(3) var1(eating*) indep1(age*) model1(cnorm) max1(90) order1(0 1 2) var2(restraint*) indep2(age*) model2(cnorm) max2(210) order2(0 1 2) var3(BMI*) indep3(age*) model3(cnorm) max3(100) order3(0 1 1)
I would appreciate your help very much, as I am quite new in working with stata.
- Badges
- Science topic
- Similar topics
- Mathematics
- Statistics
- Sampling Design
- Sample Size