gistfile1.txt

library(lme4)
library(rethinking)
library(foreign)
library(merTools)
library(dplyr)
set.seed(1234)

setwd("~/Desktop/")
BIGN <- 10000
df <- data_frame(x=rep(0:1, BIGN/2), id=(1:BIGN)%%50) %>% 
  group_by(id) %>% 
  mutate(u=rnorm(1)) %>% ungroup() %>% 
  mutate(
    y=rnorm(BIGN)+x+u,
    y_noreffs=rnorm(BIGN)+x+rnorm(BIGN)
  )

write.dta(df, "test.dta")
df %>% head
df %>% summary


m1 <- lmer(y_noreffs~x+(1|id), data=df)
m2 <- lmer(y~x+(1|id), data=df)


newd <- data_frame(x=0:1, id=1)


# raw means
df %>% group_by(x) %>% do(as.data.frame(t(quantile(.$y, probs=c(.5, .975,  .025)))))


# compare with  different predictions from stata
predictInterval(m1, n.sims=1e4, newd, level=.95)

# > predictInterval(m1, n.sims=1e4, newd, level=.95)
# fit      upr       lwr
# 1 -0.03564956 2.757970 -2.803637
# 2  0.93923557 3.669108 -1.729639



# use test
# mixed y i.x||id:
# margins x

# |            Delta-method
# |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
# -------------+----------------------------------------------------------------
#  x |
#  0  |   .0096216   .0141188     0.68   0.496    -.0180507    .0372938
#  1  |   .9752824   .0141188    69.08   0.000     .9476101    1.002955
# ------------------------------------------------------------------------------