test_that("Heteroscedastic model works", {
  mod <- galamm(
    formula = y ~ x + (1 | id),
    weights = ~ (1 | item),
    data = hsced
  )

  expect_snapshot(print(summary(mod), digits = 3))
  expect_equal(mod$model$loglik, -2058.14021326104)
  expect_equal(
    summary(mod)$AICtab,
    c(
      AIC = 4126.28042652208, BIC = 4151.73081070096,
      logLik = -2058.14021326104,
      deviance = 4116.28042652208, df.resid = 1195
    )
  )
  expect_equal(summary(mod)$Lambda, NULL)

  expect_equal(
    residuals(mod)[c(4, 8, 11)],
    c(0.879617156429014, -1.83921105781445, 0.769190482966503)
  )

  expect_equal(vcov(mod, parm = "weights"),
    structure(0.00245961267414418, dim = c(1L, 1L)),
    tolerance = 1e-4
  )

  expect_equal(vcov(mod, parm = 1L),
    structure(0.00543386086733993, dim = c(1L, 1L)),
    tolerance = 1e-4
  )

  expect_error(vcov(mod, parm = 5L), "out of bounds")
  expect_error(vcov(mod, parm = "phi"), "Parameter not found")

  expect_error(
    galamm(
      formula = y ~ x + (1 | id),
      weights = list(~ (1 | item)),
      data = hsced
    ),
    "weights must be a formula"
  )


  # Now use initial values
  mod_start <- galamm(
    formula = y ~ x + (1 | id),
    weights = ~ (1 | item),
    data = hsced,
    start = list(beta = c(.13, .70), theta = 1.01, weights = .501)
  )

  expect_equal(
    mod_start$parameters$parameter_estimates,
    mod$parameters$parameter_estimates,
    tolerance = .01
  )
})

test_that("Heteroscedastic model works with more than one group", {
  set.seed(123)
  hsced2 <- hsced
  hsced2$item[hsced2$item == 2] <- sample(c(2, 3), sum(hsced2$item == 2), replace = TRUE)

  mod <- galamm(
    formula = y ~ x + (1 | id),
    weights = ~ (1 | item),
    data = hsced2
  )

  expect_snapshot(print(summary(mod), digits = 3))
  expect_equal(mod$model$loglik, -2058.13575138403, tolerance = 1e-4)
  expect_equal(
    summary(mod)$AICtab,
    c(AIC = 4128.27150276806, BIC = 4158.81196378272,
      logLik = -2058.13575138403,
      deviance = 4116.27150276806, df.resid = 1194),
    tolerance = 1e-4
  )
  expect_equal(summary(mod)$Lambda, NULL)

  expect_equal(
    residuals(mod)[c(4, 8, 11)],
    c(0.881516323889496, -1.83825643847652, 0.77018964982384)
  )

  expect_equal(vcov(mod, parm = "weights"),
               structure(c(0.00341412313076528, 0.00142678601570138, 0.00142678601570138,
                           0.00357562468779565), dim = c(2L, 2L)),
               tolerance = 1e-4
  )

  expect_equal(vcov(mod, parm = 1L),
               structure(0.00543350591426884, dim = c(1L, 1L)),
               tolerance = 1e-4
  )

  expect_error(vcov(mod, parm = 6L), "out of bounds")
  expect_error(vcov(mod, parm = "phi"), "Parameter not found")

  # Now use initial values
  mod_start <- galamm(
    formula = y ~ x + (1 | id),
    weights = ~ (1 | item),
    data = hsced2,
    start = list(beta = c(.13, .70), theta = 1.01, weights = c(.501, .44))
  )

  expect_equal(
    mod_start$parameters$parameter_estimates,
    mod$parameters$parameter_estimates,
    tolerance = .01
  )
})