Multilevel regression using lmer function in R and Stata I have a large dataset and have performed a multilevel regression in Stata, the model is the following:
xtmixed dependent independen1 independent2 independent3  independent4    || independendt5:
So there is one grouping factor: independent5
In R I did the following:
lmer(dependent ~ independent1 + indepdendent2 + independent3 + independent4 + 1 | independent5, REML=TRUE) 

A few questions: is this identical?
Stata output gives me number of groups, 100, while R gives number of groups as 99.
Furthermore, the variances and standard deviations are not the same.
Also I would like to know how to obtain p values and coefficients from the R ouput.
I have done fixef(model) and ranef(model) but when I do coef(model) it says: Error in coef(model) : unable to align random and fixed effects
Also Stata only gives me 1 coefficient for each predictor while the R fixef(model) gives one for each group. 
So someone familiar with both Stata and R could help me with this.
 A: Number of groups differ in Stata and R
Regarding your "99 vs. 100 groups problem": Are you really sure that your R and Stata dataset are identical? In Stata run summarize, in R run summary(yourDataFrame) and compare the results. 
Fitting varying intercept/slope models in Stata and R
@Jens already has pointed out how to write Stata's xtmixed model in R. Some comment on the difference between 1 + (1|independent5) and (1|independent5): There is no difference. 1 + (...) means that an intercept is included in the model; but this is the default.  
Gelman offers a nice overview of how to fit models in R and Stata, see this page and then The (current version of) the relevant pages in the book (see Section C.4)
Please find below some sample code how to estimate simple MLMs in Stata and R. The data and some information can be found here: http://dss.princeton.edu/training/Multilevel101.pdf (here is the dataset: http://dss.princeton.edu/training/schools.dta)
Stata:
clear
use c:/tmp/schools.dta

## varying intercept model
xtmixed y x1 || school:, reml

## varying intercept and slope (x1) model
xtmixed y x1 || school:x1, reml

R: 
library(foreign)
library(lme4)
dfr <- read.dta(file="c:/tmp/schools.dta")
head(dfr)

## varying intercept model
lmer(y ~ x1 + (1|school), data=dfr)

## varying intercept and slope (x1) model
lmer(y ~ x1 + (x1|school), data=dfr)

R's lmer does not report p-values
I am not a statistician and I cannot comment on lmer's behavior. You might want to read this post by D Bates. You also might be interested in this blog post: Linear mixed-effects regression p-values in R: A likelihood ratio test function 
Hopefully @Ben Bolker will comment on this...
Extracting information from R's *mer-class-objects
To extract 


*

*the fixed effect(s): fixef(myLmerObject)

*the random effect(s): VarCor(myLmerObject) (see str(VarCorr(la)) to better understand the internal structure, e.g. VarCorr(myLmerObject)[[1]][1])

A: I think the correct formula should be
model1 <- lmer(dependent ~ independent1 + indepdendent2 + independent3 + independent4 + (1|independent5), REML=TRUE)
the (1|independent5) makes the grouping factor. you can change the '1' to a factor that stores the order of your records.
hope that helps?
