Get estimates of group levels in a multilevel model in Stata

Appendix C in Data Analysis Using Regression and Multilevel/Hierarchical Models of professor Gelman and Hill has documented commands and codes for fitting multilevel models in various software (R, SAS, SPSS, STATA and ADMB).

I once took a class called analysis of time series cross sectional data. In class, we have a demo of how to fit a linear mixed model (varying intercepts and slopes). The command is xtmixed. I remember instructor told us one pity thing about stata's xtmixed model is that it cannot give us the estimates for group levels. He thus suggests us an alternative package in R. The package of which was used a lot in the aforementioned book is lme4. I took the instructor's word for granted.

But one day, I tried to challenge his word and did a research on xtmixed. Apparently, he was wrong. We can get the estimates for random effects in stata. Here is the sample code (the fake data is here):

* fit a varying intercept and slope model
xtmixed y x, || group: x

* get estimates for the group level
predict u0 u1, reffects
list u0 u1

Use Stat/Transfer inside of Stata

As a stata user, you may find out there are limited types of data stata can read in directly. Normally, we use Stat/Transfer to transfer data into stata format. Roger Newson wrote a stata package --stcmd-- that can call Stat/Transfer inside of stata.

To install the package, simply type:

ssc install stcmd

To use stcmd, you need to set up the program path for Stat/Transfer:

global StatTransfer_path `"c:\Program Files\StatTransfer9\st.exe"'

Here is the sample code to read in an Excel data with extension xls:

inputst test.xls

To export data from stata to R inside of stata:

outputst test.Rdata

Computing Sample mean, standard errors & design effects in Stata

Stata has a package for survey sampling: svy. Here are steps to get sample mean, standard errors and design effects in Stata.


1. Declare sample design (weight, primary sampling units, strata, etc)

*stata 8.2
svyset [pweight=factor], psu(psuid) strata(strataid)

*stata 9 or 10
svyset upm [pweight=factor], strata(strataid)

2. Compute sample mean, standard errors & design effects

*stata 8.2
svymean income
gen sigma = sqrt(e(N) * el(e(V_srs),1,1))
svymean income, by(sex)
gen sigma1 = sqrt(el(e(_N),1,1) * el(e(V_srs),1,1))
gen sigma2 = sqrt(el(e(_N),1,2) * el(e(V_srs),2,2))

*version 9
svy: mean income
estat sd
estat effect *design effect
svy: mean income, over(sex)
estat sd
estat effect

The design effect is correct in stata 8.2 using svymean. But it is wrong in stata 9 or 10 using svy: mean. We have to use estat effect to get correct design effect in stata 9 or 10.

Using WinBUGS from Stata

Stata outperforms R in its ability to handle large datasets. So when we want to run WinBUGS through R to deal with a large dataset, we have to reply on Stata to help us get coda statistics. Of course, if you know how to run WinBUGS through other software that can also deal with large datasets, you don't need to learn this. Skip this post.

I will use 8 schools as an example from Andrew Gelman et al. 2004. Bayesian Data Analysis. Here are the steps: (The data, the do file and the BUGS model)

1. Install WinBUGS from Stata package.

Type the following commands in the Stata command windows:

net from http://www.hs.le.ac.uk/research/HCG/winbugsfromstata
net describe winbugsfromstata
net install winbugsfromstata

Type help winbugs in the command windows. If you install the package successfully, you will see the help window of WinBUGS from Stata.

2. Read the data into Stata.

use d:\temp\schools.dta, clear

3. Check summary statistics of the data and some data management.

sum
scalar J=8
gen y=estime

4. Convert stata data to BUGS data

There are 5 commands to do this: wbarray, wbdata, wbscalar, wbstructure, wbvector.

wbdata, c(wbscalar, sc(J) f(%3.1f) + wbvector y sd, f(%3.1f)) noprint saving(d:\temp\data.txt,replace)

(NOTE: There is no line break if you type in the above commands)

You will see a file called data.txt in your working directory.

5. Prepare the initial values for 3 chains

* gen first inits
set seed 339487731
gen theta=invnormal(uniform())
scalar mu=invnormal(uniform())
scalar sigma=uniform()
wbdata, c(wbscalar, sc(mu sigma) f(%3.1f) + wbvector theta, f(%3.1f)) nop sav(d:\temp\inits1.txt, replace)

* gen second inits
set seed 3241341
replace theta=invnormal(uniform())
scalar mu=invnormal(uniform())
scalar sigma=uniform()
wbdata, c(wbscalar, sc(mu sigma) f(%3.1f) + wbvector theta, f(%3.1f)) nop sav(d:\temp\inits2.txt, replace)

* gen third inits
set seed 12341324
replace theta=invnormal(uniform())
scalar mu=invnormal(uniform())
scalar sigma=uniform()
wbdata, c(wbscalar, sc(mu sigma) f(%3.1f) + wbvector theta, f(%3.1f)) nop sav(d:\temp\inits3.txt, replace)

You will have three files contains initial values in your working directory: inits1.txt, inits2.txt, inits3.txt.

6. Prepare a WinBUGS script.

wbscript, sav(d:\temp\script.txt, replace) model(d:\temp\model.txt) data(d:\temp\data.txt) inits(d:\temp\inits.txt) coda(d:\temp\coda) burn(1000) update(5000) chain(3) set(mu sigma theta) dic log(d:\temp\log.txt) quit thin.updater(5) stats(*)

You will have a file called script.txt in your working directory. You may want to add thin.updater(5) stats(*) into your script.txt. Please check script.txt for the right position to put these two commands. Adding these two commands will help you do the analysis in R.

7. Run WinBUGS from Stata.

wbrun, sc(d:\temp\script.txt) winbugs(c:\Program Files\WinBUGS14\winbugs14.exe)

You will have coda1.txt, coda2.txt, coda3.txt and codaIndex.txt in your working directory. You can do the analysis in R using these coda statistics. Or do the simple analysis in Stata.

8. Read in the coda statistics and do the analysis.

wbcoda, root("d:/temp/coda") chain(3) multichain clear
wbstats mu sigma theta_1 theta_2 theta_3 theta_4 theta_5 theta_7 theta_8

There are other commands you can use to do the analysis. Type help winbugs in the command windows for details.

You can convert the coda files into bugs object in R. And do all the post analysis in R.

Create lag varaible in Stata

Q. How to create a lag variable in Stata?

# if dataset is not of time series, which is wierd, but anyway
gen varlag=var[_n-1]

# if dataset is of time series or cross-sectional or panel ("L." is the operator)

tsset timevar, panelid

gen varlag=l.var

or

by panelid: gen varlag=var[_n-1]