source("Zhang_simfuns.R")

library(lme4)
glmer_4.LL_fun <- function(...) {
  s <- simfun(...)
  s$int <- 1
  g0 <- glmer(y~x+(x|id),data=s,family=binomial)
  g1 <- glmer(y~0+x+(x|id)+offset(int),data=s,family=binomial)
  g2 <- glmer(y~1+(x|id)+offset(x),data=s,family=binomial)
  p1 <- anova(g0,g1)[["Pr(>Chisq)"]][2]  ## test of b1=1
  p2 <- anova(g0,g2)[["Pr(>Chisq)"]][2]  ## test of b1=1
  ## FIXME: do I trust these/are they right???
  coef0 <- fixef(g0)
  stder0 <- sqrt(diag(vcov(g0)))
  ## p1w <- 2*pnorm(abs((fixef(g0)[1]-1)/sqrt(vcov(g0)[1,1])),lower.tail=FALSE)
  ## p2w <- 2*pnorm(abs((fixef(g0)[2]-1)/sqrt(vcov(g0)[2,2])),lower.tail=FALSE)
  p1w <- 2*pnorm(abs((coef0[1]-1)/stder0[1]),lower.tail=FALSE)
  p2w <- 2*pnorm(abs((coef0[2]-1)/stder0[2]),lower.tail=FALSE)
  c(b1_LRT=p1,b2_LRT=p2,b1_Wald=unname(p1w),b2_Wald=unname(p2w))
}


nMCMC <- 500
set.seed(1001)
fn <- "Zhang_glmersim2.RData"
reslist <- list()
reslist$res_glmer_4.LL <- sfun3(glmer_4.LL_fun,nMCMC=nMCMC,savename="Zhang_glmersim2_save")
save("reslist",file=fn)