[Joseph Coveney]

What did I wrong ... I did define the scalar, didn't I? In a wrong way?
Did I refer wrong in the scalar definition?

Yes.

You cannot define a temporary name as

tempname div_tasks div_tasks_se

and then try to refer to it as

scalar define `div_tasks_1' = _b[1.div_tasks]
scalar define `div_tasks_1_se' = _se[1.div_tasks]

[Ariane Arbol]

Ah, thank you, now I got it! And it worked! Thanks for your help and patience!

[Joseph Coveney]

You're welcome, and thank you for the closure.

[Ariane Arbol]

Hello again,

now there's a further question, I'm kind of at a loss here: The code now works fine (see above), but the estimated confidence intervals do not include the estimated coefficients. In fact, it looks like this:

This is the code:

program define chi2em, rclass
    args x cmdline ll0 coef level
    
    constraint define 1 _b[`coef'] = `x'
    
    local command = scalar(`cmdline')
    if strpos("`command'", ",") > 0 `command' constraints(1)
    else `command', constraints(1)
    
    if "`level'" == "" local level `c(level)'
    
    tempname X
    scalar define `X' = 2 * abs(`ll0' - e(ll))
    return scalar fx = (`X' - invchi2(1, `level'/100))^2
end

firthlogit care_benefits c.health_lim i.female_r c.age_r i.educ_r_cat c.hh_income_1000 c.assets_htsd i.div_tasks i.freq_help_all        ///
    if valid==1 & onecareperson==1    
tempname ll0 cmdline health_lim health_lim_se female_r female_r_se age_r age_r_se educ_r_cat educ_r_cat_se hh_income_1000 hh_income_1000_se assets_htsd assets_htsd_se div_tasks_1 div_tasks_1_se div_tasks_2 div_tasks_2_se div_tasks_3 div_tasks_3_se freq_help_all_1 freq_help_all_1_se freq_help_all_2 freq_help_all_2_se freq_help_all_3 freq_help_all_3_se freq_help_all_4 freq_help_all_4_se cons cons_se

scalar define `ll0' = e(ll)
scalar define `cmdline' = "firthlogit care_benefits c.health_lim i.female_r c.age_r i.educ_r_cat c.hh_income_1000 c.assets_htsd i.div_tasks i.freq_help_all    if valid==1 & onecareperson==1"

scalar define `health_lim' = _b[health_lim]
scalar define `health_lim_se' = _se[health_lim]

scalar define `female_r' = _b[1.female_r]
scalar define `female_r_se' = _se[1.female_r]

scalar define `age_r' = _b[age_r]
scalar define `age_r_se' = _se[age_r]

scalar define `educ_r_cat' = _b[1.educ_r_cat]
scalar define `educ_r_cat_se' = _se[1.educ_r_cat]

scalar define `hh_income_1000' = _b[hh_income_1000]
scalar define `hh_income_1000_se' = _se[hh_income_1000]

scalar define `assets_htsd' = _b[assets_htsd]
scalar define `assets_htsd_se' = _se[assets_htsd]

scalar define `div_tasks_1' = _b[1.div_tasks]
scalar define `div_tasks_1_se' = _se[1.div_tasks]

scalar define `div_tasks_2' = _b[2.div_tasks]
scalar define `div_tasks_2_se' = _se[2.div_tasks]

scalar define `div_tasks_3' = _b[3.div_tasks]
scalar define `div_tasks_3_se' = _se[3.div_tasks]

scalar define `freq_help_all_1' = _b[1.freq_help_all]
scalar define `freq_help_all_1_se' = _se[1.freq_help_all]

scalar define `freq_help_all_2' = _b[2.freq_help_all]
scalar define `freq_help_all_2_se' = _se[2.freq_help_all]

scalar define `freq_help_all_3' = _b[3.freq_help_all]
scalar define `freq_help_all_3_se' = _se[3.freq_help_all]

scalar define `freq_help_all_4' = _b[4.freq_help_all]
scalar define `freq_help_all_4_se' = _se[4.freq_help_all]

scalar define `cons' = _b[_cons]
scalar define `cons_se' = _se[_cons]


* health_lim
// Lower
minbound chi2em, range(`=`health_lim' - 3 * `health_lim_se'', `=`health_lim'') arguments(`cmdline' `ll0' health_lim 95)
display in smcl as text exp(r(x))
// Upper
minbound chi2em, range(`=`health_lim'', `=`health_lim' + 3 * `health_lim_se'') arguments(`cmdline' `ll0' health_lim 95)
display in smcl as text exp(r(x))

* female_r
// Lower
minbound chi2em, range(`=`female_r' - 3 * `female_r_se'', `=`female_r'') arguments(`cmdline' `ll0' 1.female_r 95)
display in smcl as text exp(r(x))
// Upper
minbound chi2em, range(`=`female_r'', `=`female_r' + 3 * `female_r_se'') arguments(`cmdline' `ll0' 1.female_r 95)
display in smcl as text exp(r(x))

* age_r
// Lower
minbound chi2em, range(`=`age_r' - 3 * `age_r_se'', `=`age_r'') arguments(`cmdline' `ll0' age_r 95)
display in smcl as text exp(r(x))
// Upper
minbound chi2em, range(`=`age_r'', `=`age_r' + 3 * `age_r_se'') arguments(`cmdline' `ll0' age_r 95)
display in smcl as text exp(r(x))

* educ_r_cat
// Lower
minbound chi2em, range(`=`educ_r_cat' - 3 * `educ_r_cat_se'', `=`educ_r_cat'') arguments(`cmdline' `ll0' 1.educ_r_cat 95)
display in smcl as text exp(r(x))
// Upper
minbound chi2em, range(`=`educ_r_cat'', `=`educ_r_cat' + 3 * `educ_r_cat_se'') arguments(`cmdline' `ll0' 1.educ_r_cat 95)
display in smcl as text exp(r(x))

* hh_income_1000
// Lower
minbound chi2em, range(`=`hh_income_1000' - 3 * `hh_income_1000_se'', `=`hh_income_1000'') arguments(`cmdline' `ll0' hh_income_1000 95)
display in smcl as text exp(r(x))
// Upper
minbound chi2em, range(`=`hh_income_1000'', `=`hh_income_1000' + 3 * `hh_income_1000_se'') arguments(`cmdline' `ll0' hh_income_1000 95)
display in smcl as text exp(r(x))

* assets_htsd
// Lower
minbound chi2em, range(`=`assets_htsd' - 3 * `assets_htsd_se'', `=`assets_htsd'') arguments(`cmdline' `ll0' assets_htsd 95)
display in smcl as text exp(r(x))
// Upper
minbound chi2em, range(`=`assets_htsd'', `=`assets_htsd' + 3 * `assets_htsd_se'') arguments(`cmdline' `ll0' assets_htsd 95)
display in smcl as text exp(r(x))

* div_tasks_1
// Lower
minbound chi2em, range(`=`div_tasks_1' - 3 * `div_tasks_1_se'', `=`div_tasks_1'') arguments(`cmdline' `ll0' 1.div_tasks 95)
display in smcl as text exp(r(x))
// Upper
minbound chi2em, range(`=`div_tasks_1'', `=`div_tasks_1' + 3 * `div_tasks_1_se'') arguments(`cmdline' `ll0' 1.div_tasks 95)
display in smcl as text exp(r(x))

* div_tasks_2
// Lower
minbound chi2em, range(`=`div_tasks_2' - 3 * `div_tasks_2_se'', `=`div_tasks_2'') arguments(`cmdline' `ll0' 2.div_tasks 95)
display in smcl as text exp(r(x))
// Upper
minbound chi2em, range(`=`div_tasks_2'', `=`div_tasks_2' + 3 * `div_tasks_2_se'') arguments(`cmdline' `ll0' 2.div_tasks 95)
display in smcl as text exp(r(x))

* div_tasks_3
// Lower
minbound chi2em, range(`=`div_tasks_3' - 3 * `div_tasks_3_se'', `=`div_tasks_3'') arguments(`cmdline' `ll0' 3.div_tasks 95)
display in smcl as text exp(r(x))
// Upper
minbound chi2em, range(`=`div_tasks_3'', `=`div_tasks_3' + 3 * `div_tasks_3_se'') arguments(`cmdline' `ll0' 3.div_tasks 95)
display in smcl as text exp(r(x))

* freq_help_all_1
// Lower
minbound chi2em, range(`=`freq_help_all_1' - 3 * `freq_help_all_1_se'', `=`freq_help_all_1'') arguments(`cmdline' `ll0' 1.freq_help_all 95)
display in smcl as text exp(r(x))
// Upper
minbound chi2em, range(`=`freq_help_all_1'', `=`freq_help_all_1' + 3 * `freq_help_all_1_se'') arguments(`cmdline' `ll0' 1.freq_help_all 95)
display in smcl as text exp(r(x))

* freq_help_all_2
// Lower
minbound chi2em, range(`=`freq_help_all_2' - 3 * `freq_help_all_2_se'', `=`freq_help_all_2'') arguments(`cmdline' `ll0' 2.freq_help_all 95)
display in smcl as text exp(r(x))
// Upper
minbound chi2em, range(`=`freq_help_all_2'', `=`freq_help_all_2' + 3 * `freq_help_all_2_se'') arguments(`cmdline' `ll0' 2.freq_help_all 95)
display in smcl as text exp(r(x))

* freq_help_all_3
// Lower
minbound chi2em, range(`=`freq_help_all_3' - 3 * `freq_help_all_3_se'', `=`freq_help_all_3'') arguments(`cmdline' `ll0' 3.freq_help_all 95)
display in smcl as text exp(r(x))
// Upper
minbound chi2em, range(`=`freq_help_all_3'', `=`freq_help_all_3' + 3 * `freq_help_all_3_se'') arguments(`cmdline' `ll0' 3.freq_help_all 95)
display in smcl as text exp(r(x))

* freq_help_all_4
// Lower
minbound chi2em, range(`=`freq_help_all_4' - 3 * `freq_help_all_4_se'', `=`freq_help_all_4'') arguments(`cmdline' `ll0' 4.freq_help_all 95)
display in smcl as text exp(r(x))
// Upper
minbound chi2em, range(`=`freq_help_all_4'', `=`freq_help_all_4' + 3 * `freq_help_all_4_se'') arguments(`cmdline' `ll0' 4.freq_help_all 95)
display in smcl as text exp(r(x))

* Intercept
// Lower
minbound chi2em, range(`=`cons' - 3 * `cons_se'', `=`cons'') arguments(`cmdline' `ll0' _cons 95)
display in smcl as text exp(r(x))
// Upper
minbound chi2em, range(`=`cons'', `=`cons' + 3 * `cons_se'') arguments(`cmdline' `ll0' _cons 95)
display in smcl as text exp(r(x))

This is the regression output:

. firthlogit care_benefits c.health_lim i.female_r c.age_r i.educ_r_cat c.hh_inc
> ome_1000 c.assets_htsd i.div_tasks i.freq_help_all        ///
>         if valid==1 & onecareperson==1  

initial:       penalized log likelihood = -123.39006
rescale:       penalized log likelihood = -123.39006
Iteration 0:   penalized log likelihood = -123.39006  (not concave)
Iteration 1:   penalized log likelihood = -110.65803  (not concave)
Iteration 2:   penalized log likelihood = -109.63226  (not concave)
Iteration 3:   penalized log likelihood = -90.107891  
Iteration 4:   penalized log likelihood =  -89.89721  
Iteration 5:   penalized log likelihood = -81.594757  
Iteration 6:   penalized log likelihood = -79.708637  
Iteration 7:   penalized log likelihood = -79.600919  
Iteration 8:   penalized log likelihood = -79.599539  
Iteration 9:   penalized log likelihood = -79.599539  

                                                Number of obs     =      2,640
                                                Wald chi2(13)     =      90.08
Penalized log likelihood = -79.599539           Prob > chi2       =     0.0000

-------------------------------------------------------------------------------
care_benefits |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------+----------------------------------------------------------------
   health_lim |   .3397899   .0615335     5.52   0.000     .2191864    .4603934
              |
     female_r |
      1. Yes  |   .4880948   .4577666     1.07   0.286    -.4091113    1.385301
        age_r |   .0633089   .0228712     2.77   0.006     .0184822    .1081356
 1.educ_r_cat |  -.8716873   .6089216    -1.43   0.152    -2.065152    .3217771
hh_incom~1000 |   .0048384   .0023402     2.07   0.039     .0002517    .0094251
  assets_htsd |   .0099856   .0060733     1.64   0.100    -.0019179    .0218891
              |
    div_tasks |
           1  |  -1.520905    1.73532    -0.88   0.381     -4.92207     1.88026
           2  |  -3.212227   2.300189    -1.40   0.163    -7.720514    1.296061
           3  |  -2.901894   2.290494    -1.27   0.205     -7.39118    1.587392
              |
freq_help_all |
1. less th..  |   .1663776   2.360168     0.07   0.944    -4.459466    4.792222
2. about e..  |   1.684237   2.227026     0.76   0.449    -2.680654    6.049128
3. about e..  |    2.80995   1.877108     1.50   0.134    -.8691144    6.489015
4. about d..  |   2.715191   1.727395     1.57   0.116    -.6704422    6.100823
              |
        _cons |  -9.767166   1.727155    -5.66   0.000    -13.15233   -6.382004
-------------------------------------------------------------------------------

And these are the estimated confidence intervals:

* health_lim
/* 1.2396904 - 1.5860521 */

* female_r
/* .67546004 - 4.1477109 */

* age_r
/* 1.0193941 - 1.115659 */

* educ_r_cat
/* .10617454 - 1.2379033 */

* hh_income_1000
/* .99782017 - 1.0091268 */

* assets_htsd
/* .99179947 - 1.0198916 */

* div_tasks_1
/* .00875299 - 34.486815 */

* div_tasks_2
/* .00015273 - 9.9531263 */

* div_tasks_3
/* .00020653 - 13.444368 */

* freq_help_all_1
/* .00099365 - 17.87657 */

* freq_help_all_2
/* .47018018 - 64.584904 */

* freq_help_all_3
/* .09203498 - 158.48521 */

* freq_help_all_4
/* .09631961 - 117.65046 */

* Intercept
/* 1.635e-06 - .00148425 */

There is something wrong, isn't there?

Thanks for any hint!

[Joseph Coveney]

The code now works fine (see above), but the estimated confidence intervals do not include the estimated coefficients.
There is something wrong, isn't there?

No, I don't think so.

The coefficients that are reported in the regression table are in the original logit metric. Confidence bounds shown by the likelihood profile code are exponentiated. You can use the display-as-odds-ratio option when fitting the model using -firthlogit-.

firthlogit . . . , or

[Ariane Arbol]

Thank you, Joseph, you're absolutely right!