st: Difficult ML

 From "miguel foguel" To statalist@hsphsun2.harvard.edu Subject st: Difficult ML Date Mon, 01 Mar 2004 23:01:33 +0000

Dear all,

I'm trying to estimate a competing risks model (2 risks) based on a MPH specification for each risk. Time is discrete and the observation window is 13 months. The baseline hazards are assumed piece-wise constant, so there are 13 parameters for each risk. The systematic part is modelled as exp(.), and each risk has its own continuous covariate (cov1 and cov2, respectively). As for the joint distribution of the unobservables I'm assuming two points of support for each unobservable, which gives 7 extra parameters to estimated (4 points of support + 3 probabilities). In fact, the probabilities have been transformed into a multinomial framework (see program below).

I'm using -ml- routine with method d0. The problem is that maximization is converging to a flat part of the LF. I have used -ml check- and -ml search-, and also tried interactivelly several different initial values. From -ml plot- I noticed that the parameters of the 13th dummy associated with the piece-wise baseline hazard for each risk become flat very quickly and stay there for all values I experimented with. I then normalized those parameters to be zero because I guessed that there may be a 'collinearity' problem with the means of the unobservables.

Below there is copy of the code I'm using to estimate the model (For readability I'm also attaching the code in a txt file). There are 3 components to the LF: contributions from the 1st risk (dc==1), from the 2nd risk (dc==2), and from right-censoring (dc==3). Though it is a discrete time setting with piece-wise constant baseline hazards, I arranged the data set so that each row contains only one individual. Apart from the right-censoring case, the contributions are essentially composed by: the hazard of risk j times the survivor of risk j times the survivor of risk k, where j is different from k. Each contribution is obtained by using the correspondent probability of the unobservables.

cap program drop comprisk
program define comprisk
version 8.0
args todo b lnf

tempvar lnfj
/* Parameters for the baseline hazards: ht* 1st risk; he* 2nd risk */
tempname he1 he2 he3 he4 he5 he6 he7 he8 he9 he10 he11 he12
tempname ht1 ht2 ht3 ht4 ht5 ht6 ht7 ht8 ht9 ht10 ht11 ht12
/* Parameters for mass points of the joint distributions of unobservables */
tempname v1 v2 v3 v4
/* Parameters for the probabilities of the joint distribution */
tempname q13 q14 q23 p13 p14 p23 p24
/* Parameters for the systematic part */
tempname b1 b2

/* Normalizing one of the parameters of the baseline hazards */
local he13 = 0
local ht13 = 0

mleval `he1' = `b', eq(1) scalar
mleval `he2' = `b', eq(2) scalar
mleval `he3' = `b', eq(3) scalar
mleval `he4' = `b', eq(4) scalar
mleval `he5' = `b', eq(5) scalar
mleval `he6' = `b', eq(6) scalar
mleval `he7' = `b', eq(7) scalar
mleval `he8' = `b', eq(8) scalar
mleval `he9' = `b', eq(9) scalar
mleval `he10' = `b', eq(10) scalar
mleval `he11' = `b', eq(11) scalar
mleval `he12' = `b', eq(12) scalar
mleval `ht1' = `b', eq(13) scalar
mleval `ht2' = `b', eq(14) scalar
mleval `ht3' = `b', eq(15) scalar
mleval `ht4' = `b', eq(16) scalar
mleval `ht5' = `b', eq(17) scalar
mleval `ht6' = `b', eq(18) scalar
mleval `ht7' = `b', eq(19) scalar
mleval `ht8' = `b', eq(20) scalar
mleval `ht9' = `b', eq(21) scalar
mleval `ht10' = `b', eq(22) scalar
mleval `ht11' = `b', eq(23) scalar
mleval `ht12' = `b', eq(24) scalar
mleval `v1' = `b', eq(25) scalar
mleval `v2' = `b', eq(26) scalar
mleval `v3' = `b', eq(27) scalar
mleval `v4' = `b', eq(28) scalar
mleval `q13' = `b', eq(29) scalar
mleval `q14' = `b', eq(30) scalar
mleval `q23' = `b', eq(31) scalar
mleval `b1' = `b', eq(32) scalar
mleval `b2' = `b', eq(33) scalar

/* Defining X'b for the systematic part */
tempname X1 X2 H1 S1 H2 S2
qui gen double `X1' = cov1 * `b1'
qui gen double `X2' = cov2 * `b2'

/* Defining the "X'b" for the hazard (H1) and survivor (S1) of 1st risk: demp* and surem* are dummies */
qui gen double `H1' = demp1 * `he1' + demp2 * `he2' + demp3 * `he3' + demp4 * `he4' + demp5 * `he5' + demp6 * `he6' + demp7 * `he7' + demp8 * `he8' + demp9 * `he9' + demp10 * `he10' + demp11 * `he11' + demp12 * `he12' + demp13 * `he13'
qui gen double `S1' = surem1 * `he1' + surem2 * `he2' + surem3 * `he3' + surem4 * `he4' + surem5 * `he5' + surem6 * `he6' + surem7 * `he7' + surem8 * `he8' + surem9 * `he9' + surem10 * `he10' + surem11 * `he11' + surem12 * `he12'

/* Defining the "X'b" for the hazard (H2) and survivor (S2) of 2nd risk: dtr* and surtr* are dummies */
qui gen double `H2' = dtr1 * `ht1' + dtr2 * `ht2' + dtr3 * `ht3' + dtr4 * `ht4' + dtr5 * `ht5' + dtr6 * `ht6' + dtr7 * `ht7' + dtr8 * `ht8' + dtr9 * `ht9' + dtr10 * `ht10' + dtr11 * `ht11' + dtr12 * `ht12' + dtr13 * `ht13'
qui gen double `S2' = surtr1 * `ht1' + surtr2 * `ht2' + surtr3 * `ht3' + surtr4 * `ht4' + surtr5 * `ht5' + surtr6 * `ht6' + surtr7 * `ht7' + surtr8 * `ht8' + surtr9 * `ht9' + surtr10 * `ht10' + surtr11 * `ht11' + surtr12 * `ht12'

/* Transforming probs into unrestricted multinomial probs */

qui {
scalar `p13' = exp(`q13') / (exp(`q13') + exp(`q14') + exp(`q23') + 1)
scalar `p14' = exp(`q14') / (exp(`q13') + exp(`q14') + exp(`q23') + 1)
scalar `p23' = exp(`q23') / (exp(`q13') + exp(`q14') + exp(`q23') + 1)
scalar `p24' = 1 / (exp(`q13') + exp(`q14') + exp(`q23') + 1)
}

/* Contribution to LF of 1st risk */
qui replace lnfj = ln( (`p13' * (1 - exp(-exp(`X1') * exp(`v1') * `H1')) * (exp(-exp(`X1') * exp(`v1') * `S1')) * (exp(-exp(`X2') * exp(`v3') * `S2')) ) + /*
*/ (`p14' * (1 - exp(-exp(`X1') * exp(`v2') * `H1')) * (exp(-exp(`X1') * exp(`v2') * `S1')) * (exp(-exp(`X2') * exp(`v3') * `S2')) ) + /*
*/ (`p23' * (1 - exp(-exp(`X1') * exp(`v1') * `H1')) * (exp(-exp(`X1') * exp(`v1') * `S1')) * (exp(-exp(`X2') * exp(`v4') * `S2')) ) + /* */ (`p24' * (1 - exp(-exp(`X1') * exp(`v2') * `H1')) * (exp(-exp(`X1') * exp(`v2') * `S1')) * (exp(-exp(`X2') * exp(`v4') * `S2')) ) ) /* */ if dc==1

/* Contribution to LF of 2nd risk */
qui gen double lnfj = ln( (`p13' * (1 - exp(-exp(`X2') * exp(`v3') * `H2')) * (exp(-exp(`X2') * exp(`v3') * `S2')) * (exp(-exp(`X1') * exp(`v1') * `S1')) ) + /*
*/ (`p14' * (1 - exp(-exp(`X2') * exp(`v3') * `H2')) * (exp(-exp(`X2') * exp(`v3') * `S2')) * (exp(-exp(`X1') * exp(`v2') * `S1')) ) + /*
*/ (`p23' * (1 - exp(-exp(`X2') * exp(`v4') * `H2')) * (exp(-exp(`X2') * exp(`v4') * `S2')) * (exp(-exp(`X1') * exp(`v1') * `S1')) ) + /*
*/ (`p24' * (1 - exp(-exp(`X2') * exp(`v4') * `H2')) * (exp(-exp(`X2') * exp(`v4') * `S2')) * (exp(-exp(`X1') * exp(`v2') * `S1')) ) ) /*
*/ if dc==2

/* Contribution to LF of right-censoring */
qui replace lnfj = ln( (`p13' * (exp(-exp(`X1') * exp(`v1') * `S1')) * (exp(-exp(`X2') * exp(`v3') * `S2')) ) + /*
*/ (`p14' * (exp(-exp(`X1') * exp(`v2') * `S1')) * (exp(-exp(`X2') * exp(`v3') * `S2')) ) + /*
*/ (`p23' * (exp(-exp(`X1') * exp(`v1') * `S1')) * (exp(-exp(`X2') * exp(`v4') * `S2')) ) + /*
*/ (`p24' * (exp(-exp(`X1') * exp(`v2') * `S1')) * (exp(-exp(`X2') * exp(`v4') * `S2')) ) ) /*
*/ if dc==3

mlsum `lnf' = `lnfj'

end

ml model d0 comprisk /he1 /he2 /he3 /he4 /he5 /he6 /he7 /he8 /he9 /he10 /he11 /he12 /ht1 /ht2 /ht3 /ht4 /ht5 /ht6 /ht7 /ht8 /ht9 /ht10 /ht11 /ht12 /v1 /v2 /v3 /v4 /q13 /q14 /q23 /b1 /b2

Any comments will be highly appreciated. I thank you very much in advance.

Best regards,

Miguel Foguel

_________________________________________________________________

cap program drop comprisk
program define comprisk
version 8.0
args todo b lnf

tempvar lnfj

/* Parameters for the baseline hazards: ht* 1st risk; he* 2nd risk */
tempname he1 he2 he3 he4 he5 he6 he7 he8 he9 he10 he11 he12
tempname ht1 ht2 ht3 ht4 ht5 ht6 ht7 ht8 ht9 ht10 ht11 ht12

/* Parameters for mass points of the joint distributions of unobservables */
tempname v1 v2 v3 v4

/* Parameters for the probabilities of the joint distribution */
tempname q13 q14 q23 p13 p14 p23 p24

/* Parameters for the systematic part */
tempname b1 b2

/* Normalizing one of the parameters of the baseline hazards */
local he13 = 0
local ht13 = 0

mleval `he1' = `b', eq(1) scalar
mleval `he2' = `b', eq(2) scalar
mleval `he3' = `b', eq(3) scalar
mleval `he4' = `b', eq(4) scalar
mleval `he5' = `b', eq(5) scalar
mleval `he6' = `b', eq(6) scalar
mleval `he7' = `b', eq(7) scalar
mleval `he8' = `b', eq(8) scalar
mleval `he9' = `b', eq(9) scalar
mleval `he10' = `b', eq(10) scalar
mleval `he11' = `b', eq(11) scalar
mleval `he12' = `b', eq(12) scalar
mleval `ht1' = `b', eq(13) scalar
mleval `ht2' = `b', eq(14) scalar
mleval `ht3' = `b', eq(15) scalar
mleval `ht4' = `b', eq(16) scalar
mleval `ht5' = `b', eq(17) scalar
mleval `ht6' = `b', eq(18) scalar
mleval `ht7' = `b', eq(19) scalar
mleval `ht8' = `b', eq(20) scalar
mleval `ht9' = `b', eq(21) scalar
mleval `ht10' = `b', eq(22) scalar
mleval `ht11' = `b', eq(23) scalar
mleval `ht12' = `b', eq(24) scalar
mleval `v1' = `b', eq(25) scalar
mleval `v2' = `b', eq(26) scalar
mleval `v3' = `b', eq(27) scalar
mleval `v4' = `b', eq(28) scalar
mleval `q13' = `b', eq(29) scalar
mleval `q14' = `b', eq(30) scalar
mleval `q23' = `b', eq(31) scalar
mleval `b1' = `b', eq(32) scalar
mleval `b2' = `b', eq(33) scalar

/* Defining X'b for the systematic part */
tempname X1 X2 H1 S1 H2 S2
qui gen double `X1' = cov1 * `b1'
qui gen double `X2' = cov2 * `b2'

/* Defining the "X'b" for the hazard (H1) and survivor (S1) of 1st risk: demp* and surem* are dummies */
qui gen double `H1' = demp1 * `he1' + demp2 * `he2' + demp3 * `he3' + demp4 * `he4' + demp5 * `he5' + demp6 * `he6' + demp7 * `he7' + demp8 * `he8' + demp9 * `he9' + demp10 * `he10' + demp11 * `he11' + demp12 * `he12' + demp13 * `he13'
qui gen double `S1' = surem1 * `he1' + surem2 * `he2' + surem3 * `he3' + surem4 * `he4' + surem5 * `he5' + surem6 * `he6' + surem7 * `he7' + surem8 * `he8' + surem9 * `he9' + surem10 * `he10' + surem11 * `he11' + surem12 * `he12'

/* Defining the "X'b" for the hazard (H2) and survivor (S2) of 2nd risk: dtr* and surtr* are dummies */
qui gen double `H2' = dtr1 * `ht1' + dtr2 * `ht2' + dtr3 * `ht3' + dtr4 * `ht4' + dtr5 * `ht5' + dtr6 * `ht6' + dtr7 * `ht7' + dtr8 * `ht8' + dtr9 * `ht9' + dtr10 * `ht10' + dtr11 * `ht11' + dtr12 * `ht12' + dtr13 * `ht13'
qui gen double `S2' = surtr1 * `ht1' + surtr2 * `ht2' + surtr3 * `ht3' + surtr4 * `ht4' + surtr5 * `ht5' + surtr6 * `ht6' + surtr7 * `ht7' + surtr8 * `ht8' + surtr9 * `ht9' + surtr10 * `ht10' + surtr11 * `ht11' + surtr12 * `ht12'

/* Transforming probs into unrestricted multinomial probs */

qui {
scalar `p13' = exp(`q13') / (exp(`q13') + exp(`q14') + exp(`q23') + 1)
scalar `p14' = exp(`q14') / (exp(`q13') + exp(`q14') + exp(`q23') + 1)
scalar `p23' = exp(`q23') / (exp(`q13') + exp(`q14') + exp(`q23') + 1)
scalar `p24' = 1 / (exp(`q13') + exp(`q14') + exp(`q23') + 1)
}

/* Contribution to LF of 1st risk */
qui replace lnfj = ln( (`p13' * (1 - exp(-exp(`X1') * exp(`v1') * `H1')) * (exp(-exp(`X1') * exp(`v1') * `S1')) * (exp(-exp(`X2') * exp(`v3') * `S2')) ) + /*
*/ (`p14' * (1 - exp(-exp(`X1') * exp(`v2') * `H1')) * (exp(-exp(`X1') * exp(`v2') * `S1')) * (exp(-exp(`X2') * exp(`v3') * `S2')) ) + /*
*/ (`p23' * (1 - exp(-exp(`X1') * exp(`v1') * `H1')) * (exp(-exp(`X1') * exp(`v1') * `S1')) * (exp(-exp(`X2') * exp(`v4') * `S2')) ) + /*
*/ (`p24' * (1 - exp(-exp(`X1') * exp(`v2') * `H1')) * (exp(-exp(`X1') * exp(`v2') * `S1')) * (exp(-exp(`X2') * exp(`v4') * `S2')) ) ) /*
*/ if dc==1

/* Contribution to LF of 2nd risk */
qui gen double lnfj = ln( (`p13' * (1 - exp(-exp(`X2') * exp(`v3') * `H2')) * (exp(-exp(`X2') * exp(`v3') * `S2')) * (exp(-exp(`X1') * exp(`v1') * `S1')) ) + /*
*/ (`p14' * (1 - exp(-exp(`X2') * exp(`v3') * `H2')) * (exp(-exp(`X2') * exp(`v3') * `S2')) * (exp(-exp(`X1') * exp(`v2') * `S1')) ) + /*
*/ (`p23' * (1 - exp(-exp(`X2') * exp(`v4') * `H2')) * (exp(-exp(`X2') * exp(`v4') * `S2')) * (exp(-exp(`X1') * exp(`v1') * `S1')) ) + /*
*/ (`p24' * (1 - exp(-exp(`X2') * exp(`v4') * `H2')) * (exp(-exp(`X2') * exp(`v4') * `S2')) * (exp(-exp(`X1') * exp(`v2') * `S1')) ) ) /*
*/ if dc==2

/* Contribution to LF of right-censoring */
qui replace lnfj = ln( (`p13' * (exp(-exp(`X1') * exp(`v1') * `S1')) * (exp(-exp(`X2') * exp(`v3') * `S2')) ) + /*
*/ (`p14' * (exp(-exp(`X1') * exp(`v2') * `S1')) * (exp(-exp(`X2') * exp(`v3') * `S2')) ) + /*
*/ (`p23' * (exp(-exp(`X1') * exp(`v1') * `S1')) * (exp(-exp(`X2') * exp(`v4') * `S2')) ) + /*
*/ (`p24' * (exp(-exp(`X1') * exp(`v2') * `S1')) * (exp(-exp(`X2') * exp(`v4') * `S2')) ) ) /*
*/ if dc==3

mlsum `lnf' = `lnfj'

end

ml model d0 comprisk /he1 /he2 /he3 /he4 /he5 /he6 /he7 /he8 /he9 /he10 /he11 /he12 /ht1 /ht2 /ht3 /ht4 /ht5 /ht6 /ht7 /ht8 /ht9 /ht10 /ht11 /ht12 /v1 /v2 /v3 /v4 /q13 /q14 /q23 /b1 /b2