test_that("cross-sectional binomial", {
	skip_on_cran()

	plinkfiles <- strsplit(system.file("extdata", "geno.bed", package = "GMMAT"), ".bed", fixed = TRUE)[[1]]
	gdsfile <- system.file("extdata", "geno.gds", package = "GMMAT")
	txtfile <- system.file("extdata", "geno.txt", package = "GMMAT")
	txtfile1 <- system.file("extdata", "geno.txt.gz", package = "GMMAT")
	txtfile2 <- system.file("extdata", "geno.txt.bz2", package = "GMMAT")
	data(example)
	suppressWarnings(RNGversion("3.5.0"))
	set.seed(123)
	snps <- c("SNP10", "SNP25", "SNP1", "SNP0")
	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), ]
	pheno2 <- pheno[pheno$id <= 400, ]
	kins <- diag(500)
	kins[1:400,1:400] <- example$GRM
	rownames(kins) <- colnames(kins) <- 1:500

	# REML Brent
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	out1.bed <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(signif(out1.bed$PVAL), signif(c(0.4303091, 0.8368987, 0.5435321, NA)))
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out1.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	out1.gds <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out1.bed$PVAL, out1.gds$PVAL)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out1.gds, tmpout)
	}
	
	out1.txt <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.bed$PVAL, out1.txt$PVAL)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(out1.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out1.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out1.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out1.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)

	# ML Brent
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	out2.bed <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(signif(out2.bed$PVAL), signif(c(0.3575308, 0.7640698, 0.4627067, NA)))
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out2.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	out2.gds <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out2.bed$PVAL, out2.gds$PVAL)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out2.gds, tmpout)
	}
	
	out2.txt <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.bed$PVAL, out2.txt$PVAL)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(out2.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out2.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out2.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out2.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)

	# REML Nelder-Mead
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	out3.bed <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(signif(out3.bed$PVAL), signif(c(0.4365740, 0.8411170, 0.5471632, NA)))
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out3.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	out3.gds <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out3.bed$PVAL, out3.gds$PVAL)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out3.gds, tmpout)
	}
	
	out3.txt <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.bed$PVAL, out3.txt$PVAL)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out3.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out3.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out3.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out3.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)

	# ML Nelder-Mead
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	out4.bed <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(signif(out4.bed$PVAL), signif(c(0.3575308, 0.7640698, 0.4627067, NA)))
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out4.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	out4.gds <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out4.bed$PVAL, out4.gds$PVAL)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out4.gds, tmpout)
	}
	
	out4.txt <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.bed$PVAL, out4.txt$PVAL)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out4.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out4.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out4.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out4.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)

	# REML AI
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	out5.bed <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(signif(out5.bed$PVAL), signif(c(0.4303094, 0.8368988, 0.5435322, NA)))
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out5.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	out5.gds <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out5.bed$PVAL, out5.gds$PVAL)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out5.gds, tmpout)
	}
	
	out5.txt <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.bed$PVAL, out5.txt$PVAL)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out5.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out5.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out5.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out5.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)

	# kins = NULL
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	out6.bed <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(signif(out6.bed$PVAL), signif(c(0.3575250, 0.7640674, 0.4627067, NA)))
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out6.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	out6.gds <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out6.bed$PVAL, out6.gds$PVAL)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out6.gds, tmpout)
	}
	
	out6.txt <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.bed$PVAL, out6.txt$PVAL)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out6.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out6.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out6.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out6.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)

	idx <- sample(nrow(pheno))
	pheno <- pheno[idx, ]
	idx <- sample(nrow(pheno2))
	pheno2 <- pheno2[idx, ]
	# REML Brent
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(out1.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out1.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out1.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out1.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(out1.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out1.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out1.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out1.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)

	# ML Brent
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(out2.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out2.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out2.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out2.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(out2.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out2.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out2.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out2.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)

	# REML Nelder-Mead
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out3.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out3.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out3.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out3.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out3.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out3.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out3.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out3.gds, tmpout)
	}
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)

	# ML Nelder-Mead
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out4.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out4.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out4.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out4.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out4.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out4.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out4.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out4.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)

	# REML AI
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out5.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out5.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out5.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out5.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out5.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out5.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out5.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out5.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)

	# kins = NULL
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out6.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out6.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out6.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out6.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out6.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out6.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out6.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out6.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)

	idx <- sample(nrow(kins))
	kins <- kins[idx, idx]
	# REML Brent
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(out1.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out1.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out1.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out1.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(out1.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out1.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out1.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out1.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)

	# ML Brent
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(out2.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out2.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out2.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out2.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(out2.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out2.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out2.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out2.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)

	# REML Nelder-Mead
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out3.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out3.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out3.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out3.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out3.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out3.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out3.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out3.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)

	# ML Nelder-Mead
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out4.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out4.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out4.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out4.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out4.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out4.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out4.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out4.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)

	# REML AI
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out5.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out5.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out5.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out5.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out5.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out5.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out5.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out5.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)

	# kins = NULL
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out6.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out6.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out6.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out6.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno2, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out6.bed, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out6.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out6.gds, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out6.gds, tmpout)
	}
	
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
})

test_that("cross-sectional gaussian", {
	skip_on_cran()

	plinkfiles <- strsplit(system.file("extdata", "geno.bed", package = "GMMAT"), ".bed", fixed = TRUE)[[1]]
	gdsfile <- system.file("extdata", "geno.gds", package = "GMMAT")
	txtfile <- system.file("extdata", "geno.txt", package = "GMMAT")
	txtfile1 <- system.file("extdata", "geno.txt.gz", package = "GMMAT")
	txtfile2 <- system.file("extdata", "geno.txt.bz2", package = "GMMAT")
	data(example)
	suppressWarnings(RNGversion("3.5.0"))
	set.seed(123)
	snps <- c("SNP10", "SNP25", "SNP1", "SNP0")
	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), ]
	pheno2 <- pheno[pheno$id <= 400, ]
	kins <- diag(500)
	kins[1:400,1:400] <- example$GRM
	rownames(kins) <- colnames(kins) <- 1:500

	# REML Brent
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	out1.bed <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(signif(out1.bed$PVAL, digits = 5), signif(c(0.3550174, 0.5265875, 0.6089051, NA), digits = 5))
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out1.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	out1.gds <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out1.bed$PVAL, out1.gds$PVAL)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out1.gds, tmpout)
	}
	
	out1.txt <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.bed$PVAL, out1.txt$PVAL)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(out1.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out1.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out1.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out1.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)

	# ML Brent
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	out2.bed <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(signif(out2.bed$PVAL), signif(c(0.3652024, 0.5347624, 0.5717466, NA)))
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out2.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	out2.gds <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out2.bed$PVAL, out2.gds$PVAL)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out2.gds, tmpout)
	}
	
	out2.txt <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.bed$PVAL, out2.txt$PVAL)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(out2.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out2.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out2.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out2.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)

	# REML Nelder-Mead
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	out3.bed <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(signif(out3.bed$PVAL, digits = 5), signif(c(0.3485413, 0.5254271, 0.6174975, NA), digits = 5))
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out3.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	out3.gds <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out3.bed$PVAL, out3.gds$PVAL)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out3.gds, tmpout)
	}
	
	out3.txt <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.bed$PVAL, out3.txt$PVAL)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out3.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out3.bed, tmpout)

	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out3.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out3.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)

	# ML Nelder-Mead
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	out4.bed <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(signif(out4.bed$PVAL), signif(c(0.3650012, 0.5376902, 0.5786932, NA)))
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out4.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	out4.gds <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out4.bed$PVAL, out4.gds$PVAL)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out4.gds, tmpout)
	}
	
	out4.txt <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.bed$PVAL, out4.txt$PVAL)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out4.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out4.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out4.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out4.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)

	# REML AI
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	out5.bed <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(signif(out5.bed$PVAL), signif(c(0.3550172, 0.5265874, 0.6089053, NA)))
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out5.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	out5.gds <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out5.bed$PVAL, out5.gds$PVAL)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out5.gds, tmpout)
	}
	
	out5.txt <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.bed$PVAL, out5.txt$PVAL)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out5.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out5.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out5.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out5.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)

	# kins = NULL
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	out6.bed <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(signif(out6.bed$PVAL), signif(c(0.9928327, 0.9166500, 0.1103733, NA)))
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out6.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	out6.gds <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out6.bed$PVAL, out6.gds$PVAL)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out6.gds, tmpout)
	}
	
	out6.txt <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.bed$PVAL, out6.txt$PVAL)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out6.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out6.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out6.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out6.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)

	idx <- sample(nrow(pheno))
	pheno <- pheno[idx, ]
	idx <- sample(nrow(pheno2))
	pheno2 <- pheno2[idx, ]
	# REML Brent
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(out1.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out1.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out1.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out1.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(signif(out1.bed$PVAL, digits = 5), signif(tmpout$PVAL, digits = 5))
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(signif(out1.bed$PVAL, digits = 5), signif(tmpout$PVAL, digits = 5))
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(signif(out1.gds$PVAL, digits = 5), signif(tmpout$PVAL, digits = 5))
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(signif(out1.gds$PVAL, digits = 5), signif(tmpout$PVAL, digits = 5))
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(signif(out1.txt$PVAL, digits = 5), signif(tmpout$PVAL, digits = 5))
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(signif(out1.txt$PVAL, digits = 5), signif(tmpout$PVAL, digits = 5))
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(signif(out1.txt$PVAL, digits = 5), signif(tmpout$PVAL, digits = 5))

	# ML Brent
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(out2.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out2.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out2.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out2.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(out2.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out2.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out2.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out2.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)

	# REML Nelder-Mead
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out3.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out3.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out3.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out3.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out3.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out3.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out3.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out3.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)

	# ML Nelder-Mead
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out4.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out4.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out4.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out4.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out4.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out4.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out4.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out4.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)

	# REML AI
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out5.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out5.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out5.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out5.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out5.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out5.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out5.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out5.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)

	# kins = NULL
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out6.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out6.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out6.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out6.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out6.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out6.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out6.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out6.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)

	idx <- sample(nrow(kins))
	kins <- kins[idx, idx]
	# REML Brent
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(out1.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out1.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out1.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out1.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out1.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(signif(out1.bed$PVAL, digits = 5), signif(tmpout$PVAL, digits = 5))
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(signif(out1.bed$PVAL, digits = 5), signif(tmpout$PVAL, digits = 5))
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(signif(out1.gds$PVAL, digits = 5), signif(tmpout$PVAL, digits = 5))
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(signif(out1.gds$PVAL, digits = 5), signif(tmpout$PVAL, digits = 5))
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(signif(out1.txt$PVAL, digits = 5), signif(tmpout$PVAL, digits = 5))
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(signif(out1.txt$PVAL, digits = 5), signif(tmpout$PVAL, digits = 5))
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(signif(out1.txt$PVAL, digits = 5), signif(tmpout$PVAL, digits = 5))

	# ML Brent
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(out2.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out2.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out2.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out2.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps)
	expect_equal(out2.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out2.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps)
	expect_equal(out2.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=gdsfile,snps=snps,select=select)
	expect_equal(out2.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Brent", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out2.txt, tmpout)

	# REML Nelder-Mead
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out3.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out3.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out3.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out3.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out3.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out3.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out3.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out3.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out3.txt, tmpout)

	# ML Nelder-Mead
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out4.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out4.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out4.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out4.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps)
	expect_equal(out4.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out4.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps)
	expect_equal(out4.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=gdsfile,snps=snps,select=select)
	expect_equal(out4.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "ML", method.optim = "Nelder-Mead", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out4.txt, tmpout)

	# REML AI
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out5.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out5.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out5.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out5.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out5.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out5.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out5.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out5.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)

	# kins = NULL
	select <- match(1:400, unique(pheno2[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out6.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out6.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out6.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out6.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno2, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)

	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out6.bed, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out6.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out6.gds, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out6.gds, tmpout)
	}
	
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
})

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

	plinkfiles <- strsplit(system.file("extdata", "geno.bed", package = "GMMAT"), ".bed", fixed = TRUE)[[1]]
	gdsfile <- system.file("extdata", "geno.gds", package = "GMMAT")
	txtfile <- system.file("extdata", "geno.txt", package = "GMMAT")
	txtfile1 <- system.file("extdata", "geno.txt.gz", package = "GMMAT")
	txtfile2 <- system.file("extdata", "geno.txt.bz2", package = "GMMAT")
	data(example)
	suppressWarnings(RNGversion("3.5.0"))
	set.seed(123)
	snps <- c("SNP10", "SNP25", "SNP1", "SNP0")
        pheno <- example$pheno2
        kins <- example$GRM

	# REML AI
	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	out5.bed <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(signif(out5.bed$PVAL), signif(c(0.0730724, 0.9451040, 0.3083600, NA)))
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out5.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	out5.gds <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out5.bed$PVAL, out5.gds$PVAL)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out5.gds, tmpout)
	}
	
	out5.txt <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.bed$PVAL, out5.txt$PVAL)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)

	# kins = NULL
	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	out6.bed <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(signif(out6.bed$PVAL), signif(c(0.01714963, 0.48756105, 0.96724149, NA)))
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out6.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	out6.gds <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out6.bed$PVAL, out6.gds$PVAL)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out6.gds, tmpout)
	}
	
	out6.txt <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.bed$PVAL, out6.txt$PVAL)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)

	idx <- sample(nrow(pheno))
	pheno <- pheno[idx, ]
	# REML AI
	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out5.bed, tmpout)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out5.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out5.gds, tmpout)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out5.gds, tmpout)
	}
	
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)

	# kins = NULL
	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out6.bed, tmpout)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out6.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out6.gds, tmpout)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out6.gds, tmpout)
	}
	
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)

	idx <- sample(nrow(kins))
	kins <- kins[idx, idx]
	# REML AI
	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out5.bed, tmpout)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out5.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out5.gds, tmpout)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out5.gds, tmpout)
	}
	
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = kins, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)

	# kins = NULL
	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out6.bed, tmpout)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out6.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out6.gds, tmpout)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out6.gds, tmpout)
	}
	
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(y.repeated ~ sex, data = pheno, kins = NULL, id = "id", random.slope = NULL, family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
})

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

	plinkfiles <- strsplit(system.file("extdata", "geno.bed", package = "GMMAT"), ".bed", fixed = TRUE)[[1]]
	gdsfile <- system.file("extdata", "geno.gds", package = "GMMAT")
	txtfile <- system.file("extdata", "geno.txt", package = "GMMAT")
	txtfile1 <- system.file("extdata", "geno.txt.gz", package = "GMMAT")
	txtfile2 <- system.file("extdata", "geno.txt.bz2", package = "GMMAT")
	data(example)
	suppressWarnings(RNGversion("3.5.0"))
	set.seed(123)
	snps <- c("SNP10", "SNP25", "SNP1", "SNP0")
        pheno <- example$pheno2
        kins <- example$GRM

	# REML AI
	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	out5.bed <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(signif(out5.bed$PVAL), signif(c(0.2376449, 0.2081449, 0.1638394, NA)))
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out5.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	out5.gds <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out5.bed$PVAL, out5.gds$PVAL)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out5.gds, tmpout)
	}
	
	out5.txt <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.bed$PVAL, out5.txt$PVAL)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)

	# kins = NULL
	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	out6.bed <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(signif(out6.bed$PVAL), signif(c(0.2106602, 0.4357301, 0.2180844, NA)))
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out6.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	out6.gds <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out6.bed$PVAL, out6.gds$PVAL)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out6.gds, tmpout)
	}
	
	out6.txt <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.bed$PVAL, out6.txt$PVAL)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)

	idx <- sample(nrow(pheno))
	pheno <- pheno[idx, ]
	# REML AI
	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out5.bed, tmpout)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out5.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out5.gds, tmpout)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out5.gds, tmpout)
	}
	
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)

	# kins = NULL
	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out6.bed, tmpout)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out6.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out6.gds, tmpout)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out6.gds, tmpout)
	}
	
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)

	idx <- sample(nrow(kins))
	kins <- kins[idx, idx]
	# REML AI
	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out5.bed, tmpout)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out5.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out5.gds, tmpout)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out5.gds, tmpout)
	}
	
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = kins, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out5.txt, tmpout)

	# kins = NULL
	select <- match(1:400, unique(pheno[,"id"]))
	select[is.na(select)] <- 0
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps)
	expect_equal(out6.bed, tmpout)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=plinkfiles,snps=snps,select=select)
	expect_equal(out6.bed, tmpout)
	
	if(requireNamespace("SeqArray", quietly = TRUE) && requireNamespace("SeqVarTools", quietly = TRUE)) {
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps)
	expect_equal(out6.gds, tmpout)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=gdsfile,snps=snps,select=select)
	expect_equal(out6.gds, tmpout)
	}
	
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile1,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
	tmpout <- glmm.wald(y.trend ~ sex + time, data = pheno, kins = NULL, id = "id", random.slope = "time", family = gaussian(link = "identity"), method = "REML", method.optim = "AI", infile=txtfile2,snps=snps,infile.nrow.skip = 5, infile.ncol.skip = 3, infile.ncol.print = 1:3,select=select,infile.header.print = c("SNP", "Allele1", "Allele2"))
	expect_equal(out6.txt, tmpout)
})