"Margins" after xt or mixed

[David Crow]

Dear Statalisters:

I am trying to graph random intercepts and slopes--BLUPs--for each of 50 subjects, observed on four periods each (N=50, T=4), after estimating a longitudinal growth model with "xtmixed." I would like all to appear on the same graph. (I used "xtmixed" because I'm still on Stata12.) Here's the code:

*Define Panel and Time Variables
xtset id time

*Longitudinal Growth Model
xtmixed readit time || id:time

*Graph Intercepts and Slopes
margins, at(time=(0(1)3)) predict(xbu) over(id)
marginsplot, x(time) plotdim(id)

Apparently, "margins" doesn't recover BLUPs for the random effects, as gllapred did with the "gllamm" package (for which quadrature is really slow, and maybe inaccurate, for moderately complex problems) or "ranef" does in the R package "lme4"--unless "margins" incorporates recovery of random effects in Stata13. I get this error message:

. margins, at(time=(0(1)3)) predict(xbu) over(id)
option xbu not allowed
r(198);

The graph I want would look something like the one I get with this code:

*OLS (Unpooled Model)
*Separate Intercepts and Slopes for Each Respondent
reg readit c.time##i.id

*Plots
quietly margins, at(time=(0(1)3)) over(id)

*Plot All Slopes on Single Graph
marginsplot, x(time) plotdim(id)

Any ideas on how I can do this with any of the "xt" or the "mixed" commands? Suggestions would be greatly appreciated.

Thanks,
David

[Clyde Schechter]

I think that to get what you are looking for you will have to get the estimates of the intercepts and slopes using -predict, reffects-. Then you can generate id-specific linear predictions at times 0 and 3, and then use -xtline, overlay- to get your plot.

But may I suggest a different approach. The kind of plot you're looking for is, to me, essentially unreadable, and I find it hard to extract useful information from it. I prefer to make a caterpillar plot. It's easy to do, and, more important, easy to read and understand:

predict intercept slope, reffects
predict intercept_se slope_se, reses
egen flag = tag(id)
foreach v of varlist intercept slope {
    egen `v'_rank  = rank(`v'), unique
    sort `v'_rank
    serrbar `v' `v'_se `v'_rank if flag, name(`v', replace)
}

Or if you are looking to see correlations between intercept and slope (and your time variable has been suitably centered for your purposes) you can just do a scatter plot of intercept vs slope.

[David Crow]

Thanks, Clyde, for all your suggestions. I wanted the graph I suggested as an exploratory technique (suggested by Singer and Willett in "Applied Longitudinal Data Analysis") to see if a random effects model is warranted, or if a pooled fixed effects model reasonably describes the data. The graph makes clear that there's considerable variation between individuals and that some sort of model to account for unit-level heterogeneity is called for.

I agree that the caterpillar plot is easier to read and understand (though it takes two graphs) and that a there are clearer ways than my graph of examining correlations between intercept and slope--so, thanks again.

I did find a really cheesy, low tech, poor man's solution to my problem:

*Model
xtmixed readit time || id:time

*Get Values Fitted with Random Intercepts and Slopes
predict yhat, fitted
sort id time
twoway (line yhat time, connect(ascending) lc(gray))

The trick is to connect the lines with the "connect(ascending)" option.

Best,
David

[Nate Breznau]

In order to get appropriate point estimates after an xtmixed command, I use predict and predict, reffects and then do the math by hand. Is there an easy way to get confidence intervals for these point estimates?