3
$\begingroup$

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.

Improve this question
$\endgroup$
4
  • $\begingroup$ suppose you've named your model g. fixef(g) should give you the fixed effects coefficients. Also if you type summary(g) it should give you the fixed effects estimates along with their asymptotic standard errors and p-values. To test the random effect variances, you'll have to "make your own" hypothesis testing. I suggest the likelihood ratio test (which is known to be conservative when testing the null hypothesis that a variance parameter is 0). $\endgroup$
    – Macro
    Commented Jul 6, 2011 at 14:06
  • $\begingroup$ Are you sure you're coding correctly the factor variable? I assume you tried to check that in fact there are 100 levels in R? $\endgroup$
    – Manoel Galdino
    Commented Jul 6, 2011 at 14:10
  • $\begingroup$ @Macro : yes fixef, coef, etc. works but summary(g) doesn't give me the p values only T values. I would like to have the wald Z and the corresponding p values if possible? $\endgroup$
    – DBR
    Commented Jul 6, 2011 at 14:46
  • $\begingroup$ @Manoel Galdino : the factor variable is subject id, and there are x observations from each subject id, so yes there are 100 and stata sees it as 100 $\endgroup$
    – DBR
    Commented Jul 6, 2011 at 14:47

2 Answers 2

Reset to default
6
$\begingroup$

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])
Improve this answer
$\endgroup$
6
  • 1
    $\begingroup$ @Daniel Ad 1. What do you mean by "simply | independent5"? Please check ?lmer or the forthcoming lme4-book. Ad 2. This is not possible if your data is identical in R and Stata. Did you check the frequencies? Can you provide a subset of your data? Ad 3. Sure, you could use lme(), e.g. summary(lme(y ~ x1, random = ~ 1|school, data=dfr)). $\endgroup$
    – Bernd Weiss
    Commented Jul 6, 2011 at 20:25
  • 1
    $\begingroup$ Ad 2. Thanks for the data, that really helps. See my results here. As you can see in Stata as well as in R: 99 groups. Ad 3. The R package is called nlme, the function you are interested in is called lme() (linear mixed-effects models). Or are you interested in nonlinear mixed-effects models? $\endgroup$
    – Bernd Weiss
    Commented Jul 6, 2011 at 20:58
  • 1
    $\begingroup$ @ Daniel: As you can see, the test statistics are identical, e.g. for videocond: $-1.88$ but R assumes a t- and Stata a stand. normal distribution. I would probably go with the (more conservative) t-distribution... but that's just me :-) To obtain the (two-tailed) p-values in R, just run (1-pnorm(abs(-1.88)))*2 which gives $p=0.06010808$ (z-values, i.e. assuming a stand. normal distribution). $\endgroup$
    – Bernd Weiss
    Commented Jul 6, 2011 at 21:33
  • $\begingroup$ this is extremely helpful thank you! one final question: when specifying further levels do I use the notation 1 | independen5 / independent6 / etc) where 5 is nested within 6 or is it the other way round? this is my final question and I'd like to thank you again for all the help $\endgroup$
    – DBR
    Commented Jul 6, 2011 at 22:11
  • 1
    $\begingroup$ @ Daniel: My last comment :-) It's sth like (1|level3/level2) or (1|level3:level2) + (1|level3). You might want to read Examples from Multilevel Software Comparative Reviews (esp. sec 3) $\endgroup$
    – Bernd Weiss
    Commented Jul 6, 2011 at 23:22
2
$\begingroup$

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?

Improve this answer
$\endgroup$
6
  • $\begingroup$ thanks for your help, yes this did change the model a little bit, but the number of groups is still 99 and it seems to handle constrast codes (coded as 0,1) differently, because those predictions are bit different. do you know why this might be? $\endgroup$
    – DBR
    Commented Jul 6, 2011 at 13:43
  • $\begingroup$ thanks for your help, yes this did change the model a little bit, but the number of groups is still 99 and it seems to handle constrast codes (coded as 0,1) differently, because those predictions are bit different. do you know why this might be? $\endgroup$
    – DBR
    Commented Jul 6, 2011 at 14:07
  • $\begingroup$ no I am sorry, I am not sure here. :( If I find something I´ll let you know $\endgroup$
    – Jens
    Commented Jul 6, 2011 at 14:24
  • $\begingroup$ thanks! also could you tell me exactly the differences between (1|independent5), 1+(1|independent5) and just simple |independent5? also if I want to include multiple levels, say in what order do I specify them? thank you $\endgroup$
    – DBR
    Commented Jul 6, 2011 at 14:31
  • $\begingroup$ if you have (1|independent5) and 1+(1|independent5) then you put a 1 into the equation. why? I do not know. what does happen if you do that? will the 1 be another parameter? I do not thnk so. does the intercept change? if you want to include multiple nested levels you can do (1|independent5/factor6) then I think 5 will be nested in 6 or the other way around. $\endgroup$
    – Jens
    Commented Jul 6, 2011 at 14:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.