require(testthat)

# for sleepstudy, calls in here to lmer
library(lme4)
# for grad function

options(width = 500)
options(useFancyQuotes = FALSE)

### Data Used: sleepstudy
tolerance <- 1E-3
data("sleepstudy")
### Unweigted =====================================
sleepstudyU <- sleepstudy
sleepstudyU$weight1L1 <- 1
sleepstudyU$weight1L2 <- 1

context("The model runs")
test_that("The model runs", {
  wm0 <- mix(Reaction ~ Days + (1|Subject), data=sleepstudyU, weights=c("weight1L1", "weight1L2"))
  lme1 <- lmer(Reaction ~ Days + (1|Subject), data=sleepstudyU, REML=FALSE)
  expect_equal(wm0$lnl, -897.039321502613, tolerance=tolerance*897) # value from  lmer(Reaction ~ Days + (1 | Subject), data=sleepstudy, REML=FALSE)
  expect_equal(unname(c(wm0$coef)), unname(fixef(lme1)), tolerance=tolerance)
  expect_equal(unname(wm0$theta), unname(lme1@theta), tolerance=tolerance)
  expect_equal(coef(wm0),
               expected = getME(lme1, "fixef"),
               tolerance = tolerance)
  lmevars1 <- data.frame(summary(lme1)$varcor)$sdcor
  expect_equal(unname(wm0$vars),
               expected = unname(lmevars1)^2,
               tolerance = tolerance)
  # agrees with GLLAMM
  #source: logit_command.do
  gllamm_model1 <- c("(Intercept)" = 251.4051,
                     "Days"         = 10.46729,
                     "Residual"     = 1296.8768,
                     "Subject"      = 954.52789,
                     "lnl"          = -897.03932)
  expect_equal(unname(coef(wm0)),
               expected = unname(gllamm_model1[1:2]),
               tolerance = abs(tolerance))
  expect_equal(unname(wm0$vars),
               expected = unname(gllamm_model1[3:4]),
               tolerance = tolerance)
  expect_equal(unname(wm0$lnl),
               expected=unname(gllamm_model1[5]),
               tolerance=abs(tolerance))
  
  sleepstudy2 <- sleepstudy
  sleepstudy2$block <- rep(1:6,each=30)
  sleepstudy2$weight1L1 <- sleepstudy2$weight1L2 <- sleepstudy2$weight1L3 <- 1
  # change Reaction by the artifical block we just added
  sleepstudy2$Reaction <- sleepstudy2$Reaction + sleepstudy2$block*20
  names(sleepstudy2) <- gsub("weight1L","pwt",names(sleepstudy2))
  lmr <- lmer(Reaction ~ Days + (1|Subject) + (1|block), data=sleepstudy2, REML=FALSE)
  mm <- mix(Reaction ~ Days + (1|Subject) + (1|block), data=sleepstudy2, weights=c("pwt1", "pwt2", "pwt3"))
  expect_equal(unname(mm$lnl),
               expected=unname(logLik(lmr)[[1]]),
               tolerance=tolerance)
})

context("Factor binomial")
test_that("Factor binomial", {
  sleepstudyM <- sleepstudyU
  sleepstudyM$highR <- factor(ifelse(sleepstudyM$Reaction > 340, 2, 1), levels=c(1,2), labels=c("L","H"))
  sleepstudyM$Sub <- factor(sleepstudyM$Subject, levels=300:400) 
  m1 <- mix(highR ~ Days + (1|Subject), data=sleepstudyM, weights=c("weight1L1", "weight1L2"), 
            family="binomial")
  expect_equal(coef(m1), c(`(Intercept)` = -8.69336570860094, Days = 1.17236933934874), tol=1e-3)
  summaryREF <- c("Call:",
                  "mix(formula = highR ~ Days + (1 | Subject), data = sleepstudyM, ", 
                  "    weights = c(\"weight1L1\", \"weight1L2\"), family = \"binomial\")", 
                  "",
                  "Variance terms:",
                  " Level   Group        Name Variance Std. Error Std.Dev.", 
                  "     2 Subject (Intercept)     7.99       5.14     2.83", "Groups:", 
                  " Level   Group n size mean wgt sum wgt",
                  "     2 Subject     18        1      18",
                  "     1     Obs    180        1     180",
                  "",
                  "Fixed Effects:",
                  "            Estimate Std. Error t value",
                  "(Intercept)   -8.694      2.073   -4.19", 
                  "Days           1.172      0.278    4.21",
                  "",
                  "lnl= -53.90 ")
  withr::with_options(list(digits=2),
                      co <- capture.output(summary(m1))
  )
  expect_equal(co, summaryREF)
})

context("Unweighted model with 2 random effects")
test_that("Agrees with lme4 3,handles missing data", {
  sleepstudyM <- sleepstudyU
  sleepstudyM <- sleepstudyM[order(runif(nrow(sleepstudyM))), ]
  #introduce a missing value 
  sleepstudyM$Days[3] <- NA
  sleepstudyM$Subject <- ifelse(sleepstudyM$Subject %in% c("310"),"2",as.character(sleepstudyM$Subject))
  sleepstudyM$Subject <- ifelse(sleepstudyM$Subject %in% c("309"),"1",as.character(sleepstudyM$Subject))
  sleepstudyM$Subject <- ifelse(sleepstudyM$Subject %in% c("330"),"3",as.character(sleepstudyM$Subject))

  lme2 <- lmer(Reaction ~ Days + (1 | Subject) + (0 + Days | Subject), data=sleepstudyM, REML=FALSE)
  
  # the dropped row should cause a warning
  expect_warning(wm2 <- mix(Reaction ~ Days + (1 | Subject) + (0 + Days | Subject), data=sleepstudyM, weights=c("weight1L1", "weight1L2")),"with missing data")
  expect_equal(wm2$lnl, as.numeric(logLik(lme2)), tol=2E-7*abs(as.numeric(logLik(lme2))))
  
  # check coef
  expect_equal(coef(wm2),
               expected = getME(lme2, "fixef"),
               tolerance = tolerance)
  # check vars
  lmewm2vars <- data.frame(summary(lme2)$varcor)$sdcor
  expect_equal(unname(wm2$vars),
               expected = unname(lmewm2vars)^2,
               tolerance = tolerance)
  lme4ranef <- ranef(lme2)$Subject
  # not in WeMix output
  attr(lme4ranef, "postVar") <- NULL
  expect_equal(lme4ranef, wm2$ranefMat$Subject, 0.01)
  
  expect_equal(predict(lme2), predict(wm2), tol=1e-4)
  sleepstudyM2 <- sleepstudyM[order(runif(nrow(sleepstudyM))), ]
  sleepstudyM2$Subject <- ifelse(sleepstudyM2$Subject %in% c("332"),"4",as.character(sleepstudyM$Subject))
  sleepstudyM2$Subject <- ifelse(sleepstudyM2$Subject %in% c("308"),"1",as.character(sleepstudyM$Subject))
  sleepstudyM2$Subject <- ifelse(sleepstudyM2$Subject %in% c("351"),"5",as.character(sleepstudyM$Subject))

  expect_equal(predict(lme2,sleepstudyM2, allow.new.levels=TRUE), predict(wm2,sleepstudyM2, allow.new.levels=TRUE), tol=1e-4)

})

context("Mean Centering Matches HLM results")
test_that("Mean Centering Matches HLM results", {
  wm1 <- mix(Reaction ~ Days + (1|Subject), data=sleepstudyU, weights=c("weight1L1", "weight1L2"), nQuad=13,center_group=list("Subject"= as.formula(~Days)), verbose=FALSE)
  expect_equal(unname(wm1$coef), c(298.507892, 10.467286), tolerance=1E-1)
  expect_equal(wm1$lnl, -897.0393, tolerance=tolerance*897) # value from  lmer(Reaction ~ Days + (1 | Subject), data=sleepstudy, REML=FALSE)
})

ss <- sleepstudy
ss1 <- ss
ss2 <- ss
doubles <- c(308, 309, 310) # subject with double obs

ss2 <- rbind(ss, subset(ss, Subject %in% doubles))

ss$w <- ifelse(ss$Subject %in% doubles, 2, 1)
contrasts(ss2$Subject) <- "contr.sum"
ss1$W1 <- ifelse(ss1$Subject %in% doubles, 2, 1)
ss1$W2 <- 1
ss1$bin <- ifelse(sleepstudy$Reaction<300,0,1) #for the binomial test

ss2$W2 <- ss2$W1 <- 1

doubles <- c(308, 309, 310) # subject with double obs
ss30 <- subset(ss, Subject %in% doubles)
ss30$Subject <- as.numeric(as.character(ss30$Subject)) + 1000
ss0 <- ss
ss0$Subject <- as.numeric(as.character(ss$Subject))
ss3 <- rbind(ss0, ss30)
ss3$Subject <- as.factor(ss3$Subject)

ss3$W2 <- 1
ss3$W1 <- 1

ss4 <- ss
ss4$W2 <- ifelse(ss4$Subject %in% doubles, 2, 1)
ss4$W1 <- ifelse(ss4$Subject %in% doubles, 2, 1) # make unconditional weight

context("GLM works: Binomial")
test_that("GLM works: Binomial", {
  skip_on_cran()
  #full test for binomial 
  bi_1 <- mix(bin~Days + (1|Subject), data=ss1, family=binomial(link="logit"), verbose=FALSE,
              weights=c("W1", "W2"), nQuad=13)
  expect_equal(unname(bi_1$coef), c(-3.3448,.5928), tolerance=1E-3)
  expect_equal(as.numeric(bi_1$lnl), -93.751679, tolerance=1E-5)
  sum_bi <-  summary(bi_1)
  expect_is(summary(bi_1), "summaryWeMixResults")
})

context("Repeating is the same as weighting: L1 replicate vs weighting")
test_that("Repeating is the same as weighting: L1 replicate vs weighting", {
  # mix for L1, weighted
  wmeL1W <- mix(formula=Reaction ~ Days + (1 | Subject), data=ss1,
                weights=c("W1", "W2"))
  
  # mix for L1, duplicated
  system.time(wmeL1D <- mix(formula=Reaction ~ Days + (1 | Subject), data=ss2,
                            weights=c("W1", "W2")))
  
  # check weighted agrees with duplicated lme4 results
  expect_equal(wmeL1W$lnl, -1048.34318418762, tolerance=1050*tolerance)
  
  # check duplicated agrees with duplicated lme4 results
  expect_equal(wmeL1D$lnl, -1048.34318418762, tolerance=1050*tolerance)
  
})

context("grouping factor not sorted")
test_that("grouping factor not sorted", {
  skip_on_cran()
  ss1_mixed <- ss1[c(125:180,1,100,2:99,101:124),]
  row.names(ss1_mixed) <- NULL
  # mix for L1, weighted
  wmeL1W <- mix(formula=Reaction ~ Days + (1 | Subject), data=ss1_mixed,
                weights=c("W1", "W2"), nQuad=13, run=TRUE,  verbose=FALSE)
  
  # mix for L1, duplicated
  system.time(wmeL1D <- mix(formula=Reaction ~ Days + (1 | Subject), data=ss2,
                            weights=c("W1", "W2"), nQuad=13, run=FALSE,  verbose=FALSE))
  #statares <- c(251.4619, 10.40726, 1000.7466, 1338.0865) # not used
  
  # check weighted agrees with duplicated lme4 results
  expect_equal(wmeL1W$lnl, -1048.34318418762, tolerance=tolerance)
  expect_equal(wmeL1D$lnl, -1048.34318418762, tolerance=tolerance)
  
  # check final results
  suppressWarnings(mix1 <-  mix(formula=Reaction ~ Days + (1 | Subject), data=ss1_mixed,
                                weights=c("W1", "W2")))
  suppressWarnings(mix1REF <-  mix(formula=Reaction ~ Days + (1 | Subject), data=ss1,
                                   weights=c("W1", "W2")))
  expect_equal(mix1$coef, mix1REF$coef, tolerance=1e3)
  expect_equal(mix1$vars, mix1REF$vars, tolerance=1e-3)
  expect_equal(mix1$lnl, mix1REF$lnl, tolerance=1e-3)
  #check  weighted fixed effects variances
  expect_equal(unname(sqrt(diag(mix1$cov_mat))), unname(sqrt(diag(mix1REF$cov_mat))),tolerance = tolerance)
  
})

context("Weighted three level model unsorted")
test_that("Weighted three level model unsorted", {
  skip_on_cran()
  sleepstudy2 <- sleepstudy
  sleepstudy2$Group <- 1
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
  sleepstudy2$Group <- factor(sleepstudy2$Group)
  ss2 <- sleepstudy2
  w1c <- w1 <- rep(1,180)
  w2c <- w2 <- rep(1,18)
  w3c <- w3 <- rep(1,4)
  # unbalanced (non-identical within a Subject), level-1 (obs level) weights
  w1c[1:4] <- w1[1:4] <- 2
  uR <- sleepstudy2[1:4,]
  sleepstudy2 <- rbind(sleepstudy2, uR)
  # level-2 weights
  w2c[1] <- w2[1] <- 2
  w1[ss2$Subject=="308"] <- 2*w1[ss2$Subject=="308"]
  sR <- subset(sleepstudy2, Subject == "308")
  sR$Subject <- "S308"
  sleepstudy2 <- rbind(sleepstudy2, sR)
  # level-3 weights
  w3c[1] <- w3[1] <- 2
  w2[c(1,7,12,13,16,17)] <- 2*w2[c(1,7,12,13,16,17)]
  w1[ss2$Group==1] <- 2*w1[ss2$Group==1]
  
  gR <- subset(sleepstudy2, Group %in% 1)
  gR$Subject <- paste0("G", gR$Subject)
  gR$Group <- "11"
  sleepstudy2 <- rbind(sleepstudy2, gR)
  
  # lmr for reference
  lmr <- lmer(Reaction ~ Days + (1|Subject) + (1|Group), 
              data=sleepstudy2,  control=lmerControl(optimizer="bobyqa"),
              REML=FALSE)
  
  ss2$w1 <- w1
  ss2$w2 <- rep(w2,each=10)
  ss2$w3 <- ifelse(ss2$Subject %in% c("308", "333", "350", "351", "370", "371"),2,1)
  
  ss3 <- ss2[sample(row.names(ss2), size=nrow(ss2)), ]
  wm0 <- mix(Reaction ~ Days + (1|Subject) + (1|Group), data=ss3, 
             weights=c("w1", "w2", "w3"))
  
  # check coef
  expect_equal(coef(wm0),
               expected = getME(lmr, "fixef"),
               tolerance = tolerance)

  lmeRanef <- ranef(lmr)$Subject
  # not in WeMix output
  attr(lmeRanef, "postVar") <- NULL
  # only first 18 of lmeRanef are unique/comparable 
  expect_equal(lmeRanef[1:18,,drop=FALSE], wm0$ranefMat$Subject, 20*sqrt(.Machine$double.eps))
  
  # check vars
  lmewm2vars <- data.frame(summary(lmr)$varcor)$sdcor
  expect_equal(unname(wm0$vars),
               expected = unname(lmewm2vars)^2,
               tolerance = tolerance)
  
  # var beta not expected to be equal  
})

context("repeating is the same as weighting: L2 replicate vs weighting")
test_that("L2 replicate vs weighting", {
  # mix for L2, duplicated
  system.time(wmeL2D <- mix(formula=Reaction ~ Days + (1 | Subject),
                            data=ss3, weights=c("W1", "W2")))
  
  # mix for L2, weighted
  system.time(wmeL2W <- mix(formula=Reaction ~ Days + (1 | Subject), data=ss4,
                            weights=c("W1", "W2")))
  
  expect_equal(wmeL2W$lnl, -1055.34690957995, tolerance=2E-7)
  expect_equal(wmeL2D$lnl, -1055.34690957995, tolerance=2E-7)
})

context("Repeating is the same as weighting: L1 replicate vs weighting, 2 REs")
test_that("Repeating is the same as weighting: L1 replicate vs weighting, 2 REs", {
  # mix for L1, weighted, 2 REs
  wmeL1WRE2 <- mix(formula=Reaction ~ Days + (1 | Subject) + (0+Days|Subject),
                   data=ss1, weights=c("W1", "W2"), nQuad=13, run=FALSE, verbose=FALSE)
  
  # mix for L1, duplicated, 2 REs
  wmeL1DRE2 <- mix(formula=Reaction ~ Days + (1 | Subject) + (0+Days|Subject),
                   data=ss2, weights=c("W1", "W2"),nQuad=13, run=FALSE, verbose=FALSE)
  
  expect_equal(wmeL1WRE2$lnl, -1018.29298875158, tolerance=1050*2E-7)
  
  expect_equal(wmeL1DRE2$lnl, -1018.29298875158, tolerance=1050*2E-7)
})

ssB <- sleepstudy
set.seed(2)
ssB$Reaction <- ssB$Days * 3.141 + rnorm(nrow(ssB))
ssB$W2 <- 1
ssB$W1 <- 1

context("Zero variance estimate")
test_that("Zero variance estimate", {
  skip_on_cran()
  # this has 0 variance estimate in lmer
  lmeB <- lmer(Reaction ~ Days + (1|Subject), data=ssB, REML=FALSE)
  mixB <- mix(formula=Reaction ~ Days + (1 | Subject), data=ssB,
              weights=c("W1", "W2"),  nQuad=13, run=TRUE,   verbose=FALSE)
  expect_equal(mixB$lnl, as.numeric(logLik(lmeB)), tol=1e-5)
  expect_equal(coef(mixB), fixef(lmeB), tol=1e-5)
  expect_equal(unname(mixB$vars[length(mixB$vars)]), unname(lmeB@devcomp$cmp["sigmaML"]^2), tol=1e-5)
  
  ss1 <- sleepstudy
  #add group variables for 3 level model 
  ss1$Group <- 1
  ss1$Group <- ifelse(ss1$Subject %in% c(349,335,330, 352, 337, 369), 2, ss1$Group)
  
  # Create weights
  ss1$W1 <- ifelse(ss1$Subject %in% c(308, 309, 310), 2, 1)
  ss1$W2 <- 1
  ss1$W3 <- ifelse(ss1$Group == 2, 2, 1)
  
  #Run three level model with random slope and intercept. 
  three_level <- mix(Reaction ~ Days + (1|Subject) + (1+Days|Group), data=ss1, weights = c("W1","W2","W3"))
  expect_is(three_level, "WeMixResults")
})

context("Unweighted three level model")
test_that("Unweighted three level model", {
  sleepstudy2 <- sleepstudy
  sleepstudy2$Group <- 1
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
  sleepstudy2$Group <- factor(sleepstudy2$Group)
  sleepstudy2$w1 <- 1 
  sleepstudy2$w2 <- 1
  sleepstudy2$w3 <- 1
  wm0 <- mix(Reaction ~ Days + (1|Subject) + (0+Days|Subject) + (1 | Group), data=sleepstudy2, weights=c("w1", "w2","w3"), verbose=FALSE, run=TRUE)
  lm0 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject) + (1 | Group), data=sleepstudy2,REML=FALSE)
  wm0RE <- wm0$ranefMat
  lm0RE <- ranef(lm0)
  expect_equal(names(wm0RE), names(lm0RE))
  attr(lm0RE$Subject, "postVar") <- NULL
  expect_equal(wm0RE$Subject, lm0RE$Subject, tol=.Machine$double.eps^0.25)
  attr(lm0RE$Group, "postVar") <- NULL
  expect_equal(wm0RE$Group, lm0RE$Group, tol=.Machine$double.eps^0.25)

  wm0RE <- wm0$ranefMat
  lm0RE <- ranef(lm0)
  attr(lm0RE$Subject, "postVar") <- NULL
  expect_equal(wm0RE$Subject, lm0RE$Subject, tol=.Machine$double.eps^0.25)
  
  # check vars
  lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
  expect_equal(unname(wm0$vars),
               expected = unname(lmevars1)^2,
               tolerance = 1e-3)
  expect_equal(wm0$lnl, as.numeric(logLik(lm0)), tol=1e-3)
  expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
  expect_equal(unname(wm0$vars[length(wm0$vars)]), unname(lm0@devcomp$cmp["sigmaML"]^2), tol=1e-4)
  ranef <- ranef(lm0)
  # WeMix does not have postVar attribute
  attr(ranef$Subject, "postVar") <- NULL
  attr(ranef$Group, "postVar") <- NULL
  expect_equal(ranef$Subject, wm0$ranefMat$Subject, (.Machine$double.eps)^0.25)
  expect_equal(ranef$Group, wm0$ranefMat$Group, (.Machine$double.eps)^0.25)
})

context("Unweighted three level model, binomial")
test_that("Unweighted three level model, binomial", {
  sleepstudy2 <- sleepstudy
  set.seed(564)
  sleepstudy2$rand_slope <- rnorm(180,mean=3.25,sd=1.79)
  sleepstudy2$highR <- ifelse(sleepstudy2$Reaction > 340, 1, 0)
  sleepstudy2$Group <- 1
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
  sleepstudy2$Group <- factor(sleepstudy2$Group)
  sleepstudy2$w1 <- 1 
  sleepstudy2$w2 <- 1
  sleepstudy2$w3 <- 1
  sleepstudy2 <- sleepstudy2[sample(1:nrow(sleepstudy2), nrow(sleepstudy2)),]
  # up the max_iterations for more complex model structure
  suppressWarnings(wm0 <- mix(highR ~ Days + (rand_slope|Subject) + (1 | Group), data=sleepstudy2,
             family=binomial(),weights=c("w1", "w2","w3"), verbose=FALSE, run=TRUE,
             nQuad=1,max_iteration = 20))
  lm0 <- glmer(highR ~ Days + (rand_slope|Subject) + (1 | Group), data=sleepstudy2,
               family=binomial())
  # check vars
  lmevars1 <- data.frame(summary(lm0)$varcor)$vcov
  expect_equal(unname(wm0$vars),
               expected = unname(lmevars1),
               tolerance = 1e-3)
  expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
  expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
  expect_equal(predict(wm0), predict(lm0), tol=1e-4)
  expect_equal(predict(wm0, sleepstudy2), predict(lm0, sleepstudy2), tol=1e-4)
  sleepstudy3 <- sleepstudy2[sample(1:nrow(sleepstudy2), floor(nrow(sleepstudy2)/4)),]
  expect_equal(a1 <- predict(wm0, sleepstudy3), b1 <- predict(lm0, sleepstudy3), tol=1e-4)
  sleepstudy4 <- sleepstudy2
  sleepstudy4$Subject <- ifelse(sleepstudy4$Subject == "332", "1", as.character(sleepstudy4$Subject))
  expect_equal(a2 <- predict(wm0, sleepstudy4, allow.new.levels=TRUE), b2 <- predict(lm0, sleepstudy4, allow.new.levels=TRUE), tol=1e-4)
})


context("Weighted v unweighted replicated two level model")
test_that("Weighted v unweighted replicated two level model", {
  sleepstudy2 <- sleepstudy
  sleepstudy2$w1 <- 1 
  sleepstudy2$w2 <- 1
  sleepstudy2$w3 <- 1
  sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
  g1 <- subset(sleepstudy2, sleepstudy2$Subject %in% c("310", "309", "349", "335"))
  g1$Subject <- paste0("R", g1$Subject)
  sleepstudyrep <- rbind(sleepstudy2, g1)
  sleepstudy2$Subject <- factor(sleepstudy2$Subject)
  
  sleepstudy2$w2 <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, 1)
  
  # two level
  wm0 <- mix(Reaction ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"), verbose=FALSE, run=TRUE)
  wmw <- mix(Reaction ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"), verbose=FALSE, run=TRUE, cWeights=TRUE)
  lm0 <- lmer(Reaction ~ Days + (1|Subject), data=sleepstudyrep, REML=FALSE)
  # check vars
  lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
  expect_equal(unname(wm0$vars),
               expected = unname(lmevars1)^2,
               tolerance = 1e-3)
  expect_equal(wm0$lnl, as.numeric(logLik(lm0)), tol=1e-3)
  expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
  expect_equal(unname(wm0$vars[length(wm0$vars)]), unname(lm0@devcomp$cmp["sigmaML"]^2), tol=1e-4)
  
  lm0ranef <- ranef(lm0)$Subject
  # not in WeMix output
  attr(lm0ranef, "postVar") <- NULL
  expect_equal(lm0ranef, wm0$ranefMat$Subject, 20*sqrt(.Machine$double.eps))
  
  
  expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
  expect_equal(coef(wmw), coef(wm0), tol=1e-4)
  expect_equal(wmw$vars, wm0$vars, tol=1e-4)
  
  # should be the same if not conditional as well
  sleepstudy2$w1 <- sleepstudy2$w2
  wmw <- mix(Reaction ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"), verbose=FALSE, run=TRUE, cWeights=FALSE)
  expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
  expect_equal(coef(wmw), coef(wm0), tol=1e-4)
  expect_equal(wmw$vars, wm0$vars, tol=1e-4)
  
  # sub 1 conditional weights work
  sleepstudy2 <- sleepstudy
  sleepstudy2$w1 <- 1 
  sleepstudy2$w2 <- 1
  sleepstudy2$w3 <- 1
  sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
  g1 <- subset(sleepstudy2, sleepstudy2$Subject %in% c("310", "309", "349", "335"))
  # do not renaem subject
  sleepstudyrep <- rbind(sleepstudy2, g1)
  sleepstudyrep$w1 <- ifelse(sleepstudyrep$Subject %in% c("310", "309", "349", "335"), 0.5, 1)
  
  # two level
  wm0 <- mix(Reaction ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"), verbose=FALSE, run=TRUE, cWeights=TRUE)
  wmw <- mix(Reaction ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"), verbose=FALSE, run=TRUE)
  lm0 <- lmer(Reaction ~ Days + (1|Subject), data=sleepstudy2, REML=FALSE)
  # check vars
  lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
  expect_equal(unname(wm0$vars),
               expected = unname(lmevars1)^2,
               tolerance = 1e-3)
  expect_equal(wm0$lnl, as.numeric(logLik(lm0)), tol=1e-3)
  expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
  expect_equal(unname(wm0$vars[length(wm0$vars)]), unname(lm0@devcomp$cmp["sigmaML"]^2), tol=1e-4)
  
  expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
  expect_equal(coef(wmw), coef(wm0), tol=1e-4)
  expect_equal(wmw$vars, wm0$vars, tol=1e-4)
  
  # should be the same if not conditional as well
  wm0 <- mix(Reaction ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"), verbose=FALSE, run=TRUE)
  expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
  expect_equal(coef(wmw), coef(wm0), tol=1e-4)
  expect_equal(wmw$vars, wm0$vars, tol=1e-4)
})

context("Weighted v unweighted replicated three level model")
test_that("Weighted v unweighted replicated three level model", {
  sleepstudy2 <- sleepstudy
  sleepstudy2$Group <- 1
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
  #sleepstudy2$Group <- factor(sleepstudy2$Group)
  sleepstudy2$w1 <- 1 
  sleepstudy2$w2 <- 1
  sleepstudy2$w3 <- 1
  sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
  g1 <- subset(sleepstudy2, Group == 1)
  g1$Group <- 5
  g1$Subject <- paste0("5", g1$Subject)
  sleepstudyrep <- rbind(sleepstudy2, g1)
  s308 <- subset(sleepstudy2, Subject %in% "308")
  s308$Subject <- "R308"
  sleepstudyrep <- rbind(sleepstudyrep, s308)
  s308$Group <- 5
  s308$Subject <- "5R308"
  sleepstudyrep <- rbind(sleepstudyrep, s308)
  sleepstudy2$w1u <- ifelse(sleepstudy2$Group == 1, 2, 1)
  s308unit <- subset(sleepstudyrep, Subject %in% c("308", "5308", "R308", "5R308") & Days %in% c(0,3))
  sleepstudyrep <- rbind(sleepstudyrep, s308unit)
  
  sleepstudy2$Group <- factor(sleepstudy2$Group)
  sleepstudy2$Subject <- factor(sleepstudy2$Subject)
  
  sleepstudy2$w1 <- ifelse(sleepstudy2$Subject %in%  c("308", "R308", "5R308") & sleepstudy2$Days %in% c(0,3), 2, 1)
  sleepstudy2$w2 <- ifelse(sleepstudy2$Subject == "308", 2, 1)
  sleepstudy2$w3 <- ifelse(sleepstudy2$Group == 1, 2, 1)
  sleepstudy2$w2u <- ifelse(sleepstudy2$Group == 1, 2, 1)
  sleepstudy2$w2u <- ifelse(sleepstudy2$Subject == "308", 4, sleepstudy2$w2u)
  sleepstudy2$w1u <- ifelse(sleepstudy2$Group == 1, 2, 1)
  sleepstudy2$w1u <- ifelse(sleepstudy2$Subject == "308", 4, sleepstudy2$w1u)
  sleepstudy2$w1u <- ifelse(sleepstudy2$Subject == "308" & sleepstudy2$Days %in% c(0,3), 8, sleepstudy2$w1u)
  
  # replication
  wm0 <- mix(Reaction ~ Days + (1|Subject) + (1 | Group), data=sleepstudyrep, weights=c("w1", "w2","w3"), verbose=FALSE, run=TRUE)
  # conditional weights
  wmw <- mix(Reaction ~ Days + (1|Subject) + (1 | Group), data=sleepstudy2, weights=c("w1", "w2","w3"), verbose=FALSE, run=TRUE, cWeights=TRUE)
  expect_equal(wm0$lnl, wmw$lnl, tol=20*sqrt(.Machine$double.eps))
  expect_equal(coef(wmw), coef(wm0), tol=20*sqrt(.Machine$double.eps))
  # unconditional weights
  wmwc <- mix(Reaction ~ Days + (1|Subject) + (1 | Group), data=sleepstudy2, weights=c("w1u", "w2u","w3"), verbose=FALSE, run=TRUE, cWeights=FALSE)
  expect_equal(wm0$lnl, wmwc$lnl, tol=20*sqrt(.Machine$double.eps))
  expect_equal(coef(wmwc), coef(wm0), tol=20*sqrt(.Machine$double.eps))
  lm0 <- lmer(Reaction ~ Days + (1|Subject) + (1 | Group), data=sleepstudyrep,REML=FALSE)
  # check vars
  lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
  expect_equal(unname(wm0$vars),
               expected = unname(lmevars1)^2,
               tolerance = 1e-3)
  expect_equal(wm0$lnl, as.numeric(logLik(lm0)), tol=1e-3)
  expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
  expect_equal(unname(wm0$vars[length(wm0$vars)]), unname(lm0@devcomp$cmp["sigmaML"]^2), tol=1e-4)
  
  expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
  expect_equal(coef(wmw), coef(wm0), tol=1e-4)
  expect_equal(wmw$vars, wm0$vars, tol=1e-4)
  
  expect_equal(wm0$ranefMat$Subject[1:18,,drop=FALSE], wmw$ranefMat$Subject, 20*sqrt(.Machine$double.eps))
  expect_equal(wm0$ranefMat$Group[1:4,,drop=FALSE], wmw$ranefMat$Group, 20*sqrt(.Machine$double.eps))
})

context("Weighted v unweighted replicated two level model, binomial")
test_that("Weighted v unweighted replicated two level model, binomial", {
  sleepstudy2 <- sleepstudy
  sleepstudy2$w1 <- 1 
  sleepstudy2$w2 <- 1
  sleepstudy2$w3 <- 1
  sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
  sleepstudy2$highR <- ifelse(sleepstudy2$Reaction > 340, 1, 0)
  g1 <- subset(sleepstudy2, sleepstudy2$Subject %in% c("310", "309", "349", "335"))
  g1$Subject <- paste0("R", g1$Subject)
  sleepstudyrep <- rbind(sleepstudy2, g1)
  sleepstudy2$Subject <- factor(sleepstudy2$Subject)
  
  sleepstudy2$w2 <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, 1)
  
  # two level
  wm0 <- mix(highR ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"), 
             family=binomial(), verbose=FALSE)
  wmw <- mix(highR ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"), 
             family=binomial(), verbose=FALSE, cWeights=TRUE,
             start=c(wm0$coef,wm0$vars))
  lm0 <- glmer(highR ~ Days + (1|Subject), data=sleepstudyrep, family=binomial(),
               nAGQ = 13)
  # check vars
  lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
  expect_equal(unname(wm0$vars),
               expected = unname(lmevars1)^2,
               tolerance = 1e-3)
  expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
  expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
  
  expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
  expect_equal(coef(wmw), coef(wm0), tol=1e-4)
  expect_equal(wmw$vars, wm0$vars, tol=1e-4)
  
  # should be the same if not conditional as well
  sleepstudy2$w1 <- sleepstudy2$w2
  wmw <- mix(highR ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"), 
             family=binomial(),
             cWeights=FALSE)
  expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
  expect_equal(coef(wmw), coef(wm0), tol=1e-4)
  expect_equal(wmw$vars, wm0$vars, tol=1e-4)
  
  # sub 1 conditional weights work
  sleepstudy2 <- sleepstudy
  sleepstudy2$highR <- ifelse(sleepstudy2$Reaction > 340, 1, 0)
  sleepstudy2$w1 <- 1 
  sleepstudy2$w2 <- 1
  sleepstudy2$w3 <- 1
  sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
  g1 <- subset(sleepstudy2, sleepstudy2$Subject %in% c("310", "309", "349", "335"))
  # do not renaem subject
  sleepstudyrep <- rbind(sleepstudy2, g1)
  sleepstudyrep$w1 <- ifelse(sleepstudyrep$Subject %in% c("310", "309", "349", "335"), 0.5, 1)
  
  # two level
  wm0 <- mix(highR ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"), 
             family=binomial(),verbose=FALSE, run=TRUE, cWeights=TRUE)
  wmw <- mix(highR ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"), 
             family=binomial(),verbose=FALSE, run=TRUE)
  lm0 <- glmer(highR ~ Days + (1|Subject), data=sleepstudy2, family=binomial(),
               nAGQ=13)
  # check vars
  lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
  expect_equal(unname(wm0$vars),
               expected = unname(lmevars1)^2,
               tolerance = 1e-3)
  expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
  expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
  
  expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
  expect_equal(coef(wmw), coef(wm0), tol=1e-4)
  expect_equal(wmw$vars, wm0$vars, tol=1e-4)
  
  # should be the same if not conditional as well
  wm0 <- mix(highR ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"), 
             family=binomial(),verbose=FALSE, run=TRUE)
  expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
  expect_equal(coef(wmw), coef(wm0), tol=1e-4)
  expect_equal(wmw$vars, wm0$vars, tol=1e-4)
})

context("Weighted v unweighted replicated two level model, binomial with probit link")
test_that("Weighted v unweighted replicated two level model, binomial with probit link", {
  sleepstudy2 <- sleepstudy
  sleepstudy2$w1 <- 1 
  sleepstudy2$w2 <- 1
  sleepstudy2$w3 <- 1
  sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
  sleepstudy2$highR <- ifelse(sleepstudy2$Reaction > 340, 1, 0)
  g1 <- subset(sleepstudy2, sleepstudy2$Subject %in% c("310", "309", "349", "335"))
  g1$Subject <- paste0("R", g1$Subject)
  sleepstudyrep <- rbind(sleepstudy2, g1)
  sleepstudy2$Subject <- factor(sleepstudy2$Subject)
  
  sleepstudy2$w2 <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, 1)
  
  # two level
  suppressWarnings(wm0 <- mix(highR ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"), 
             family=binomial(link="probit"), verbose=FALSE))
  suppressWarnings(wmw <- mix(highR ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"), 
             family=binomial(link="probit"), verbose=FALSE, cWeights=TRUE,
             start=c(wm0$coef,wm0$vars)))
  lm0 <- glmer(highR ~ Days + (1|Subject), data=sleepstudyrep, family=binomial(link="probit"),
               nAGQ = 13)
  # check vars
  lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
  expect_equal(unname(wm0$vars),
               expected = unname(lmevars1)^2,
               tolerance = 1e-3)
  expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
  expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
  
  expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
  expect_equal(coef(wmw), coef(wm0), tol=1e-4)
  expect_equal(wmw$vars, wm0$vars, tol=1e-4)
  
  # should be the same if not conditional as well
  sleepstudy2$w1 <- sleepstudy2$w2
  suppressWarnings(wmw <- mix(highR ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"), 
             family=binomial(link="probit"),
             cWeights=FALSE))
  expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
  expect_equal(coef(wmw), coef(wm0), tol=1e-4)
  expect_equal(wmw$vars, wm0$vars, tol=1e-4)
  
  # sub 1 conditional weights work
  sleepstudy2 <- sleepstudy
  sleepstudy2$highR <- ifelse(sleepstudy2$Reaction > 340, 1, 0)
  sleepstudy2$w1 <- 1 
  sleepstudy2$w2 <- 1
  sleepstudy2$w3 <- 1
  sleepstudy2$Subject <- as.character(sleepstudy2$Subject)
  g1 <- subset(sleepstudy2, sleepstudy2$Subject %in% c("310", "309", "349", "335"))
  # do not renaem subject
  sleepstudyrep <- rbind(sleepstudy2, g1)
  sleepstudyrep$w1 <- ifelse(sleepstudyrep$Subject %in% c("310", "309", "349", "335"), 0.5, 1)
  
  # two level
  suppressWarnings(wm0 <- mix(highR ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"), 
             family=binomial(link="probit"),verbose=FALSE, run=TRUE, cWeights=TRUE))
  suppressWarnings(wmw <- mix(highR ~ Days + (1|Subject), data=sleepstudy2, weights=c("w1", "w2"), 
             family=binomial(link="probit"),verbose=FALSE, run=TRUE))
  lm0 <- glmer(highR ~ Days + (1|Subject), data=sleepstudy2, family=binomial(link="probit"),
               nAGQ=13)
  # check vars
  lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
  expect_equal(unname(wm0$vars),
               expected = unname(lmevars1)^2,
               tolerance = 1e-3)
  expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
  expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
  
  expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
  expect_equal(coef(wmw), coef(wm0), tol=1e-4)
  expect_equal(wmw$vars, wm0$vars, tol=1e-4)
  
  # should be the same if not conditional as well
  suppressWarnings(wm0 <- mix(highR ~ Days + (1|Subject), data=sleepstudyrep, weights=c("w1", "w2"), 
             family=binomial(link="probit"),verbose=FALSE, run=TRUE))
  expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
  expect_equal(coef(wmw), coef(wm0), tol=1e-4)
  expect_equal(wmw$vars, wm0$vars, tol=1e-4)
})

context("Weighted v unweighted replicated two level model, poisson")
test_that("Weighted v unweighted replicated two level model, poisson", {
  skip_on_cran()
  skip_if_not_installed("glmmTMB")
  require(glmmTMB)
  owls2 <- Owls
  owls2$w1 <- 1 
  owls2$w2 <- 1
  owls2$w3 <- 1

  g1 <- subset(owls2, owls2$Nest %in% c("Jeuss", "Lully", "Seiry", "Oleyes"))
  g1$Nest <- paste0("R", g1$Nest)
  owlsrep <- rbind(owls2, g1)
  
  owls2$w2 <- ifelse(owls2$Nest %in% c("Jeuss", "Lully", "Seiry", "Oleyes"), 2, 1)
  
  # glmer seems to have a bug with nAGQ > 1 for poisson models
  # (see https://stat.ethz.ch/pipermail/r-sig-mixed-models/2020q2/028606.html),
  # so can only compare with nQuad = 1
  suppressWarnings(wm0 <- mix(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owlsrep, weights=c("w1", "w2"), 
             family=poisson(), verbose=FALSE,nQuad=1))
  suppressWarnings(wmw <- mix(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owls2, weights=c("w1", "w2"), nQuad = 1,
             family=poisson(), verbose=FALSE, cWeights=TRUE))
  lm0 <- glmer(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owlsrep, 
               family="poisson")
  # check vars
  lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor

  expect_equal(unname(wm0$vars),
               expected = unname(lmevars1)^2,
               tolerance = 1e-3)
  expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
  expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
  
  expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
  expect_equal(coef(wmw), coef(wm0), tol=1e-4)
  expect_equal(wmw$vars, wm0$vars, tol=1e-4)
  
  # should be the same if not conditional as well
  owls2$w1 <- owls2$w2
  suppressWarnings(wmw <- mix(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owls2, weights=c("w1", "w2"), 
             nQuad = 1,family=poisson(), verbose=FALSE, cWeights=FALSE))
  expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
  expect_equal(coef(wmw), coef(wm0), tol=1e-4)
  expect_equal(wmw$vars, wm0$vars, tol=1e-4)
  
  # sub 1 conditional weights work
  owls2 <- Owls
  owls2$w1 <- 1 
  owls2$w2 <- 1
  owls2$w3 <- 1
  owls2$Nest <- as.character(owls2$Nest)
  g1 <- subset(owls2, owls2$Nest %in% c("Jeuss", "Lully", "Seiry", "Oleyes"))
  # do not renaem subject
  owlsrep <- rbind(owls2, g1)
  owlsrep$w1 <- ifelse(owlsrep$Nest %in% c("Jeuss", "Lully", "Seiry", "Oleyes"), 0.5, 1)
  
  # two level
  suppressWarnings(wm0 <- mix(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owlsrep, weights=c("w1", "w2"), 
             family=poisson(), verbose=FALSE,cWeights=TRUE,nQuad=1))
  suppressWarnings(wmw <- mix(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owls2, weights=c("w1", "w2"), nQuad = 1,
             family=poisson(), verbose=FALSE))
  lm0 <- glmer(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owls2, 
               family="poisson")
  # check vars
  lmevars1 <- data.frame(summary(lm0)$varcor)$sdcor
  expect_equal(unname(wm0$vars),
               expected = unname(lmevars1)^2,
               tolerance = 1e-3)
  expect_equal(as.numeric(wm0$lnl), as.numeric(logLik(lm0)), tol=1e-3)
  expect_equal(coef(wm0), fixef(lm0), tol=1e-4)
  
  expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
  expect_equal(coef(wmw), coef(wm0), tol=1e-4)
  expect_equal(wmw$vars, wm0$vars, tol=1e-4)
  
  # should be the same if not conditional as well
  suppressWarnings(wm0 <- mix(SiblingNegotiation ~ ArrivalTime + (1|Nest), data=owlsrep, weights=c("w1", "w2"), 
             family=poisson(), verbose=FALSE,nQuad=1))
  expect_equal(wmw$lnl, wm0$lnl, tol=1e-3)
  expect_equal(coef(wmw), coef(wm0), tol=1e-4)
  expect_equal(wmw$vars, wm0$vars, tol=1e-4)
})

context("Quadrature works, 3 RE at level 2")
test_that("Quadrature works, 3 RE at level 2", {
  skip_on_cran()
  skip_if_not_installed("glmmTMB")
  require(glmmTMB)
  owls2 <- Owls
  owls2$w1 <- 1 
  owls2$w2 <- 1
  owls2$w3 <- 1
  
  suppressWarnings(wm0 <- mix(SiblingNegotiation ~ ArrivalTime + (ArrivalTime + NegPerChick|Nest), 
                              data=owls2, weights=c("w1", "w2"), 
                              family=poisson(),verbose=FALSE,nQuad=7))
  expect_is(wm0, "WeMixResults")
  
})

context("Three level model slash and colon")
test_that("Three level model slash and colon", {
  sleepstudy2 <- sleepstudy
  sleepstudy2$Group <- 1
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
  sleepstudy2$Group <- factor(sleepstudy2$Group)
  # make subj that restarts at 1 per group, so there are four Subjects with subj=1, one per group.
  sleepstudy2 <- sleepstudy2[order(sleepstudy2$Group, sleepstudy2$Subject),]
  sleepstudy2$subj <- factor(rep(c(1:6,1:4,1:4,1:4),each=10))
  # table(sleepstudy2$subj, sleepstudy2$Group) # shows each group has a subject 1:4 and group 1 as a 5 and 6 too
  ss2 <- sleepstudy2
  ss2$w1 <- w1c <- w1 <- rep(1,180)
  w2c <- w2 <- rep(1,18)
  w3c <- w3 <- rep(1,4)
  ss2$w2 <- rep(1, 180)
  ss2$w3 <- rep(1, 180)
  
  lmr <- lmer(Reaction ~ Days + (1|Group) + (1|Subject), data=sleepstudy2, REML=FALSE)
  
  # next line is a bad specification: confounds group=1, subj=1 person with group=2, subj=2 person
  # lmr <- lmer(Reaction ~ Days + (1|Group) + (1|subj), data=sleepstudy2, REML=FALSE)
  wm0 <- mix(Reaction ~ Days + (1|Group:subj) + (1|Group) , data=ss2, weights=c("w1", "w2","w3"))
  expect_equal(wm0$lnl, as.numeric(logLik(lmr)), tol=1e-6)
  expect_equal(coef(wm0), fixef(lmr), tol=1e-6)
  expect_equal(unname(wm0$vars[length(wm0$vars)]), unname(lmr@devcomp$cmp["sigmaML"]^2), tol=1e-4)
  
  #group mean center Days to test group mean centering as well
  # we know the mean is 4.5 for all groups because all individuals were tested on all days
  sleepstudy2$gmc_days <- sleepstudy2$Days-4.5
  
  lmr <- lmer(Reaction ~ gmc_days + (1|Group/subj), data=sleepstudy2, REML=FALSE)
  # above is same as:
  # lmr <- lmer(Reaction ~ Days + (1|Group) + (1|Group:subj), data=sleepstudy2, REML=FALSE)
  wm0 <- mix(Reaction ~ Days + (1|Group/subj), data=ss2, weights=c("w1", "w2","w3"),center_group = list("Group/subj" = ~Days))
  expect_equal(wm0$lnl, as.numeric(logLik(lmr)), tol=1e-6)
  expect_equal(unname(coef(wm0)), unname(fixef(lmr)), tol=1e-6)
  expect_equal(unname(wm0$vars[length(wm0$vars)]), unname(lmr@devcomp$cmp["sigmaML"]^2), tol=1e-4)
  
})

# check the format of summary output
context("Summary output format")
test_that("summary output format", {
  skip_on_cran()
  sleepstudy2 <- sleepstudy
  sleepstudy2$Group <- 1
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
  sleepstudy2$Group <- factor(sleepstudy2$Group)
  # make subj that restarts at 1 per group, so there are four Subjects with subj=1, one per group.
  sleepstudy2 <- sleepstudy2[order(sleepstudy2$Group, sleepstudy2$Subject),]
  sleepstudy2$subj <- factor(rep(c(1:6,1:4,1:4,1:4),each=10))
  sleepstudy2$carPr <- pnorm(sleepstudy2$Reaction -285,sd=50)
  set.seed(2) 
  sleepstudy2$car <- runif(180) < sleepstudy2$carPr
  # table(sleepstudy2$subj, sleepstudy2$Group) # shows each group has a subject 1:4 and group 1 as a 5 and 6 too
  ss2 <- sleepstudy2
  ss2$w1 <- w1c <- w1 <- rep(1,180)
  w2c <- w2 <- rep(1,18)
  w3c <- w3 <- rep(1,4)
  ss2$w2 <- rep(1, 180)
  ss2$w3 <- rep(1, 180)
  
  co0 <- c("Call:",
           "mix(formula = Reaction ~ Days + (Days + car | Subject), data = ss2, ", 
           "    weights = c(\"w1\", \"w2\"))", "",
           "Variance terms:",
           " Level    Group        Name Variance Std. Error Std.Dev. Corr1 Corr2", 
           "     2  Subject (Intercept)    629.4      254.1    25.09            ", 
           "     2  Subject        Days     36.7       12.2     6.06  0.17      ", 
           "     2  Subject     carTRUE    648.5      309.3    25.46 -0.42 -0.39", 
           "     1 Residual                528.8      145.7    23.00            ", 
           "Groups:",
           " Level   Group n size mean wgt sum wgt",
           "     2 Subject     18        1      18",
           "     1     Obs    180        1     180",
           "",
           "Fixed Effects:",
           "            Estimate Std. Error t value", 
           "(Intercept)   252.72       6.53   38.70",
           "Days           11.19       1.50    7.44",
           "",
           "lnl= -868.13 ",
           "Intraclass Correlation= 0.713 ")
  
  wm0 <- mix(Reaction ~ Days + (Days+car|Subject), data=ss2, weights=c("w1", "w2"))
  withr::with_options(list(digits=2),
                      co <- capture.output(summary(wm0))
  )
  expect_equal(co, co0)
  # this model is not accuractly estimated so is not a good test of output format. 
  # nevertheless, it uses data from this block, so appears here
  expect_warning(wm0 <- mix(Reaction ~ Days + (car||Subject), data=ss2, weights=c("w1", "w2")), "standard error of variance terms")
  
  varsref <- structure(list(Variance    = c(698.125, 213.50 , 1618.92, 808.25  ),
                            `Std. Error` = c(385.8  , 232    , 850.5  , 162.0972),
                            Std.Dev.    = c(26.4220, 14.6116, 40.2358, 28.42977)),
                       class = "data.frame",
                       row.names = c("Subject.(Intercept)", "Subject.carFALSE", "Subject.carTRUE", "Residual"))
  coefref <- structure(c(236.842700778804, 9.59410909072588, 6.7061583949811, 
                         1.56079000196884, 35.317194558968, 6.14695704010372), .Dim = 2:3,
                       .Dimnames = list(c("(Intercept)", "Days"),
                                        c("Estimate", "Std. Error", "t value")))
  sum0 <- summary(wm0)
  expect_equal(sum0$coef, coefref, tolerance=1e-4)
  expect_equal(sum0$vars, varsref, tolerance=4e-3)
})

context("Weighted three level model")
test_that("Weighted three level model", {
  sleepstudy2 <- sleepstudy
  sleepstudy2$Group <- 1
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
  sleepstudy2$Group <- factor(sleepstudy2$Group)
  
  #sleepstudy2$NewVar <- sleepstudy2$Reaction/100 + rnorm(180,10,5)
  sleepstudy2$NewVar <- c(9.43457771946133, 12.7183652886321, 11.242490364576, 8.32049769489214, 21.4769136458094, 12.801691804759, 11.2944480319612, 23.2877967612343, 17.6554807720438, 15.6839576150421,
                          12.0415732799454, 12.6326929342708, 10.7497191412678, 8.67408880103278, 16.7872963953428, 14.6248001125594, 14.3553882674867, 11.8275049579491, 18.7525279773655, 5.57650843995549, 
                          19.7623973348143, 19.1858339464439, 16.9845555399241, 5.79214586016569, 13.6760775160985, 8.20870365747105, 10.469622698121, 9.5294364080553, 0.474272552642729, 9.90998530696626,
                          7.85598926755497, 13.6908324775899, 8.363578606421, 12.4585008913942, 20.8371544197279, 6.95572306644593, 4.49104017816579, 13.529348727477, 8.31255622061886, 8.78452997885101, 
                          17.216378291436, 13.6427335980228, 19.498592420836, 7.77906722601152, 15.1266474871889, 12.8030124499139, 6.62054433117513, 8.66974559938289, 11.9540398510289, 14.4271543736958,
                          10.0496588773239, 15.5177485124749, 9.12750122498998, 10.1801282848309, 14.2562515194416, 19.3249322019648, 10.7021090784807, 19.5336266451282, 7.58133671725254, 19.889529248361, 
                          5.32630146259616, 13.6044477039579, 4.98570236262749, 23.8801945918399, 14.447634009947, 15.6043025989833, 22.411260279495, 16.66695588071, 13.2723030644539, 17.6752877679585,
                          13.4050771535649, 16.8476633007641, 14.4754165647234, 13.2597805356783, 12.2720489082021, 16.3515751846435, 18.2589186962967, 13.5608781488381, 22.128982372665, 14.0340540905538,
                          2.08454544532827, 10.8863647522786, 21.0103204678449, 18.4595691535641, 11.9914499297614, 9.96791953740012, 20.0728999616675, 10.7237225063556, 18.6117287544885, 15.6290844263721,
                          13.584186380461, 12.8188973592142, 8.73588740878364, 22.1049445716224, 8.82190889109436, 16.7501795316447, 11.7894955061621, 8.67751729426691, 8.28560204932531, 25.3725616103908,
                          14.0620375346182, 14.9882805901214, 17.181272042065, 19.6581190406994, 11.6757977643058, 16.638959090019, 9.68436466184703, 23.1282879805077, 20.7690736926973, 14.3401954001056,
                          16.2747349705209, 8.61220295948687, 19.869302293393, 5.78478905807136, 6.59988762289215, 20.489844013367, 27.239468499143, 20.5797203990527, 19.4544929129352, 12.5320006698243, 
                          18.2664525619468, 16.6022500953422, 13.0216995414837, 12.4573456587016, 18.6922622366627, 7.37338081469625, 18.2016743123929, 12.9216600044986, 17.2543801572655, 20.1850819959469,
                          9.46342758905614, 16.1505649461445, 11.644647125078, 17.6992809329738, 13.5726089392145, 11.9006738563255, 9.11529530502925, 11.9473408244809, 13.7037965361481, 21.052192466607, 
                          8.94718178506257, 15.165337112518, 8.23261412105177, 12.2188157585207, 14.1455895553682, 4.94479286361234, 16.9113105275149, 12.7165809453375, 5.58774136619933, 8.31670629685991,
                          6.67728085523989, 11.7741431392146, 7.91590113917014, 9.13428547474213, 18.5689544076891, 18.1995195536489, 13.3510327364336, 21.4047919488063, 17.7821335119331, 9.19515401988151, 
                          10.7918240635605, 6.63952336871624, 15.3733302621749, 15.5262403960572, 13.4421915705244, 11.5980148647886, 17.3920393426159, 6.6883267354037, 10.6865661488425, 15.1104585129107,
                          9.99190979270527, 16.54240123493, 13.5687018760124, 13.5592897858269, 19.4019579232211, 15.0561754718535, 9.36036662501835, 4.27173388144395, 19.2040511963462, 6.36268714321505)
  
  ss2 <- sleepstudy2
  w1c <- w1 <- rep(1,180)
  w2c <- w2 <- rep(1,18)
  w3c <- w3 <- rep(1,4)
  # unbalanced (non-identical within a Subject), level-1 (obs level) weights
  w1c[1:4] <- w1[1:4] <- 2
  uR <- sleepstudy2[1:4,]
  sleepstudy2 <- rbind(sleepstudy2, uR)
  # level-2 weights
  w2c[1] <- w2[1] <- 2
  w1[ss2$Subject=="308"] <- 2*w1[ss2$Subject=="308"]
  sR <- subset(sleepstudy2, Subject == "308")
  sR$Subject <- "S308"
  sleepstudy2 <- rbind(sleepstudy2, sR)
  # level-3 weights
  w3c[1] <- w3[1] <- 2
  w2[c(1,7,12,13,16,17)] <- 2*w2[c(1,7,12,13,16,17)]
  w1[ss2$Group==1] <- 2*w1[ss2$Group==1]
  
  gR <- subset(sleepstudy2, Group %in% 1)
  gR$Subject <- paste0("G", gR$Subject)
  gR$Group <- "11"
  sleepstudy2 <- rbind(sleepstudy2, gR)
  
  ss2$w1 <- w1
  ss2$w2 <- rep(w2,each=10)
  ss2$w3 <- ifelse(ss2$Subject %in% c("308", "333", "350", "351", "370", "371"),2,1)
  
  # lmr for reference
  lmr <- lmer(Reaction ~ Days + (1|Subject) + (1|Group), data=sleepstudy2, REML=FALSE)
  
  #these are the conditional weights, used in the stata comparison
  ss2$c1 <- w1c
  ss2$c2 <- rep(w2c,each=10)
  ss2$c3 <- ifelse(ss2$Subject %in% c("308", "333", "350", "351", "370", "371"),2,1)
  
  #compare against mixed stata resutls
  wm0 <- mix(Reaction ~ Days + (1|Subject) + (1|Group), data=ss2, weights=c("w1", "w2","w3"))
  expect_equal(wm0$lnl, -1389.0983, tolerance=1e-3)
  expect_equal(unname(wm0$coef), c(243.6831,12.67954),tolerance = 1e-3)
  expect_equal(unname(wm0$vars), c(360.97 ,756.5857, 1153.704 ),tolerance = 1e-3)
  expect_equal(unname( wm0$SE), c( 11.68881 , 2.468474  ),tolerance = 1e-3)
  
  wm1 <- mix(Reaction ~ Days + (1 |Subject) + (1|Group)+ (0+Days|Group), data=ss2, weights=c("w1", "w2","w3"))
  expect_equal(wm1$lnl,  -1377.6876  , tolerance=1e-3)
  expect_equal(unname(wm1$coef), c(247.6234 , 11.71211 ),tolerance = 1e-3)
  expect_equal(unname(wm1$vars), c(397.6301 , 220.4182  , 16.86835 , 1022.88  ),tolerance = 1e-3)
  expect_equal(unname( wm1$SE), c(9.97378  , 2.60348  ),tolerance = 1e-3)
  
})

context("Wald Test")
test_that("Test for Wald Tests using  ", {
  wm1<- mix(Reaction ~ Days + (0+Days|Subject), data=sleepstudyU, weights=c("weight1L1", "weight1L2"))
  
  beta_1 <- waldTest(wm1, type="beta", coefs="Days", hypothesis=10.4073- 1.5458 *1.96)
  #wald test and T test are the same for individual betas 
  expect_equal(beta_1$p,.05,.01)
  
  beta_2<- waldTest(wm1, type="beta", coefs="(Intercept)", hypothesis=251.4051- 6.8246 *1.96)
  #wald test and T test are the same for individual betas 
  expect_equal(beta_2$p,.05,.01)
  
  #  Ensure that tests for Lambda run
  wm2<- mix(Reaction ~ Days + (1+Days|Subject), data=sleepstudyU, weights=c("weight1L1", "weight1L2"))
  waldTest(wm2,type="Lambda")
  waldTest(wm2,type="Lambda",coefs = c("Subject.Days.(Intercept)","Subject.Days" ))
})

context("Complex weighted three level model")
test_that("Complex weighted three level model", {
  skip_on_cran()
  sleepstudy2 <- sleepstudy
  sleepstudy2$Group <- 1
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
  sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
  sleepstudy2$Group <- factor(sleepstudy2$Group)
  sleepstudy2$NewVar <- c(9.43457771946133, 12.7183652886321, 11.242490364576, 8.32049769489214, 21.4769136458094, 12.801691804759, 11.2944480319612, 23.2877967612343, 17.6554807720438, 15.6839576150421,
                          12.0415732799454, 12.6326929342708, 10.7497191412678, 8.67408880103278, 16.7872963953428, 14.6248001125594, 14.3553882674867, 11.8275049579491, 18.7525279773655, 5.57650843995549, 
                          19.7623973348143, 19.1858339464439, 16.9845555399241, 5.79214586016569, 13.6760775160985, 8.20870365747105, 10.469622698121, 9.5294364080553, 0.474272552642729, 9.90998530696626,
                          7.85598926755497, 13.6908324775899, 8.363578606421, 12.4585008913942, 20.8371544197279, 6.95572306644593, 4.49104017816579, 13.529348727477, 8.31255622061886, 8.78452997885101, 
                          17.216378291436, 13.6427335980228, 19.498592420836, 7.77906722601152, 15.1266474871889, 12.8030124499139, 6.62054433117513, 8.66974559938289, 11.9540398510289, 14.4271543736958,
                          10.0496588773239, 15.5177485124749, 9.12750122498998, 10.1801282848309, 14.2562515194416, 19.3249322019648, 10.7021090784807, 19.5336266451282, 7.58133671725254, 19.889529248361, 
                          5.32630146259616, 13.6044477039579, 4.98570236262749, 23.8801945918399, 14.447634009947, 15.6043025989833, 22.411260279495, 16.66695588071, 13.2723030644539, 17.6752877679585,
                          13.4050771535649, 16.8476633007641, 14.4754165647234, 13.2597805356783, 12.2720489082021, 16.3515751846435, 18.2589186962967, 13.5608781488381, 22.128982372665, 14.0340540905538,
                          2.08454544532827, 10.8863647522786, 21.0103204678449, 18.4595691535641, 11.9914499297614, 9.96791953740012, 20.0728999616675, 10.7237225063556, 18.6117287544885, 15.6290844263721,
                          13.584186380461, 12.8188973592142, 8.73588740878364, 22.1049445716224, 8.82190889109436, 16.7501795316447, 11.7894955061621, 8.67751729426691, 8.28560204932531, 25.3725616103908,
                          14.0620375346182, 14.9882805901214, 17.181272042065, 19.6581190406994, 11.6757977643058, 16.638959090019, 9.68436466184703, 23.1282879805077, 20.7690736926973, 14.3401954001056,
                          16.2747349705209, 8.61220295948687, 19.869302293393, 5.78478905807136, 6.59988762289215, 20.489844013367, 27.239468499143, 20.5797203990527, 19.4544929129352, 12.5320006698243, 
                          18.2664525619468, 16.6022500953422, 13.0216995414837, 12.4573456587016, 18.6922622366627, 7.37338081469625, 18.2016743123929, 12.9216600044986, 17.2543801572655, 20.1850819959469,
                          9.46342758905614, 16.1505649461445, 11.644647125078, 17.6992809329738, 13.5726089392145, 11.9006738563255, 9.11529530502925, 11.9473408244809, 13.7037965361481, 21.052192466607, 
                          8.94718178506257, 15.165337112518, 8.23261412105177, 12.2188157585207, 14.1455895553682, 4.94479286361234, 16.9113105275149, 12.7165809453375, 5.58774136619933, 8.31670629685991,
                          6.67728085523989, 11.7741431392146, 7.91590113917014, 9.13428547474213, 18.5689544076891, 18.1995195536489, 13.3510327364336, 21.4047919488063, 17.7821335119331, 9.19515401988151, 
                          10.7918240635605, 6.63952336871624, 15.3733302621749, 15.5262403960572, 13.4421915705244, 11.5980148647886, 17.3920393426159, 6.6883267354037, 10.6865661488425, 15.1104585129107,
                          9.99190979270527, 16.54240123493, 13.5687018760124, 13.5592897858269, 19.4019579232211, 15.0561754718535, 9.36036662501835, 4.27173388144395, 19.2040511963462, 6.36268714321505)
  
  ss2 <- sleepstudy2
  w1c <- w1 <- rep(1,180)
  w2c <- w2 <- rep(1,18)
  w3c <- w3 <- rep(1,4)
  # unbalanced (non-identical within a Subject), level-1 (obs level) weights
  w1c[1:4] <- w1[1:4] <- 2
  uR <- sleepstudy2[1:4,]
  sleepstudy2 <- rbind(sleepstudy2, uR)
  # level-2 weights
  w2c[1] <- w2[1] <- 2
  w1[ss2$Subject=="308"] <- 2*w1[ss2$Subject=="308"]
  sR <- subset(sleepstudy2, Subject == "308")
  sR$Subject <- "S308"
  sleepstudy2 <- rbind(sleepstudy2, sR)
  # level-3 weights
  w3c[1] <- w3[1] <- 2
  w2[c(1,7,12,13,16,17)] <- 2*w2[c(1,7,12,13,16,17)]
  w1[ss2$Group==1] <- 2*w1[ss2$Group==1]
  
  gR <- subset(sleepstudy2, Group %in% 1)
  gR$Subject <- paste0("G", gR$Subject)
  gR$Group <- "11"
  sleepstudy2 <- rbind(sleepstudy2, gR)
  
  
  ss2$w1 <- w1
  ss2$w2 <- rep(w2,each=10)
  ss2$w3 <- ifelse(ss2$Subject %in% c("308", "333", "350", "351", "370", "371"),2,1)
  
  ss2$n1 <- ss2$n2 <- ss2$n3 <- 1
  #lme reference with duplicated data
  lmr <- lmer(Reaction ~ Days + (1|Subject) + (1+Days|Group), data=sleepstudy2, REML=FALSE)
  
  #Does WeMix match to duplicated lmr with one+two random effects
  wmr <- mix(Reaction ~ Days + (1|Subject) + (1+Days|Group), data=ss2,weights=c("w1","w2","w3"))
  expect_equal(wmr$lnl, as.numeric(logLik(lmr)), tol=1e-3)
  expect_equal(coef(wmr), fixef(lmr), tol=1e-4)
  lmevars1 <- data.frame(summary(lmr)$varcor)
  vars <- lmevars1[is.na(lmevars1$var2),"vcov"]
  expect_equal(unname(wmr$vars),
               expected = unname(vars),
               tolerance = 1e-3)
  expect_equal(coef(wmr), fixef(lmr), tol=1e-4)
  expect_equal(unname(wmr$vars[length(wmr$vars)]), unname(lmr@devcomp$cmp["sigmaML"]^2), tol=1e-4)
  
  #lme reference with duplicated data
  lmr2 <- lmer(Reaction ~ Days + (1 + NewVar |Subject) + (1+Days|Group), 
               data=sleepstudy2, control=lmerControl(optimizer="bobyqa"), REML=FALSE)
  #wemix with weights
  suppressWarnings(wmr2 <- mix(Reaction ~ Days + (1+ NewVar |Subject) + (1+Days|Group), data=ss2,weights=c("w1","w2","w3")))
  #Does WeMix match to duplicated lmr with two correlated random effects at two levesl 
  expect_equal(wmr2$lnl, as.numeric(logLik(lmr2)), tol=1e-3)
  expect_equal(coef(wmr2), fixef(lmr2), tol=1e-4)
  
  lmevars2 <- data.frame(summary(lmr2)$varcor)
  vars <- lmevars2[is.na(lmevars2$var2),"vcov"]
  expect_equal(unname(wmr2$vars),
               expected = unname(vars),
               tolerance = 1e-3)
  
  expect_equal(coef(wmr2), fixef(lmr2), tol=1e-4)
  expect_equal(unname(wmr2$vars[length(wmr2$vars)]), unname(lmr2@devcomp$cmp["sigmaML"]^2), tol=1e-4)
}) 

context("Complex weighted three level model, binomial")
test_that("Complex weighted three level model, binomial", {
    skip_on_cran()
    sleepstudy2 <- sleepstudy
    sleepstudy2$highR <- ifelse(sleepstudy2$Reaction > 340, 1, 0)
    sleepstudy2$Group <- 1
    sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("310", "309", "349", "335"), 2, sleepstudy2$Group)
    sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("330", "352", "337", "369"), 3, sleepstudy2$Group)
    sleepstudy2$Group <- ifelse(sleepstudy2$Subject %in% c("331", "332", "334", "372"), 4, sleepstudy2$Group)
    sleepstudy2$Group <- factor(sleepstudy2$Group)
    sleepstudy2$NewVar <- c(9.43457771946133, 12.7183652886321, 11.242490364576, 8.32049769489214, 21.4769136458094, 12.801691804759, 11.2944480319612, 23.2877967612343, 17.6554807720438, 15.6839576150421,
                            12.0415732799454, 12.6326929342708, 10.7497191412678, 8.67408880103278, 16.7872963953428, 14.6248001125594, 14.3553882674867, 11.8275049579491, 18.7525279773655, 5.57650843995549, 
                            19.7623973348143, 19.1858339464439, 16.9845555399241, 5.79214586016569, 13.6760775160985, 8.20870365747105, 10.469622698121, 9.5294364080553, 0.474272552642729, 9.90998530696626,
                            7.85598926755497, 13.6908324775899, 8.363578606421, 12.4585008913942, 20.8371544197279, 6.95572306644593, 4.49104017816579, 13.529348727477, 8.31255622061886, 8.78452997885101, 
                            17.216378291436, 13.6427335980228, 19.498592420836, 7.77906722601152, 15.1266474871889, 12.8030124499139, 6.62054433117513, 8.66974559938289, 11.9540398510289, 14.4271543736958,
                            10.0496588773239, 15.5177485124749, 9.12750122498998, 10.1801282848309, 14.2562515194416, 19.3249322019648, 10.7021090784807, 19.5336266451282, 7.58133671725254, 19.889529248361, 
                            5.32630146259616, 13.6044477039579, 4.98570236262749, 23.8801945918399, 14.447634009947, 15.6043025989833, 22.411260279495, 16.66695588071, 13.2723030644539, 17.6752877679585,
                            13.4050771535649, 16.8476633007641, 14.4754165647234, 13.2597805356783, 12.2720489082021, 16.3515751846435, 18.2589186962967, 13.5608781488381, 22.128982372665, 14.0340540905538,
                            2.08454544532827, 10.8863647522786, 21.0103204678449, 18.4595691535641, 11.9914499297614, 9.96791953740012, 20.0728999616675, 10.7237225063556, 18.6117287544885, 15.6290844263721,
                            13.584186380461, 12.8188973592142, 8.73588740878364, 22.1049445716224, 8.82190889109436, 16.7501795316447, 11.7894955061621, 8.67751729426691, 8.28560204932531, 25.3725616103908,
                            14.0620375346182, 14.9882805901214, 17.181272042065, 19.6581190406994, 11.6757977643058, 16.638959090019, 9.68436466184703, 23.1282879805077, 20.7690736926973, 14.3401954001056,
                            16.2747349705209, 8.61220295948687, 19.869302293393, 5.78478905807136, 6.59988762289215, 20.489844013367, 27.239468499143, 20.5797203990527, 19.4544929129352, 12.5320006698243, 
                            18.2664525619468, 16.6022500953422, 13.0216995414837, 12.4573456587016, 18.6922622366627, 7.37338081469625, 18.2016743123929, 12.9216600044986, 17.2543801572655, 20.1850819959469,
                            9.46342758905614, 16.1505649461445, 11.644647125078, 17.6992809329738, 13.5726089392145, 11.9006738563255, 9.11529530502925, 11.9473408244809, 13.7037965361481, 21.052192466607, 
                            8.94718178506257, 15.165337112518, 8.23261412105177, 12.2188157585207, 14.1455895553682, 4.94479286361234, 16.9113105275149, 12.7165809453375, 5.58774136619933, 8.31670629685991,
                            6.67728085523989, 11.7741431392146, 7.91590113917014, 9.13428547474213, 18.5689544076891, 18.1995195536489, 13.3510327364336, 21.4047919488063, 17.7821335119331, 9.19515401988151, 
                            10.7918240635605, 6.63952336871624, 15.3733302621749, 15.5262403960572, 13.4421915705244, 11.5980148647886, 17.3920393426159, 6.6883267354037, 10.6865661488425, 15.1104585129107,
                            9.99190979270527, 16.54240123493, 13.5687018760124, 13.5592897858269, 19.4019579232211, 15.0561754718535, 9.36036662501835, 4.27173388144395, 19.2040511963462, 6.36268714321505)
    
    ss2 <- sleepstudy2
    w1c <- w1 <- rep(1,180)
    w2c <- w2 <- rep(1,18)
    w3c <- w3 <- rep(1,4)
    # unbalanced (non-identical within a Subject), level-1 (obs level) weights
    w1c[1:4] <- w1[1:4] <- 2
    uR <- sleepstudy2[1:4,]
    sleepstudy2 <- rbind(sleepstudy2, uR)
    # level-2 weights
    w2c[1] <- w2[1] <- 2
    w1[ss2$Subject=="308"] <- 2*w1[ss2$Subject=="308"]
    sR <- subset(sleepstudy2, Subject == "308")
    sR$Subject <- "S308"
    sleepstudy2 <- rbind(sleepstudy2, sR)
    # level-3 weights
    w3c[1] <- w3[1] <- 2
    w2[c(1,7,12,13,16,17)] <- 2*w2[c(1,7,12,13,16,17)]
    w1[ss2$Group==1] <- 2*w1[ss2$Group==1]
    
    gR <- subset(sleepstudy2, Group %in% 1)
    gR$Subject <- paste0("G", gR$Subject)
    gR$Group <- "11"
    sleepstudy2 <- rbind(sleepstudy2, gR)
    
    
    ss2$w1 <- w1
    ss2$w2 <- rep(w2,each=10)
    ss2$w3 <- ifelse(ss2$Subject %in% c("308", "333", "350", "351", "370", "371"),2,1)
    
    ss2$n1 <- ss2$n2 <- ss2$n3 <- 1
    #lme reference with duplicated data
    glmr <- glmer(highR ~ Days + (1|Subject) + (1+Days|Group), data=sleepstudy2, 
                  family=binomial())
    
    #Does WeMix match to duplicated lmr with one+two random effects
    suppressWarnings(wmr <- mix(highR ~ Days + (1|Subject) + (Days|Group), data=ss2,
               weights=c("w1","w2","w3"),
               family=binomial(),nQuad=1,verbose=FALSE))
    
    expect_equal(as.numeric(wmr$lnl), as.numeric(logLik(glmr)), tol=1e-3)
    expect_equal(coef(wmr), fixef(glmr), tol=1e-4)
    lmevars1 <- data.frame(summary(glmr)$varcor)
    vars <- lmevars1$vcov
    expect_equal(unname(wmr$vars),
                 expected = unname(vars),
                 tolerance = 1e-3)
    
    #lme reference with duplicated data
    glmr2 <- glmer(highR ~ Days + (1 + NewVar|Subject) + (1+Days|Group), 
                 data=sleepstudy2,
                 family=binomial(),control=glmerControl(optimizer="bobyqa"))
    #wemix with weights
    suppressWarnings(wmr2 <- mix(highR ~ Days + (1+ NewVar |Subject) + (1+Days|Group), data=ss2,
                weights=c("w1","w2","w3"), family=binomial(),nQuad=1,
                max_iteration = 20,verbose=FALSE))
    #Does WeMix match to duplicated lmr with two correlated random effects at two levesl 
    expect_equal(as.numeric(wmr2$lnl), as.numeric(logLik(glmr2)), tol=1e-3)
    expect_equal(coef(wmr2), fixef(glmr2), tol=1e-4)
    
    lmevars2 <- data.frame(summary(glmr2)$varcor)
    vars <- lmevars2$vcov
    expect_equal(unname(wmr2$vars),
                 expected = unname(vars),
                 tolerance = 1e-3)
    
  }) 

context("Model Matrix has a hard time with")
test_that("Model Matrix has a hard time with", {
  skip_on_cran()
  require(EdSurvey)
  sdf <- readNAEP(system.file("extdata/data", "M36NT2PM.dat", package = "NAEPprimer"))
  gg <- EdSurvey::getData(varnames=c("composite", "dsex", "b017451", "scrpsu", "origwt", "smsrswt"), data=sdf, returnJKreplicates=FALSE)
  gg2 <- gg[gg$origwt > 0 & gg$smsrswt > 0,]
  suppressMessages(m4 <- mix(mrpcm2 ~ dsex + b017451 + (1|scrpsu), data=gg2, weights=c("origwt", "smsrswt")))
  expect_equal(m4$lnl, -81882.3634148408, tol=1e-5)
  expect_equal(m4$coef, c(`(Intercept)` = 270.571817625202, dsexFemale = -2.14907551600309, `b017451Once every few weeks` = 3.84433128452533, `b017451About once a week` = 9.19954009166631, `b0174512 or 3 times a week` = 12.8701977366809, `b017451Every day` = 6.29635843233831))
  expect_equal(m4$SE, c(`(Intercept)` = 1.14859092700564, dsexFemale = 0.636207386283148, `b017451Once every few weeks` = 1.05190060207538, `b017451About once a week` = 0.986062663493941, `b0174512 or 3 times a week` = 1.00366545196335, `b017451Every day` = 1.10858750669702))
  m4varDF <- structure(list(grp = c("scrpsu", "Residual"),
                            var1 = c("(Intercept)", NA),
                            var2 = c(NA_character_, NA_character_),
                            vcov = c(300.695102536901, 969.501918646696),
                            ngrp = c(672, 16331),
                            level = c(2, 1),
                            SEvcov = c(25.0004604746421, 16.2937862202021),
                            fullGroup = c("scrpsu.(Intercept)", "Residual")),
                       row.names = c(NA, -2L),
                       class = "data.frame")
  expect_equal(m4$varDF, m4varDF, tol=1e-5)
})



context("examples run")
test_that("examples run", {
  skip_on_cran()
  ss1 <- sleepstudy
  
  # Create weights
  ss1$W1 <- ifelse(ss1$Subject %in% c(308, 309, 310), 2, 1)
  ss1$W2 <- 1
  
  # Run random-intercept 2-level model 
  two_level <- mix(Reaction ~ Days + (1|Subject), data=ss1, weights=c("W1", "W2"))
  expect_is(two_level, "WeMixResults")
  
  #Run random-intercept 2-level model with group-mean centering
  grp_centered <- mix(Reaction ~ Days + (1|Subject), data=ss1,
                      weights = c("W1", "W2"),
                      center_group = list("Subject" = ~Days))
  expect_is(grp_centered, "WeMixResults")
  
  #Run three level model with random slope and intercept. 
  #add group variables for 3 level model 
  ss1$Group <- 3
  ss1$Group <- ifelse(as.numeric(ss1$Subject) %% 10 < 7, 2, ss1$Group)
  ss1$Group <- ifelse(as.numeric(ss1$Subject) %% 10 < 4, 1, ss1$Group)
  # level-3 weights
  ss1$W3 <- ifelse(ss1$Group == 2, 2, 1)
  
  three_level <- mix(Reaction ~ Days + (1|Subject) + (1+Days|Group), data=ss1, 
                     weights=c("W1", "W2", "W3"))
  expect_is(three_level, "WeMixResults")
  
})