Firthlogit, marginal effects and interactions

Ariane Arbol

Join Date: May 2021
Posts: 36

Firthlogit, marginal effects and interactions

21 May 2021, 04:36

Dear all,

I tried to run a logit model with an interaction effect which unfortunately does not work due to separation:

Code:

. logit care_benefits c.health_lim##i.onecareperson_3 if valid==1 & onecarepersononly==1

note: 2.onecareperson_3 != 0 predicts failure perfectly
      2.onecareperson_3 dropped and 65 obs not used

note: 3.onecareperson_3#c.health_lim != 0 predicts failure perfectly
      3.onecareperson_3#c.health_lim dropped and 32 obs not used

note: 2.onecareperson_3#c.health_lim omitted because of collinearity
Iteration 0:   log likelihood = -47.883416  
Iteration 1:   log likelihood = -39.596133  
Iteration 2:   log likelihood = -35.881327  
Iteration 3:   log likelihood = -35.798874  
Iteration 4:   log likelihood = -35.798488  
Iteration 5:   log likelihood = -35.798488  

Logistic regression                             Number of obs     =        320
                                                LR chi2(2)        =      24.17
                                                Prob > chi2       =     0.0000
Log likelihood = -35.798488                     Pseudo R2         =     0.2524

----------------------------------------------------------------------------------------------
               care_benefits |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-----------------------------+----------------------------------------------------------------
                  health_lim |   .3637094    .085982     4.23   0.000     .1951877     .532231
                             |
             onecareperson_3 |
               wider family  |          0  (empty)
                 not family  |  -.3931463   1.182315    -0.33   0.739    -2.710441    1.924148
                             |
onecareperson_3#c.health_lim |
               wider family  |          0  (empty)
                 not family  |          0  (omitted)
                             |
                       _cons |  -4.360444   .6238969    -6.99   0.000    -5.583259   -3.137628
----------------------------------------------------------------------------------------------

I tried to fix that with firthlogit, what actually worked:

Code:

. firthlogit care_benefits c.health_lim##i.onecareperson_3 if valid==1 & onecarepersononly==1

initial:       penalized log likelihood = -45.352147
rescale:       penalized log likelihood = -45.352147
Iteration 0:   penalized log likelihood = -45.352147  
Iteration 1:   penalized log likelihood = -32.617258  
Iteration 2:   penalized log likelihood = -31.545871  
Iteration 3:   penalized log likelihood = -31.462159  
Iteration 4:   penalized log likelihood =   -31.4621  
Iteration 5:   penalized log likelihood =   -31.4621  

                                                Number of obs     =        417
                                                Wald chi2(5)      =      27.73
Penalized log likelihood =   -31.4621           Prob > chi2       =     0.0000

----------------------------------------------------------------------------------------------
               care_benefits |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-----------------------------+----------------------------------------------------------------
                  health_lim |   .3502754   .0817553     4.28   0.000      .190038    .5105129
                             |
             onecareperson_3 |
               wider family  |  -.6839636   1.635612    -0.42   0.676    -3.889704    2.521777
                 not family  |  -.4077488   1.045765    -0.39   0.697    -2.457411    1.641913
                             |
onecareperson_3#c.health_lim |
               wider family  |  -.0106092   .3302318    -0.03   0.974    -.6578516    .6366332
                 not family  |  -.0372747   .3046006    -0.12   0.903     -.634281    .5597315
                             |
                       _cons |  -4.206488   .5869453    -7.17   0.000    -5.356879   -3.056096
----------------------------------------------------------------------------------------------

However, to be able to interpret the effects more clearly, I would like to generate the marginal effects from this. I therefore tried the margins postestimation command with the option expression(invlogit(predict(xb))) which I found (among others) in this thread. The result is the following:

Code:

. margins, dydx(*) expression(invlogit(predict(xb)))

Average marginal effects                        Number of obs     =        417
Model VCE    : OIM

Expression   : invlogit(predict(xb))
dy/dx w.r.t. : health_lim 2.onecareperson_3 3.onecareperson_3

---------------------------------------------------------------------------------
                |            Delta-method
                |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
----------------+----------------------------------------------------------------
     health_lim |   .0091291    .002512     3.63   0.000     .0042057    .0140525
                |
onecareperson_3 |
  wider family  |  -.0174202   .0227347    -0.77   0.444    -.0619794    .0271389
    not family  |   -.014628   .0252153    -0.58   0.562    -.0640492    .0347931
---------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.

I'm a little bit confused now about the lack of the interaction effects. How does the command work correctly? Is it correct at all? Which approach do you recommend here?

Thanks for any help!

Tags: firthlogit, interaction, interactions, logit, margins

Joseph Coveney

Join Date: Apr 2014

Posts: 4406
#2

21 May 2021, 05:26

Originally posted by Ariane Arbol View Post

I'm a little bit confused now about the lack of the interaction effects. How does the command work correctly?

Does the following help?

.ÿ
.ÿversionÿ17.0

.ÿ
.ÿclearÿ*

.ÿ
.ÿsetÿseedÿ`=strreverse("1610971")'

.ÿ
.ÿquietlyÿsetÿobsÿ10000

.ÿ
.ÿgenerateÿbyteÿcatÿ=ÿruniform()ÿ<ÿ0.5

.ÿgenerateÿdoubleÿconÿ=ÿruniform(-0.5,ÿ0.5)

.ÿ
.ÿgenerateÿdoubleÿoutÿ=ÿ100ÿ+ÿcond(cat,ÿ-con,ÿcon)ÿ+ÿrnormal(0,ÿ0.25)

.ÿ
.ÿregressÿoutÿc.con##i.cat,

ÿÿÿÿÿÿSourceÿ|ÿÿÿÿÿÿÿSSÿÿÿÿÿÿÿÿÿÿÿdfÿÿÿÿÿÿÿMSÿÿÿÿÿÿNumberÿofÿobsÿÿÿ=ÿÿÿÿ10,000
-------------+----------------------------------ÿÿÿF(3,ÿ9996)ÿÿÿÿÿÿ=ÿÿÿ4359.33
ÿÿÿÿÿÿÿModelÿ|ÿÿ826.676254ÿÿÿÿÿÿÿÿÿ3ÿÿ275.558751ÿÿÿProbÿ>ÿFÿÿÿÿÿÿÿÿ=ÿÿÿÿ0.0000
ÿÿÿÿResidualÿ|ÿÿ631.859149ÿÿÿÿÿ9,996ÿÿ.063211199ÿÿÿR-squaredÿÿÿÿÿÿÿ=ÿÿÿÿ0.5668
-------------+----------------------------------ÿÿÿAdjÿR-squaredÿÿÿ=ÿÿÿÿ0.5667
ÿÿÿÿÿÿÿTotalÿ|ÿÿÿ1458.5354ÿÿÿÿÿ9,999ÿÿ.145868127ÿÿÿRootÿMSEÿÿÿÿÿÿÿÿ=ÿÿÿÿ.25142

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿoutÿ|ÿCoefficientÿÿStd.ÿerr.ÿÿÿÿÿÿtÿÿÿÿP>|t|ÿÿÿÿÿ[95%ÿconf.ÿinterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿconÿ|ÿÿÿ.9808643ÿÿÿ.0123234ÿÿÿÿ79.59ÿÿÿ0.000ÿÿÿÿÿ.9567079ÿÿÿÿ1.005021
ÿÿÿÿÿÿÿ1.catÿ|ÿÿÿ.0023312ÿÿÿ.0050294ÿÿÿÿÿ0.46ÿÿÿ0.643ÿÿÿÿ-.0075274ÿÿÿÿ.0121898
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿcat#c.conÿ|
ÿÿÿÿÿÿÿÿÿÿ1ÿÿ|ÿÿ-1.990248ÿÿÿ.0174072ÿÿ-114.33ÿÿÿ0.000ÿÿÿÿ-2.024369ÿÿÿ-1.956126
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿ_consÿ|ÿÿÿ99.99942ÿÿÿ.0035544ÿÿ2.8e+04ÿÿÿ0.000ÿÿÿÿÿ99.99245ÿÿÿÿ100.0064
------------------------------------------------------------------------------

.ÿ
.ÿmarginsÿcat,ÿdydx(con)

AverageÿmarginalÿeffectsÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿ=ÿ10,000
ModelÿVCE:ÿOLS

Expression:ÿLinearÿprediction,ÿpredict()
dy/dxÿwrt:ÿÿcon

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿÿÿÿÿÿDelta-method
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿdy/dxÿÿÿstd.ÿerr.ÿÿÿÿÿÿtÿÿÿÿP>|t|ÿÿÿÿÿ[95%ÿconf.ÿinterval]
-------------+----------------------------------------------------------------
conÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿcatÿ|
ÿÿÿÿÿÿÿÿÿÿ0ÿÿ|ÿÿÿ.9808643ÿÿÿ.0123234ÿÿÿÿ79.59ÿÿÿ0.000ÿÿÿÿÿ.9567079ÿÿÿÿ1.005021
ÿÿÿÿÿÿÿÿÿÿ1ÿÿ|ÿÿ-1.009383ÿÿÿ.0122941ÿÿÿ-82.10ÿÿÿ0.000ÿÿÿÿ-1.033482ÿÿÿ-.9852843
------------------------------------------------------------------------------

.ÿ
.ÿexit

endÿofÿdo-file

.
Comment
Ariane Arbol

Join Date: May 2021

Posts: 36
#3

22 May 2021, 10:06

Thank you, Joseph!

I think I do understand that the given effects are the resulting interaction effects, right? My output therefore would be interpreted as:
"The marginal effect of health_lim for respondents with a careperson from the close family [reference category] is 0,009.
The marginal effect of health_lim for respondents with a careperson from the wider family is -0,017.
The marginal effect of health_lim for respondents with a careperson not from the family is -0,015."

Is that correct?

Would the following code produce a correct graphical output of the marginal effects?:

Code:

margins, at(health_lim=(0(1)15) onecareperson=(1(1)3)) expression(invlogit(predict(xb))) marginsplot

I would appreciate some feedback!
Comment
Joseph Coveney

Join Date: Apr 2014

Posts: 4406
#4

22 May 2021, 19:57

Originally posted by Ariane Arbol View Post

I think I do understand that the given effects are the resulting interaction effects, right?

Just from inspection of the regression table output shown above in #1 (the Wald test results for the two components of the interaction term have P = 1.0 and P = 0.9), I would say that there are no "interaction effects" worth concerning oneself over.
Comment
Ariane Arbol

Join Date: May 2021

Posts: 36
#5

07 Jun 2021, 10:13

Joseph, thank you again for your answer. (Sorry for the late reply, unfortunately I got ill and was barely able to think ...)

Originally posted by Joseph Coveney View Post

Just from inspection of the regression table output shown above in #1 (the Wald test results for the two components of the interaction term have P = 1.0 and P = 0.9), I would say that there are no "interaction effects" worth concerning oneself over.

Right, there are no significant interaction effects. But just for understanding, would my interpretation be correct?
Comment
Joseph Coveney

Join Date: Apr 2014

Posts: 4406
#6

07 Jun 2021, 22:20

Originally posted by Ariane Arbol View Post

would my interpretation be correct?

Sorry, I don't use -margins- after -firthlogit-.
Comment

Ariane Arbol

Join Date: May 2021
Posts: 36

12 Jun 2021, 10:15

Thank you, Joseph, either way!

If anyone is interested:
To get the average marginal effects of (in my case) health_lim for the different careperson-groups, you would have to use:

Code:

margins, dydx(health_lim) at(onecareperson_3=1) at(onecareperson_3=2) expression(invlogit(predict(xb)))

Code:

margins onecareperson_3, dydx(health_lim) expression(invlogit(predict(xb)))

(same results).

To get the 'marginal interaction effects', you would type:

Code:

margins onecareperson_3, dydx(health_lim) pwcompare expression(invlogit(predict(xb)))

(or just calculate the differences between the categories of the output from the first command (see above)).

Sources:
https://www.stata.com/stata-news/news31-3/spotlight/
https://www.statalist.org/forums/for...argins-command

___

Joseph, may I ask a further question?:
I was wondering how to get the Profile Likelihood Confidence Intervals recommended by Heinze & Schemper (2002) for Firth regression. For this I found this thread and tried the code you provided:

Code:

program define chi2em, rclass
    args x cmdline 110 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(`110' - e(11))
    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 110 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 div_tasks_se freq_help_all freq_help_all_se cons cons_se

scalar define `110' = e(11)
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' = _b[1.div_tasks]
scalar define `div_tasks_se' = _se[1.div_tasks]

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

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

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

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

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

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

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

Unfortunately this did not work, I got the following error message:

Code:

. scalar define `110' = e(11)
11 invalid name
r(198);

I certainly didn't understand everything you wrote down in your code. Did I do anything wrong?

Very thankful for hints and corrections!

Comment

Joseph Coveney

Join Date: Apr 2014

Posts: 4406
#8

12 Jun 2021, 17:10

For the text editor that you're using to write code, I recommend using a typeface that makes it easy to distinguish between Arabic numeral one and lower-case letter el. I think that StataCorp uses Consolas by default now for its do-file editor, but you can experiment to find something that works best for you. (Obviously, the Courier New that this forum's software uses for its code blocks doesn't cut it.)

In your case, change

Code:

scalar define `110' = e(11)

to

Code:

scalar define `ll0’ = e(ll)

which returns the log-likelihood in an estimation (ereturn) scalar.
Comment

Ariane Arbol

Join Date: May 2021
Posts: 36

13 Jun 2021, 08:09

Ah, thank you again, Joseph!

I corrected the mistyping and the code given above is working fine now. I therefore added the next steps:

Code:

* health_lim
// Lower
minbound chi2em, range(`=`health_lim' - 3 * `health_lim_se'', `=`health_lim'') arguments(`cmdline' `ll0' health_lim 95)
// Upper
minbound chi2em, range(`=`health_lim' + 3 * `health_lim_se'', `=`health_lim'') arguments(`cmdline' `ll0' health_lim 95)

* female_r
// Lower
minbound chi2em, range(`=`female_r' - 3 * `female_r_se'', `=`female_r'') arguments(`cmdline' `ll0' 1.female_r 95)
// Upper
minbound chi2em, range(`=`female_r' + 3 * `female_r_se'', `=`female_r'') arguments(`cmdline' `ll0' 1.female_r 95)

* age_r
// Lower
minbound chi2em, range(`=`age_r' - 3 * `age_r_se'', `=`age_r'') arguments(`cmdline' `ll0' age_r 95)
// Upper
minbound chi2em, range(`=`age_r' + 3 * `age_r_se'', `=`age_r'') arguments(`cmdline' `ll0' age_r 95)

* 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)
// Upper
minbound chi2em, range(`=`educ_r_cat' + 3 * `educ_r_cat_se'', `=`educ_r_cat'') arguments(`cmdline' `ll0' 1.educ_r_cat 95)

* 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)
// Upper
minbound chi2em, range(`=`hh_income_1000' + 3 * `hh_income_1000_se'', `=`hh_income_1000'') arguments(`cmdline' `ll0' hh_income_1000 95)

* assets_htsd
// Lower
minbound chi2em, range(`=`assets_htsd' - 3 * `assets_htsd_se'', `=`assets_htsd'') arguments(`cmdline' `ll0' assets_htsd 95)
// Upper
minbound chi2em, range(`=`assets_htsd' + 3 * `assets_htsd_se'', `=`assets_htsd'') arguments(`cmdline' `ll0' assets_htsd 95)

* 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)
// Upper
minbound chi2em, range(`=`div_tasks_1' + 3 * `div_tasks_1_se'', `=`div_tasks_1'') arguments(`cmdline' `ll0' 1.div_tasks 95)

* 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)
// Upper
minbound chi2em, range(`=`div_tasks_2' + 3 * `div_tasks_2_se'', `=`div_tasks_2'') arguments(`cmdline' `ll0' 2.div_tasks 95)

* 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)
// Upper
minbound chi2em, range(`=`div_tasks_3' + 3 * `div_tasks_3_se'', `=`div_tasks_3'') arguments(`cmdline' `ll0' 3.div_tasks 95)

* 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)
// Upper
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)

* 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)
// Upper
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)

* 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)
// Upper
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)

* 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)
// Upper
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)

and immediately got the following error message:

Code:

. * health_lim
. // Lower
. minbound chi2em, range(`=`health_lim' - 3 * `health_lim_se'', `=`health_lim'')
>  arguments(`cmdline' `ll0' health_lim 95)

. // Upper
. minbound chi2em, range(`=`health_lim' + 3 * `health_lim_se'', `=`health_lim'')
>  arguments(`cmdline' `ll0' health_lim 95)

range() invalid
lowerbound larger than upperbound
r(198);

Argh! Any idea what to do with that? Did I do something wrong again?

Comment

Joseph Coveney

Join Date: Apr 2014

Posts: 4406
#10

13 Jun 2021, 17:29

Originally posted by Ariane Arbol View Post

range() invalid
lowerbound larger than upperbound
r(198);
Argh! Any idea what to do with that? Did I do something wrong again?

Try swapping the two arguments in the -range()- option.

When your regression coefficient is negative, your starting values in the search need to be swapped so that the most-negative value is the left argument and the least negative value is on the righthand side.

In case you weren't aware of it, there's a user-written command that automates profile likelihood searches.

Code:

search pllf

will lead you to where you can download and install it.

I've never used it and so I don't know whether it will work with -firthlogit-, but if you're having trouble following what I showed to profile the likelihood manually, then you might want to give it a try.
Comment

Ariane Arbol

Join Date: May 2021
Posts: 36

#11

14 Jun 2021, 13:33

Originally posted by Joseph Coveney View Post

Try swapping the two arguments in the -range()- option.

When your regression coefficient is negative, your starting values in the search need to be swapped so that the most-negative value is the left argument and the least negative value is on the righthand side.

Hm, the regression coefficient of health_lim is not negative. Actually, the regression output looks like this:

Code:

. 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
-------------------------------------------------------------------------------

Following your advice, I swapped the two arguments in the -range()- option every time the above error occurred (changes in bold):

Code:

* 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' + 3 * `div_tasks_1_se'', `=`div_tasks_1'') 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' + 3 * `div_tasks_2_se'', `=`div_tasks_2'') 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' + 3 * `div_tasks_3_se'', `=`div_tasks_3'') 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' + 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))

* 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' + 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))

* 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' + 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))

* 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' + 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))

* 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' + 3 * `cons_se'', `=`cons'') arguments(`cmdline' `ll0' _cons 95)
display in smcl as text exp(r(x))

Is that correct up to this point? I am kind of worried about producing artifacts here due to the error focused procedure ...

And then, next question and problem, with 'div_tasks_1' I got an error message again:

Code:

. . * 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)
invalid syntax
invalid syntax
invalid syntax
r(198);

This error has appeared before, namely when I do not start the 'minbound' command series from the very top (at 'health_lim'), but only execute single 'minbound' commands. I thought therefore that the command series must be executed as a whole in order to work correctly.
Now I'm wondering if the calculation at 'assets_htsd' takes so long that Stata stops after that (I don't know if that's really a coherent explanation, since I only marginally understand the background of the commands). Or whether I've put in a mistake somewhere again that I do not find. Can you see or guess what could be causing this error message?

Originally posted by Joseph Coveney View Post

In case you weren't aware of it, there's a user-written command that automates profile likelihood searches.

Yes, I am aware of the pllf command, but it didn't work with firthlogit when I tried it. So I used your code. :-)

Thank you again for help!

Comment

Joseph Coveney

Join Date: Apr 2014

Posts: 4406
#12

14 Jun 2021, 17:42

You're better off, I think, breaking the problem down, taking one thing at a time, and debug problems using tracing (set trace on) until each piece works.

Once you get everything ironed out, you can then re-assemble and run the entire list of coefficients to the end in one fell swoop in your final do-file.
Comment

Ariane Arbol

Join Date: May 2021
Posts: 36

#13

15 Jun 2021, 16:45

Okay, I executed the whole code (from 'program' to the end of the 'minbound' series) with the preceding 'set trace on' command. That took a very long time and in the end, I got the following (only end of the output where the error happens):

Code:

. * 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)
invalid syntax
invalid syntax
  ----------------------------------------------------------- begin minbound ---
  - version 9
  - syntax anything(name=func id=program) , RANge(str) [ FROM(str) ARGuments(str
> ) TOLerance(real 1e-5) ITERate(int 100) MISSing LOg TRace ]
  - if "`log'" != "" {
  = if "" != "" {
    local trace trace
    }
  - tempname cgold zeps
  - scalar `cgold' = 0.3819660
  = scalar __00000K = 0.3819660
  - scalar `zeps' = sqrt(c(epsdouble))
  = scalar __00000L = sqrt(c(epsdouble))
  - tempname a ax b bx d e etemp p q r tol1 tol2
  - tempname u v w x xmid fa fb fu fv fw fx
  - gettoken toka range: range, match(lmac) parse(", ")
  - gettoken tokb range: range, match(lmac) parse(", ")
  - if `"`tokb'"'=="," {
  = if `""'=="," {
    gettoken tokb range: range, match(lmac) parse(", ")
    }
  - if `"`range'"'!="" {
  = if `""'!="" {
    dis as err "range() invalid"
    dis as err "unexpected text found: " `"`range'"'
    exit 198
    }
  - scalar `ax' = `toka'
  = scalar __00000N = ,
invalid syntax
  ------------------------------------------------------------- end minbound ---
r(198);

I tried to compare this with previous output, but unfortunately that output is too long (for the output window) or too large (for the log file), so I can't currently access it. However, the three bold lines somehow indicate to me that the error might be there? But why should 'range' suddenly be unknown, it worked before, didn't it?

(I know that's quite a lot to ask to find the error. If you don't see it right away, please don't waste your time. However, I am grateful for any hint that may get me closer to the solution. Any idea?)

Comment

Joseph Coveney

Join Date: Apr 2014

Posts: 4406
#14

15 Jun 2021, 19:04

It's not the lines that you bolded, but the scalar assignment afterward. The righthand side is empty, which means that the argument that you supplied -range()- is empty or undefined.

Just before the -minbound- command, put in a line

Code:

display in smcl as text "Lower = " as result `=`div_tasks_1' - 3 * `div_tasks_1_se''

and see whether there is anything after the "Lower = ".
Comment

Ariane Arbol

Join Date: May 2021
Posts: 36

#15

16 Jun 2021, 16:05

Originally posted by Joseph Coveney View Post

Just before the -minbound- command, put in a line

Code:

display in smcl as text "Lower = " as result `=`div_tasks_1' - 3 * `div_tasks_1_se''

and see whether there is anything after the "Lower = ".

I did and got this:

Code:

. . * div_tasks_1
. // Lower
. display in smcl as text "Lower = " as result `=`div_tasks_1' - 3 * `div_tasks_
> 1_se''
invalid syntax
Lower =

. minbound chi2em, range(`=`div_tasks_1' - 3 * `div_tasks_1_se'', `=`div_tasks_1
> '') arguments(`cmdline' `ll0' 1.div_tasks 95)
invalid syntax
invalid syntax
  ----------------------------------------------------------- begin minbound ---
  - version 9
  - syntax anything(name=func id=program) , RANge(str) [ FROM(str) ARGuments(str
> ) TOLerance(real 1e-5) ITERate(int 100) MISSing LOg TRace ]
  - if "`log'" != "" {
  = if "" != "" {
    local trace trace
    }
  - tempname cgold zeps
  - scalar `cgold' = 0.3819660
  = scalar __00000K = 0.3819660
  - scalar `zeps' = sqrt(c(epsdouble))
  = scalar __00000L = sqrt(c(epsdouble))
  - tempname a ax b bx d e etemp p q r tol1 tol2
  - tempname u v w x xmid fa fb fu fv fw fx
  - gettoken toka range: range, match(lmac) parse(", ")
  - gettoken tokb range: range, match(lmac) parse(", ")
  - if `"`tokb'"'=="," {
  = if `""'=="," {
    gettoken tokb range: range, match(lmac) parse(", ")
    }
  - if `"`range'"'!="" {
  = if `""'!="" {
    dis as err "range() invalid"
    dis as err "unexpected text found: " `"`range'"'
    exit 198
    }
  - scalar `ax' = `toka'
  = scalar __00000N = ,
invalid syntax
  ------------------------------------------------------------- end minbound ---
r(198);

Looks like there is a line with "Lower = " and nothing else. So what now? How to fix that? What did I wrong ... I did define the scalar, didn't I? In a wrong way? The entire command sequence currently looks like this:

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 div_tasks_se freq_help_all freq_help_all_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
display in smcl as text "Lower = " as result `=`div_tasks_1' - 3 * `div_tasks_1_se''
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' + 3 * `div_tasks_1_se'', `=`div_tasks_1'') 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' + 3 * `div_tasks_2_se'', `=`div_tasks_2'') 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' + 3 * `div_tasks_3_se'', `=`div_tasks_3'') 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' + 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))

* 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' + 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))

* 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' + 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))

* 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' + 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))

* 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' + 3 * `cons_se'', `=`cons'') arguments(`cmdline' `ll0' _cons 95)
display in smcl as text exp(r(x))

Did I refer wrong in the scalar definition?

Announcement

Firthlogit, marginal effects and interactions

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment