test_that("cross-sectional gaussian", { skip_if_not_installed("SeqArray") skip_if_not_installed("SeqVarTools") gdsfile <- system.file("extdata", "geno.gds", package = "MAGEE") bgenfile <- system.file("extdata", "geno.bgen", package = "MAGEE") samplefile <- system.file("extdata", "geno.sample", package = "MAGEE") group.file <- system.file("extdata", "SetID.singlevariant.txt", package = "MAGEE") data(example) set.seed(123) pheno <- rbind(example$pheno, example$pheno[1:100,]) pheno$id <- 1:500 pheno$disease[sample(1:500,20)] <- NA pheno$age[sample(1:500,20)] <- NA pheno$sex[sample(1:500,20)] <- NA pheno <- pheno[sample(1:500,450), ] pheno <- pheno[pheno$id <= 400, ] pheno$obese <- rbinom(nrow(pheno), 1, 0.3) kins <- example$GRM ### single thread with kins and one interaction variable obj1 <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity")) out1 <- MAGEE(null.obj=obj1, interaction="sex", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj1.outfile.gds <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=gdsfile, outfile=obj1.outfile.gds) obj1.gds <- read.table(obj1.outfile.gds, header = T, as.is = T) expect_equal(signif(range(obj1.gds$P_Value_Marginal)), signif(c( 0.007431214, 0.993737081))) expect_equal(signif(range(obj1.gds$P_Value_Interaction)), signif(c( 0.02633466, 0.99253758))) expect_equal(signif(range(obj1.gds$P_Value_Joint)), signif(c(0.00626638, 0.98730215))) expect_equal(signif(out1$MV.pval,4), signif(obj1.gds$P_Value_Marginal,4)) expect_equal(signif(out1$MF.pval), signif(obj1.gds$P_Value_Marginal)) expect_equal(signif(out1$IV.pval,4), signif(obj1.gds$P_Value_Interaction,4)) expect_equal(signif(out1$IF.pval), signif(obj1.gds$P_Value_Interaction)) expect_equal(signif(pchisq(qchisq(out1$MV.pval,2,lower=F)+qchisq(out1$IV.pval,2,lower=F),4,lower=F)), signif(out1$JV.pval)) expect_equal(signif(out1$JV.pval,3), signif(out1$JF.pval,3)) expect_equal(signif(out1$JV.pval,3), signif(out1$JD.pval,3)) expect_equal(signif(pchisq(qchisq(obj1.gds$P_Value_Marginal,1,lower=F)+qchisq(obj1.gds$P_Value_Interaction,1,lower=F),2,lower=F)), signif(obj1.gds$P_Value_Joint)) obj1.outfile.bgen <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen) obj1.bgen <- read.table(obj1.outfile.bgen, header = T, as.is = T) expect_equal(obj1.gds, obj1.bgen[,-2],tolerance=.0001) unlink(c(obj1.outfile.gds, obj1.outfile.bgen)) ### adding one interaction covariate out1.cov <- MAGEE(null.obj=obj1, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj1.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, outfile=obj1.outfile.gds.cov) obj1.gds.cov <- read.table(obj1.outfile.gds.cov, header = T, as.is = T) expect_equal(signif(range(obj1.gds.cov$P_Value_Marginal)), signif(c(0.007431214, 0.993737081))) expect_equal(signif(range(obj1.gds.cov$P_Value_Interaction)), signif(c(0.0172582, 0.9990436))) expect_equal(signif(range(obj1.gds.cov$P_Value_Joint)), signif(c(0.006160225, 0.994098424))) expect_equal(signif(out1.cov$MV.pval,4), signif(obj1.gds.cov$P_Value_Marginal,4)) expect_equal(signif(out1.cov$MF.pval), signif(obj1.gds.cov$P_Value_Marginal)) expect_equal(signif(out1.cov$IV.pval,4), signif(obj1.gds.cov$P_Value_Interaction,4)) expect_equal(signif(out1.cov$IF.pval), signif(obj1.gds.cov$P_Value_Interaction)) expect_equal(signif(pchisq(qchisq(pchisq(qchisq(obj1.gds.cov$P_Value_Joint, df=2, lower=F)-qchisq(obj1.gds.cov$P_Value_Interaction, df=1, lower=F), df=1, lower=F),2,lower=F)+qchisq(out1.cov$IV.pval,2,lower=F),4,lower=F),4), signif(out1.cov$JV.pval,4)) expect_equal(signif(out1.cov$JV.pval,3), signif(out1.cov$JF.pval,3)) expect_equal(signif(out1.cov$JV.pval,3), signif(out1.cov$JD.pval,3)) obj1.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen.cov) obj1.bgen.cov <- read.table(obj1.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj1.gds.cov, obj1.bgen.cov[,-2],tolerance=.0001) unlink(c(obj1.outfile.gds.cov, obj1.outfile.bgen.cov)) obj1.glmm.score.outfile.gds <- tempfile() glmm.score(obj=obj1, infile=gdsfile, outfile=obj1.glmm.score.outfile.gds) obj1.glmm.score.gds <- read.table(obj1.glmm.score.outfile.gds,header=T,as.is=T) expect_equal(obj1.gds$P_Value_Marginal, obj1.glmm.score.gds$PVAL) obj1.glmm.score.outfile.bgen <- tempfile() glmm.score(obj=obj1, infile=bgenfile, BGEN.samplefile = samplefile, outfile=obj1.glmm.score.outfile.bgen) obj1.glmm.score.bgen <- read.table(obj1.glmm.score.outfile.bgen,header=T,as.is=T) expect_equal(signif(obj1.gds$P_Value_Marginal), signif(obj1.glmm.score.bgen$PVAL)) unlink(c(obj1.glmm.score.outfile.gds, obj1.glmm.score.outfile.bgen)) ### single thread without kins obj2 <- glmmkin(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity")) out2 <- MAGEE(null.obj=obj2, interaction="sex", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj2.outfile.gds <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=gdsfile, outfile=obj2.outfile.gds) obj2.gds <- read.table(obj2.outfile.gds, header = T, as.is = T) expect_equal(signif(range(obj2.gds$P_Value_Marginal)), signif(c( 0.0002174546, 0.9996967239))) expect_equal(signif(range(obj2.gds$P_Value_Interaction)), signif(c(0.02536793, 0.99887262))) expect_equal(signif(range(obj2.gds$P_Value_Joint)), signif(c( 0.0006177245, 0.9979718241))) expect_equal(signif(out2$MV.pval,3), signif(obj2.gds$P_Value_Marginal,3)) expect_equal(signif(out2$MF.pval), signif(obj2.gds$P_Value_Marginal)) expect_equal(signif(out2$IV.pval,3), signif(obj2.gds$P_Value_Interaction,3)) expect_equal(signif(out2$IF.pval), signif(obj2.gds$P_Value_Interaction)) expect_equal(signif(pchisq(qchisq(out2$MV.pval,2,lower=F)+qchisq(out2$IV.pval,2,lower=F),4,lower=F)), signif(out2$JV.pval)) expect_equal(signif(out2$JV.pval,4), signif(out2$JF.pval,4)) expect_equal(signif(out2$JV.pval,4), signif(out2$JD.pval,4)) expect_equal(signif(pchisq(qchisq(obj2.gds$P_Value_Marginal,1,lower=F)+qchisq(obj2.gds$P_Value_Interaction,1,lower=F),2,lower=F)), signif(obj2.gds$P_Value_Joint)) obj2.outfile.bgen <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen) obj2.bgen <- read.table(obj2.outfile.bgen, header = T, as.is = T) expect_equal(obj2.gds, obj2.bgen[,-2],tolerance=.0001) unlink(c(obj2.outfile.gds, obj2.outfile.bgen)) ### adding one interaction covariate out2.cov <- MAGEE(null.obj=obj2, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj2.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, outfile=obj2.outfile.gds.cov) obj2.gds.cov <- read.table(obj2.outfile.gds.cov, header = T, as.is = T) expect_equal(signif(range(obj2.gds.cov$P_Value_Marginal)), signif(c(0.0002174546, 0.9996967239))) expect_equal(signif(range(obj2.gds.cov$P_Value_Interaction)), signif(c(0.01803992, 0.98953447))) expect_equal(signif(range(obj2.gds.cov$P_Value_Joint)), signif(c(0.0004792136, 0.9856456114))) expect_equal(signif(out2.cov$MV.pval,3), signif(obj2.gds.cov$P_Value_Marginal,3)) expect_equal(signif(out2.cov$MF.pval), signif(obj2.gds.cov$P_Value_Marginal)) expect_equal(signif(out2.cov$IV.pval,4), signif(obj2.gds.cov$P_Value_Interaction,4)) expect_equal(signif(out2.cov$IF.pval), signif(obj2.gds.cov$P_Value_Interaction)) expect_equal(signif(pchisq(qchisq(pchisq(qchisq(obj2.gds.cov$P_Value_Joint, df=2, lower=F)-qchisq(obj2.gds.cov$P_Value_Interaction, df=1, lower=F), df=1, lower=F),2,lower=F)+qchisq(out2.cov$IV.pval,2,lower=F),4,lower=F),4), signif(out2.cov$JV.pval,4)) expect_equal(signif(out2.cov$JV.pval,3), signif(out2.cov$JF.pval,3)) expect_equal(signif(out2.cov$JV.pval,3), signif(out2.cov$JD.pval,3)) obj2.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen.cov) obj2.bgen.cov <- read.table(obj2.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj2.gds.cov, obj2.bgen.cov[,-2],tolerance=.0001) unlink(c(obj2.outfile.gds.cov, obj2.outfile.bgen.cov)) obj2.glmm.score.outfile.gds <- tempfile() glmm.score(obj=obj2, infile=gdsfile, outfile=obj2.glmm.score.outfile.gds) obj2.glmm.score.gds <- read.table(obj2.glmm.score.outfile.gds,header=T,as.is=T) expect_equal(obj2.gds$P_Value_Marginal, obj2.glmm.score.gds$PVAL) obj2.glmm.score.outfile.bgen <- tempfile() glmm.score(obj=obj2, infile=bgenfile, BGEN.samplefile = samplefile, outfile=obj2.glmm.score.outfile.bgen) obj2.glmm.score.bgen <- read.table(obj2.glmm.score.outfile.bgen,header=T,as.is=T) expect_equal(signif(obj2.gds$P_Value_Marginal), signif(obj2.glmm.score.bgen$PVAL)) unlink(c(obj2.glmm.score.outfile.gds, obj2.glmm.score.outfile.bgen)) ###single thread with kins and two interaction terms obj3 <- glmmkin(trait ~ age + sex + obese, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity")) out3 <- MAGEE(null.obj=obj3, interaction=c("sex","obese"), geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj3.outfile.gds <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj3.outfile.gds) obj3.gds <- read.table(obj3.outfile.gds, header = T, as.is = T) expect_equal(signif(range(obj3.gds$P_Value_Marginal)), signif(c( 0.006841261, 0.995945684))) expect_equal(signif(range(obj3.gds$P_Value_Interaction)), signif(c(0.03765037, 0.97863054))) expect_equal(signif(range(obj3.gds$P_Value_Joint)), signif(c( 0.009402029, 0.990308837))) expect_equal(signif(out3$MV.pval,4), signif(obj3.gds$P_Value_Marginal,4)) expect_equal(signif(out3$MF.pval), signif(obj3.gds$P_Value_Marginal)) expect_equal(signif(pchisq(qchisq(out3$MV.pval, df =2 ,lower=F)+qchisq(out3$IV.pval, df=2, lower=F), df=4, lower=F),4 ),signif(out3$JV.pval,4 )) expect_equal(signif(pchisq(qchisq(out3$MF.pval, df =2 ,lower=F)+qchisq(out3$IF.pval, df=2, lower=F), df=4, lower=F),4 ),signif(out3$JD.pval,4 )) expect_equal(signif(pchisq(qchisq(obj3.gds$P_Value_Marginal[-83],1,lower=F)+qchisq(obj3.gds$P_Value_Interaction[-83],2,lower=F),3,lower=F)), signif(obj3.gds$P_Value_Joint[-83])) obj3.outfile.bgen <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen) obj3.bgen <- read.table(obj3.outfile.bgen, header = T, as.is = T) expect_equal(obj3.gds[-83,], obj3.bgen[-83,-2],tolerance=.0001) unlink(c(obj3.outfile.gds, obj3.outfile.bgen)) ### adding one interaction covariate out3.cov <- MAGEE(null.obj=obj3, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj3.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, outfile=obj3.outfile.gds.cov) obj3.gds.cov <- read.table(obj3.outfile.gds.cov, header = T, as.is = T) expect_equal(signif(range(obj3.gds.cov$P_Value_Marginal)), signif(c( 0.006841261, 0.995945684))) expect_equal(signif(range(obj3.gds.cov$P_Value_Interaction)), signif(c(0.02681612, 0.98444877))) expect_equal(signif(range(obj3.gds.cov$P_Value_Joint)), signif(c(0.009402014, 0.981364975))) expect_equal(signif(out3.cov$MV.pval,4), signif(obj3.gds.cov$P_Value_Marginal,4)) expect_equal(signif(out3.cov$MF.pval), signif(obj3.gds.cov$P_Value_Marginal)) expect_equal(signif(pchisq(qchisq(pchisq(qchisq(obj3.gds.cov$P_Value_Joint[c(-64,-83)], df=3, lower=F)-qchisq(obj3.gds.cov$P_Value_Interaction[c(-64,-83)], df=2, lower=F), df=1, lower=F), df =2 ,lower=F)+qchisq(out3.cov$IV.pval[c(-64,-83)], df=2, lower=F), df=4, lower=F),4 ),signif(out3.cov$JV.pval[c(-64,-83)],4 )) expect_equal(signif(pchisq(qchisq(pchisq(qchisq(obj3.gds.cov$P_Value_Joint[c(-64,-83)], df=3, lower=F)-qchisq(obj3.gds.cov$P_Value_Interaction[c(-64,-83)], df=2, lower=F), df=1, lower=F), df =2 ,lower=F)+qchisq(out3.cov$IF.pval[c(-64,-83)], df=2, lower=F), df=4, lower=F),4 ),signif(out3.cov$JD.pval[c(-64,-83)],4 )) obj3.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen.cov) obj3.bgen.cov <- read.table(obj3.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj3.gds.cov[c(-64,-83),], obj3.bgen.cov[c(-64,-83),-2],tolerance=.0001) unlink(c(obj3.outfile.gds.cov, obj3.outfile.bgen.cov)) obj3.glmm.score.outfile.gds <- tempfile() glmm.score(obj=obj3, infile=gdsfile, outfile=obj3.glmm.score.outfile.gds) obj3.glmm.score.gds <- read.table(obj3.glmm.score.outfile.gds,header=T,as.is=T) expect_equal(obj3.gds$P_Value_Marginal, obj3.glmm.score.gds$PVAL) obj3.glmm.score.outfile.bgen <- tempfile() glmm.score(obj=obj3, infile=bgenfile, BGEN.samplefile = samplefile, outfile=obj3.glmm.score.outfile.bgen) obj3.glmm.score.bgen <- read.table(obj3.glmm.score.outfile.bgen,header=T,as.is=T) expect_equal(signif(obj3.gds$P_Value_Marginal), signif(obj3.glmm.score.bgen$PVAL)) unlink(c(obj3.glmm.score.outfile.gds, obj3.glmm.score.outfile.bgen)) ###single thread without kins and with two interaction terms obj4 <- glmmkin(trait ~ age + sex + obese, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity")) out4 <- MAGEE(null.obj=obj4, interaction=c("sex","obese"), geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj4.outfile.gds <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj4.outfile.gds) obj4.gds <- read.table(obj4.outfile.gds, header = T, as.is = T) expect_equal(signif(range(obj4.gds$P_Value_Marginal)), signif(c(0.000172054, 0.998167132 ))) expect_equal(signif(range(obj4.gds$P_Value_Interaction)), signif(c(0.0428305, 0.9992470))) expect_equal(signif(range(obj4.gds$P_Value_Joint)), signif(c( 0.001257493, 0.999691629 ))) expect_equal(signif(out4$MV.pval,3), signif(obj4.gds$P_Value_Marginal,3)) expect_equal(signif(out4$MF.pval), signif(obj4.gds$P_Value_Marginal)) expect_equal(signif(pchisq(qchisq(out4$MV.pval, df =2 ,lower=F)+qchisq(out4$IV.pval, df=2, lower=F), df=4, lower=F),4 ),signif(out4$JV.pval,4 )) expect_equal(signif(pchisq(qchisq(out4$MF.pval, df =2 ,lower=F)+qchisq(out4$IF.pval, df=2, lower=F), df=4, lower=F),4 ),signif(out4$JD.pval,4 )) expect_equal(signif(pchisq(qchisq(obj4.gds$P_Value_Marginal[-83],1,lower=F)+qchisq(obj4.gds$P_Value_Interaction[-83],2,lower=F),3,lower=F)), signif(obj4.gds$P_Value_Joint[-83])) obj4.outfile.bgen <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen) obj4.bgen <- read.table(obj4.outfile.bgen, header = T, as.is = T) expect_equal(obj4.gds[-83,], obj4.bgen[-83,-2],tolerance=.0001) unlink(c(obj4.outfile.gds, obj4.outfile.bgen)) ### adding one interaction covariate out4.cov <- MAGEE(null.obj=obj4, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj4.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, outfile=obj4.outfile.gds.cov) obj4.gds.cov <- read.table(obj4.outfile.gds.cov, header = T, as.is = T) expect_equal(signif(range(obj4.gds.cov$P_Value_Marginal)), signif(c( 0.000172054, 0.998167132))) expect_equal(signif(range(obj4.gds.cov$P_Value_Interaction)), signif(c(0.03288822, 0.99213068))) expect_equal(signif(range(obj4.gds.cov$P_Value_Joint)), signif(c(0.0009759281, 0.9808774302))) expect_equal(signif(out4.cov$MV.pval,3), signif(obj4.gds.cov$P_Value_Marginal,3)) expect_equal(signif(out4.cov$MF.pval), signif(obj4.gds.cov$P_Value_Marginal)) expect_equal(signif(pchisq(qchisq(pchisq(qchisq(obj4.gds.cov$P_Value_Joint[c(-64,-83)], df=3, lower=F)-qchisq(obj4.gds.cov$P_Value_Interaction[c(-64,-83)], df=2, lower=F), df=1, lower=F), df =2 ,lower=F)+qchisq(out4.cov$IV.pval[c(-64,-83)], df=2, lower=F), df=4, lower=F),4 ),signif(out4.cov$JV.pval[c(-64,-83)],4 )) expect_equal(signif(pchisq(qchisq(pchisq(qchisq(obj4.gds.cov$P_Value_Joint[c(-64,-83)], df=3, lower=F)-qchisq(obj4.gds.cov$P_Value_Interaction[c(-64,-83)], df=2, lower=F), df=1, lower=F), df =2 ,lower=F)+qchisq(out4.cov$IF.pval[c(-64,-83)], df=2, lower=F), df=4, lower=F),4 ),signif(out4.cov$JD.pval[c(-64,-83)],4 )) obj4.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen.cov) obj4.bgen.cov <- read.table(obj4.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj4.gds.cov[c(-64,-83),], obj4.bgen.cov[c(-64,-83),-2],tolerance=.0001) unlink(c(obj4.outfile.gds.cov, obj4.outfile.bgen.cov)) obj4.glmm.score.outfile.gds <- tempfile() glmm.score(obj=obj4, infile=gdsfile, outfile=obj4.glmm.score.outfile.gds) obj4.glmm.score.gds <- read.table(obj4.glmm.score.outfile.gds,header=T,as.is=T) expect_equal(obj4.gds$P_Value_Marginal, obj4.glmm.score.gds$PVAL) obj4.glmm.score.outfile.bgen <- tempfile() glmm.score(obj=obj4, infile=bgenfile, BGEN.samplefile = samplefile, outfile=obj4.glmm.score.outfile.bgen) obj4.glmm.score.bgen <- read.table(obj4.glmm.score.outfile.bgen,header=T,as.is=T) expect_equal(signif(obj4.gds$P_Value_Marginal), signif(obj4.glmm.score.bgen$PVAL)) unlink(c(obj4.glmm.score.outfile.gds, obj4.glmm.score.outfile.bgen)) ### multi-thread ### obj1 skip_on_cran() obj1.outfile.gds.tmp <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=gdsfile, outfile=obj1.outfile.gds.tmp, ncores=2) obj1.gds.tmp <- read.table(obj1.outfile.gds.tmp, header = T, as.is = T) expect_equal(obj1.gds, obj1.gds.tmp) obj1.outfile.bgen.tmp <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen.tmp, ncores=2) obj1.bgen.tmp <- read.table(obj1.outfile.bgen.tmp, header = T, as.is = T) expect_equal(obj1.bgen, obj1.bgen.tmp) unlink(c(obj1.outfile.gds.tmp, obj1.outfile.bgen.tmp)) ### adding one interaction covariate for obj1 obj1.outfile.gds.tmp.cov <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", interaction.covariates ="age",geno.file=gdsfile, outfile=obj1.outfile.gds.tmp.cov, ncores=2) obj1.gds.tmp.cov <- read.table(obj1.outfile.gds.tmp.cov, header = T, as.is = T) expect_equal(obj1.gds.cov, obj1.gds.tmp.cov) obj1.outfile.bgen.tmp.cov <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen.tmp.cov, ncores=2) obj1.bgen.tmp.cov <- read.table(obj1.outfile.bgen.tmp.cov, header = T, as.is = T) expect_equal(obj1.bgen.cov, obj1.bgen.tmp.cov) unlink(c(obj1.outfile.gds.tmp.cov, obj1.outfile.bgen.tmp.cov)) ### obj2 obj2.outfile.gds.tmp <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=gdsfile, outfile=obj2.outfile.gds.tmp, ncores=2) obj2.gds.tmp <- read.table(obj2.outfile.gds.tmp, header = T, as.is = T) expect_equal(obj2.gds, obj2.gds.tmp) obj2.outfile.bgen.tmp <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen.tmp, ncores=2) obj2.bgen.tmp <- read.table(obj2.outfile.bgen.tmp, header = T, as.is = T) expect_equal(obj2.bgen, obj2.bgen.tmp) unlink(c(obj2.outfile.gds.tmp, obj2.outfile.bgen.tmp)) ### adding one interaction covariate for obj2 obj2.outfile.gds.tmp.cov <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", interaction.covariates ="age",geno.file=gdsfile, outfile=obj2.outfile.gds.tmp.cov, ncores=2) obj2.gds.tmp.cov <- read.table(obj2.outfile.gds.tmp.cov, header = T, as.is = T) expect_equal(obj2.gds.cov, obj2.gds.tmp.cov) obj2.outfile.bgen.tmp.cov <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen.tmp.cov, ncores=2) obj2.bgen.tmp.cov <- read.table(obj2.outfile.bgen.tmp.cov, header = T, as.is = T) expect_equal(obj2.bgen.cov, obj2.bgen.tmp.cov) unlink(c(obj2.outfile.gds.tmp.cov, obj2.outfile.bgen.tmp.cov)) ### multi-thread with two interaction terms ### obj3 obj3.outfile.gds.tmp <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj3.outfile.gds.tmp, ncores=2) obj3.gds.tmp <- read.table(obj3.outfile.gds.tmp, header = T, as.is = T) expect_equal(obj3.gds, obj3.gds.tmp) obj3.outfile.bgen.tmp <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen.tmp, ncores=2) obj3.bgen.tmp <- read.table(obj3.outfile.bgen.tmp, header = T, as.is = T) expect_equal(obj3.bgen[-83,], obj3.bgen.tmp[-83,]) unlink(c(obj3.outfile.gds.tmp, obj3.outfile.bgen.tmp)) ### adding one interaction covariate for obj3 obj3.outfile.gds.tmp.cov <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, outfile=obj3.outfile.gds.tmp.cov, ncores=2) obj3.gds.tmp.cov <- read.table(obj3.outfile.gds.tmp.cov, header = T, as.is = T) expect_equal(obj3.gds.cov, obj3.gds.tmp.cov) obj3.outfile.bgen.tmp.cov <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"),interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen.tmp.cov, ncores=2) obj3.bgen.tmp.cov <- read.table(obj3.outfile.bgen.tmp.cov, header = T, as.is = T) expect_equal(obj3.bgen.cov[c(-64,-83),], obj3.bgen.tmp.cov[c(-64,-83),]) unlink(c(obj3.outfile.gds.tmp.cov, obj3.outfile.bgen.tmp.cov)) ### obj4 obj4.outfile.gds.tmp <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj4.outfile.gds.tmp, ncores=2) obj4.gds.tmp <- read.table(obj4.outfile.gds.tmp, header = T, as.is = T) expect_equal(obj4.gds, obj4.gds.tmp) obj4.outfile.bgen.tmp <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen.tmp, ncores=2) obj4.bgen.tmp <- read.table(obj4.outfile.bgen.tmp, header = T, as.is = T) expect_equal(obj4.bgen[c(-64,-83),], obj4.bgen.tmp[c(-64,-83),]) unlink(c(obj4.outfile.gds.tmp, obj4.outfile.bgen.tmp)) ### adding one interaction covariate for obj4 obj4.outfile.gds.tmp.cov <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, outfile=obj4.outfile.gds.tmp.cov, ncores=2) obj4.gds.tmp.cov <- read.table(obj4.outfile.gds.tmp.cov, header = T, as.is = T) expect_equal(obj4.gds.cov, obj4.gds.tmp.cov) obj4.outfile.bgen.tmp.cov <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"),interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen.tmp.cov, ncores=2) obj4.bgen.tmp.cov <- read.table(obj4.outfile.bgen.tmp.cov, header = T, as.is = T) expect_equal(obj4.bgen.cov[c(-64,-83),], obj4.bgen.tmp.cov[c(-64,-83),]) unlink(c(obj4.outfile.gds.tmp.cov, obj4.outfile.bgen.tmp.cov)) ### re-order pheno id ### obj1 #skip_on_cran() idx <- sample(nrow(pheno)) pheno <- pheno[idx, ] obj1 <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity")) obj1.outfile.gds <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=gdsfile, outfile=obj1.outfile.gds) tmpout <- read.table(obj1.outfile.gds, header = T, as.is = T) expect_equal(obj1.gds, tmpout) obj1.outfile.bgen <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen) tmpout <- read.table(obj1.outfile.bgen, header = T, as.is = T) expect_equal(obj1.bgen, tmpout) unlink(c(obj1.outfile.gds, obj1.outfile.bgen)) ### adding one interaction covariate for obj1 obj1.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, outfile=obj1.outfile.gds.cov) tmpout.cov <- read.table(obj1.outfile.gds.cov, header = T, as.is = T) expect_equal(obj1.gds.cov, tmpout.cov) obj1.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj1, interaction="sex",interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen.cov) tmpout.cov <- read.table(obj1.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj1.bgen.cov, tmpout.cov) unlink(c(obj1.outfile.gds.cov, obj1.outfile.bgen.cov)) ### obj2 obj2 <- glmmkin(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity")) obj2.outfile.gds <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=gdsfile, outfile=obj2.outfile.gds) tmpout <- read.table(obj2.outfile.gds, header = T, as.is = T) expect_equal(obj2.gds, tmpout) obj2.outfile.bgen <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen) tmpout <- read.table(obj2.outfile.bgen, header = T, as.is = T) expect_equal(obj2.bgen, tmpout) unlink(c(obj2.outfile.gds, obj2.outfile.bgen)) ### adding one interaction covariate for obj2 obj2.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, outfile=obj2.outfile.gds.cov) tmpout.cov <- read.table(obj2.outfile.gds.cov, header = T, as.is = T) expect_equal(obj2.gds.cov, tmpout.cov) obj2.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj2, interaction="sex",interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen.cov) tmpout.cov <- read.table(obj2.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj2.bgen.cov, tmpout.cov) unlink(c(obj2.outfile.gds.cov, obj2.outfile.bgen.cov)) ### re-order pheno id with two interaction terms ### obj3 obj3 <- glmmkin(trait ~ age + sex + obese, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity")) obj3.outfile.gds <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj3.outfile.gds) tmpout <- read.table(obj3.outfile.gds, header = T, as.is = T) expect_equal(obj3.gds, tmpout) obj3.outfile.bgen <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen,nperbatch=2) tmpout <- read.table(obj3.outfile.bgen, header = T, as.is = T) expect_equal(obj3.bgen[-83,], tmpout[-83,]) unlink(c(obj3.outfile.gds, obj3.outfile.bgen)) ### adding one interaction covariate for obj3 obj3.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, outfile=obj3.outfile.gds.cov) tmpout.cov <- read.table(obj3.outfile.gds.cov, header = T, as.is = T) expect_equal(obj3.gds.cov, tmpout.cov) obj3.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"),interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen.cov) tmpout.cov <- read.table(obj3.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj3.bgen.cov[c(-64,-83),], tmpout.cov[c(-64,-83),]) unlink(c(obj3.outfile.gds.cov, obj3.outfile.bgen.cov)) ### obj4 obj4 <- glmmkin(trait ~ age + sex + obese, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity")) obj4.outfile.gds <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj4.outfile.gds) tmpout <- read.table(obj4.outfile.gds, header = T, as.is = T) expect_equal(obj4.gds, tmpout) obj4.outfile.bgen <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen) tmpout <- read.table(obj4.outfile.bgen, header = T, as.is = T) expect_equal(obj4.bgen[-83,], tmpout[-83,]) unlink(c(obj4.outfile.gds, obj4.outfile.bgen)) ### adding one interaction covariate for obj4 obj4.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, outfile=obj4.outfile.gds.cov) tmpout.cov <- read.table(obj4.outfile.gds.cov, header = T, as.is = T) expect_equal(obj4.gds.cov, tmpout.cov) obj4.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"),interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen.cov) tmpout.cov <- read.table(obj4.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj4.bgen.cov[c(-64,-83),], tmpout.cov[c(-64,-83),]) unlink(c(obj4.outfile.gds.cov, obj4.outfile.bgen.cov)) ### re-order kins id ### obj1 idx <- sample(nrow(kins)) kins <- kins[idx, idx] obj1 <- glmmkin(trait ~ age + sex, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity")) obj1.outfile.gds <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=gdsfile, outfile=obj1.outfile.gds) tmpout <- read.table(obj1.outfile.gds, header = T, as.is = T) expect_equal(obj1.gds, tmpout) obj1.outfile.bgen <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen) tmpout <- read.table(obj1.outfile.bgen, header = T, as.is = T) expect_equal(obj1.bgen, tmpout) unlink(c(obj1.outfile.gds, obj1.outfile.bgen)) ### adding one interaction covariate for obj1 obj1.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, outfile=obj1.outfile.gds.cov) tmpout.cov <- read.table(obj1.outfile.gds.cov, header = T, as.is = T) expect_equal(obj1.gds.cov, tmpout.cov) obj1.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj1, interaction="sex",interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen.cov) tmpout.cov <- read.table(obj1.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj1.bgen.cov, tmpout.cov) unlink(c(obj1.outfile.gds.cov, obj1.outfile.bgen.cov)) ### obj2 obj2 <- glmmkin(trait ~ age + sex, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity")) obj2.outfile.gds <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=gdsfile, outfile=obj2.outfile.gds) tmpout <- read.table(obj2.outfile.gds, header = T, as.is = T) expect_equal(obj2.gds, tmpout) obj2.outfile.bgen <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen) tmpout <- read.table(obj2.outfile.bgen, header = T, as.is = T) expect_equal(obj2.bgen, tmpout) unlink(c(obj2.outfile.gds, obj2.outfile.bgen)) ### adding one interaction covariate for obj2 obj2.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, outfile=obj2.outfile.gds.cov) tmpout.cov <- read.table(obj2.outfile.gds.cov, header = T, as.is = T) expect_equal(obj2.gds.cov, tmpout.cov) obj2.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj2, interaction="sex",interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen.cov) tmpout.cov <- read.table(obj2.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj2.bgen.cov, tmpout.cov) unlink(c(obj2.outfile.gds.cov, obj2.outfile.bgen.cov)) ### re-order kins id with two interaction terms ### obj3 obj3 <- glmmkin(trait ~ age + sex + obese, data = pheno, kins = kins, id = "id", family = gaussian(link = "identity")) obj3.outfile.gds <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj3.outfile.gds) tmpout <- read.table(obj3.outfile.gds, header = T, as.is = T) expect_equal(obj3.gds, tmpout) obj3.outfile.bgen <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen,nperbatch=2) tmpout <- read.table(obj3.outfile.bgen, header = T, as.is = T) expect_equal(obj3.bgen[-83,], tmpout[-83,]) unlink(c(obj3.outfile.gds, obj3.outfile.bgen)) ### adding one interaction covariate for obj3 obj3.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, outfile=obj3.outfile.gds.cov) tmpout.cov <- read.table(obj3.outfile.gds.cov, header = T, as.is = T) expect_equal(obj3.gds.cov, tmpout.cov) obj3.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"),interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen.cov) tmpout.cov <- read.table(obj3.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj3.bgen.cov[c(-64,-83),], tmpout.cov[c(-64,-83),]) unlink(c(obj3.outfile.gds.cov, obj3.outfile.bgen.cov)) ### obj4 obj4 <- glmmkin(trait ~ age + sex + obese, data = pheno, kins = NULL, id = "id", family = gaussian(link = "identity")) obj4.outfile.gds <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj4.outfile.gds) tmpout <- read.table(obj4.outfile.gds, header = T, as.is = T) expect_equal(obj4.gds, tmpout) obj4.outfile.bgen <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen) tmpout <- read.table(obj4.outfile.bgen, header = T, as.is = T) expect_equal(obj4.bgen[-83,], tmpout[-83,]) unlink(c(obj4.outfile.gds, obj4.outfile.bgen)) ### adding one interaction covariate for obj4 obj4.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, outfile=obj4.outfile.gds.cov) tmpout.cov <- read.table(obj4.outfile.gds.cov, header = T, as.is = T) expect_equal(obj4.gds.cov, tmpout.cov) obj4.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"),interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen.cov) tmpout.cov <- read.table(obj4.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj4.bgen.cov[c(-64,-83),], tmpout.cov[c(-64,-83),]) unlink(c(obj4.outfile.gds.cov, obj4.outfile.bgen.cov)) }) test_that("cross-sectional binomial", { skip_on_cran() skip_if_not_installed("SeqArray") skip_if_not_installed("SeqVarTools") gdsfile <- system.file("extdata", "geno.gds", package = "MAGEE") bgenfile <- system.file("extdata", "geno.bgen", package = "MAGEE") samplefile <- system.file("extdata", "geno.sample", package = "MAGEE") group.file <- system.file("extdata", "SetID.singlevariant.txt", package = "MAGEE") data(example) set.seed(123) pheno <- rbind(example$pheno, example$pheno[1:100, ]) pheno$id <- 1:500 pheno$disease[sample(1:500,20)] <- NA pheno$age[sample(1:500,20)] <- NA pheno$sex[sample(1:500,20)] <- NA pheno <- pheno[sample(1:500,450), ] pheno <- pheno[pheno$id <= 400, ] pheno$obese <- rbinom(nrow(pheno), 1, 0.3) kins <- example$GRM ### single thread with kins obj1 <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit")) out1 <- MAGEE(null.obj=obj1, interaction="sex", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj1.outfile.gds <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=gdsfile, outfile=obj1.outfile.gds) obj1.gds <- read.table(obj1.outfile.gds, header = T, as.is = T) expect_equal(signif(range(obj1.gds$P_Value_Marginal)), signif(c(0.008599802, 0.993475648))) expect_equal(signif(range(obj1.gds$P_Value_Interaction)), signif(c(0.01813711, 0.99990766))) expect_equal(signif(range(obj1.gds$P_Value_Joint)), signif(c(0.007484391, 0.998834862))) expect_equal(signif(out1$MV.pval,5), signif(obj1.gds$P_Value_Marginal,5)) expect_equal(signif(out1$MF.pval), signif(obj1.gds$P_Value_Marginal)) expect_equal(signif(out1$IV.pval,5), signif(obj1.gds$P_Value_Interaction,5)) expect_equal(signif(out1$IF.pval), signif(obj1.gds$P_Value_Interaction)) expect_equal(signif(pchisq(qchisq(out1$MV.pval,2,lower=F)+qchisq(out1$IV.pval,2,lower=F),4,lower=F)), signif(out1$JV.pval)) expect_equal(signif(out1$JV.pval,4), signif(out1$JF.pval,4)) expect_equal(signif(out1$JV.pval,4), signif(out1$JD.pval,4)) expect_equal(signif(pchisq(qchisq(obj1.gds$P_Value_Marginal,1,lower=F)+qchisq(obj1.gds$P_Value_Interaction,1,lower=F),2,lower=F)), signif(obj1.gds$P_Value_Joint)) obj1.outfile.bgen <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen) obj1.bgen <- read.table(obj1.outfile.bgen, header = T, as.is = T) expect_equal(obj1.gds, obj1.bgen[,-2],tolerance=.0001) unlink(c(obj1.outfile.gds, obj1.outfile.bgen)) ### adding one interaction covariate out1.cov <- MAGEE(null.obj=obj1, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj1.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, outfile=obj1.outfile.gds.cov) obj1.gds.cov <- read.table(obj1.outfile.gds.cov, header = T, as.is = T) expect_equal(signif(range(obj1.gds.cov$P_Value_Marginal)), signif(c(0.008599802, 0.993475648))) expect_equal(signif(range(obj1.gds.cov$P_Value_Interaction)), signif(c(0.02126704, 0.99966930))) expect_equal(signif(range(obj1.gds.cov$P_Value_Joint)), signif(c(0.008000314, 0.999792779))) expect_equal(signif(out1.cov$MV.pval,4), signif(obj1.gds.cov$P_Value_Marginal,4)) expect_equal(signif(out1.cov$MF.pval), signif(obj1.gds.cov$P_Value_Marginal)) expect_equal(signif(out1.cov$IV.pval,4), signif(obj1.gds.cov$P_Value_Interaction,4)) expect_equal(signif(out1.cov$IF.pval), signif(obj1.gds.cov$P_Value_Interaction)) expect_equal(signif(pchisq(qchisq(pchisq(qchisq(obj1.gds.cov$P_Value_Joint, df=2, lower=F)-qchisq(obj1.gds.cov$P_Value_Interaction, df=1, lower=F), df=1, lower=F),2,lower=F)+qchisq(out1.cov$IV.pval,2,lower=F),4,lower=F),4), signif(out1.cov$JV.pval,4)) expect_equal(signif(out1.cov$JV.pval,3), signif(out1.cov$JF.pval,3)) expect_equal(signif(out1.cov$JV.pval,3), signif(out1.cov$JD.pval,3)) obj1.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen.cov) obj1.bgen.cov <- read.table(obj1.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj1.gds.cov, obj1.bgen.cov[,-2],tolerance=.0001) unlink(c(obj1.outfile.gds.cov, obj1.outfile.bgen.cov)) obj1.glmm.score.outfile.gds <- tempfile() glmm.score(obj=obj1, infile=gdsfile, outfile=obj1.glmm.score.outfile.gds) obj1.glmm.score.gds <- read.table(obj1.glmm.score.outfile.gds,header=T,as.is=T) expect_equal(obj1.gds$P_Value_Marginal, obj1.glmm.score.gds$PVAL) obj1.glmm.score.outfile.bgen <- tempfile() glmm.score(obj=obj1, infile=bgenfile, BGEN.samplefile = samplefile, outfile=obj1.glmm.score.outfile.bgen) obj1.glmm.score.bgen <- read.table(obj1.glmm.score.outfile.bgen,header=T,as.is=T) expect_equal(signif(obj1.gds$P_Value_Marginal), signif(obj1.glmm.score.bgen$PVAL)) unlink(c(obj1.glmm.score.outfile.gds, obj1.glmm.score.outfile.bgen)) ### single thread without kins obj2 <- glmmkin(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit")) out2 <- MAGEE(null.obj=obj2, interaction="sex", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj2.outfile.gds <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=gdsfile, outfile=obj2.outfile.gds) obj2.gds <- read.table(obj2.outfile.gds, header = T, as.is = T) expect_equal(signif(range(obj2.gds$P_Value_Marginal)), signif(c(0.003470366, 0.973400784))) expect_equal(signif(range(obj2.gds$P_Value_Interaction)), signif(c(0.01136652, 0.98224264))) expect_equal(signif(range(obj2.gds$P_Value_Joint)), signif(c(0.003101226, 0.993232042))) expect_equal(signif(out2$MV.pval,4), signif(obj2.gds$P_Value_Marginal,4)) expect_equal(signif(out2$MF.pval,4), signif(obj2.gds$P_Value_Marginal,4)) expect_equal(signif(out2$IV.pval,3), signif(obj2.gds$P_Value_Interaction,3)) expect_equal(signif(out2$IF.pval,3), signif(obj2.gds$P_Value_Interaction,3)) expect_equal(signif(pchisq(qchisq(out2$MV.pval,2,lower=F)+qchisq(out2$IV.pval,2,lower=F),4,lower=F)), signif(out2$JV.pval)) expect_equal(signif(out2$JV.pval,3), signif(out2$JF.pval,3)) expect_equal(signif(out2$JV.pval,3), signif(out2$JD.pval,3)) expect_equal(signif(pchisq(qchisq(obj2.gds$P_Value_Marginal,1,lower=F)+qchisq(obj2.gds$P_Value_Interaction,1,lower=F),2,lower=F)), signif(obj2.gds$P_Value_Joint)) obj2.outfile.bgen <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen) obj2.bgen <- read.table(obj2.outfile.bgen, header = T, as.is = T) expect_equal(obj2.gds, obj2.bgen[,-2],tolerance=.0001) unlink(c(obj2.outfile.gds, obj2.outfile.bgen)) ### adding one interaction covariate out2.cov <- MAGEE(null.obj=obj2, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj2.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, outfile=obj2.outfile.gds.cov) obj2.gds.cov <- read.table(obj2.outfile.gds.cov, header = T, as.is = T) expect_equal(signif(range(obj2.gds.cov$P_Value_Marginal)), signif(c(0.003470366, 0.973400784))) expect_equal(signif(range(obj2.gds.cov$P_Value_Interaction)), signif(c(0.01338328, 0.98135298))) expect_equal(signif(range(obj2.gds.cov$P_Value_Joint)), signif(c(0.003427858, 0.998889784))) expect_equal(signif(out2.cov$MV.pval,3), signif(obj2.gds.cov$P_Value_Marginal,3)) expect_equal(signif(out2.cov$MF.pval,4), signif(obj2.gds.cov$P_Value_Marginal,4)) expect_equal(signif(out2.cov$IV.pval,4), signif(obj2.gds.cov$P_Value_Interaction,4)) expect_equal(signif(out2.cov$IF.pval,4), signif(obj2.gds.cov$P_Value_Interaction,4)) expect_equal(signif(pchisq(qchisq(pchisq(qchisq(obj2.gds.cov$P_Value_Joint, df=2, lower=F)-qchisq(obj2.gds.cov$P_Value_Interaction, df=1, lower=F), df=1, lower=F),2,lower=F)+qchisq(out2.cov$IV.pval,2,lower=F),4,lower=F),3), signif(out2.cov$JV.pval,3)) expect_equal(signif(out2.cov$JV.pval,3), signif(out2.cov$JF.pval,3)) expect_equal(signif(out2.cov$JV.pval,3), signif(out2.cov$JD.pval,3)) obj2.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen.cov) obj2.bgen.cov <- read.table(obj2.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj2.gds.cov, obj2.bgen.cov[,-2],tolerance=.0001) unlink(c(obj2.outfile.gds.cov, obj2.outfile.bgen.cov)) obj2.glmm.score.outfile.gds <- tempfile() glmm.score(obj=obj2, infile=gdsfile, outfile=obj2.glmm.score.outfile.gds) obj2.glmm.score.gds <- read.table(obj2.glmm.score.outfile.gds,header=T,as.is=T) expect_equal(obj2.gds$P_Value_Marginal, obj2.glmm.score.gds$PVAL) obj2.glmm.score.outfile.bgen <- tempfile() glmm.score(obj=obj2, infile=bgenfile, BGEN.samplefile = samplefile, outfile=obj2.glmm.score.outfile.bgen) obj2.glmm.score.bgen <- read.table(obj2.glmm.score.outfile.bgen,header=T,as.is=T) expect_equal(signif(obj2.gds$P_Value_Marginal), signif(obj2.glmm.score.bgen$PVAL)) unlink(c(obj2.glmm.score.outfile.gds, obj2.glmm.score.outfile.bgen)) ### single thread with kins and two interaction terms obj3 <- glmmkin(disease ~ age + sex + obese, data = pheno, kins = kins, id = "id", family = binomial(link = "logit")) out3 <- MAGEE(null.obj=obj3, interaction=c("sex","obese"), geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj3.outfile.gds <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj3.outfile.gds) obj3.gds <- read.table(obj3.outfile.gds, header = T, as.is = T) expect_equal(signif(range(obj3.gds$P_Value_Marginal)), signif(c(0.009277771, 0.983908371))) expect_equal(signif(range(obj3.gds$P_Value_Interaction)), signif(c(0.0287344, 0.9999111))) expect_equal(signif(range(obj3.gds$P_Value_Joint)), signif(c(0.02165704, 0.95845023))) expect_equal(signif(out3$MV.pval,4), signif(obj3.gds$P_Value_Marginal,4)) expect_equal(signif(out3$MF.pval), signif(obj3.gds$P_Value_Marginal)) expect_equal(signif(pchisq(qchisq(out3$MV.pval, df =2 ,lower=F)+qchisq(out3$IV.pval, df=2, lower=F), df=4, lower=F),4 ),signif(out3$JV.pval,4 )) expect_equal(signif(pchisq(qchisq(out3$MF.pval, df =2 ,lower=F)+qchisq(out3$IF.pval, df=2, lower=F), df=4, lower=F),4 ),signif(out3$JD.pval,4 )) expect_equal(signif(pchisq(qchisq(obj3.gds$P_Value_Marginal[-83],1,lower=F)+qchisq(obj3.gds$P_Value_Interaction[-83],2,lower=F),3,lower=F)), signif(obj3.gds$P_Value_Joint[-83])) obj3.outfile.bgen <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen) obj3.bgen <- read.table(obj3.outfile.bgen, header = T, as.is = T) expect_equal(obj3.gds[-83,], obj3.bgen[-83,-2],tolerance=.0001) unlink(c(obj3.outfile.gds, obj3.outfile.bgen)) ### adding one interaction covariate out3.cov <- MAGEE(null.obj=obj3, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj3.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, outfile=obj3.outfile.gds.cov) obj3.gds.cov <- read.table(obj3.outfile.gds.cov, header = T, as.is = T) expect_equal(signif(range(obj3.gds.cov$P_Value_Marginal)), signif(c(0.009277771, 0.983908371))) expect_equal(signif(range(obj3.gds.cov$P_Value_Interaction)), signif(c(0.03262695, 0.99999310))) expect_equal(signif(range(obj3.gds.cov$P_Value_Joint),7), signif(c(0.02305297, 0.96998850),7)) expect_equal(signif(out3.cov$MV.pval,4), signif(obj3.gds.cov$P_Value_Marginal,4)) expect_equal(signif(out3.cov$MF.pval), signif(obj3.gds.cov$P_Value_Marginal)) expect_equal(signif(pchisq(qchisq(pchisq(qchisq(obj3.gds.cov$P_Value_Joint[c(-64,-83)], df=3, lower=F)-qchisq(obj3.gds.cov$P_Value_Interaction[c(-64,-83)], df=2, lower=F), df=1, lower=F), df =2 ,lower=F)+qchisq(out3.cov$IV.pval[c(-64,-83)], df=2, lower=F), df=4, lower=F),4 ),signif(out3.cov$JV.pval[c(-64,-83)],4 )) expect_equal(signif(pchisq(qchisq(pchisq(qchisq(obj3.gds.cov$P_Value_Joint[c(-64,-83)], df=3, lower=F)-qchisq(obj3.gds.cov$P_Value_Interaction[c(-64,-83)], df=2, lower=F), df=1, lower=F), df =2 ,lower=F)+qchisq(out3.cov$IF.pval[c(-64,-83)], df=2, lower=F), df=4, lower=F),4 ),signif(out3.cov$JD.pval[c(-64,-83)],4 )) obj3.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen.cov) obj3.bgen.cov <- read.table(obj3.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj3.gds.cov[c(-64,-83),], obj3.bgen.cov[c(-64,-83),-2],tolerance=.0001) unlink(c(obj3.outfile.gds.cov, obj3.outfile.bgen.cov)) obj3.glmm.score.outfile.gds <- tempfile() glmm.score(obj=obj3, infile=gdsfile, outfile=obj3.glmm.score.outfile.gds) obj3.glmm.score.gds <- read.table(obj3.glmm.score.outfile.gds,header=T,as.is=T) expect_equal(obj3.gds$P_Value_Marginal, obj3.glmm.score.gds$PVAL) obj3.glmm.score.outfile.bgen <- tempfile() glmm.score(obj=obj3, infile=bgenfile, BGEN.samplefile = samplefile, outfile=obj3.glmm.score.outfile.bgen) obj3.glmm.score.bgen <- read.table(obj3.glmm.score.outfile.bgen,header=T,as.is=T) expect_equal(signif(obj3.gds$P_Value_Marginal), signif(obj3.glmm.score.bgen$PVAL)) unlink(c(obj3.glmm.score.outfile.gds, obj3.glmm.score.outfile.bgen)) ### single thread without kins and with two interaction terms obj4 <- glmmkin(disease ~ age + sex + obese, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit")) out4 <- MAGEE(null.obj=obj4, interaction=c("sex","obese"), geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj4.outfile.gds <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj4.outfile.gds) obj4.gds <- read.table(obj4.outfile.gds, header = T, as.is = T) expect_equal(signif(range(obj4.gds$P_Value_Marginal)), signif(c(0.003684986, 0.997820900))) expect_equal(signif(range(obj4.gds$P_Value_Interaction)), signif(c(0.01649446, 0.99966122))) expect_equal(signif(range(obj4.gds$P_Value_Joint)), signif(c(0.009495391, 0.966728864))) expect_equal(signif(out4$MV.pval,2), signif(obj4.gds$P_Value_Marginal,2)) expect_equal(signif(out4$MF.pval,2), signif(obj4.gds$P_Value_Marginal,2)) expect_equal(signif(pchisq(qchisq(out4$MV.pval, df =2 ,lower=F)+qchisq(out4$IV.pval, df=2, lower=F), df=4, lower=F),4 ),signif(out4$JV.pval,4 )) expect_equal(signif(pchisq(qchisq(out4$MF.pval, df =2 ,lower=F)+qchisq(out4$IF.pval, df=2, lower=F), df=4, lower=F),4 ),signif(out4$JD.pval,4 )) expect_equal(signif(pchisq(qchisq(obj4.gds$P_Value_Marginal[-83],1,lower=F)+qchisq(obj4.gds$P_Value_Interaction[-83],2,lower=F),3,lower=F)), signif(obj4.gds$P_Value_Joint[-83])) obj4.outfile.bgen <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen) obj4.bgen <- read.table(obj4.outfile.bgen, header = T, as.is = T) expect_equal(obj4.gds[-83,], obj4.bgen[-83,-2],tolerance=.0001) unlink(c(obj4.outfile.gds, obj4.outfile.bgen)) ### adding one interaction covariate out4.cov <- MAGEE(null.obj=obj4, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj4.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, outfile=obj4.outfile.gds.cov) obj4.gds.cov <- read.table(obj4.outfile.gds.cov, header = T, as.is = T) expect_equal(signif(range(obj4.gds.cov$P_Value_Marginal)), signif(c(0.003684986, 0.997820900))) expect_equal(signif(range(obj4.gds.cov$P_Value_Interaction)), signif(c(0.01882705, 0.99947529))) expect_equal(signif(range(obj4.gds.cov$P_Value_Joint),5), signif(c(0.0104292, 0.9645215),5)) expect_equal(signif(out4.cov$MV.pval,2), signif(obj4.gds.cov$P_Value_Marginal,2)) expect_equal(signif(out4.cov$MF.pval,2), signif(obj4.gds.cov$P_Value_Marginal,2)) expect_equal(signif(pchisq(qchisq(pchisq(qchisq(obj4.gds.cov$P_Value_Joint[c(-64,-83)], df=3, lower=F)-qchisq(obj4.gds.cov$P_Value_Interaction[c(-64,-83)], df=2, lower=F), df=1, lower=F), df =2 ,lower=F)+qchisq(out4.cov$IV.pval[c(-64,-83)], df=2, lower=F), df=4, lower=F),3),signif(out4.cov$JV.pval[c(-64,-83)],3)) expect_equal(signif(pchisq(qchisq(pchisq(qchisq(obj4.gds.cov$P_Value_Joint[c(-64,-83)], df=3, lower=F)-qchisq(obj4.gds.cov$P_Value_Interaction[c(-64,-83)], df=2, lower=F), df=1, lower=F), df =2 ,lower=F)+qchisq(out4.cov$IF.pval[c(-64,-83)], df=2, lower=F), df=4, lower=F),2),signif(out4.cov$JD.pval[c(-64,-83)],2)) obj4.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen.cov) obj4.bgen.cov <- read.table(obj4.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj4.gds.cov[c(-64,-83),], obj4.bgen.cov[c(-64,-83),-2],tolerance=.0001) unlink(c(obj4.outfile.gds.cov, obj4.outfile.bgen.cov)) obj4.glmm.score.outfile.gds <- tempfile() glmm.score(obj=obj4, infile=gdsfile, outfile=obj4.glmm.score.outfile.gds) obj4.glmm.score.gds <- read.table(obj4.glmm.score.outfile.gds,header=T,as.is=T) expect_equal(obj4.gds$P_Value_Marginal, obj4.glmm.score.gds$PVAL) obj4.glmm.score.outfile.bgen <- tempfile() glmm.score(obj=obj4, infile=bgenfile, BGEN.samplefile = samplefile, outfile=obj4.glmm.score.outfile.bgen) obj4.glmm.score.bgen <- read.table(obj4.glmm.score.outfile.bgen,header=T,as.is=T) expect_equal(signif(obj4.gds$P_Value_Marginal), signif(obj4.glmm.score.bgen$PVAL)) unlink(c(obj4.glmm.score.outfile.gds, obj4.glmm.score.outfile.bgen)) ### multi-thread ### obj1 obj1.outfile.gds.tmp <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=gdsfile, outfile=obj1.outfile.gds.tmp, ncores=2) obj1.gds.tmp <- read.table(obj1.outfile.gds.tmp, header = T, as.is = T) expect_equal(obj1.gds, obj1.gds.tmp) obj1.outfile.bgen.tmp <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen.tmp, ncores=2) obj1.bgen.tmp <- read.table(obj1.outfile.bgen.tmp, header = T, as.is = T) expect_equal(obj1.bgen, obj1.bgen.tmp) unlink(c(obj1.outfile.gds.tmp, obj1.outfile.bgen.tmp)) ### adding one interaction covariate for obj1 obj1.outfile.gds.tmp.cov <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", interaction.covariates ="age",geno.file=gdsfile, outfile=obj1.outfile.gds.tmp.cov, ncores=2) obj1.gds.tmp.cov <- read.table(obj1.outfile.gds.tmp.cov, header = T, as.is = T) expect_equal(obj1.gds.cov, obj1.gds.tmp.cov) obj1.outfile.bgen.tmp.cov <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen.tmp.cov, ncores=2) obj1.bgen.tmp.cov <- read.table(obj1.outfile.bgen.tmp.cov, header = T, as.is = T) expect_equal(obj1.bgen.cov, obj1.bgen.tmp.cov) unlink(c(obj1.outfile.gds.tmp.cov, obj1.outfile.bgen.tmp.cov)) ### obj2 obj2.outfile.gds.tmp <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=gdsfile, outfile=obj2.outfile.gds.tmp, ncores=2) obj2.gds.tmp <- read.table(obj2.outfile.gds.tmp, header = T, as.is = T) expect_equal(obj2.gds, obj2.gds.tmp) obj2.outfile.bgen.tmp <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen.tmp, ncores=2) obj2.bgen.tmp <- read.table(obj2.outfile.bgen.tmp, header = T, as.is = T) expect_equal(obj2.bgen, obj2.bgen.tmp) unlink(c(obj2.outfile.gds.tmp, obj2.outfile.bgen.tmp)) ### adding one interaction covariate for obj2 obj2.outfile.gds.tmp.cov <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", interaction.covariates ="age",geno.file=gdsfile, outfile=obj2.outfile.gds.tmp.cov, ncores=2) obj2.gds.tmp.cov <- read.table(obj2.outfile.gds.tmp.cov, header = T, as.is = T) expect_equal(obj2.gds.cov, obj2.gds.tmp.cov) obj2.outfile.bgen.tmp.cov <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen.tmp.cov, ncores=2) obj2.bgen.tmp.cov <- read.table(obj2.outfile.bgen.tmp.cov, header = T, as.is = T) expect_equal(obj2.bgen.cov, obj2.bgen.tmp.cov) unlink(c(obj2.outfile.gds.tmp.cov, obj2.outfile.bgen.tmp.cov)) ### multi-thread with two interaction terms ### obj3 obj3.outfile.gds.tmp <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj3.outfile.gds.tmp, ncores=2) obj3.gds.tmp <- read.table(obj3.outfile.gds.tmp, header = T, as.is = T) expect_equal(obj3.gds, obj3.gds.tmp) obj3.outfile.bgen.tmp <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen.tmp, ncores=2) obj3.bgen.tmp <- read.table(obj3.outfile.bgen.tmp, header = T, as.is = T) expect_equal(obj3.bgen[-83,], obj3.bgen.tmp[-83,]) unlink(c(obj3.outfile.gds.tmp, obj3.outfile.bgen.tmp)) ### adding one interaction covariate for obj3 obj3.outfile.gds.tmp.cov <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, outfile=obj3.outfile.gds.tmp.cov, ncores=2) obj3.gds.tmp.cov <- read.table(obj3.outfile.gds.tmp.cov, header = T, as.is = T) expect_equal(obj3.gds.cov, obj3.gds.tmp.cov) obj3.outfile.bgen.tmp.cov <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"),interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen.tmp.cov, ncores=2) obj3.bgen.tmp.cov <- read.table(obj3.outfile.bgen.tmp.cov, header = T, as.is = T) expect_equal(obj3.bgen.cov[c(-64,-83),], obj3.bgen.tmp.cov[c(-64,-83),]) unlink(c(obj3.outfile.gds.tmp.cov, obj3.outfile.bgen.tmp.cov)) ### obj4 obj4.outfile.gds.tmp <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj4.outfile.gds.tmp, ncores=2) obj4.gds.tmp <- read.table(obj4.outfile.gds.tmp, header = T, as.is = T) expect_equal(obj4.gds, obj4.gds.tmp) obj4.outfile.bgen.tmp <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen.tmp, ncores=2) obj4.bgen.tmp <- read.table(obj4.outfile.bgen.tmp, header = T, as.is = T) expect_equal(obj4.bgen[-83,], obj4.bgen.tmp[-83,]) unlink(c(obj4.outfile.gds.tmp, obj4.outfile.bgen.tmp)) ### adding one interaction covariate for obj4 obj4.outfile.gds.tmp.cov <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, outfile=obj4.outfile.gds.tmp.cov, ncores=2) obj4.gds.tmp.cov <- read.table(obj4.outfile.gds.tmp.cov, header = T, as.is = T) expect_equal(obj4.gds.cov, obj4.gds.tmp.cov,tolerance=.000001) obj4.outfile.bgen.tmp.cov <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"),interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen.tmp.cov, ncores=2) obj4.bgen.tmp.cov <- read.table(obj4.outfile.bgen.tmp.cov, header = T, as.is = T) expect_equal(obj4.bgen.cov[c(-64,-83),], obj4.bgen.tmp.cov[c(-64,-83),]) expect_equal(obj4.bgen.cov[-83,], obj4.bgen.tmp.cov[-83,]) unlink(c(obj4.outfile.gds.tmp.cov, obj4.outfile.bgen.tmp.cov)) ### re-order pheno id ### obj1 idx <- sample(nrow(pheno)) pheno <- pheno[idx, ] obj1 <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit")) obj1.outfile.gds <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=gdsfile, outfile=obj1.outfile.gds) tmpout <- read.table(obj1.outfile.gds, header = T, as.is = T) expect_equal(obj1.gds, tmpout) obj1.outfile.bgen <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen) tmpout <- read.table(obj1.outfile.bgen, header = T, as.is = T) expect_equal(obj1.bgen, tmpout) unlink(c(obj1.outfile.gds, obj1.outfile.bgen)) ### adding one interaction covariate for obj1 obj1.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, outfile=obj1.outfile.gds.cov) tmpout.cov <- read.table(obj1.outfile.gds.cov, header = T, as.is = T) expect_equal(obj1.gds.cov, tmpout.cov) obj1.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj1, interaction="sex",interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen.cov) tmpout.cov <- read.table(obj1.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj1.bgen.cov, tmpout.cov) unlink(c(obj1.outfile.gds.cov, obj1.outfile.bgen.cov)) ### obj2 obj2 <- glmmkin(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit")) obj2.outfile.gds <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=gdsfile, outfile=obj2.outfile.gds) tmpout <- read.table(obj2.outfile.gds, header = T, as.is = T) expect_equal(obj2.gds, tmpout) obj2.outfile.bgen <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen) tmpout <- read.table(obj2.outfile.bgen, header = T, as.is = T) expect_equal(obj2.bgen, tmpout) unlink(c(obj2.outfile.gds, obj2.outfile.bgen)) ### adding one interaction covariate for obj2 obj2.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, outfile=obj2.outfile.gds.cov) tmpout.cov <- read.table(obj2.outfile.gds.cov, header = T, as.is = T) expect_equal(obj2.gds.cov, tmpout.cov) obj2.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj2, interaction="sex",interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen.cov) tmpout.cov <- read.table(obj2.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj2.bgen.cov, tmpout.cov) unlink(c(obj2.outfile.gds.cov, obj2.outfile.bgen.cov)) ### re-order pheno id with two interaction terms ### obj3 obj3 <- glmmkin(disease ~ age + sex + obese, data = pheno, kins = kins, id = "id", family = binomial(link = "logit")) obj3.outfile.gds <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj3.outfile.gds) tmpout <- read.table(obj3.outfile.gds, header = T, as.is = T) expect_equal(obj3.gds, tmpout) obj3.outfile.bgen <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen) tmpout <- read.table(obj3.outfile.bgen, header = T, as.is = T) expect_equal(obj3.bgen[-83,], tmpout[-83,]) unlink(c(obj3.outfile.gds, obj3.outfile.bgen)) ### adding one interaction covariate for obj3 obj3.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, outfile=obj3.outfile.gds.cov) tmpout.cov <- read.table(obj3.outfile.gds.cov, header = T, as.is = T) expect_equal(obj3.gds.cov, tmpout.cov) obj3.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"),interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen.cov) tmpout.cov <- read.table(obj3.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj3.bgen.cov[c(-64,-83),], tmpout.cov[c(-64,-83),]) unlink(c(obj3.outfile.gds.cov, obj3.outfile.bgen.cov)) ### obj4 obj4 <- glmmkin(disease ~ age + sex + obese, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit")) obj4.outfile.gds <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj4.outfile.gds) tmpout <- read.table(obj4.outfile.gds, header = T, as.is = T) expect_equal(obj4.gds, tmpout,tolerance=.000001) obj4.outfile.bgen <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen) tmpout <- read.table(obj4.outfile.bgen, header = T, as.is = T) expect_equal(obj4.bgen[-83,], tmpout[-83,]) unlink(c(obj4.outfile.gds, obj4.outfile.bgen)) ### adding one interaction covariate for obj4 obj4.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, outfile=obj4.outfile.gds.cov) tmpout.cov <- read.table(obj4.outfile.gds.cov, header = T, as.is = T) expect_equal(obj4.gds.cov, tmpout.cov,tolerance=.000001) obj4.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"),interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen.cov) tmpout.cov <- read.table(obj4.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj4.bgen.cov[c(-64,-83),], tmpout.cov[c(-64,-83),]) unlink(c(obj4.outfile.gds.cov, obj4.outfile.bgen.cov)) ### re-order kins id ### obj1 idx <- sample(nrow(kins)) kins <- kins[idx, idx] obj1 <- glmmkin(disease ~ age + sex, data = pheno, kins = kins, id = "id", family = binomial(link = "logit")) obj1.outfile.gds <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=gdsfile, outfile=obj1.outfile.gds) tmpout <- read.table(obj1.outfile.gds, header = T, as.is = T) expect_equal(obj1.gds, tmpout) obj1.outfile.bgen <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen) tmpout <- read.table(obj1.outfile.bgen, header = T, as.is = T) expect_equal(obj1.bgen, tmpout) unlink(c(obj1.outfile.gds, obj1.outfile.bgen)) ### adding one interaction covariate for obj1 obj1.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, outfile=obj1.outfile.gds.cov) tmpout.cov <- read.table(obj1.outfile.gds.cov, header = T, as.is = T) expect_equal(obj1.gds.cov, tmpout.cov) obj1.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj1, interaction="sex",interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen.cov) tmpout.cov <- read.table(obj1.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj1.bgen.cov, tmpout.cov) unlink(c(obj1.outfile.gds.cov, obj1.outfile.bgen.cov)) ### obj2 obj2 <- glmmkin(disease ~ age + sex, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit")) obj2.outfile.gds <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=gdsfile, outfile=obj2.outfile.gds) tmpout <- read.table(obj2.outfile.gds, header = T, as.is = T) expect_equal(obj2.gds, tmpout) obj2.outfile.bgen <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen) tmpout <- read.table(obj2.outfile.bgen, header = T, as.is = T) expect_equal(obj2.bgen, tmpout) unlink(c(obj2.outfile.gds, obj2.outfile.bgen)) ### adding one interaction covariate for obj2 obj2.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, outfile=obj2.outfile.gds.cov) tmpout.cov <- read.table(obj2.outfile.gds.cov, header = T, as.is = T) expect_equal(obj2.gds.cov, tmpout.cov) obj2.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj2, interaction="sex",interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen.cov) tmpout.cov <- read.table(obj2.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj2.bgen.cov, tmpout.cov) unlink(c(obj2.outfile.gds.cov, obj2.outfile.bgen.cov)) ### re-order kins id with two interaction terms ### obj3 obj3 <- glmmkin(disease ~ age + sex + obese, data = pheno, kins = kins, id = "id", family = binomial(link = "logit")) obj3.outfile.gds <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj3.outfile.gds) tmpout <- read.table(obj3.outfile.gds, header = T, as.is = T) expect_equal(obj3.gds, tmpout) obj3.outfile.bgen <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen) tmpout <- read.table(obj3.outfile.bgen, header = T, as.is = T) expect_equal(obj3.bgen[-83,], tmpout[-83,]) unlink(c(obj3.outfile.gds, obj3.outfile.bgen)) ### adding one interaction covariate for obj3 obj3.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, outfile=obj3.outfile.gds.cov) tmpout.cov <- read.table(obj3.outfile.gds.cov, header = T, as.is = T) expect_equal(obj3.gds.cov, tmpout.cov) obj3.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"),interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen.cov) tmpout.cov <- read.table(obj3.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj3.bgen.cov[c(-64,-83),], tmpout.cov[c(-64,-83),]) unlink(c(obj3.outfile.gds.cov, obj3.outfile.bgen.cov)) ### obj4 obj4 <- glmmkin(disease ~ age + sex + obese, data = pheno, kins = NULL, id = "id", family = binomial(link = "logit")) obj4.outfile.gds <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj4.outfile.gds) tmpout <- read.table(obj4.outfile.gds, header = T, as.is = T) expect_equal(obj4.gds, tmpout,tolerance=.000001) obj4.outfile.bgen <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen) tmpout <- read.table(obj4.outfile.bgen, header = T, as.is = T) expect_equal(obj4.bgen[-83,], tmpout[-83,]) unlink(c(obj4.outfile.gds, obj4.outfile.bgen)) ### adding one interaction covariate for obj4 obj4.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, outfile=obj4.outfile.gds.cov) tmpout.cov <- read.table(obj4.outfile.gds.cov, header = T, as.is = T) expect_equal(obj4.gds.cov, tmpout.cov,tolerance=.000001) obj4.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"),interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen.cov) tmpout.cov <- read.table(obj4.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj4.bgen.cov[c(-64,-83),], tmpout.cov[c(-64,-83),]) unlink(c(obj4.outfile.gds.cov, obj4.outfile.bgen.cov)) }) test_that("longitudinal random time trend gaussian", { skip_on_cran() skip_if_not_installed("SeqArray") skip_if_not_installed("SeqVarTools") gdsfile <- system.file("extdata", "geno.gds", package = "MAGEE") bgenfile <- system.file("extdata", "geno.bgen", package = "MAGEE") samplefile <- system.file("extdata", "geno.sample", package = "MAGEE") group.file <- system.file("extdata", "SetID.singlevariant.txt", package = "MAGEE") data(example) set.seed(123) pheno <- example$pheno2 pheno$y.trend[sample(1:2000,20)] <- NA pheno$sex[sample(1:2000,20)] <- NA pheno$time[sample(1:2000,20)] <- NA pheno <- pheno[sample(1:2000,1950), ] pheno$obese <- rbinom(nrow(pheno), 1, 0.3) pheno$age<-rbinom(nrow(pheno),80,0.5) kins <- example$GRM ### single thread with kins obj1 <- glmmkin(y.trend ~ sex + time + age, data = pheno, kins = kins, id = "id",random.slope = "time", family = gaussian(link = "identity")) out1 <- MAGEE(null.obj=obj1, interaction="sex", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj1.outfile.gds <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=gdsfile, outfile=obj1.outfile.gds) obj1.gds <- read.table(obj1.outfile.gds, header = T, as.is = T) expect_equal(signif(range(obj1.gds$P_Value_Marginal)), signif(c(0.01776223, 0.99662787))) expect_equal(signif(range(obj1.gds$P_Value_Interaction)), signif(c(0.007458706, 0.980579058))) expect_equal(signif(range(obj1.gds$P_Value_Joint)), signif(c(0.01283483, 0.99295392))) expect_equal(signif(out1$MV.pval,5), signif(obj1.gds$P_Value_Marginal,5)) expect_equal(signif(out1$MF.pval), signif(obj1.gds$P_Value_Marginal)) expect_equal(signif(out1$IV.pval,4), signif(obj1.gds$P_Value_Interaction,4)) expect_equal(signif(out1$IF.pval), signif(obj1.gds$P_Value_Interaction)) expect_equal(signif(pchisq(qchisq(out1$MV.pval,2,lower=F)+qchisq(out1$IV.pval,2,lower=F),4,lower=F)), signif(out1$JV.pval)) expect_equal(signif(out1$JV.pval,3), signif(out1$JF.pval,3)) expect_equal(signif(out1$JV.pval,3), signif(out1$JD.pval,3)) expect_equal(signif(pchisq(qchisq(obj1.gds$P_Value_Marginal,1,lower=F)+qchisq(obj1.gds$P_Value_Interaction,1,lower=F),2,lower=F)), signif(obj1.gds$P_Value_Joint)) obj1.outfile.bgen <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen) obj1.bgen <- read.table(obj1.outfile.bgen, header = T, as.is = T) expect_equal(obj1.gds, obj1.bgen[,-2],tolerance=.0001) unlink(c(obj1.outfile.gds, obj1.outfile.bgen)) ### adding one interaction covariate out1.cov <- MAGEE(null.obj=obj1, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj1.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, outfile=obj1.outfile.gds.cov) obj1.gds.cov <- read.table(obj1.outfile.gds.cov, header = T, as.is = T) expect_equal(signif(range(obj1.gds.cov$P_Value_Marginal)), signif(c(0.01776223, 0.99662787))) expect_equal(signif(range(obj1.gds.cov$P_Value_Interaction)), signif(c(0.006627583, 0.987577820))) expect_equal(signif(range(obj1.gds.cov$P_Value_Joint)), signif(c(0.01240654, 0.99399595))) expect_equal(signif(out1.cov$MV.pval,4), signif(obj1.gds.cov$P_Value_Marginal,4)) expect_equal(signif(out1.cov$MF.pval), signif(obj1.gds.cov$P_Value_Marginal)) expect_equal(signif(out1.cov$IV.pval,4), signif(obj1.gds.cov$P_Value_Interaction,4)) expect_equal(signif(out1.cov$IF.pval), signif(obj1.gds.cov$P_Value_Interaction)) expect_equal(signif(pchisq(qchisq(pchisq(qchisq(obj1.gds.cov$P_Value_Joint, df=2, lower=F)-qchisq(obj1.gds.cov$P_Value_Interaction, df=1, lower=F), df=1, lower=F),2,lower=F)+qchisq(out1.cov$IV.pval,2,lower=F),4,lower=F),4), signif(out1.cov$JV.pval,4)) expect_equal(signif(out1.cov$JV.pval,3), signif(out1.cov$JF.pval,3)) expect_equal(signif(out1.cov$JV.pval,3), signif(out1.cov$JD.pval,3)) obj1.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen.cov) obj1.bgen.cov <- read.table(obj1.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj1.gds.cov, obj1.bgen.cov[,-2],tolerance=.0001) unlink(c(obj1.outfile.gds.cov, obj1.outfile.bgen.cov)) obj1.glmm.score.outfile.gds <- tempfile() glmm.score(obj=obj1, infile=gdsfile, outfile=obj1.glmm.score.outfile.gds) obj1.glmm.score.gds <- read.table(obj1.glmm.score.outfile.gds,header=T,as.is=T) expect_equal(obj1.gds$P_Value_Marginal, obj1.glmm.score.gds$PVAL) obj1.glmm.score.outfile.bgen <- tempfile() glmm.score(obj=obj1, infile=bgenfile, BGEN.samplefile = samplefile, outfile=obj1.glmm.score.outfile.bgen) obj1.glmm.score.bgen <- read.table(obj1.glmm.score.outfile.bgen,header=T,as.is=T) expect_equal(signif(obj1.gds$P_Value_Marginal), signif(obj1.glmm.score.bgen$PVAL)) unlink(c(obj1.glmm.score.outfile.gds, obj1.glmm.score.outfile.bgen)) out1.time <- MAGEE(null.obj=obj1, interaction="time", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj1.outfile.gds <- tempfile() glmm.gei(null.obj=obj1, interaction="time", geno.file=gdsfile, outfile=obj1.outfile.gds) obj1.time.gds <- read.table(obj1.outfile.gds, header = T, as.is = T) expect_equal(obj1.gds$P_Value_Marginal, obj1.time.gds$P_Value_Marginal) expect_equal(signif(range(obj1.time.gds$P_Value_Interaction)), signif(c(0.01464601, 0.98577260))) expect_equal(signif(range(obj1.time.gds$P_Value_Joint)), signif(c(0.0282489, 0.9875144))) expect_equal(out1$MV.pval, out1.time$MV.pval) expect_equal(out1$MF.pval, out1.time$MF.pval) expect_equal(signif(out1.time$IV.pval,3), signif(obj1.time.gds$P_Value_Interaction,3)) expect_equal(signif(out1.time$IF.pval), signif(obj1.time.gds$P_Value_Interaction)) expect_equal(signif(pchisq(qchisq(out1.time$MV.pval,2,lower=F)+qchisq(out1.time$IV.pval,2,lower=F),4,lower=F)), signif(out1.time$JV.pval)) expect_equal(signif(out1.time$JV.pval,3), signif(out1.time$JF.pval,3)) expect_equal(signif(out1.time$JV.pval,3), signif(out1.time$JD.pval,3)) expect_equal(signif(pchisq(qchisq(obj1.time.gds$P_Value_Marginal,1,lower=F)+qchisq(obj1.time.gds$P_Value_Interaction,1,lower=F),2,lower=F)), signif(obj1.time.gds$P_Value_Joint)) obj1.outfile.bgen <- tempfile() glmm.gei(null.obj=obj1, interaction="time", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen) obj1.time.bgen <- read.table(obj1.outfile.bgen, header = T, as.is = T) expect_equal(obj1.time.gds, obj1.time.bgen[,-2],tolerance=.0001) unlink(c(obj1.outfile.gds, obj1.outfile.bgen)) ### single thread without kins obj2 <- glmmkin(y.trend ~ sex + time + age, data = pheno, kins = NULL, id = "id",random.slope = "time", family = gaussian(link = "identity")) out2 <- MAGEE(null.obj=obj2, interaction="sex", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj2.outfile.gds <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=gdsfile, outfile=obj2.outfile.gds) obj2.gds <- read.table(obj2.outfile.gds, header = T, as.is = T) expect_equal(signif(range(obj2.gds$P_Value_Marginal)), signif(c(0.009131402, 0.999837921))) expect_equal(signif(range(obj2.gds$P_Value_Interaction)), signif(c(0.01209104, 0.98826124))) expect_equal(signif(range(obj2.gds$P_Value_Joint)), signif(c(0.01220809, 0.99232575))) expect_equal(signif(out2$MV.pval,4), signif(obj2.gds$P_Value_Marginal,4)) expect_equal(signif(out2$MF.pval), signif(obj2.gds$P_Value_Marginal)) expect_equal(signif(out2$IV.pval,4), signif(obj2.gds$P_Value_Interaction,4)) expect_equal(signif(out2$IF.pval), signif(obj2.gds$P_Value_Interaction)) expect_equal(signif(pchisq(qchisq(out2$MV.pval,2,lower=F)+qchisq(out2$IV.pval,2,lower=F),4,lower=F)), signif(out2$JV.pval)) expect_equal(signif(out2$JV.pval,4), signif(out2$JF.pval,4)) expect_equal(signif(out2$JV.pval,4), signif(out2$JD.pval,4)) expect_equal(signif(pchisq(qchisq(obj2.gds$P_Value_Marginal,1,lower=F)+qchisq(obj2.gds$P_Value_Interaction,1,lower=F),2,lower=F)), signif(obj2.gds$P_Value_Joint)) obj2.outfile.bgen <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen) obj2.bgen <- read.table(obj2.outfile.bgen, header = T, as.is = T) expect_equal(obj2.gds, obj2.bgen[,-2],tolerance=.0001) unlink(c(obj2.outfile.gds, obj2.outfile.bgen)) ### adding one interaction covariate out2.cov <- MAGEE(null.obj=obj2, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj2.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", interaction.covariates ="age", geno.file=gdsfile, outfile=obj2.outfile.gds.cov) obj2.gds.cov <- read.table(obj2.outfile.gds.cov, header = T, as.is = T) expect_equal(signif(range(obj2.gds.cov$P_Value_Marginal),7), signif(c(0.009131402, 0.999837921),7)) expect_equal(signif(range(obj2.gds.cov$P_Value_Interaction)), signif(c(0.01089319, 0.98903876))) expect_equal(signif(range(obj2.gds.cov$P_Value_Joint)), signif(c(0.01220468, 0.99843355))) expect_equal(signif(out2.cov$MV.pval,3), signif(obj2.gds.cov$P_Value_Marginal,3)) expect_equal(signif(out2.cov$MF.pval,4), signif(obj2.gds.cov$P_Value_Marginal,4)) expect_equal(signif(out2.cov$IV.pval,4), signif(obj2.gds.cov$P_Value_Interaction,4)) expect_equal(signif(out2.cov$IF.pval,4), signif(obj2.gds.cov$P_Value_Interaction,4)) expect_equal(signif(pchisq(qchisq(pchisq(qchisq(obj2.gds.cov$P_Value_Joint, df=2, lower=F)-qchisq(obj2.gds.cov$P_Value_Interaction, df=1, lower=F), df=1, lower=F),2,lower=F)+qchisq(out2.cov$IV.pval,2,lower=F),4,lower=F),3), signif(out2.cov$JV.pval,3)) expect_equal(signif(out2.cov$JV.pval,3), signif(out2.cov$JF.pval,3)) expect_equal(signif(out2.cov$JV.pval,3), signif(out2.cov$JD.pval,3)) obj2.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen.cov) obj2.bgen.cov <- read.table(obj2.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj2.gds.cov, obj2.bgen.cov[,-2],tolerance=.0001) unlink(c(obj2.outfile.gds.cov, obj2.outfile.bgen.cov)) obj2.glmm.score.outfile.gds <- tempfile() glmm.score(obj=obj2, infile=gdsfile, outfile=obj2.glmm.score.outfile.gds) obj2.glmm.score.gds <- read.table(obj2.glmm.score.outfile.gds,header=T,as.is=T) expect_equal(obj2.gds$P_Value_Marginal, obj2.glmm.score.gds$PVAL) obj2.glmm.score.outfile.bgen <- tempfile() glmm.score(obj=obj2, infile=bgenfile, BGEN.samplefile = samplefile, outfile=obj2.glmm.score.outfile.bgen) obj2.glmm.score.bgen <- read.table(obj2.glmm.score.outfile.bgen,header=T,as.is=T) expect_equal(signif(obj2.gds$P_Value_Marginal), signif(obj2.glmm.score.bgen$PVAL)) unlink(c(obj2.glmm.score.outfile.gds, obj2.glmm.score.outfile.bgen)) out2.time <- MAGEE(null.obj=obj2, interaction="time", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj2.outfile.gds <- tempfile() glmm.gei(null.obj=obj2, interaction="time", geno.file=gdsfile, outfile=obj2.outfile.gds) obj2.time.gds <- read.table(obj2.outfile.gds, header = T, as.is = T) expect_equal(obj2.gds$P_Value_Marginal, obj2.time.gds$P_Value_Marginal) expect_equal(signif(range(obj2.time.gds$P_Value_Interaction)), signif(c(0.006924606, 0.992354799))) expect_equal(signif(range(obj2.time.gds$P_Value_Joint),7), signif(c(0.0139802, 0.9988815),7)) expect_equal(out2$MV.pval, out2.time$MV.pval) expect_equal(out2$MF.pval, out2.time$MF.pval) expect_equal(signif(out2.time$IV.pval,4), signif(obj2.time.gds$P_Value_Interaction,4)) expect_equal(signif(out2.time$IF.pval), signif(obj2.time.gds$P_Value_Interaction)) expect_equal(signif(pchisq(qchisq(out2.time$MV.pval,2,lower=F)+qchisq(out2.time$IV.pval,2,lower=F),4,lower=F)), signif(out2.time$JV.pval)) expect_equal(signif(out2.time$JV.pval,3), signif(out2.time$JF.pval,3)) expect_equal(signif(out2.time$JV.pval,3), signif(out2.time$JD.pval,3)) expect_equal(signif(pchisq(qchisq(obj2.time.gds$P_Value_Marginal,1,lower=F)+qchisq(obj2.time.gds$P_Value_Interaction,1,lower=F),2,lower=F)), signif(obj2.time.gds$P_Value_Joint)) obj2.outfile.bgen <- tempfile() glmm.gei(null.obj=obj2, interaction="time", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen) obj2.time.bgen <- read.table(obj2.outfile.bgen, header = T, as.is = T) expect_equal(obj2.time.gds, obj2.time.bgen[,-2],tolerance=.0001) unlink(c(obj2.outfile.gds, obj2.outfile.bgen)) ### single thread with kins and with two interaction variables obj3 <- glmmkin(y.trend ~ sex + time + obese +age, data = pheno, kins = kins, id = "id",random.slope = "time", family = gaussian(link = "identity")) out3 <- MAGEE(null.obj=obj3, interaction=c("sex","obese"), geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj3.outfile.gds <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj3.outfile.gds) obj3.gds <- read.table(obj3.outfile.gds, header = T, as.is = T) expect_equal(signif(range(obj3.gds$P_Value_Marginal)), signif(c(0.01603888, 0.99980279))) expect_equal(signif(range(obj3.gds$P_Value_Interaction),7), signif(c( 0.002690285, 0.987263619),7)) expect_equal(signif(range(obj3.gds$P_Value_Joint)), signif(c(0.00789263, 0.99718114))) expect_equal(signif(out3$MV.pval,5), signif(obj3.gds$P_Value_Marginal,5)) expect_equal(signif(out3$MF.pval), signif(obj3.gds$P_Value_Marginal)) expect_equal(signif(pchisq(qchisq(out3$MV.pval, df =2 ,lower=F)+qchisq(out3$IV.pval, df=2, lower=F), df=4, lower=F),4 ),signif(out3$JV.pval,4 )) expect_equal(signif(pchisq(qchisq(out3$MF.pval, df =2 ,lower=F)+qchisq(out3$IF.pval, df=2, lower=F), df=4, lower=F),4 ),signif(out3$JD.pval,4 )) expect_equal(signif(pchisq(qchisq(obj3.gds$P_Value_Marginal,1,lower=F)+qchisq(obj3.gds$P_Value_Interaction,2,lower=F),3,lower=F)), signif(obj3.gds$P_Value_Joint)) obj3.outfile.bgen <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen) obj3.bgen <- read.table(obj3.outfile.bgen, header = T, as.is = T) expect_equal(obj3.gds, obj3.bgen[,-2],tolerance=.0001) unlink(c(obj3.outfile.gds, obj3.outfile.bgen)) ### adding one interaction covariate out3.cov <- MAGEE(null.obj=obj3, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj3.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, outfile=obj3.outfile.gds.cov) obj3.gds.cov <- read.table(obj3.outfile.gds.cov, header = T, as.is = T) expect_equal(signif(range(obj3.gds.cov$P_Value_Marginal)), signif(c(0.01603888, 0.99980279))) expect_equal(signif(range(obj3.gds.cov$P_Value_Interaction)), signif(c(0.00267669, 0.99726949))) expect_equal(signif(range(obj3.gds.cov$P_Value_Joint)), signif(c(0.007747924, 0.999369233))) expect_equal(signif(out3.cov$MV.pval,4), signif(obj3.gds.cov$P_Value_Marginal,4)) expect_equal(signif(out3.cov$MF.pval), signif(obj3.gds.cov$P_Value_Marginal)) expect_equal(signif(pchisq(qchisq(pchisq(qchisq(obj3.gds.cov$P_Value_Joint[c(-64,-83)], df=3, lower=F)-qchisq(obj3.gds.cov$P_Value_Interaction[c(-64,-83)], df=2, lower=F), df=1, lower=F), df =2 ,lower=F)+qchisq(out3.cov$IV.pval[c(-64,-83)], df=2, lower=F), df=4, lower=F),4 ),signif(out3.cov$JV.pval[c(-64,-83)],4 )) expect_equal(signif(pchisq(qchisq(pchisq(qchisq(obj3.gds.cov$P_Value_Joint[c(-64,-83)], df=3, lower=F)-qchisq(obj3.gds.cov$P_Value_Interaction[c(-64,-83)], df=2, lower=F), df=1, lower=F), df =2 ,lower=F)+qchisq(out3.cov$IF.pval[c(-64,-83)], df=2, lower=F), df=4, lower=F),4 ),signif(out3.cov$JD.pval[c(-64,-83)],4 )) obj3.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen.cov) obj3.bgen.cov <- read.table(obj3.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj3.gds.cov[c(-64,-83),], obj3.bgen.cov[c(-64,-83),-2],tolerance=.0001) unlink(c(obj3.outfile.gds.cov, obj3.outfile.bgen.cov)) obj3.glmm.score.outfile.gds <- tempfile() glmm.score(obj=obj3, infile=gdsfile, outfile=obj3.glmm.score.outfile.gds) obj3.glmm.score.gds <- read.table(obj3.glmm.score.outfile.gds,header=T,as.is=T) expect_equal(obj3.gds$P_Value_Marginal, obj3.glmm.score.gds$PVAL) obj3.glmm.score.outfile.bgen <- tempfile() glmm.score(obj=obj3, infile=bgenfile, BGEN.samplefile = samplefile, outfile=obj3.glmm.score.outfile.bgen) obj3.glmm.score.bgen <- read.table(obj3.glmm.score.outfile.bgen,header=T,as.is=T) expect_equal(signif(obj3.gds$P_Value_Marginal), signif(obj3.glmm.score.bgen$PVAL)) unlink(c(obj3.glmm.score.outfile.gds, obj3.glmm.score.outfile.bgen)) out3.time <- MAGEE(null.obj=obj3, interaction=c("time","sex"), geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj3.outfile.gds <- tempfile() glmm.gei(null.obj=obj3, interaction=c("time","sex"), geno.file=gdsfile, outfile=obj3.outfile.gds) obj3.time.gds <- read.table(obj3.outfile.gds, header = T, as.is = T) expect_equal(obj3.gds$P_Value_Marginal, obj3.time.gds$P_Value_Marginal) expect_equal(signif(range(obj3.time.gds$P_Value_Interaction)), signif(c(0.01105401, 0.97060979))) expect_equal(signif(range(obj3.time.gds$P_Value_Joint)), signif(c(0.01922166, 0.98723879))) expect_equal(out3$MV.pval, out3.time$MV.pval) expect_equal(out3$MF.pval, out3.time$MF.pval) expect_equal(signif(pchisq(qchisq(out3.time$MV.pval, df =2 ,lower=F)+qchisq(out3.time$IV.pval, df=2, lower=F), df=4, lower=F),4 ),signif(out3.time$JV.pval,4 )) expect_equal(signif(pchisq(qchisq(out3.time$MF.pval, df =2 ,lower=F)+qchisq(out3.time$IF.pval, df=2, lower=F), df=4, lower=F),4 ),signif(out3.time$JD.pval,4 )) expect_equal(signif(pchisq(qchisq(obj3.time.gds$P_Value_Marginal,1,lower=F)+qchisq(obj3.time.gds$P_Value_Interaction,2,lower=F),3,lower=F)), signif(obj3.time.gds$P_Value_Joint)) obj3.outfile.bgen <- tempfile() glmm.gei(null.obj=obj3, interaction=c("time","sex"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen) obj3.time.bgen <- read.table(obj3.outfile.bgen, header = T, as.is = T) expect_equal(obj3.time.gds, obj3.time.bgen[,-2],tolerance=.0001) unlink(c(obj3.outfile.gds, obj3.outfile.bgen)) ### single thread without kins and with two interaction variables obj4 <- glmmkin(y.trend ~ sex + time + obese + age, data = pheno, kins = NULL, id = "id",random.slope = "time", family = gaussian(link = "identity")) out4 <- MAGEE(null.obj=obj4, interaction=c("sex","obese"), geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj4.outfile.gds <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj4.outfile.gds) obj4.gds <- read.table(obj4.outfile.gds, header = T, as.is = T) expect_equal(signif(range(obj4.gds$P_Value_Marginal)), signif(c(0.008155485, 0.990967332))) expect_equal(signif(range(obj4.gds$P_Value_Interaction)), signif(c(0.002348382, 0.989754496))) expect_equal(signif(range(obj4.gds$P_Value_Joint)), signif(c(0.007019809, 0.997505989))) expect_equal(signif(out4$MV.pval,3), signif(obj4.gds$P_Value_Marginal,3)) expect_equal(signif(out4$MF.pval), signif(obj4.gds$P_Value_Marginal)) expect_equal(signif(pchisq(qchisq(out4$MV.pval, df =2 ,lower=F)+qchisq(out4$IV.pval, df=2, lower=F), df=4, lower=F),4 ),signif(out4$JV.pval,4 )) expect_equal(signif(pchisq(qchisq(out4$MF.pval, df =2 ,lower=F)+qchisq(out4$IF.pval, df=2, lower=F), df=4, lower=F),4 ),signif(out4$JD.pval,4 )) expect_equal(signif(pchisq(qchisq(obj4.gds$P_Value_Marginal,1,lower=F)+qchisq(obj4.gds$P_Value_Interaction,2,lower=F),3,lower=F)), signif(obj4.gds$P_Value_Joint)) obj4.outfile.bgen <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen) obj4.bgen <- read.table(obj4.outfile.bgen, header = T, as.is = T) expect_equal(obj4.gds, obj4.bgen[,-2],tolerance=.0001) unlink(c(obj4.outfile.gds, obj4.outfile.bgen)) ### adding one interaction covariate out4.cov <- MAGEE(null.obj=obj4, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj4.outfile.gds.cov <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=gdsfile, outfile=obj4.outfile.gds.cov) obj4.gds.cov <- read.table(obj4.outfile.gds.cov, header = T, as.is = T) expect_equal(signif(range(obj4.gds.cov$P_Value_Marginal)), signif(c(0.008155485, 0.990967332))) expect_equal(signif(range(obj4.gds.cov$P_Value_Interaction)), signif(c(0.00233817, 0.99006853))) expect_equal(signif(range(obj4.gds.cov$P_Value_Joint)), signif(c(0.006957274, 0.998853901))) expect_equal(signif(out4.cov$MV.pval,2), signif(obj4.gds.cov$P_Value_Marginal,2)) expect_equal(signif(out4.cov$MF.pval,2), signif(obj4.gds.cov$P_Value_Marginal,2)) expect_equal(signif(pchisq(qchisq(pchisq(qchisq(obj4.gds.cov$P_Value_Joint[c(-64,-83)], df=3, lower=F)-qchisq(obj4.gds.cov$P_Value_Interaction[c(-64,-83)], df=2, lower=F), df=1, lower=F), df =2 ,lower=F)+qchisq(out4.cov$IV.pval[c(-64,-83)], df=2, lower=F), df=4, lower=F),3),signif(out4.cov$JV.pval[c(-64,-83)],3)) expect_equal(signif(pchisq(qchisq(pchisq(qchisq(obj4.gds.cov$P_Value_Joint[c(-64,-83)], df=3, lower=F)-qchisq(obj4.gds.cov$P_Value_Interaction[c(-64,-83)], df=2, lower=F), df=1, lower=F), df =2 ,lower=F)+qchisq(out4.cov$IF.pval[c(-64,-83)], df=2, lower=F), df=4, lower=F),4 ),signif(out4.cov$JD.pval[c(-64,-83)],4 )) obj4.outfile.bgen.cov <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), interaction.covariates ="age", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen.cov) obj4.bgen.cov <- read.table(obj4.outfile.bgen.cov, header = T, as.is = T) expect_equal(obj4.gds.cov[c(-64,-83),], obj4.bgen.cov[c(-64,-83),-2],tolerance=.0001) unlink(c(obj4.outfile.gds.cov, obj4.outfile.bgen.cov)) obj4.glmm.score.outfile.gds <- tempfile() glmm.score(obj=obj4, infile=gdsfile, outfile=obj4.glmm.score.outfile.gds) obj4.glmm.score.gds <- read.table(obj4.glmm.score.outfile.gds,header=T,as.is=T) expect_equal(obj4.gds$P_Value_Marginal, obj4.glmm.score.gds$PVAL) obj4.glmm.score.outfile.bgen <- tempfile() glmm.score(obj=obj4, infile=bgenfile, BGEN.samplefile = samplefile, outfile=obj4.glmm.score.outfile.bgen) obj4.glmm.score.bgen <- read.table(obj4.glmm.score.outfile.bgen,header=T,as.is=T) expect_equal(signif(obj4.gds$P_Value_Marginal), signif(obj4.glmm.score.bgen$PVAL)) unlink(c(obj4.glmm.score.outfile.gds, obj4.glmm.score.outfile.bgen)) out4.time <- MAGEE(null.obj=obj4, interaction=c("time","sex"), geno.file=gdsfile, group.file=group.file, tests = c("JV", "JF", "JD"), use.minor.allele = T, auto.flip = T) obj4.outfile.gds <- tempfile() glmm.gei(null.obj=obj4, interaction=c("time","sex"), geno.file=gdsfile, outfile=obj4.outfile.gds) obj4.time.gds <- read.table(obj4.outfile.gds, header = T, as.is = T) expect_equal(obj4.gds$P_Value_Marginal, obj4.time.gds$P_Value_Marginal) expect_equal(signif(range(obj4.time.gds$P_Value_Interaction)), signif(c(0.005718586, 0.988761356))) expect_equal(signif(range(obj4.time.gds$P_Value_Joint)), signif(c(0.00902537, 0.98766190))) expect_equal(out4$MV.pval, out4.time$MV.pval) expect_equal(out4$MF.pval, out4.time$MF.pval) expect_equal(signif(pchisq(qchisq(out4.time$MV.pval, df =2 ,lower=F)+qchisq(out4.time$IV.pval, df=2, lower=F), df=4, lower=F),4 ),signif(out4.time$JV.pval,4 )) expect_equal(signif(pchisq(qchisq(out4.time$MF.pval, df =2 ,lower=F)+qchisq(out4.time$IF.pval, df=2, lower=F), df=4, lower=F),4 ),signif(out4.time$JD.pval,4 )) expect_equal(signif(pchisq(qchisq(obj4.time.gds$P_Value_Marginal,1,lower=F)+qchisq(obj4.time.gds$P_Value_Interaction,2,lower=F),3,lower=F)), signif(obj4.time.gds$P_Value_Joint)) obj4.outfile.bgen <- tempfile() glmm.gei(null.obj=obj4, interaction=c("time","sex"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen) obj4.time.bgen <- read.table(obj4.outfile.bgen, header = T, as.is = T) expect_equal(obj4.time.gds, obj4.time.bgen[,-2],tolerance=.0001) unlink(c(obj4.outfile.gds, obj4.outfile.bgen)) ### multi-thread obj1.outfile.gds.tmp <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=gdsfile, outfile=obj1.outfile.gds.tmp, ncores=2) obj1.gds.tmp <- read.table(obj1.outfile.gds.tmp, header = T, as.is = T) expect_equal(obj1.gds, obj1.gds.tmp) obj1.outfile.bgen.tmp <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen.tmp, ncores=2) obj1.bgen.tmp <- read.table(obj1.outfile.bgen.tmp, header = T, as.is = T) expect_equal(obj1.bgen, obj1.bgen.tmp) unlink(c(obj1.outfile.gds.tmp, obj1.outfile.bgen.tmp)) obj1.outfile.gds.tmp <- tempfile() glmm.gei(null.obj=obj1, interaction="time", geno.file=gdsfile, outfile=obj1.outfile.gds.tmp, ncores=2) obj1.time.gds.tmp <- read.table(obj1.outfile.gds.tmp, header = T, as.is = T) expect_equal(obj1.time.gds, obj1.time.gds.tmp) obj1.outfile.bgen.tmp <- tempfile() glmm.gei(null.obj=obj1, interaction="time", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen.tmp, ncores=2) obj1.time.bgen.tmp <- read.table(obj1.outfile.bgen.tmp, header = T, as.is = T) expect_equal(obj1.time.bgen, obj1.time.bgen.tmp) unlink(c(obj1.outfile.gds.tmp, obj1.outfile.bgen.tmp)) obj2.outfile.gds.tmp <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=gdsfile, outfile=obj2.outfile.gds.tmp, ncores=2) obj2.gds.tmp <- read.table(obj2.outfile.gds.tmp, header = T, as.is = T) expect_equal(obj2.gds, obj2.gds.tmp) obj2.outfile.bgen.tmp <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen.tmp, ncores=2) obj2.bgen.tmp <- read.table(obj2.outfile.bgen.tmp, header = T, as.is = T) expect_equal(obj2.bgen, obj2.bgen.tmp) unlink(c(obj2.outfile.gds.tmp, obj2.outfile.bgen.tmp)) obj2.outfile.gds.tmp <- tempfile() glmm.gei(null.obj=obj2, interaction="time", geno.file=gdsfile, outfile=obj2.outfile.gds.tmp, ncores=2) obj2.time.gds.tmp <- read.table(obj2.outfile.gds.tmp, header = T, as.is = T) expect_equal(obj2.time.gds, obj2.time.gds.tmp) obj2.outfile.bgen.tmp <- tempfile() glmm.gei(null.obj=obj2, interaction="time", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen.tmp, ncores=2) obj2.time.bgen.tmp <- read.table(obj2.outfile.bgen.tmp, header = T, as.is = T) expect_equal(obj2.time.bgen, obj2.time.bgen.tmp) unlink(c(obj2.outfile.gds.tmp, obj2.outfile.bgen.tmp)) ### multi-thread with two interaction variables obj3.outfile.gds.tmp <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj3.outfile.gds.tmp, ncores=2) obj3.gds.tmp <- read.table(obj3.outfile.gds.tmp, header = T, as.is = T) expect_equal(obj3.gds, obj3.gds.tmp) obj3.outfile.bgen.tmp <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen.tmp, ncores=2) obj3.bgen.tmp <- read.table(obj3.outfile.bgen.tmp, header = T, as.is = T) expect_equal(obj3.bgen, obj3.bgen.tmp) unlink(c(obj3.outfile.gds.tmp, obj3.outfile.bgen.tmp)) obj3.outfile.gds.tmp <- tempfile() glmm.gei(null.obj=obj3, interaction=c("time","sex"), geno.file=gdsfile, outfile=obj3.outfile.gds.tmp, ncores=2) obj3.time.gds.tmp <- read.table(obj3.outfile.gds.tmp, header = T, as.is = T) expect_equal(obj3.time.gds, obj3.time.gds.tmp) obj3.outfile.bgen.tmp <- tempfile() glmm.gei(null.obj=obj3, interaction=c("time","sex"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen.tmp, ncores=2) obj3.time.bgen.tmp <- read.table(obj3.outfile.bgen.tmp, header = T, as.is = T) expect_equal(obj3.time.bgen, obj3.time.bgen.tmp) unlink(c(obj3.outfile.gds.tmp, obj3.outfile.bgen.tmp)) obj4.outfile.gds.tmp <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj4.outfile.gds.tmp, ncores=2) obj4.gds.tmp <- read.table(obj4.outfile.gds.tmp, header = T, as.is = T) expect_equal(obj4.gds, obj4.gds.tmp) obj4.outfile.bgen.tmp <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen.tmp, ncores=2) obj4.bgen.tmp <- read.table(obj4.outfile.bgen.tmp, header = T, as.is = T) expect_equal(obj4.bgen, obj4.bgen.tmp) unlink(c(obj4.outfile.gds.tmp, obj4.outfile.bgen.tmp)) obj4.outfile.gds.tmp <- tempfile() glmm.gei(null.obj=obj4, interaction=c("time","sex"), geno.file=gdsfile, outfile=obj4.outfile.gds.tmp, ncores=2) obj4.time.gds.tmp <- read.table(obj4.outfile.gds.tmp, header = T, as.is = T) expect_equal(obj4.time.gds, obj4.time.gds.tmp) obj4.outfile.bgen.tmp <- tempfile() glmm.gei(null.obj=obj4, interaction=c("time","sex"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen.tmp, ncores=2) obj4.time.bgen.tmp <- read.table(obj4.outfile.bgen.tmp, header = T, as.is = T) expect_equal(obj4.time.bgen, obj4.time.bgen.tmp) unlink(c(obj4.outfile.gds.tmp, obj4.outfile.bgen.tmp)) ### re-order pheno id idx <- sample(nrow(pheno)) pheno <- pheno[idx, ] obj1 <- glmmkin(y.trend ~ sex + time +age, data = pheno, kins = kins, id = "id",random.slope = "time", family = gaussian(link = "identity")) obj1.outfile.gds <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=gdsfile, outfile=obj1.outfile.gds) tmpout <- read.table(obj1.outfile.gds, header = T, as.is = T) expect_equal(obj1.gds[,-(10:13)], tmpout[,-(10:13)]) obj1.outfile.bgen <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen) tmpout <- read.table(obj1.outfile.bgen, header = T, as.is = T) expect_equal(obj1.bgen[,-(11:14)], tmpout[,-(11:14)]) unlink(c(obj1.outfile.gds, obj1.outfile.bgen)) obj1.outfile.gds <- tempfile() glmm.gei(null.obj=obj1, interaction="time", geno.file=gdsfile, outfile=obj1.outfile.gds) tmpout <- read.table(obj1.outfile.gds, header = T, as.is = T) expect_equal(obj1.time.gds[,-(10:19)], tmpout[,-(10:19)]) obj1.outfile.bgen <- tempfile() glmm.gei(null.obj=obj1, interaction="time", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen) tmpout <- read.table(obj1.outfile.bgen, header = T, as.is = T) expect_equal(obj1.time.bgen[,-(11:20)], tmpout[,-(11:20)]) unlink(c(obj1.outfile.gds, obj1.outfile.bgen)) obj2 <- glmmkin(y.trend ~ sex + time + age, data = pheno, kins = NULL, id = "id",random.slope = "time", family = gaussian(link = "identity")) obj2.outfile.gds <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=gdsfile, outfile=obj2.outfile.gds) tmpout <- read.table(obj2.outfile.gds, header = T, as.is = T) expect_equal(obj2.gds[,-(10:13)], tmpout[,-(10:13)]) obj2.outfile.bgen <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen) tmpout <- read.table(obj2.outfile.bgen, header = T, as.is = T) expect_equal(obj2.bgen[,-(11:14)], tmpout[,-(11:14)]) unlink(c(obj2.outfile.gds, obj2.outfile.bgen)) obj2.outfile.gds <- tempfile() glmm.gei(null.obj=obj2, interaction="time", geno.file=gdsfile, outfile=obj2.outfile.gds) tmpout <- read.table(obj2.outfile.gds, header = T, as.is = T) expect_equal(obj2.time.gds[,-(10:19)], tmpout[,-(10:19)]) obj2.outfile.bgen <- tempfile() glmm.gei(null.obj=obj2, interaction="time", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen) tmpout <- read.table(obj2.outfile.bgen, header = T, as.is = T) expect_equal(obj2.time.bgen[,-(11:20)], tmpout[,-(11:20)]) unlink(c(obj2.outfile.gds, obj2.outfile.bgen)) ### re-order pheno id with two interaction terms obj3 <- glmmkin(y.trend ~ sex + time+ obese + age, data = pheno, kins = kins, id = "id",random.slope = "time", family = gaussian(link = "identity")) obj3.outfile.gds <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj3.outfile.gds) tmpout <- read.table(obj3.outfile.gds, header = T, as.is = T) expect_equal(obj3.gds[,-(10:17)], tmpout[,-(10:17)]) obj3.outfile.bgen <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen) tmpout <- read.table(obj3.outfile.bgen, header = T, as.is = T) expect_equal(obj3.bgen[,-(11:18)], tmpout[,-(11:18)]) unlink(c(obj3.outfile.gds, obj3.outfile.bgen)) obj3.outfile.gds <- tempfile() glmm.gei(null.obj=obj3, interaction=c("time","sex"), geno.file=gdsfile, outfile=obj3.outfile.gds) tmpout <- read.table(obj3.outfile.gds, header = T, as.is = T) expect_equal(obj3.time.gds[,-(10:29)], tmpout[,-(10:29)]) obj3.outfile.bgen <- tempfile() glmm.gei(null.obj=obj3, interaction=c("time","sex"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen) tmpout <- read.table(obj3.outfile.bgen, header = T, as.is = T) expect_equal(obj3.time.bgen[,-(11:30)], tmpout[,-(11:30)]) unlink(c(obj3.outfile.gds, obj3.outfile.bgen)) obj4 <- glmmkin(y.trend ~ sex + time + obese + age, data = pheno, kins = NULL, id = "id",random.slope = "time", family = gaussian(link = "identity")) obj4.outfile.gds <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj4.outfile.gds) tmpout <- read.table(obj4.outfile.gds, header = T, as.is = T) expect_equal(obj4.gds[,-(10:17)], tmpout[,-(10:17)]) obj4.outfile.bgen <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen) tmpout <- read.table(obj4.outfile.bgen, header = T, as.is = T) expect_equal(obj4.bgen[,-(11:18)], tmpout[,-(11:18)]) unlink(c(obj4.outfile.gds, obj4.outfile.bgen)) obj4.outfile.gds <- tempfile() glmm.gei(null.obj=obj4, interaction=c("time","sex"), geno.file=gdsfile, outfile=obj4.outfile.gds) tmpout <- read.table(obj4.outfile.gds, header = T, as.is = T) expect_equal(obj4.time.gds[,-(10:29)], tmpout[,-(10:29)]) obj4.outfile.bgen <- tempfile() glmm.gei(null.obj=obj4, interaction=c("time","sex"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen) tmpout <- read.table(obj4.outfile.bgen, header = T, as.is = T) expect_equal(obj4.time.bgen[,-(11:30)], tmpout[,-(11:30)]) unlink(c(obj4.outfile.gds, obj4.outfile.bgen)) ### re-order id idx <- sample(nrow(kins)) kins <- kins[idx, idx] obj1 <- glmmkin(y.trend ~ sex + time + age, data = pheno, kins = kins, id = "id",random.slope = "time", family = gaussian(link = "identity")) obj1.outfile.gds <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=gdsfile, outfile=obj1.outfile.gds) tmpout <- read.table(obj1.outfile.gds, header = T, as.is = T) expect_equal(obj1.gds[,-(10:13)], tmpout[,-(10:13)]) obj1.outfile.bgen <- tempfile() glmm.gei(null.obj=obj1, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen) tmpout <- read.table(obj1.outfile.bgen, header = T, as.is = T) expect_equal(obj1.bgen[,-(11:14)], tmpout[,-(11:14)]) unlink(c(obj1.outfile.gds, obj1.outfile.bgen)) obj1.outfile.gds <- tempfile() glmm.gei(null.obj=obj1, interaction="time", geno.file=gdsfile, outfile=obj1.outfile.gds) tmpout <- read.table(obj1.outfile.gds, header = T, as.is = T) expect_equal(obj1.time.gds[,-(10:19)], tmpout[,-(10:19)]) obj1.outfile.bgen <- tempfile() glmm.gei(null.obj=obj1, interaction="time", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj1.outfile.bgen) tmpout <- read.table(obj1.outfile.bgen, header = T, as.is = T) expect_equal(obj1.time.bgen[,-(11:20)], tmpout[,-(11:20)]) unlink(c(obj1.outfile.gds, obj1.outfile.bgen)) obj2 <- glmmkin(y.trend ~ sex + time + age, data = pheno, kins = NULL, id = "id",random.slope = "time", family = gaussian(link = "identity")) obj2.outfile.gds <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=gdsfile, outfile=obj2.outfile.gds) tmpout <- read.table(obj2.outfile.gds, header = T, as.is = T) expect_equal(obj2.gds[,-(10:13)], tmpout[,-(10:13)]) obj2.outfile.bgen <- tempfile() glmm.gei(null.obj=obj2, interaction="sex", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen) tmpout <- read.table(obj2.outfile.bgen, header = T, as.is = T) expect_equal(obj2.bgen[,-(11:14)], tmpout[,-(11:14)]) unlink(c(obj2.outfile.gds, obj2.outfile.bgen)) obj2.outfile.gds <- tempfile() glmm.gei(null.obj=obj2, interaction="time", geno.file=gdsfile, outfile=obj2.outfile.gds) tmpout <- read.table(obj2.outfile.gds, header = T, as.is = T) expect_equal(obj2.time.gds[,-(10:19)], tmpout[,-(10:19)]) obj2.outfile.bgen <- tempfile() glmm.gei(null.obj=obj2, interaction="time", geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj2.outfile.bgen) tmpout <- read.table(obj2.outfile.bgen, header = T, as.is = T) expect_equal(obj2.time.bgen[,-(11:20)], tmpout[,-(11:20)]) unlink(c(obj2.outfile.gds, obj2.outfile.bgen)) ### re-order kins with two interaction terms obj3 <- glmmkin(y.trend ~ sex + time +obese +age, data = pheno, kins = kins, id = "id",random.slope = "time", family = gaussian(link = "identity")) obj3.outfile.gds <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj3.outfile.gds) tmpout <- read.table(obj3.outfile.gds, header = T, as.is = T) expect_equal(obj3.gds[,-(10:17)], tmpout[,-(10:17)]) obj3.outfile.bgen <- tempfile() glmm.gei(null.obj=obj3, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen) tmpout <- read.table(obj3.outfile.bgen, header = T, as.is = T) expect_equal(obj3.bgen[,-(11:18)], tmpout[,-(11:18)]) unlink(c(obj3.outfile.gds, obj3.outfile.bgen)) obj3.outfile.gds <- tempfile() glmm.gei(null.obj=obj3, interaction=c("time","sex"), geno.file=gdsfile, outfile=obj3.outfile.gds) tmpout <- read.table(obj3.outfile.gds, header = T, as.is = T) expect_equal(obj3.time.gds[,-(10:29)], tmpout[,-(10:29)]) obj3.outfile.bgen <- tempfile() glmm.gei(null.obj=obj3, interaction=c("time","sex"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj3.outfile.bgen) tmpout <- read.table(obj3.outfile.bgen, header = T, as.is = T) expect_equal(obj3.time.bgen[,-(11:30)], tmpout[,-(11:30)]) unlink(c(obj3.outfile.gds, obj3.outfile.bgen)) obj4 <- glmmkin(y.trend ~ sex + time +obese +age, data = pheno, kins = NULL, id = "id",random.slope = "time", family = gaussian(link = "identity")) obj4.outfile.gds <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=gdsfile, outfile=obj4.outfile.gds) tmpout <- read.table(obj4.outfile.gds, header = T, as.is = T) expect_equal(obj4.gds[,-(10:17)], tmpout[,-(10:17)]) obj4.outfile.bgen <- tempfile() glmm.gei(null.obj=obj4, interaction=c("sex","obese"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen) tmpout <- read.table(obj4.outfile.bgen, header = T, as.is = T) expect_equal(obj4.bgen[,-(11:18)], tmpout[,-(11:18)]) unlink(c(obj4.outfile.gds, obj4.outfile.bgen)) obj4.outfile.gds <- tempfile() glmm.gei(null.obj=obj4, interaction=c("time","sex"), geno.file=gdsfile, outfile=obj4.outfile.gds) tmpout <- read.table(obj4.outfile.gds, header = T, as.is = T) expect_equal(obj4.time.gds[,-(10:29)], tmpout[,-(10:29)]) obj4.outfile.bgen <- tempfile() glmm.gei(null.obj=obj4, interaction=c("time","sex"), geno.file=bgenfile, bgen.samplefile = samplefile, outfile=obj4.outfile.bgen) tmpout <- read.table(obj4.outfile.bgen, header = T, as.is = T) expect_equal(obj4.time.bgen[,-(11:30)], tmpout[,-(11:30)]) unlink(c(obj4.outfile.gds, obj4.outfile.bgen)) })