Giving R the strengths of Stata

Useful information.

Generate predicted confidence intervals from logistic regression

using Stata.

Performance difference between Stata and R

With respect to multinomial logit model, the performance difference between the two packages are quite large, based on this post.

GEE using Stata vs. R

I am running GEE logistic regression model for my fetal loss paper. As usual, I compare results between Stata and R and make sure they are consistent. To my surprise, the models assuming independent correlation structure give similar results but the models assuming exchangeable correlation structure give drastically different results.

It turns out that there is only one woman in my sample who reported a total number of eleven pregnancies (all others reported ten or less) and the presence of this single observation had huge influence on the algorithm used in R but not the one used in Stata. After excluding this single observation, the two sets of results look identical.

How to get marginal effects with "mvrs"

I found this very helpful post by Ryan sometime ago about how to get marginal effect from estimated models with "mvrs" command, very helpful.

Stata 12

There are two things I like Stata 12 most: (1) automatic memory management, and (2) a decent program editor for Linux (I know the editor has been available on the Windows platform since version 11, but I am a Linux person). The updated "margin" commands also look good.

Replicating multilevel analysis Mplus results

Since I am relatively unfamiliar with Mplus, The first step toward building a complicated model is to replicate some simpler results using statistical packages I have more knowledge on such Stata and R, then build it up piece by piece. 

As the effort to understand the multilevel mixture modeling capacity of Mplus, I have decided to use the Mplus sample data set, "ex10.1.dat", to try alternative model specifications in both Mplus, Stata and R. I am going to post a series of code and results. 

As the first step, I try this really simple Mplus code: 

-----------------------------------------------------------

DATA: FILE IS ex10.1.dat;

VARIABLE: NAMES ARE y x1 x2 w class clus;

USEVARIABLES = y x1 x2;

WITHIN = x1 x2;

CLUSTER = clus;

ANALYSIS: TYPE = TWOLEVEL;

STARTS = 0;

MODEL:

%WITHIN%

y ON x1 x2;

%BETWEEN%

y ON ;

-----------------------------------------------------------

which gives exactly the same results as the following Stata code: 

-----------------------------------------------------------

xtreg y x1 x2, i(clus) mle

-----------------------------------------------------------

How can I graph the results of the margins command?

A useful post from UCLA ATS Computing explaining how to plot the results from the "margins" command. This one explains how to explain interaction results. More generally, the "Stata Frequently Asked Questions" section is also helpful.

Post-estimation results manipulation: Stata vs. R

Stata's new "margins" command is as powerful as "Zelig" package for R for the purpose of conducting post-estimation manipulation of the results (e.g. producing predicted probabilities, confidence intervals, etc.). While "margins" uses delta method, "Zelig" uses predictive simulation. The best thing is that the results produced by one package can be checked against the other to ensure the consistency.

Labelbook

The Stata command "labelbook" can be used to create output that shows both the raw values and value labels side by side, a feature very helpful dealing with second-hand data sources.

Estimation of quantile treatment effects with Stata

Here is a very useful article published on Stata Journal.

Handbook on Impact Evaluation: Quantitative Methods and Practices

Here is a the book, did not receive very positive review from the Stata community though.

Some Fortran code

Here are some Fortran code for Stata-related operations.

Stata 11

Got my Stata 11 today, was a little disappointed because the new program editor feature is not available on Linux platform.

I like the new multiple imputation feature though.

Generalized linear mixed effect model problem

I am trying to compare cohort difference in infant mortality using generalized linear mixed model. I first estimated the model in Stata:

xi:xtlogit inftmort i.cohort, i(code)

which converged nicely:


Fitting comparison model:


Iteration 0:   log likelihood = -1754.4476
Iteration 1:   log likelihood = -1749.3366
Iteration 2:   log likelihood = -1749.2491
Iteration 3:   log likelihood = -1749.2491


Fitting full model:


tau =  0.0     log likelihood = -1749.2491
tau =  0.1     log likelihood = -1743.8418
tau =  0.2     log likelihood = -1739.0769
tau =  0.3     log likelihood = -1736.4914
tau =  0.4     log likelihood = -1739.5415


Iteration 0:   log likelihood = -1736.4914  
Iteration 1:   log likelihood = -1722.6629  
Iteration 2:   log likelihood = -1694.9114  
Iteration 3:   log likelihood = -1694.6509  
Iteration 4:   log likelihood =  -1694.649  
Iteration 5:   log likelihood =  -1694.649  


Random-effects logistic regression              Number of obs      =     21694
Group variable: code                            Number of groups   =     10789


Random effects u_i ~ Gaussian                   Obs per group: min =         1
                                                               avg =       2.0
                                                               max =         9


                                                Wald chi2(2)       =      8.05
Log likelihood  =  -1694.649                    Prob > chi2        =    0.0178


------------------------------------------------------------------------------
    inftmort |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
  _Icohort_2 |  -.5246846   .1850328    -2.84   0.005    -.8873422   -.1620269
  _Icohort_3 |  -.1424331    .140369    -1.01   0.310    -.4175513     .132685
       _cons |  -5.214642   .1839703   -28.35   0.000    -5.575217   -4.854067
-------------+----------------------------------------------------------------
    /lnsig2u |   .9232684   .1416214                      .6456956    1.200841
-------------+----------------------------------------------------------------
     sigma_u |   1.586665   .1123528                      1.381055    1.822885
         rho |   .4335015   .0347791                      .3669899    .5024984
------------------------------------------------------------------------------
Likelihood-ratio test of rho=0: chibar2(01) =   109.20 Prob >= chibar2 = 0.000

Then I tried the same model using lme4:

m2 <- glmer(inftmort ~ (1|code) + as.factor(cohort), family=binomial, data=d)

I got:


Warning message:
In mer_finalize(ans) : false convergence (8)

And the results are quite different from what I got from Stata. I tried to estimate the model using "glmmML" and also got into trouble. This time the error message is:


Warning message:
In glmmML.fit(X, Y, weights, cluster.weights, start.coef, start.sigma,  :
  Hessian non-positive definite. No variance!

I then tried MCMCglmm:


prior=list(R=list(V=1, nu=0, fix=1), G=list(G1=list(V=1, nu=0)))
m2 <- MCMCglmm(inftmort ~ 1 + as.factor(cohort), random=~code, family = "categorical", data=d, prior=prior)

It seems to be working and produced estimates that are comparable to what Stata produced (not identical, of course):


Iterations = 3001:12991
Thinning interval = 10 
Number of chains = 1 
Sample size per chain = 1000 


1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:


                      Mean     SD Naive SE Time-series SE
(Intercept)        -5.7145 0.1805 0.005708        0.02295
as.factor(cohort)2 -0.5633 0.1788 0.005653        0.02194
as.factor(cohort)3 -0.1888 0.1471 0.004653        0.01912


2. Quantiles for each variable:


                      2.5%     25%     50%      75%   97.5%
(Intercept)        -6.0464 -5.8324 -5.7251 -5.61120 -5.2817
as.factor(cohort)2 -0.8977 -0.6985 -0.5551 -0.43084 -0.2371
as.factor(cohort)3 -0.4614 -0.2947 -0.1867 -0.09644  0.1232 

This is puzzling.

Event History Analysis with Stata

Book web site here.

Some Stata tips

Here are some useful Stata tips.

Create LaTeX table for descriptive statistics using Estout

I want to create a table of descriptive statistics. To be more specific, The table should contains 6 rows and 4 columns, something like this:

Descriptive Statistics

-----------------------------------------------------------

Cohort A Cohort B Cohort C

-----------------------------------------------------------

Urban Residence 0.263 0.313 0.264

Years of Schooling 5.585 7.280 7.592

Ethnic Majority 0.928 0.918 0.915

Want Children 0.990 0.992 0.991

Ever Pregnant 0.991 0.987 0.992

-----------------------------------------------------------

Observations 3553 2450 4694

-----------------------------------------------------------

The official document for Estout provides some clues but does not make it crystal clear. Here are the steps I followed:

1. Generate an arbitrary variable;
2. Run regression, separately for each cohort;
3. Collect the results;
4. Make the table

Here are the code:

eststo clear

gen y=uniform()

quietly regress y urban eduy han want_birth preg if cohort==1, noconstant

estadd mean

eststo d1

quietly regress y urban eduy han want_birth preg if cohort==2, noconstant

estadd mean

eststo d2

quietly regress y urban eduy han want_birth preg if cohort==3, noconstant

estadd mean

eststo d3

esttab d1 d2 d3, cells(mean(fmt(%8.3f))) label nonotes nonumber nodepvars mtitles("Cohort A" "Cohort B" "Cohort C") width(1\hsize) title(Descriptive Statistics)

Using Stata with a make file

This blog describes how to run Stata via make files.

This one describes how to use make with LaTeX.

Using the "foreign" package for data conversion

I was in a rush to convert a SPSS data into Stata format. Somehow my Stattransfer v.8 for Linux was lost and I did not want pause my work and go back to Windows just to get this one file converted. So fire Emacs+ESS+R, load the "foreign" package, did the file conversion, and it worked!