Thanks very much, Kit and Maarten. This gives me a lot to think about.
I hadn't missed Kit's first hints about lincom, I was just unsure about
what to do with them. :)
-JW
-----Original Message-----
From: [email protected]
[mailto:[email protected]] On Behalf Of Maarten buis
Sent: Thursday, October 11, 2007 8:09 AM
To: [email protected]
Subject: Re: st: RE: re: missing dummy variable
--- Maarten buis <[email protected]> wrote:
> I think that the display of this many parameters would call for a
> graph, so in the example below I show how I would add such
> information together into one graph. A lot of the tricks I use are
> explained here:
> http://home.fsw.vu.nl/m.buis/wp/pvalue.html
I have added some extra bells and wistles to the graph. I have probably
added too many bells and wistles, but
a) you can pick from this code whatever element you like, and
b) it nicely illustrates that models focusing on the indidual means
of
the different groups or on deviatons from a refence group, or a
grand mean are just the same model.
Under normal surcomstances I would leave out the right y-axis and
probably the confidence interval around the grand mean, but that would
completely defeat the purpose of the original question.
*------------------ begin example ---------------
sysuse auto, clear
gen domestic = !foreign
rename make desc
gen make = word( desc,1)
tab make
xi3: regress mpg e.make
unab varlist : _I*
local test : subinstr local varlist " _" " + _", all
lincom -1*(`test')
sort make
by make : keep if _n == 1
keep make foreign
mata
b = st_matrix("e(b)")'
// remove constant
k = rows(b) - 1
b = b[1::k,1]
// add dropped level
b = st_numscalar("r(estimate)") \b
V = st_matrix("e(V)")
se = diagonal(cholesky(diag(V)))
// remove constant
se = se[1::k,1]
// add droped level
se = st_numscalar("r(se)") \ se
df = st_numscalar("e(df_r)")
ci = b :- invttail(df,0.025):*se, b :+ invttail(df,0.025):*se
idx = st_addvar("float",("b","lb","ub"))
st_store(.,idx,(b,ci))
end
egen xaxis = axis(b), label(make)
gen lb_mean = - invttail(e(df_r), 0.025)
gen ub_mean = invttail(e(df_r), 0.025)
gen mean = 0
forvalues mpg = 0(10)40 {
local diff = `mpg' - _b[_cons]
local newsc `"`newsc' `diff' "`mpg'" "'
}
local min = -_b[_cons]
sum ub, meanonly
local max = r(max)
twoway rarea lb_mean ub_mean xaxis, color(gs14) || /*
*/ line mean xaxis, clstyle(solid) ||/*
*/ rcap ub lb xaxis || /*
*/ rcap ub lb xaxis, /*
*/ lpattern(blank) yaxis(2) || /*
*/ scatter b xaxis if foreign == 1, /*
*/ mcolor(red) msymbol(O) || /*
*/ scatter b xaxis if foreign == 0, /*
*/ mcolor(blue) msymbol(O) /*
*/ xlab(1/23,valuelabel angle(45)) /*
*/ xtitle("") /*
*/ ytitle("mpg", axis(1))/*
*/ ytitle("deviation from grand mean", axis(2)) /*
*/ yscale(range(`min' `max') axis(1)) /*
*/ yscale(range(`min' `max') axis(2)) /*
*/ ylab(`newsc', axis(1)) /*
*/ ylab(-20(10)20, axis(2)) /*
*/ xsize(7) /*
*/ legend(order(- "Grand mean" /*
*/ 1 "95%" "conf. int." /*
*/ 2 "mean" /*
*/ - " " /*
*/ - "By make" /*
*/ 3 "95%" "conf. int." /*
*/ 4 "foreign" /*
*/ 5 "domestic" ) /*
*/ position(4) cols(1))
*------------------ end example ---------------
Hope this is useful,
Maarten
-----------------------------------------
Maarten L. Buis
Department of Social Research Methodology
Vrije Universiteit Amsterdam
Boelelaan 1081
1081 HV Amsterdam
The Netherlands
visiting address:
Buitenveldertselaan 3 (Metropolitan), room Z434
+31 20 5986715
http://home.fsw.vu.nl/m.buis/
-----------------------------------------
___________________________________________________________
Want ideas for reducing your carbon footprint? Visit Yahoo! For Good
http://uk.promotions.yahoo.com/forgood/environment.html
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
