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]

From |
"Scheetz, Marc" <mschee@midwestern.edu> |

To |
"statalist@hsphsun2.harvard.edu" <statalist@hsphsun2.harvard.edu> |

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: owner-statalist@hsphsun2.harvard.edu [mailto:owner-statalist@hsphsun2.harvard.edu] On Behalf Of Steve Samuels Sent: Tuesday, October 15, 2013 9:18 PM To: statalist@hsphsun2.harvard.edu 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: owner-statalist@hsphsun2.harvard.edu [mailto:owner-statalist@hsphsun2.harvard.edu] On Behalf Of Nick Cox Sent: Tuesday, October 15, 2013 3:38 AM To: statalist@hsphsun2.harvard.edu 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 njcoxstata@gmail.com On 15 October 2013 03:37, Steve Samuels <sjsamuels@gmail.com> 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 > sjsamuels@gmail.com > > > >> >> 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/

**References**:**st: Negative probabilities after a margins command for a categorical variable (post logistic model).***From:*"Scheetz, Marc" <mschee@midwestern.edu>

**Re: st: Negative probabilities after a margins command for a categorical variable (post logistic model).***From:*Steve Samuels <sjsamuels@gmail.com>

**Re: st: Negative probabilities after a margins command for a categorical variable (post logistic model).***From:*Nick Cox <njcoxstata@gmail.com>

**RE: st: Negative probabilities after a margins command for a categorical variable (post logistic model).***From:*"Scheetz, Marc" <mschee@midwestern.edu>

**Re: st: Negative probabilities after a margins command for a categorical variable (post logistic model).***From:*Steve Samuels <sjsamuels@gmail.com>

- Prev by Date:
**st: RE: dictionary problem** - Next by Date:
**st: RE: RE: dictionary problem** - Previous by thread:
- Next by thread:
**RE: RE: st: Negative probabilities after a margins command for a categorical variable (post logistic model).** - Index(es):