Good morning to everybody I have 4 variables measured at 3 timepoints: 12 months, 18 months, and 24 months. Is the syntax for choosing the GBTM model correct? I have 106 adults(at least two measurements per outcome).
Censoring limits were defined outcome by outcome in accordance with the observed empirical range, with a small margin beyond the extremes, as no unambiguous theoretical limits were available for these standardized variables.
Based on the data structure and criteria of parsimony, stability, and interpretability, with three time points, the search was limited to polynomial forms of order 0/1, without exploring quadratic terms. Although a quadratic specification is technically possible with three surveys, it is often poorly informative and potentially unstable, especially with multiple outcomes and a small sample size. What do you think? Thanks in advanced to everybody
Code:
****************************************************
* GBTM MULTI-OUTCOME (4 outcomes, cnorm) * FINAL OPERATIONAL VERSION * CORRECT VERSION: pass uses OCC_pp, also checks TotProb and adds diagnostic entropy * Consistent with: Klijn + Nagin multitrajectory + recent review * * LOGIC: * STEP 0 = preliminary univariate exploration of individual outcomes * STEP 1 = choice of K in the multi-outcome model with equal initial order * STEP 2 = fixed K, structured comparison of all plausible 0/1 models * STEP 2B = refit/inspection of finalist models * * CRITERIA: * - BIC = primary criterion * - APPA / OCC_pp / minP / minTotProb / mismatch = adequacy/support criteria * - relative entropy = additional diagnostic of assignment clarity; NOT included in the pass * - absolute number of groups = descriptive; Does NOT qualify * - DELTABIC <= 2 = competing models * - final decision = BIC + parsimony + interpretability + classification diagnostics * - with 3 time points and MAXORDER=1, the possible orders are 0/1
****************************************************
clear all
set more off
set seed 12345
set sortseed 12345
cd "C:\Users\xxxxxxxxxxx\Desktop\LCA_prova"
global DATAFILE "databasex.dta"
global IDVAR "id"
****************************************************
* OUTCOME
****************************************************
global VAR1 "var1_12 var1_18 var1_24"
global VAR2 "var2_12 var2_18 var2_24"
global VAR3 "var3_12 var3_18 var3_24"
global VAR4 "var4_12 var4_18 var4_24"
****************************************************
* CNORM RANGE - OPTIMIZED ON EMPIRICAL DATA * Expanded outward to ensure numerical stability and avoid artificial clipping
****************************************************
global MIN1 -8
global MAX1 6
global MIN2 -9
global MAX2 13
global MIN3 -15
global MAX3 11
global MIN4 -3
global MAX4 9
****************************************************
* RICERCA
****************************************************
global MAXK 2
global MAXORDER 1
global STARTORDER 1
global NREFIT 5
****************************************************
* SOGLIE DI ADEGUATEZZA
****************************************************
global THR_MINP 0.05 // minimum assigned proportion of the group global THR_MINTOTPROB 0.05 // minimum estimated proportion from posterior probabilities global THR_MINAPP 0.70 // minimum average posterior probability global THR_MINOCC 5 // minimum OCC_pp global THR_MAXMIS 0.05 // maximum mismatch
global DELTABIC 2
****************************************************
*
PROGRAM: create times
****************************************************
capture program drop make_time
program define make_time
capture drop t1 t2 t3
gen t1 = 0
gen t2 = 1
gen t3 = 2
end
****************************************************
* PROGRAMMA: statistiche post-traj
****************************************************
capture program drop gbtm_stats
program define gbtm_stats, rclass
syntax , K(integer)
capture drop Mp countG counter APP p n d OCC TotProb mismatch d_pp OCC_pp SD_post __sdtmp
gen double Mp = 0
foreach pr of varlist _traj_ProbG* {
replace Mp = `pr' if `pr' > Mp
}
sort _traj_Group
by _traj_Group: gen countG = _N
by _traj_Group: gen counter = _n
by _traj_Group: egen double APP = mean(Mp)
gen double p = countG/_N
gen double TotProb = .
forvalues gg = 1/`k' {
quietly summarize _traj_ProbG`gg', meanonly
replace TotProb = r(mean) if _traj_Group == `gg'
}
gen double mismatch = abs(TotProb - p)
gen double OCC = .
gen double OCC_pp = .
if `k' > 1 {
gen double n = APP/(1-APP)
gen double d = p/(1-p)
replace OCC = n/d
gen double d_pp = TotProb/(1-TotProb)
replace OCC_pp = n/d_pp
}
else {
replace OCC = 999
replace OCC_pp = 999
}
* PROTEZIONE: Evita crash di Stata se un sottogruppo contiene un solo record (SD non calcolabile)
gen double SD_post = .
forvalues gg = 1/`k' {
capture by _traj_Group: egen double __sdtmp = sd(_traj_ProbG`gg') if _traj_Group == `gg'
if !_rc {
replace SD_post = __sdtmp if _traj_Group == `gg'
drop __sdtmp
}
}
* Relative entropy (0-1)
tempvar __hsum __plnp
local entropy = 1
if `k' > 1 {
gen double `__hsum' = 0
forvalues gg = 1/`k' {
gen double `__plnp' = cond(_traj_ProbG`gg' > 0, _traj_ProbG`gg' * ln(_traj_ProbG`gg'), 0)
replace `__hsum' = `__hsum' + `__plnp'
drop `__plnp'
}
quietly summarize `__hsum', meanonly
local entropy = 1 + (r(sum) / (_N * ln(`k')))
}
preserve
keep if counter == 1
quietly summarize APP, meanonly
local minAPP = r(min)
local meanAPP = r(mean)
quietly summarize p, meanonly
local minP = r(min)
quietly summarize TotProb, meanonly
local minTotProb = r(min)
quietly summarize mismatch, meanonly
local maxMismatch = r(max)
quietly summarize OCC, meanonly
local minOCC = r(min)
quietly summarize OCC_pp, meanonly
local minOCCpp = r(min)
restore
local pass = (`minP' >= $THR_MINP) & ///
(`minTotProb' >= $THR_MINTOTPROB) & ///
(`minAPP' >= $THR_MINAPP) & ///
(`minOCCpp' >= $THR_MINOCC) & ///
(`maxMismatch' <= $THR_MAXMIS)
return scalar minAPP = `minAPP'
return scalar meanAPP = `meanAPP'
return scalar minP = `minP'
return scalar minTotProb = `minTotProb'
return scalar maxMismatch = `maxMismatch'
return scalar entropy = `entropy'
return scalar minOCC = `minOCC'
return scalar minOCCpp = `minOCCpp'
return scalar pass = `pass'
end
****************************************************
* TEMPORARY FILES
****************************************************
tempfile phase0tmp step1tmp step2tmp finalists4 step2ranked
****************************************************
* PHASE 0: PRELIMINARY UNIVARIATE EXPLORATION
****************************************************
tempname h0
capture postclose `h0'
postfile `h0' str8 outcome int K str20 orders ///
double ll aic bic minAPP minOCC minOCCpp minP minTotProb maxMismatch entropy pass ///
using `phase0tmp', replace
forvalues vv = 1/4 {
forvalues k = 1/$MAXK {
use "$DATAFILE", clear
sort $IDVAR, stable
make_time
local indep t1 t2 t3
local oo ""
forvalues g = 1/`k' {
local oo "`oo' $STARTORDER"
}
local oo : list retok oo
quietly capture traj, ///
var(${VAR`vv'}) indep(`indep') order(`oo') model(cnorm) min(${MIN`vv'}) max(${MAX`vv'})
if _rc continue
quietly gbtm_stats, k(`k')
post `h0' ("VAR`vv'") (`k') ("`oo'") (e(ll)) (e(AIC)) (e(BIC_n_subjects)) ///
(r(minAPP)) (r(minOCC)) (r(minOCCpp)) ///
(r(minP)) (r(minTotProb)) (r(maxMismatch)) (r(entropy)) (r(pass))
}
}
postclose `h0'
use `phase0tmp', clear
save phase0_univariate_scan_4var.dta, replace
****************************************************
* STEP 1: choice of K in multi-outcome
****************************************************
tempname h1
capture postclose `h1'
postfile `h1' ///
str5 stage int K str20 o1 str20 o2 str20 o3 str20 o4 ///
int group nG ///
double p TotProb APP OCC OCC_pp mismatch SD_post ///
double ll aic bic minAPP meanAPP minOCC minOCCpp minP minTotProb maxMismatch entropy pass ///
using `step1tmp', replace
forvalues k = 1/$MAXK {
use "$DATAFILE", clear
sort $IDVAR, stable
make_time
local indep t1 t2 t3
local o1 ""
local o2 ""
local o3 ""
local o4 ""
forvalues g = 1/`k' {
local o1 "`o1' $STARTORDER"
local o2 "`o2' $STARTORDER"
local o3 "`o3' $STARTORDER"
local o4 "`o4' $STARTORDER"
}
local o1 : list retok o1
local o2 : list retok o2
local o3 : list retok o3
local o4 : list retok o4
quietly capture traj, multgroups(`k') ///
var1($VAR1) indep1(`indep') order1(`o1') model1(cnorm) min1($MIN1) max1($MAX1) ///
var2($VAR2) indep2(`indep') order2(`o2') model2(cnorm) min2($MIN2) max2($MAX2) ///
var3($VAR3) indep3(`indep') order3(`o3') model3(cnorm) min3($MIN3) max3($MAX3) ///
var4($VAR4) indep4(`indep') order4(`o4') model4(cnorm) min4($MIN4) max4($MAX4)
if _rc continue
quietly gbtm_stats, k(`k')
local ll = e(ll)
local aic = e(AIC)
local bic = e(BIC_n_subjects)
local minAPP = r(minAPP)
local meanAPP = r(meanAPP)
local minOCC = r(minOCC)
local minOCCpp = r(minOCCpp)
local maxMismatch = r(maxMismatch)
local entropy = r(entropy)
local minP = r(minP)
local minTotProb = r(minTotProb)
local pass = r(pass)
forvalues gg = 1/`k' {
quietly summarize countG if _traj_Group == `gg', meanonly
local nG = r(mean)
quietly summarize p if _traj_Group == `gg', meanonly
local pg = r(mean)
quietly summarize TotProb if _traj_Group == `gg', meanonly
local tpg = r(mean)
quietly summarize APP if _traj_Group == `gg', meanonly
local appg = r(mean)
quietly summarize OCC if _traj_Group == `gg', meanonly
local occg = r(mean)
quietly summarize OCC_pp if _traj_Group == `gg', meanonly
local occppg = r(mean)
quietly summarize mismatch if _traj_Group == `gg', meanonly
local misg = r(mean)
local sdg = .
quietly count if _traj_Group == `gg'
if r(N) > 1 {
quietly summarize SD_post if _traj_Group == `gg', meanonly
local sdg = r(mean)
}
post `h1' ("STEP1") (`k') ("`o1'") ("`o2'") ("`o3'") ("`o4'") ///
(`gg') (`nG') (`pg') (`tpg') (`appg') (`occg') (`occppg') (`misg') (`sdg') ///
(`ll') (`aic') (`bic') (`minAPP') (`meanAPP') (`minOCC') (`minOCCpp') ///
(`minP') (`minTotProb') (`maxMismatch') (`entropy') (`pass')
}
}
postclose `h1'
use `step1tmp', clear
egen byte tagmodel = tag(K o1 o2 o3 o4)
keep if tagmodel
drop tagmodel
save step1_kselection_4var.dta, replace
gsort -pass -bic
count if pass == 1
if r(N) > 0 {
keep if pass == 1
gsort -bic
}
else {
gsort -bic
}
quietly summarize K in 1, meanonly
local BESTK = r(min)
di as result "K selezionato = `BESTK'"
****************************************************
* STEP 2: STRUCTURED SEARCH (Safe Combinatorial Logic)
****************************************************
tempname h2
capture postclose `h2'
postfile `h2' ///
str5 stage int K str20 o1 str20 o2 str20 o3 str20 o4 ///
int group nG ///
double p TotProb APP OCC OCC_pp mismatch SD_post ///
double ll aic bic minAPP meanAPP minOCC minOCCpp minP minTotProb maxMismatch entropy pass ///
using `step2tmp', replace
local k = `BESTK'
local base = $MAXORDER + 1
local ncomb = `base'^`k'
forvalues i1 = 1/`ncomb' {
local o1 ""
forvalues g = 1/`k' {
local div = `base'^(`k' - `g')
local digit = mod(int((`i1' - 1)/`div'), `base')
local o1 "`o1' `digit'"
}
local o1 : list retok o1
forvalues i2 = 1/`ncomb' {
local o2 ""
forvalues g = 1/`k' {
local div = `base'^(`k' - `g')
local digit = mod(int((`i2' - 1)/`div'), `base')
local o2 "`o2' `digit'"
}
local o2 : list retok o2
forvalues i3 = 1/`ncomb' {
local o3 ""
forvalues g = 1/`k' {
local div = `base'^(`k' - `g')
local digit = mod(int((`i3' - 1)/`div'), `base')
local o3 "`o3' `digit'"
}
local o3 : list retok o3
forvalues i4 = 1/`ncomb' {
local o4 ""
forvalues g = 1/`k' {
local div = `base'^(`k' - `g')
local digit = mod(int((`i4' - 1)/`div'), `base')
local o4 "`o4' `digit'"
}
local o4 : list retok o4
use "$DATAFILE", clear
sort $IDVAR, stable
make_time
local indep t1 t2 t3
quietly capture traj, multgroups(`k') ///
var1($VAR1) indep1(`indep') order1(`o1') model1(cnorm) min1($MIN1) max1($MAX1) ///
var2($VAR2) indep2(`indep') order2(`o2') model2(cnorm) min2($MIN2) max2($MAX2) ///
var3($VAR3) indep3(`indep') order3(`o3') model3(cnorm) min3($MIN3) max3($MAX3) ///
var4($VAR4) indep4(`indep') order4(`o4') model4(cnorm) min4($MIN4) max4($MAX4)
if _rc continue
quietly gbtm_stats, k(`k')
local ll = e(ll)
local aic = e(AIC)
local bic = e(BIC_n_subjects)
local minAPP = r(minAPP)
local meanAPP = r(meanAPP)
local minOCC = r(minOCC)
local minOCCpp = r(minOCCpp)
local maxMismatch = r(maxMismatch)
local entropy = r(entropy)
local minP = r(minP)
local minTotProb = r(minTotProb)
local pass = r(pass)
forvalues gg = 1/`k' {
quietly summarize countG if _traj_Group == `gg', meanonly
local nG = r(mean)
quietly summarize p if _traj_Group == `gg', meanonly
local pg = r(mean)
quietly summarize TotProb if _traj_Group == `gg', meanonly
local tpg = r(mean)
quietly summarize APP if _traj_Group == `gg', meanonly
local appg = r(mean)
quietly summarize OCC if _traj_Group == `gg', meanonly
local occg = r(mean)
quietly summarize OCC_pp if _traj_Group == `gg', meanonly
local occppg = r(mean)
quietly summarize mismatch if _traj_Group == `gg', meanonly
local misg = r(mean)
local sdg = .
quietly count if _traj_Group == `gg'
if r(N) > 1 {
quietly summarize SD_post if _traj_Group == `gg', meanonly
local sdg = r(mean)
}
post `h2' ("STEP2") (`k') ("`o1'") ("`o2'") ("`o3'") ("`o4'") ///
(`gg') (`nG') (`pg') (`tpg') (`appg') (`occg') (`occppg') (`misg') (`sdg') ///
(`ll') (`aic') (`bic') (`minAPP') (`meanAPP') (`minOCC') (`minOCCpp') ///
(`minP') (`minTotProb') (`maxMismatch') (`entropy') (`pass')
}
}
}
}
}
postclose `h2'
use `step2tmp', clear
save step2_models_4var.dta, replace
egen byte tagmodel = tag(K o1 o2 o3 o4)
keep if tagmodel
drop tagmodel
keep if K == `BESTK'
count if pass == 1
if r(N) > 0 {
keep if pass == 1
}
gsort -bic
quietly summarize bic, meanonly
local bestbic = r(max)
keep if bic >= (`bestbic' - $DELTABIC)
gsort -bic -minAPP -minOCCpp maxMismatch
gen rank_finalista = _n
save `step2ranked', replace
save finalists_step2_4var.dta, replace
count
local NFINAL = r(N)
di as result "Numero modelli finalisti entro DeltaBIC = `NFINAL'"
list rank_finalista K o1 o2 o3 o4 bic minAPP minOCCpp maxMismatch entropy minP minTotProb pass, noobs
****************************************************
* STEP 2B: refit of the finalist models
****************************************************
local NINSPECT = cond(`NFINAL' < $NREFIT, `NFINAL', $NREFIT)
forvalues i = 1/`NINSPECT' {
use finalists_step2_4var.dta, clear
local CK = K[`i']
local CO1 = o1[`i']
local CO2 = o2[`i']
local CO3 = o3[`i']
local CO4 = o4[`i']
capture log close candlog
log using "candidate4_`i'_K`CK'.smcl", replace name(candlog)
use "$DATAFILE", clear
sort $IDVAR, stable
make_time
local indep t1 t2 t3
traj, multgroups(`CK') ///
var1($VAR1) indep1(`indep') order1(`CO1') model1(cnorm) min1($MIN1) max1($MAX1) ///
var2($VAR2) indep2(`indep') order2(`CO2') model2(cnorm) min2($MIN2) max2($MAX2) ///
var3($VAR3) indep3(`indep') order3(`CO3') model3(cnorm) min3($MIN3) max3($MAX3) ///
var4($VAR4) indep4(`indep') order4(`CO4') model4(cnorm) min4($MIN4) max4($MAX4)
di as result "BIC = " e(BIC_n_subjects)
di as result "AIC = " e(AIC)
di as result "LL = " e(ll)
log close candlog
}
****************************************************
* LEAD CANDIDATE ACCORDING TO PRE-SPECIFIED CRITERIA
****************************************************
use `step2ranked', clear
gen byte _pick = (_n == 1)
quietly summarize K if _pick, meanonly
local FK = r(min)
levelsof o1 if _pick, local(FO1) clean
levelsof o2 if _pick, local(FO2) clean
levelsof o3 if _pick, local(FO3) clean
levelsof o4 if _pick, local(FO4) clean
drop _pick
di as result "CANDIDATO PRINCIPALE:"
di as result "K = `FK'"
di as result "order1 = `FO1'"
di as result "order2 = `FO2'"
di as result "order3 = `FO3'"
di as result "order4 = `FO4'"
use "$DATAFILE", clear
sort $IDVAR, stable
make_time
local indep t1 t2 t3
traj, multgroups(`FK') ///
var1($VAR1) indep1(`indep') order1(`FO1') model1(cnorm) min1($MIN1) max1($MAX1) ///
var2($VAR2) indep2(`indep') order2(`FO2') model2(cnorm) min2($MIN2) max2($MAX2) ///
var3($VAR3) indep3(`indep') order3(`FO3') model3(cnorm) min3($MIN3) max3($MAX3) ///
var4($VAR4) indep4(`indep') order4(`FO4') model4(cnorm) min4($MIN4) max4($MAX4)
di as result "BIC finale = " e(BIC_n_subjects)
di as result "AIC finale = " e(AIC)
di as result "LL finale = " e(ll)
****************************************************
* FINAL STATISTICS OF THE SELECTED MODEL
****************************************************
quietly gbtm_stats, k(`FK')
di as result "minAPP finale = " r(minAPP)
di as result "meanAPP finale = " r(meanAPP)
di as result "minP finale = " r(minP)
di as result "minTotProb finale = " r(minTotProb)
di as result "minOCC finale = " r(minOCC)
di as result "minOCCpp finale = " r(minOCCpp)
di as result "maxMismatch finale = " r(maxMismatch)
di as result "entropy finale = " r(entropy)
di as result "pass finale = " r(pass)
Comment