R/order_boot_glmer.R

In fartist/orderboot:

# install.packages("lmeresampler")
library(lmeresampler)
library(lme4)
mixed <- glmer(incidence/size ~ period + (1|herd),
               weights=size, data=cbpp, family=binomial)
fixef(mixed)
FUN <- function(fit) {
  return(fixef(fit))
}
# cbpp$herd
# result <- bootMer(mixed, FUN, nsim = 100)
# 
# ?gl
# dat<-data.frame(x=runif(100,-2,2),ind=gl(n=10,k=10))
# dat$y<-1+2*dat$x+rnorm(10,0,1.2)[dat$ind]+rnorm(100,0,0.5)
# m<-lmer(y~x+(1|ind),dat)
# b_par<-bootMer(x=m,FUN=fixef,nsim=200)
# boot.ci(b_par,type="perc",index=1)
# 
# pl
# 
# 
# library(ltm)
# dsc <- descript(Abortion)
# dsc
# 
# 
# abo_long <- 
#   Abortion %>% mutate(person = 1:n()) %>% 
#   pivot_longer(cols = starts_with("Item"),
#                           names_to = "item",
#                           values_to = "score")
# 
# 
# mixed <- glmer(score ~ 0 + item + (1|person), 
#                data=abo_long, 
#                family=binomial, nAGQ = 20)
# beta_hat <- fixef(mixed)
# fit <- mixed
# r = .1
# N = dim(abo_long)[1]
# 
# beta_order <- as.double(sort(beta_hat))
# c_hat_order <- map(beta_order, ~ (1 - N ^ (r - .5)) * (.x - beta_hat))
# FUN <- function(fit) {
#   return(fixef(fit))
# }
# timestamp()
# result <- bootMer(mixed, FUN, nsim = 500)
# timestamp()
# n_boot = 500
# results <-   
#   map(1:n_boot, function(index) 
#   {map2_dbl(c_hat_order,
#             1:K,
#             ~ (.x + result$t[index, ]) %>% 
#               sort_unique %>%
#               `[`(.y))} - beta_order) %>%
#   do.call("rbind", .)
# (est  <-  sweep(-results, 2, -beta_order) %>% colMeans())
# 
# tibble(beta_order, est) %>% knitr::kable("latex")
# beta_order_boot_sample <-
#   sweep(-results, 2,-beta_order) %>% as.data.frame()
# M <-  ((((
#   beta_order_boot_sample - th
# ) >= 0) %>% colMeans()) > .05) %>% sum()
# cr  <-  map(beta_order_boot_sample, get_ci_1) %>% map2_lgl(
#   .,
#   1:K,
#   ~ is_cover_1(lower = .x[1], point = beta[.y])
# )
# 
# 
# # do it for each row of result (replace beta_hat_star with each row of beta_hat_star)
# map2_dbl(c_hat_order,
#          1:K,
#          ~ (.x + beta_hat_star) %>% 
#            sort_unique %>% 
#            `[`(.y)) - beta_order
# ) %>% do.call("rbind", .)
# 
# 
# result$t
# map2(c_hat_order,
#          1:K,
#          ~ (.x + result$t[2, ]) %>% 
#            sort_unique %>% 
#            `[`(.y)) 
# 
# 
# 
# # - beta_order
# 
# est  <-  sweep(-result$t, 2, -beta_order) %>% colMeans()
# beta_order_boot_sample <-
#   sweep(-result$t, 2,-beta_order) %>% as.data.frame()
# M <-  ((((
#   beta_order_boot_sample - th
# ) >= 0) %>% colMeans()) > .05) %>% sum()
# cr  <-  map(beta_order_boot_sample, get_ci_1) %>% map2_lgl(
#   .,
#   1:K,
#   ~ is_cover_1(lower = .x[1], point = beta[.y])
# )

library(lme4)
library(tidyverse)
# illustration of bias
J = 30
I = 10
sim_rim <- function(I = 10, J = 30, beta = rep(0, 10), u_sd = 1) {
  person = rep(1:J, each = I)
  u <- rnorm(n = J, 0, u_sd) %>% rep(., each = I)
  x = model.matrix(~ 0 + factor(1:10))
  X = do.call("rbind", replicate(J, x, simplify = FALSE))
  p <- 1 / (1 + exp(-(X %*% beta + u)))
  y <- map_dbl(p, ~ rbinom(n = 1, size = 1, .x))
  return(data.frame(item = factor(rep(1:I, J)), y, person))
}

original_mat <- matrix(nrow = 100, ncol = 10)
boot_mat <- matrix(nrow = 100, ncol = 10)
original_mat[1:20, ] %>% colMeans()
boot_mat[1:20, ] %>% colMeans()
# rbind(original_mat, boot_mat) %>% write.csv("mirt_boot_res.csv")
n_boot = 200

iter_i = 1
for (iter_i in 1:200){
  data <- sim_rim(J = 30)
  fit1 <-
    glmer(y ~ 0 + item + (1 | person), family = "binomial", data = data, nAGQ = 20)
  original_mat[iter_i,] <- fixef(fit1) %>% sort
  timestamp()
  result <- bootMer(fit1, FUN, nsim = n_boot)
  timestamp()
  
  K = 10
  beta_hat <- fixef(fit1)
  r = .01
  N = dim(data)[1]
  beta_order <- as.double(sort(beta_hat))
  c_hat_order <-
    map(beta_order, ~ (1 - N ^ (r - .5)) * (.x - beta_hat))
  results <-
    map(1:n_boot, function(index)
    {
      map2_dbl(c_hat_order,
               1:K,
               ~ (.x + result$t[index,]) %>%
                 sort_unique %>%
                 `[`(.y))
    } - beta_order) %>%
    do.call("rbind", .)
  est  <-  sweep(-results, 2,-beta_order) %>% colMeans() %>% sort
  boot_mat[iter_i,] <- est
}

res <- read_csv("C:/Users/ifeng/Google Drive/Jingshen/res.csv", 
                +     col_types = cols(X1 = col_skip()))
res$X1 <- res$X1_1

res %>% select(X1_1:X10)