Plotting the random+fixed effects of slopes after multilevel models are fit with `xtmixed`

[Louis Ashwall]

Hello all,

That's my first post in here, so let me know if I do something wrong. Well, my problem is the following.

I am fitting a linear multilevel model of the following form:

    Y = b0 * X0 + b1 * X1 + b2 * X2 + b3 * X3 + u
    b0 = e0
    b3 = e3

The units in the first level are 10000 students and the clusters of the second level are 30 schools.

So, what I am trying to do is to plot each of the 30 versions of `b3`, i.e. the random effects slope of each cluster. I also need to plot that if confidence intervals of any type.

In R, I know how to do it. Assuming the model fitted is saved in the `mymodel` object, one can get the random + fixed effects of a multilevel model in R as follows:

    #This will extract simulated random effects from the model:
    myranef <- ranef(arm::sim(mymodel, 1000))
    myranef <- as.data.frame(myranef)
    myranef <- myranef[,grep("X3",colnames(bb))]
   
    #This will extract simulated fixed effects from the model:
    myfixef <- as.data.frame(fixef(arm::sim(mymodel, 1000)))
    names(myfixef) <- names(fixef(mymodel))
    
    #Sum random and fixed effects
    for(c in 1:ncol(myranef)){
      myranef[,c] <- myranef[,c] + myfixef[,"X3"]
    }

In the end, each column of `myranef` represents one of the 30 schools and each line is one of the 1000 simulation-iterations. From then on, it's easy to make the plot using the mean or median of each col and the desired percentiles of each col.

My problem is how to achieve a similar thing in Stata. I guess I should start by doing the following:

    predict re*, reffects  //  linear predictor for the random effects
    predict fe, xb         // linear predictor for the fixed effects

However, then I don't know how to:

1) get the combined version of random and fixed effects for each of the 30 schools
2) get a version with confidence intervals of some sort


Would anyone be able to kindly shed some light on how to proceed in Stata similarly to what I have in R?

[Stephen Jenkins]

You can get BLUPs of random effects using predict. See help mixed postestimation

[Louis Ashwall]

Hi Dr. Stephen Jenkins, thanks for your quick reply.

Sure thing, the trick is in using the "reses" option. So it seems to be as easy as to do the following:
predict re*, reffects // linear predictor for the random effects predict re_se*, reses //standard errors of the random effects predict fe, xb // linear predictor for the fixed effects However, my main problem is how to get the final random + fixed effects and its SE or confidence interval of the specific variable of interest.

[Clyde Schechter]

You just add them.

predict fe, xb
predict re, reffects reses(reses)
gen predicted = fe + re
gen lower_ci = predicted - 1.96*reses
gen upper_ci = predicted + 1.96*reses

[Louis Ashwall]

Hi Dr. Schechter, many thanks for your answer.

So, that is the point I am getting trouble with. Let me make it clearer. Two things are in my way here.

First, in my case I have more than one random slope, therefore "re" has to be a list of variables, while "fe" not necessarily:

predict fe, xb
predict re*, reffects

So the issue is: how to sum fe + re* or, more precisely, for one of the variables of interest (in my example, X3?

Second, I think the following is wrong:

predict re, reffects reses(reses)

Because it gives me the error "option reses() incorrectly specified"

[Clyde Schechter]

So, I take it the issue is that you have random slopes that you want to take into account in building up a predicted value for each school. Here's a concrete example of doing this using a two-level random-slopes model fitted to one of the Stata manual data sets:

use "http://www.stata-press.com/data/r14/pig.dta", clear


// TWO LEVEL RANDOM SLOPES MODEL
mixed weight week || id: week

predict xb, xb
predict u*, reffects reses(u_se*)

des u*

gen predicted_value = xb + u2 + u1*week
gen predicted_value_se = sqrt(u_se2^2 + (u_se1*week)^2)

Note that this does not give a unique value for each unit at the second level (school in your case) because with a random slope, the predicted values also depend on the variable(s) that have random coefficients.

[Louis Ashwall]

Great, many thank Dr. Schechter. That was helpful.

Quick question, though. Does it work the same for 'xtmixed'? Because I getting the following error after running the model with xtmixed in Stata 13:

. predict u*, reffects reses(u_se*)
option reses() incorrectly specified

[Clyde Schechter]

-xtmixed- is the old name for the -mixed- command. They are exactly the same. The issue here is Stata version. I'm using 14.1, the current version. I think the syntax for getting random effects standard errors in Stata 13.1 was different, but I don't remember it. You should check -help xtmixed postestimation- on your computer, then click on the link to -predict- and see what it says there.