st: Graphing Marginal Effects after -xtprobit-

["Pavlos C. Symeou" <p.symeou@gmail.com>]

Dear Statalisters,

I am attempting to graph the marginal effects of an interaction 
term whose coefficient is estimated with the -xtprobit- command, 
inspired by the Stata code which can be found here 
http://homepages.nyu.edu/~mrg217/interaction.html. The original 
code is based on a -probit- model. Stata interrupts my do file in 
line 4 and gives me the following error:

. drawnorm MG_b1-MG_b6, n(10000) means(e(b)) cov(e(V)) clear;
6 variables specified but specified means(e(b)) is 1 x 7
expected a 6-vector
r(503);

I speculate the missing variable is the unobserved error term 
encountered in panel models. How do I account for this? Apparently, 
I will need to address this in lines 13, 14, and 15 where I give 
values to my variables (including the constant).

Could you suggest a solution to this problem?


Pavlos


**********Code_Follows**************************************************************************

1. #delimit ;
2. xtprobit str_logit   dev_log_ch_lag1   pr_lag1  
dev_log_ch_pr_lag1  log_emp roa, re;
3. preserve;
4. drawnorm MG_b1-MG_b6, n(10000) means(e(b)) cov(e(V)) clear;

5. postutil clear;
6. postfile mypost prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat 
diff_lo diff_hi
            using sim , replace;
            noisily display "start";

7. local a=0 ;
8. while `a' <= 2 { ; /*a here represents the range of values that 
pr takes 0,1,2*/

9    {;

10. scalar h_dev_log_ch_lag1=0.00000000173;
11. scalar h_log_emp=1.243789;
12. scalar h_roa=0.0228954;
13. scalar h_constant=1;

14.     generate x_betahat0 = MG_b1*h_dev_log_ch_lag1
                            + MG_b2*(`a')
                            + MG_b3*h_dev_log_ch_lag1*(`a')
                            + MG_b4*h_log_emp
                            + MG_b5*h_roa
                            + MG_b6*h_constant;


15.    generate x_betahat1 = MG_b1*(h_dev_log_ch_lag1+1.75814) 
/*1.75814 is the SD of dev_log_ch_lag1*/
                            + MG_b2*`a'
                            + MG_b3*(h_dev_log_ch_lag1+1.75814)*(`a')
                            + MG_b4*h_log_emp
                            + MG_b5*h_roa
                            + MG_b6*h_constant;

16.    gen prob0=norm(x_betahat0);
    gen prob1=norm(x_betahat1);
    gen diff=prob1-prob0;

    egen probhat0=mean(prob0);
    egen probhat1=mean(prob1);
    egen diffhat=mean(diff);

    tempname prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat diff_lo 
diff_hi ;

    _pctile prob0, p(2.5,97.5) ;
    scalar `lo0' = r(r1);
    scalar `hi0' = r(r2);

    _pctile prob1, p(2.5,97.5);
    scalar `lo1'= r(r1);
    scalar `hi1'= r(r2);

    _pctile diff, p(2.5,97.5);
    scalar `diff_lo'= r(r1);
    scalar `diff_hi'= r(r2);


    scalar `prob_hat0'=probhat0;
    scalar `prob_hat1'=probhat1;
    scalar `diff_hat'=diffhat;

    post mypost (`prob_hat0') (`lo0') (`hi0') (`prob_hat1') 
(`lo1') (`hi1')
                (`diff_hat') (`diff_lo') (`diff_hi') ;
    };
    drop x_betahat0 x_betahat1 prob0 prob1 diff probhat0 probhat1 
diffhat ;
    local a=`a'+ 1 ; /*A unit increase in pr_lag1. It can be 
whatever I think it makes sense*/
    display "." _c;
    } ;

display "";

postclose mypost;

use sim, clear;

gen MV = _n-1;

graph twoway line diff_hat MV, clwidth(medium) clcolor(blue) 
clcolor(black)
        ||   line diff_lo  MV, clpattern(dash) clwidth(thin) 
clcolor(black)
        ||   line diff_hi  MV, clpattern(dash) clwidth(thin) 
clcolor(black)
        ||  ,
            xlabel(0 1 2, labsize(3)) /*I think these should be 
the values of the modifying variable*/
            ylabel(-.1 0 0.1 0.2, labsize(3))
            yscale(noline)
            xscale(noline)
            yline(0)
            legend(off)
            title("Marginal Effect of Log_Ch on Env_Performance As 
Pressure Changes", size(4))
            subtitle(" " "Dependent Variable: 
Environmental_Performance" " ", size(3))
            xtitle(pr_lag1, size(3.5)  )
            ytitle("Marginal Effect of Log_Ch", size(3.5))
            xsca(titlegap(2))
            ysca(titlegap(2))
            scheme(s2mono) graphregion(fcolor(white));

graph export  h:\figure1.eps, replace;

translate @Graph h:\figure1.wmf;

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