Bookmark and Share

Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

RE: st: Negative probabilities after a margins command for a categorical variable (post logistic model).


From   "Scheetz, Marc" <[email protected]>
To   "[email protected]" <[email protected]>
Subject   RE: st: Negative probabilities after a margins command for a categorical variable (post logistic model).
Date   Wed, 16 Oct 2013 13:15:58 +0000

Many thanks Steve-  your comments and are salient, and amended code is helpful. 
Marc



-----Original Message-----
From: [email protected] [mailto:[email protected]] On Behalf Of Steve Samuels
Sent: Tuesday, October 15, 2013 9:18 PM
To: [email protected]
Cc: Nick Cox
Subject: Re: st: Negative probabilities after a margins command for a categorical variable (post logistic model).


Offline, Nick Cox pointed out to me that the four line mata block can be
reduced to one line,

. mata: st_replacematrix("R", invlogit(st_matrix("R")))

I  rejected this originally because I thought that three levels of nested
functions would lack clarity. I've changed my mind. For this problem, I
find the one-liner to be more readable than my proposal below.

Steve


You are welcome, Marc. Yes, you can transform A, but that produces a
matrix with much extraneous information and no row or column names.
As you must extract and format the needed elements by hand, the
chance of error is non-negligible. Here's a better approach. Note
the use of the transform operator "'" to put the results in standard
results form.

***********Code Begins*************
sysuse auto, clear
recode rep78 1/3= 1 4=2 5=3
logistic foreign i.rep78 turn

margins, at(rep78=(1(1)3)) predict(xb)
matrix list r(table)
matrix  A = r(table)
/* Matrix to Hold Results */
matrix  R = (A["b",1...]  \ A["ll",1...] \ A["ul",1...])'
mata:
C = invlogit(st_matrix("R"))
st_replacematrix("R", C)
end

mat list R , format(%6.4f)
**********CODE ENDS***************



On Oct 15, 2013, at 11:45 AM, Scheetz, Marc wrote:

Dear Nick, Steve, and all:

Many thanks-  this fixed my problem, and now my calculated probabilities are as expected.  I greatly appreciate your response/help.

Mata is a learning curve for me.  It appears that if one is only after the 95% CIs, they can be obtained by obtaining the invlogit of matrix A (since ul and ll are already part of the matrix).

sysuse auto, clear
recode rep78 1/3= 1 4=2 5=3
logistic foreign i.rep78 turn

margins, at(rep78=(1(1)3)) predict(xb)
matrix list r(table)
matrix  A = r(table)

\\ modification step
mata: invlogit(st_matrix("A"))


Once again,  many thanks for your help.

Marc

-----Original Message-----
From: [email protected] [mailto:[email protected]] On Behalf Of Nick Cox
Sent: Tuesday, October 15, 2013 3:38 AM
To: [email protected]
Subject: Re: st: Negative probabilities after a margins command for a categorical variable (post logistic model).

Stata and Mata have -invlogit()- too. So, once you understand the
logic that Steve has clearly laid out step by step, you could also do
this:

sysuse auto, clear
recode rep78 1/3= 1 4=2 5=3
logistic foreign i.rep78 turn

margins, at(rep78=(1(1)3)) predict(xb)
matrix list r(table)
matrix  A = r(table)
/* Get rows corresponding to confidence limits */
matrix  C = A["ll",1...] \ A["ul",1...]
matrix list C
mata: invlogit(st_matrix("C"))

Nick
[email protected]


On 15 October 2013 03:37, Steve Samuels <[email protected]> wrote:
> 
> Marc:
> 
> 
> Section 5 of the FAQ lists reasons why questions don't get answered, and
> none really apply to your question. Reposting, once, after a week is
> quite acceptable under the circumstances.
> 
> Of course probabilities are restricted to [0,1]. The phenomenon you
> observed occurs when a CI is based on the formula: estimate +/- 1.96 SE,
> but there are, in fact, bounds on the true value. The "illegal"
> endpoints occur more often then you would expect.
> 
> The way around this problem is to ask -margins- to operate on the logit
> scale and then to back transform the CI endpoints. (This is what -svy:
> tabulate- does, by the way.)
> 
> Luckily -margins- returns the displayed table in a Stata matrix
> r(table). The lower and upper CIs are in rows "ll" and "ul". Below is
> code to do the work. I use Mata to simplify the calculation.
> Type: "help m2_op_colon" to understand how this worked.
> 
> 
> ***********Code Begins*************
> sysuse auto, clear
> recode rep78 1/3= 1 4=2 5=3
> logistic foreign i.rep78 turn
> 
> margins, at(rep78=(1(1)3)) predict(xb)
> matrix list r(table)
> matrix  A = r(table)
> /* Get rows corresponding to confidence limits */
> matrix  C = A["ll",1...] \ A["ul",1...]
> matrix list C
> mata:
> L =st_matrix("C")
> /* Now transform to prob scale using:
>  P  = 1/(1 + exp(-xb) */
> CI = 1:/(1 :+exp(-L))
> CI
> end
> **********CODE ENDS***************
> 
> Steve
> [email protected]
> 
> 
> 
>> 
>> On Oct 14, 2013, at 10:02 AM, Scheetz, Marc wrote:
>> 
>> Dear Listserv,
>> 
>> I am reposting a question from last week in hopes of receiving a response.  This is my first content post to the listserv; I appreciate your consideration.  Please let me know if I violated any rules for posting.
>> 
>> I am wondering if anyone can help explain the scenario below to me.  I am running Stata IC v13.0.  I am using the margins command after a multivariate-logistic model with the outcome of "died".  I am attempting to characterize the probabilities of death according to each categorical increase of the variable "log2X".  The referent category below is 2^0=1.  I have modeled the variable as categorical since I lose power due to uneven sample size in some of the categories.
>> 
>> My question is that I receive 95% CIs that have negative margins in  2 of the categories (i.e. 2._at: log2X=1, 4._at:log2X= 3).
>> 
>> Perhaps this is  a rudimentary question, but I thought that probabilities calculated from Odds Ratios could not be negative.  Is this because it is a probability relative to the referent category?  Do you see other errors in my syntax (below)?  Sincerely,
>> 
>> 
>> Marc Scheetz, PharmD, MSc
>> 
>> 
>> . logistic died i.log2X a2_day0 log10_days_to_pos_cx
>> note: 4.log2X != 0 predicts failure perfectly
>>    4.log2X dropped and 5 obs not used
>> 
>> note: 5.log2X != 0 predicts failure perfectly
>>    5.log2X dropped and 3 obs not used
>> 
>> 
>> Logistic regression                               Number of obs   =         83
>>                                                LR chi2(6)      =      18.58
>>                                                Prob > chi2     =     0.0049
>> Log likelihood = -35.358908                       Pseudo R2       =     0.2081
>> 
>> --------------------------------------------------------------------------------------
>>              died | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
>> ---------------------+--------------------------------------------------
>>       log2X       |
>>                1  |   .7903086     .92746    -0.20   0.841     .0792275    7.883466
>>                2  |   6.471551   5.420137     2.23   0.026     1.253427    33.41317
>>                3  |   1.587899   1.492738     0.49   0.623     .2515548    10.02335
>>                4  |          1  (empty)
>>                5  |          1  (empty)
>>                6  |   6.542207   6.159993     1.99   0.046     1.033362    41.41868
>>                   |
>>           a2_day0 |   1.075268   .0732118     1.07   0.286     .9409374    1.228775
>> log10_days_to_pos_cx |   4.854903   3.054503     2.51   0.012      1.41462    16.66177
>>             _cons |    .012261   .0189665    -2.85   0.004     .0005913    .2542394
>> 
>> 
>> . margins, at(log2X=(0(1)6))
>> 
>> Predictive margins                                Number of obs   =         83
>> Model VCE    : OIM
>> 
>> Expression   : Pr(died), predict()
>> 
>> 1._at        : log2X     =           0
>> 
>> 2._at        : log2X     =           1
>> 
>> 3._at        : log2X     =           2
>> 
>> 4._at        : log2X     =           3
>> 
>> 5._at        : log2X     =           4
>> 
>> 6._at        : log2X     =           5
>> 
>> 7._at        : log2X     =           6
>> 
>> ------------------------------------------------------------------------------
>>           |            Delta-method
>>           |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
>> -------------+----------------------------------------------------------
>>       _at |
>>        1  |   .1518137   .0514472     2.95   0.003     .0509791    .2526484
>>        2  |   .1259233   .1108515     1.14   0.256    -.0913416    .3431883
>>        3  |   .4764486   .1462235     3.26   0.001     .1898558    .7630413
>>        4  |   .2137861   .1241819     1.72   0.085    -.0296061    .4571782
>>        5  |          .  (not estimable)
>>        6  |          .  (not estimable)
>>        7  |   .4787175   .1765856     2.71   0.007      .132616     .824819
>> ------------------------------------------------------------------------------
>> 
>> 
> 
> *
> *   For searches and help try:
> *   http://www.stata.com/help.cgi?search
> *   http://www.stata.com/support/faqs/resources/statalist-faq/
> *   http://www.ats.ucla.edu/stat/stata/
> 
> 
> *
> *   For searches and help try:
> *   http://www.stata.com/help.cgi?search
> *   http://www.stata.com/support/faqs/resources/statalist-faq/
> *   http://www.ats.ucla.edu/stat/stata/

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/


*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/


© Copyright 1996–2018 StataCorp LLC   |   Terms of use   |   Privacy   |   Contact us   |   Site index