Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.
| 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/