I've used Stata 12 to estimate a hierarchical logit model (xtmelogit) with random intercepts. The outcome is incarceration (1=yes) for a series of convicted criminal offenders. There are a slew of interaction terms and I'm trying to understand these by looking at the predicted odds at different levels of X while holding other values as observed.
The margins command won't work here since it cannot include the random effects portion of the model. So I'm doing it by hand. After setting my values where I want them I calculate the fixed estimate (the logodds for the fixed portion of the model) and add that to the estimated random intercept to get a point estimate
The code looks something like this:
predict randeffect, reffects
to fit with a random intercept and
predict randse, reses
to get the standard error for the random intercept
The following code is done within a
forvalues loop where I vary the value of some covariates and then estimate the fixed effect at each level... what I've posted is a simplification that demonstrates the principle without the complexity of the
predict fixedeffect, xb
which gets the point estimate for the fixed portion of the model (log odds) and
predict fixedse, stdp
which gets the standard error for the fixed portion of the model (log odds).
To get the point estimate at each level I add the fixed effect to the
randeffect and take the exponent of it to get the predicted odds. That is,
gen odds = exp(randeffect + fixedeffect).
And now the tricky part...
To get the confidence interval my first instinct was to do all the addition with the log odds and then take the exponent to get the upper and lower bounds. This results in very asymmetrical bounds because the exponent is (obviously) not a linear transformation of the values. The code is as follows:
gen lowerbound = exp((randeffect + fixedeffect) - (1.96 * (fixedse + randse)))
The above makes sense but the confidence intervals are not symmetrical. However, I see in the output of the logit command that all odds ratios have symmetrical confidence intervals.
Alternatively, I can calculate the confidence interval and then take the exponent and then add and subtract that to the calculated odds ratio. This gives me nice a nice symmetrical confidence band around my point estimate when I graph it. This code looks like
gen lowerbound = odds - (exp(1.96 * randeffect + fixedeffect))
So which way is right? Do I exponentiate first and then add and subtract? Or do I add and subtract and then exponentiate?
Thanks for any insight you can provide.
EDIT: I wanted to follow up with a caveat for anyone considering doing the same in Stata. First, Stata can make predictions based on the fixed part of the model which basically ignores the intercept. Since the fixed effects are the same for all clusters, it's like estimating a model where all level 1 units are in the same cluster. It's like an apples to apples comparison. I've chosen this approach to understand the fixed effects and then separately consider variation in the intercepts across my level 2 units. It's much simpler this way.
Second, from what I can tell the method described in my original post does not correctly model interaction effects which were the focus of the analysis. The problem, I think, is that when you vary the value for X2 to see how the effect of X1 on Y changes, you must also change/reset the value of X2*X1 (the interaction term). This sounds simple enough but the problem quickly becomes a thorny mess of code when your X1 variable is involved in more than one interaction (i.e. you have X2*X1 as well as X3*X1 in the model). So this is just something to be watchful for if you, after reading this, get the bright idea to take this route.
I personally gave up on it, preferring instead to use the margins command with the predict (mu fixedonly) option
margins x1, dydx(x2) predict (mu fixedonly)
But I sincerely thank those of you who contributed to the effort!