Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
puzzling_mlmpreds.R
library(parallel)
cl <- makeCluster(detectCores())
results <- parLapply(cl, 1:5000, function(placeholder) {
library(tidyverse)
library(lme4)
inv_logit <- function(x) exp(x) / (1 + exp(x))
j <- 400
id <- 1:j
g00 <- -1.5
g01 <- .2
u0sd <- 1
ij <- 3
x <- rbinom(j, 1, .5)
u0 <- rnorm(j, 0, u0sd)
yjlogit <- g00 + g01 * x + u0
dat <- data.frame(id, g00, g01, x, u0, yjlogit)
dat <- dat[rep(seq_len(nrow(dat)), each = ij), ]
dat$y <- rbinom(nrow(dat), 1, inv_logit(dat$yjlogit))
mod <- glmer(y ~ x + (1 | id), dat, binomial)
# replicate distribution of 1, 2, and 3 obs like in data
dat2 <- dat[dat$id < 41, ]
dat2 <- dat[dat$id > 40 & dat$id < 161, ] %>%
group_by(id) %>%
slice(1:2) %>%
bind_rows(dat2, .)
dat2 <- dat[dat$id > 160, ] %>%
group_by(id) %>%
slice(1) %>%
bind_rows(dat2, .)
mod2 <- glmer(y ~ x + (1 | id), dat2, binomial)
# is just 1 obs problem? just 2 and 3
dat3 <- dat[dat$id < 41, ]
dat3 <- dat[dat$id > 40, ] %>%
group_by(id) %>%
slice(1:2) %>%
bind_rows(dat3, .)
mod3 <- glmer(y ~ x + (1 | id), dat3, binomial)
# compare performance
out <- c(
# mod
u0_3obs = diff(c(u0sd, as.data.frame(VarCorr(mod))$sdcor)),
g00_3obs = diff(c(g00, fixef(mod)[[1]])),
g01_3obs = diff(c(g01, fixef(mod)[[2]])),
pred0_3obs = inv_logit(fixef(mod)[[1]]),
pred1_3obs = inv_logit(sum(fixef(mod))),
# mod2
u0_123obs = diff(c(u0sd, as.data.frame(VarCorr(mod2))$sdcor)),
g00_123obs = diff(c(g00, fixef(mod2)[[1]])),
g01_123obs = diff(c(g01, fixef(mod2)[[2]])),
pred0_123obs = inv_logit(fixef(mod2)[[1]]),
pred1_123obs = inv_logit(sum(fixef(mod2))),
# mod3
u0_23obs = diff(c(u0sd, as.data.frame(VarCorr(mod3))$sdcor)),
g00_23obs = diff(c(g00, fixef(mod3)[[1]])),
g01_23obs = diff(c(g01, fixef(mod3)[[2]])),
pred0_23obs = inv_logit(fixef(mod3)[[1]]),
pred1_23obs = inv_logit(sum(fixef(mod3)))
)
})
stopCluster(cl)
results <- as.data.frame(do.call(rbind, results))
results %>%
mutate(iter = 1:nrow(.)) %>%
gather(... = -iter) %>%
separate("key", c("param", "mod")) %>%
spread(param, value) %>%
group_by(mod) %>%
summarise_at(vars(g00:u0), mean)
results %>%
mutate(iter = 1:nrow(.)) %>%
gather(... = -iter) %>%
separate("key", c("param", "mod")) %>%
filter(param %in% c("pred0", "pred1")) %>%
ggplot(aes(x = value, fill = param)) +
geom_density(alpha = .7) +
facet_wrap(~ mod)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment