test_that("cross-sectional id le 400 binomial", {
	skip_on_cran()

	data(example)
	suppressWarnings(RNGversion("3.5.0"))
	set.seed(123)
	pheno <- rbind(example$pheno, example$pheno[1:100, ])
	pheno$id <- 1:500
	pheno$disease[sample(1:500,20)] <- NA
	pheno$age[sample(1:500,20)] <- NA
	pheno$sex[sample(1:500,20)] <- NA
	pheno <- pheno[sample(1:500,450), ]
	pheno <- pheno[pheno$id <= 400, ]
	kins <- example$GRM
	obj1 <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent")
	expect_equal(signif(as.numeric(obj1$theta)), signif(c(1, 0.1925217)))
	expect_equal(signif(as.numeric(obj1$coef)), signif(c(1.01676248, -0.01506251, -0.33240660)))
	obj2 <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent")
	expect_equal(as.numeric(obj2$theta), c(1, 0))
	expect_equal(signif(as.numeric(obj2$coef)), signif(c(0.94253959, -0.01429532, -0.32823930)))
	obj3 <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead")
	expect_equal(signif(as.numeric(obj3$theta)), c(1, 0.1925))
	expect_equal(signif(as.numeric(obj3$coef)), signif(c(1.0224921, -0.0151259, -0.3328061)))
	obj4 <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead")
	expect_equal(as.numeric(obj4$theta), c(1, 0))
	expect_equal(signif(as.numeric(obj4$coef)), signif(c(0.94253959, -0.01429532, -0.32823930)))
	obj5 <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI")
	expect_equal(signif(as.numeric(obj5$theta)), signif(c(1, 0.1925225)))
	expect_equal(signif(as.numeric(obj5$coef)), signif(c(1.01676230, -0.01506251, -0.33240659)))
	obj6 <- glmmkin(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI")
	expect_equal(c(obj6$theta), 1)
	expect_equal(signif(as.numeric(obj6$coef)), signif(c(0.94253958, -0.01429532, -0.32823930)))
	obj <- glm(disease ~ age + sex, data = pheno, family = binomial(link = "logit"))
	expect_equal(obj6$coef, obj$coef)

	idx <- sample(nrow(pheno))
	pheno <- pheno[idx, ]
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent")
	expect_equal(obj1$theta, obj$theta)
	expect_equal(obj1$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent")
	expect_equal(obj2$theta, obj$theta)
	expect_equal(obj2$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead")
	expect_equal(obj3$theta, obj$theta)
	expect_equal(obj3$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead")
	expect_equal(obj4$theta, obj$theta)
	expect_equal(obj4$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI")
	expect_equal(obj5$theta, obj$theta)
	expect_equal(obj5$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI")
	expect_equal(obj6$theta, obj$theta)
	expect_equal(obj6$coef, obj$coef)
	obj <- glm(disease ~ age + sex, data = pheno, family = binomial(link = "logit"))
	expect_equal(obj6$coef, obj$coef)

	idx <- sample(nrow(kins))
	kins <- kins[idx, idx]
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent")
	expect_equal(obj1$theta, obj$theta)
	expect_equal(obj1$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent")
	expect_equal(obj2$theta, obj$theta)
	expect_equal(obj2$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead")
	expect_equal(obj3$theta, obj$theta)
	expect_equal(obj3$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead")
	expect_equal(obj4$theta, obj$theta)
	expect_equal(obj4$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI")
	expect_equal(obj5$theta, obj$theta)
	expect_equal(obj5$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI")
	expect_equal(obj6$theta, obj$theta)
	expect_equal(obj6$coef, obj$coef)
})

test_that("cross-sectional id gt 400 binomial", {
	skip_on_cran()

	data(example)
	suppressWarnings(RNGversion("3.5.0"))
	set.seed(123)
	pheno <- rbind(example$pheno, example$pheno[1:100, ])
	pheno$id <- 1:500
	pheno$disease[sample(1:500,20)] <- NA
	pheno$age[sample(1:500,20)] <- NA
	pheno$sex[sample(1:500,20)] <- NA
	pheno <- pheno[sample(1:500,450), ]
	kins <- diag(500)
	kins[1:400, 1:400] <- example$GRM
	rownames(kins) <- colnames(kins) <- 1:500
	obj1 <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent")
	expect_equal(signif(as.numeric(obj1$theta)), signif(c(1, 0.1263608)))
	expect_equal(signif(as.numeric(obj1$coef)), signif(c(0.92300941, -0.01457307, -0.18165857)))
	obj2 <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent")
	expect_equal(as.numeric(obj2$theta), c(1, 0))
	expect_equal(signif(as.numeric(obj2$coef)), signif(c(0.86840326, -0.01402939, -0.17898775)))
	obj3 <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead")
	expect_equal(signif(as.numeric(obj3$theta)), c(1, 0.121))
	expect_equal(signif(as.numeric(obj3$coef), digits = 5), signif(c(0.92527071, -0.01459762, -0.18179150), digits = 5))
	obj4 <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead")
	expect_equal(as.numeric(obj4$theta), c(1, 0))
	expect_equal(signif(as.numeric(obj4$coef)), signif(c(0.86840326, -0.01402939, -0.17898775)))
	obj5 <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI")
	expect_equal(signif(as.numeric(obj5$theta)), signif(c(1, 0.1263611)))
	expect_equal(signif(as.numeric(obj5$coef)), signif(c(0.92300945, -0.01457307, -0.18165858)))
	obj6 <- glmmkin(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI")
	expect_equal(c(obj6$theta), 1)
	expect_equal(signif(as.numeric(obj6$coef)), signif(c(0.86840325, -0.01402939, -0.17898775)))
	obj <- glm(disease ~ age + sex, data = pheno, family = binomial(link = "logit"))
	expect_equal(obj6$coef, obj$coef)

	idx <- sample(nrow(pheno))
	pheno <- pheno[idx, ]
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent")
	expect_equal(obj1$theta, obj$theta)
	expect_equal(obj1$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent")
	expect_equal(obj2$theta, obj$theta)
	expect_equal(obj2$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead")
	expect_equal(obj3$theta, obj$theta)
	expect_equal(obj3$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead")
	expect_equal(obj4$theta, obj$theta)
	expect_equal(obj4$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI")
	expect_equal(obj5$theta, obj$theta)
	expect_equal(obj5$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI")
	expect_equal(obj6$theta, obj$theta)
	expect_equal(obj6$coef, obj$coef)
	obj <- glm(disease ~ age + sex, data = pheno, family = binomial(link = "logit"))
	expect_equal(obj6$coef, obj$coef)

	idx <- sample(nrow(kins))
	kins <- kins[idx, idx]
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent")
	expect_equal(obj1$theta, obj$theta)
	expect_equal(obj1$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent")
	expect_equal(obj2$theta, obj$theta)
	expect_equal(obj2$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead")
	expect_equal(obj3$theta, obj$theta)
	expect_equal(obj3$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead")
	expect_equal(obj4$theta, obj$theta)
	expect_equal(obj4$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI")
	expect_equal(obj5$theta, obj$theta)
	expect_equal(obj5$coef, obj$coef)
	obj <- glmmkin(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI")
	expect_equal(obj6$theta, obj$theta)
	expect_equal(obj6$coef, obj$coef)
})

test_that("cross-sectional id le 400 gaussian", {
	skip_on_cran()

	data(example)
	suppressWarnings(RNGversion("3.5.0"))
	set.seed(123)
	pheno <- rbind(example$pheno, example$pheno[1:100, ])
	pheno$id <- 1:500
	pheno$disease[sample(1:500,20)] <- NA
	pheno$age[sample(1:500,20)] <- NA
	pheno$sex[sample(1:500,20)] <- NA
	pheno <- pheno[sample(1:500,450), ]
	pheno <- pheno[pheno$id <= 400, ]
	kins <- example$GRM
	obj1 <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent")
	expect_equal(signif(as.numeric(obj1$theta)), signif(c(0.7682478, 1.3037470)))
	expect_equal(signif(as.numeric(obj1$coef)), signif(c(3.7634934, 0.0346562, 0.3062784)))
	obj2 <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent")
	expect_equal(signif(as.numeric(obj2$theta)), signif(c(0.9391236, 1.0851816)))
	expect_equal(signif(as.numeric(obj2$coef)), signif(c(3.76521261, 0.03451443, 0.30509914)))
	obj3 <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead")
	expect_equal(signif(as.numeric(obj3$theta), digits = 5), signif(c(0.7174382, 1.2196450), digits = 5))
	expect_equal(signif(as.numeric(obj3$coef)), signif(c(3.76315625, 0.03469074, 0.30656586)))
	obj4 <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead")
	expect_equal(signif(as.numeric(obj4$theta)), signif(c(0.9622361, 1.0584597)))
	expect_equal(signif(as.numeric(obj4$coef)), signif(c(3.76551377, 0.03449328, 0.30492363)))
	obj5 <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI")
	expect_equal(signif(as.numeric(obj5$theta)), signif(c(0.7682481, 1.3037467)))
	expect_equal(signif(as.numeric(obj5$coef)), signif(c(3.7634933, 0.0346562, 0.3062784)))
	obj6 <- glmmkin(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI")
	expect_equal(signif(as.numeric(obj6$theta)), signif(1.996857))
	expect_equal(signif(as.numeric(obj6$coef)), signif(c(3.89665633, 0.03156906, 0.27860778)))
	obj <- lm(trait ~ age + sex, data = pheno)
	expect_equal(c(obj6$theta), summary(obj)$sigma^2)
	expect_equal(obj6$coef, obj$coef)

	idx <- sample(nrow(pheno))
	pheno <- pheno[idx, ]
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent")
	expect_equal(obj1$theta, obj$theta)
	expect_equal(obj1$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent")
	expect_equal(obj2$theta, obj$theta)
	expect_equal(obj2$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead")
	expect_equal(obj3$theta, obj$theta)
	expect_equal(obj3$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead")
	expect_equal(obj4$theta, obj$theta)
	expect_equal(obj4$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI")
	expect_equal(obj5$theta, obj$theta)
	expect_equal(obj5$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI")
	expect_equal(obj6$theta, obj$theta)
	expect_equal(obj6$coef, obj$coef)
	obj <- lm(trait ~ age + sex, data = pheno)
	expect_equal(c(obj6$theta), summary(obj)$sigma^2)
	expect_equal(obj6$coef, obj$coef)

	idx <- sample(nrow(kins))
	kins <- kins[idx, idx]
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent")
	expect_equal(obj1$theta, obj$theta)
	expect_equal(obj1$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent")
	expect_equal(obj2$theta, obj$theta)
	expect_equal(obj2$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead")
	expect_equal(obj3$theta, obj$theta)
	expect_equal(obj3$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead")
	expect_equal(obj4$theta, obj$theta)
	expect_equal(obj4$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI")
	expect_equal(obj5$theta, obj$theta)
	expect_equal(obj5$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI")
	expect_equal(obj6$theta, obj$theta)
	expect_equal(obj6$coef, obj$coef)
})

test_that("cross-sectional id gt 400 gaussian", {
	skip_on_cran()

	data(example)
	suppressWarnings(RNGversion("3.5.0"))
	set.seed(123)
	pheno <- rbind(example$pheno, example$pheno[1:100, ])
	pheno$id <- 1:500
	pheno$disease[sample(1:500,20)] <- NA
	pheno$age[sample(1:500,20)] <- NA
	pheno$sex[sample(1:500,20)] <- NA
	pheno <- pheno[sample(1:500,450), ]
	kins <- diag(500)
	kins[1:400, 1:400] <- example$GRM
	rownames(kins) <- colnames(kins) <- 1:500
	obj1 <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent")
	expect_equal(signif(as.numeric(obj1$theta)), signif(c(0.7205609, 1.3795640)))
	expect_equal(signif(as.numeric(obj1$coef)), signif(c(3.90459575, 0.03650594, 0.37958367)))
	obj2 <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent")
	expect_equal(signif(as.numeric(obj2$theta)), signif(c(0.7485909, 1.3321148)))
	expect_equal(signif(as.numeric(obj2$coef)), signif(c(3.90396978, 0.03649397, 0.37960264)))
	obj3 <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead")
	expect_equal(signif(as.numeric(obj3$theta)), signif(c(0.6044086, 1.2692581)))
	expect_equal(signif(as.numeric(obj3$coef)), signif(c(3.90659373, 0.03654332, 0.37947361)))
	obj4 <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead")
	expect_equal(signif(as.numeric(obj4$theta)), signif(c(0.7188051, 1.2219687)))
	expect_equal(signif(as.numeric(obj4$coef)), signif(c(3.90453205, 0.03650473, 0.37958592)))
	obj5 <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI")
	expect_equal(signif(as.numeric(obj5$theta)), signif(c(0.7205609, 1.3795639)))
	expect_equal(signif(as.numeric(obj5$coef)), signif(c(3.90459568, 0.03650594, 0.37958367)))
	obj6 <- glmmkin(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI")
	expect_equal(signif(as.numeric(obj6$theta)), signif(2.051809))
	expect_equal(signif(as.numeric(obj6$coef)), signif(c(3.7836398, 0.0344577, 0.3602852)))
	obj <- lm(trait ~ age + sex, data = pheno)
	expect_equal(c(obj6$theta), summary(obj)$sigma^2)
	expect_equal(obj6$coef, obj$coef)

	idx <- sample(nrow(pheno))
	pheno <- pheno[idx, ]
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent")
	expect_equal(obj1$theta, obj$theta)
	expect_equal(obj1$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent")
	expect_equal(obj2$theta, obj$theta)
	expect_equal(obj2$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead")
	expect_equal(obj3$theta, obj$theta)
	expect_equal(obj3$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead")
	expect_equal(obj4$theta, obj$theta)
	expect_equal(obj4$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI")
	expect_equal(obj5$theta, obj$theta)
	expect_equal(obj5$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI")
	expect_equal(obj6$theta, obj$theta)
	expect_equal(obj6$coef, obj$coef)
	obj <- lm(trait ~ age + sex, data = pheno)
	expect_equal(c(obj6$theta), summary(obj)$sigma^2)
	expect_equal(obj6$coef, obj$coef)

	idx <- sample(nrow(kins))
	kins <- kins[idx, idx]
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent")
	expect_equal(obj1$theta, obj$theta)
	expect_equal(obj1$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent")
	expect_equal(obj2$theta, obj$theta)
	expect_equal(obj2$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead")
	expect_equal(obj3$theta, obj$theta)
	expect_equal(obj3$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead")
	expect_equal(obj4$theta, obj$theta)
	expect_equal(obj4$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI")
	expect_equal(obj5$theta, obj$theta)
	expect_equal(obj5$coef, obj$coef)
	obj <- glmmkin(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI")
	expect_equal(obj6$theta, obj$theta)
	expect_equal(obj6$coef, obj$coef)
})

test_that("longitudinal repeated measures gaussian", {
	skip_on_cran()

	data(example)
	suppressWarnings(RNGversion("3.5.0"))
	set.seed(123)
	pheno <- example$pheno2
	kins <- example$GRM
	obj1 <- glmmkin(y.repeated ~ sex, data = pheno, kins = kins, id = "id",random.slope = NULL, family = gaussian(link = "identity"))
	expect_equal(signif(as.numeric(obj1$theta)), signif(c(0.5089983, 0.2410866, 0.2129918)))
	expect_equal(signif(as.numeric(obj1$coef)), signif(c(0.9871441, 0.6134051)))
	obj2 <- glmmkin(y.repeated ~ sex, data = pheno, kins = NULL, id = "id",random.slope = NULL, family = gaussian(link = "identity"))
	expect_equal(signif(as.numeric(obj2$theta)), signif(c(0.5089982, 0.4374597)))
	expect_equal(signif(as.numeric(obj2$coef)), signif(c(1.0237, 0.6141)))

	idx <- sample(nrow(pheno))
	pheno <- pheno[idx, ]
	obj <- glmmkin(y.repeated ~ sex, data = pheno, kins = kins, id = "id",random.slope = NULL, family = gaussian(link = "identity"))
	expect_equal(obj1$theta, obj$theta)
	expect_equal(obj1$coef, obj$coef)
	obj <- glmmkin(y.repeated ~ sex, data = pheno, kins = NULL, id = "id",random.slope = NULL, family = gaussian(link = "identity"))
	expect_equal(obj2$theta, obj$theta)
	expect_equal(obj2$coef, obj$coef)

	idx <- sample(nrow(kins))
	kins <- kins[idx, idx]
	obj <- glmmkin(y.repeated ~ sex, data = pheno, kins = kins, id = "id",random.slope = NULL, family = gaussian(link = "identity"))
	expect_equal(obj1$theta, obj$theta)
	expect_equal(obj1$coef, obj$coef)
	obj <- glmmkin(y.repeated ~ sex, data = pheno, kins = NULL, id = "id",random.slope = NULL, family = gaussian(link = "identity"))
	expect_equal(obj2$theta, obj$theta)
	expect_equal(obj2$coef, obj$coef)
})

test_that("longitudinal random time trend gaussian", {
	skip_on_cran()

	data(example)
	suppressWarnings(RNGversion("3.5.0"))
	set.seed(123)
	pheno <- example$pheno2
	kins <- example$GRM
	obj1 <- glmmkin(y.trend ~ sex + time, data = pheno, kins = kins, id = "id",random.slope = "time", family = gaussian(link = "identity"))
	expect_equal(signif(as.numeric(obj1$theta), digits = 5), signif(c(2.0139758, 1.5056913, 0.5443502, 1.2539360, -0.2469174, 1.1690885, 0.6931906), digits = 5))
	expect_equal(signif(as.numeric(obj1$coef)), signif(c(1.1834620, 0.6320665, 1.0701016)))
	obj2 <- glmmkin(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id",random.slope = "time", family = gaussian(link = "identity"))
	expect_equal(signif(as.numeric(obj2$theta)), signif(c(2.0139758, 1.8949301, 0.8748916, 1.7403904)))
	expect_equal(signif(as.numeric(obj2$coef), digits = 5), signif(c(1.1532040, 0.6177671, 1.0189475), digits = 5))

	idx <- sample(nrow(pheno))
	pheno <- pheno[idx, ]
	obj <- glmmkin(y.trend ~ sex + time, data = pheno, kins = kins, id = "id",random.slope = "time", family = gaussian(link = "identity"))
	expect_equal(obj1$theta, obj$theta)
	expect_equal(obj1$coef, obj$coef)
	obj <- glmmkin(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id",random.slope = "time", family = gaussian(link = "identity"))
	expect_equal(obj2$theta, obj$theta)
	expect_equal(obj2$coef, obj$coef)

	idx <- sample(nrow(kins))
	kins <- kins[idx, idx]
	obj <- glmmkin(y.trend ~ sex + time, data = pheno, kins = kins, id = "id",random.slope = "time", family = gaussian(link = "identity"))
	expect_equal(obj1$theta, obj$theta)
	expect_equal(obj1$coef, obj$coef)
	obj <- glmmkin(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id",random.slope = "time", family = gaussian(link = "identity"))
	expect_equal(obj2$theta, obj$theta)
	expect_equal(obj2$coef, obj$coef)
})

test_that("multiple phenotypes gaussian", {
	skip_on_cran()

	data(example)
	suppressWarnings(RNGversion("3.5.0"))
	set.seed(103)
	kins <- example$GRM
	tau1 <- matrix(c(3,0.5,0,0.5,2.5,-0.1,0,-0.1,3),3,3)
	tau2 <- matrix(c(2.5,0.8,0.2,0.8,4.8,-0.1,0.2,-0.1,2.8),3,3)
	kins.chol <- chol(tau1 %x% kins + tau2 %x% diag(400))
	tmp <- as.vector(crossprod(kins.chol, rnorm(1200)))
	x1 <- rnorm(400)
	x2 <- rbinom(400,1,0.5)
	pheno <- data.frame(id = 1:400, x1 = x1, x2 = x2, y1 = 0.5*x1+0.8*x2+tmp[1:400], y2 = x1-0.3*x2+tmp[401:800], y3 = x2+tmp[801:1200])
	obj1 <- glmmkin(cbind(y1,y2,y3)~x1+x2, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"))
	expect_equal(signif(c(as.vector(obj1$theta[[1]]),as.vector(obj1$theta[[2]])), digits = 5), signif(c(3.13246847, 1.52966809, 0.27819033, 1.52966809, 5.21123625, -1.87847885, 0.27819033, -1.87847885, 3.04660044, 1.69332546, -0.36423974, -0.07708463, -0.36423974, 2.07095534, 1.62893810, -0.07708463, 1.62893810, 2.61803110), digits = 5))
	expect_equal(signif(c(obj1$coef)), signif(c(0.20216323, 0.30601599, 0.74140916, 0.04119201, 1.04136078, -0.40556256, 1.08891927, 0.02044031, 1.01281152)))
	obj2 <- glmmkin(cbind(y1,y2,y3)~x1+x2, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"))
	expect_equal(signif(c(obj2$theta[[1]])), signif(c(4.6708159, 1.1614428, 0.1580549, 1.1614428, 7.2537122, -0.2297423, 0.1580549, -0.2297423, 5.6175721)))
	expect_equal(signif(c(obj2$coef), digits = 5), signif(c(0.13897937, 0.31965037, 0.77867243, 0.26349000, 1.03272869, -0.46414252, 1.25165117, 0.01778239, 0.95851562), digits = 5))

	idx <- sample(nrow(pheno))
	pheno <- pheno[idx, ]
	obj <- glmmkin(cbind(y1,y2,y3)~x1+x2, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"))
	expect_equal(obj1$theta, obj$theta)
	expect_equal(obj1$coef, obj$coef)
	obj <- glmmkin(cbind(y1,y2,y3)~x1+x2, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"))
	expect_equal(obj2$theta, obj$theta)
	expect_equal(obj2$coef, obj$coef)

	idx <- sample(nrow(kins))
	kins <- kins[idx, idx]
	obj <- glmmkin(cbind(y1,y2,y3)~x1+x2, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"))
	expect_equal(obj1$theta, obj$theta)
	expect_equal(obj1$coef, obj$coef)
	obj <- glmmkin(cbind(y1,y2,y3)~x1+x2, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"))
	expect_equal(obj2$theta, obj$theta)
	expect_equal(obj2$coef, obj$coef)
})