Predicted probabilities of square term in Stata

Question

I am estimating a multivariate logit model.

I have included a square term on my main independent variable, IDV and IDV:IDV. And of course a set of control variables. The results indicate an inverted U shaped relationship; IDV is positive, while IDV:IDV is negative. I have graphed the probabilities using -postgr3-:

postgr3 IDV, asis(IDV IDV*IDV)

However, I am unsure about how to get the probabilities of different values of IDV. To be clear, I am looking for something similar to the information -clarify- gives, i.e. the probabilities for different percentiles of my IDV, while taking into account the squared term.

Any ideas? It would be very much appreciated. Let me know if I my question is not clear. I recently started using Stata.

Answer 1

I am unfortunately unfamiliar with postgr3 and -clarify-, but I have had a similar problem before and ended up employing a solution used by Andreas Wimmer & Brian Min (2006) (article content unrelated).

Basically it boils down to four steps:

1. create copies of your control variables.
2. regress with the copies.
3. fix the value of the control variable copies for all cases using whatever assumptions you deem appropriate.
4. use the regression equation to compute predicted values for all cases.

With sorted data this produces a nice distribution of predicted probabilities estimates which can be easily plotted. You can also easily get summary data or run additional statistics on the predicted probabilities.

/*logistic regression, c1, c2, and c3 are controls*/
logit y xlin xsq c1 c2 c3

/*generate graphing copies*/
    gen c1g = c1
    gen c2g = c2
    gen c3g = c3

/*regress on copies, resultis will be idential to above regression*/
    logit y xlin xsq c1g c2g c3g

/*replace copies with fixed values, in this case median values*/
    quietly sum c1, det
    replace c1g = r(p50)
    quietly sum c2, det
    replace c2g = r(p50)
    quietly sum c3, det
    replace c3g = r(p50)

/*generate predicted values*/
    sort xlin
    predict pr, p /*probability*/
    predict logodds, xb
    predict stderr, stdp
    generate lodds_lb = logodds - 1.96*stderr
    generate lodds_ub = logodds + 1.96*stderr
    generate ub_p = exp(lodds_ub)/(1+exp(lodds_ub)) /*upper confidence band*/
    generate lb_p = exp(lodds_lb)/(1+exp(lodds_lb)) /*lower confidence band*/

/* PLOT RESULTS */
        twoway (rarea lb_p ub_p xlin, bcolor(gs14)) ///
        (line pr xlin, clcolor(black) clwidth(medthick)), ///
            xline(0, lp(dash) lc(gs14) lw(thin)) ///
            ylabel(#8, labsize(small)) xlabel(#20, labsize(small)) ///
            ytitle(Predicted probability of y, size(small)) xtitle(x, size(small)) ///
            legend(order(2 "Probability of y" 1 "95% confidence interval") size(small) rows(1)) ///
            graphregion(color(white))