Joint p-value from SPA or score test?

[pottj]

Hi,

Many thanks for developing this really nice and efficient tool!

My question is about the joint p-value. The way I understood the documentation, this should be the harmonic mean of betapval and taupval, obtained from either the saddle point approximation (usespa=true, default) or score test (usespa=false). The harmonic mean should always be between those two values. However, in your example data the jointpval is sometimes smaller, and not equal to the harmonic mean of the two p-values.

I reran your example code with usespa=false, and obtain different betapval and taupval but the exactly same jointpval, which now are the harmonic mean of the score betapval and score taupval (see code and tables below).

Did I misunderstand the jointpval documentation, and it is supposed to be always the joint score test pvalue?

Thank you for your help!

# as in your example (basic usage) with usespa=true (default)
nullpath = "trajgwas.null.txt"
datadir = normpath(joinpath(dirname(pathof(TrajGWAS)), "../data/"))
pvalpath = "trajgwas.pval.txt"
trajgwas(@formula(y ~ 1 + sex + onMeds),
               @formula(y ~ 1),
               @formula(y ~ 1 + sex + onMeds),
               :id,
               datadir * "trajgwas_plinkex.csv",
               datadir * "hapmap3",
               pvalfile = pvalpath,
               nullfile = nullpath)
pretty_table(first(CSV.read("trajgwas.pval2.txt", DataFrame), 8), tf = `tf_markdown)

chr	pos	snpid	allele1	allele2	maf	hwepval	betapval	betadir	taupval	taudir	jointpval
1	554484	rs10458597	0	2	0.0	1.0	1.0	0	1.0	0	1.0
1	758311	rs12562034	1	2	0.0776398	0.409876	0.38939	-1	0.922822	-1	0.547682
1	967643	rs2710875	1	2	0.324074	4.07625e-7	5.83856e-6	1	3.35408e-6	-1	4.93598e-7
1	1168108	rs11260566	1	2	0.191589	0.128568	0.000850889	1	0.0559055	1	0.00101827
1	1375074	rs1312568	1	2	0.441358	2.5376e-19	0.000171683	-1	1.84811e-13	1	2.07091e-15
1	1588771	rs35154105	0	2	0.0	1.0	1.0	0	1.0	0	1.0
1	1789051	rs16824508	1	2	0.00462963	0.933278	0.295035	1	0.304109	1	0.299504
1	1990452	rs2678939	1	2	0.453704	5.07696e-11	1.27871e-7	1	7.73809e-9	-1	3.3986e-9

# now with usespa=false
pvalpath = "trajgwas.pval2.txt"
trajgwas(@formula(y ~ 1 + sex + onMeds),
               @formula(y ~ 1),
               @formula(y ~ 1 + sex + onMeds),
               :id,
               datadir * "trajgwas_plinkex.csv",
               datadir * "hapmap3",
               pvalfile = pvalpath,
               nullfile = nullpath,usespa=false)
pretty_table(first(CSV.read("trajgwas.pval2.txt", DataFrame), 8), tf = `tf_markdown)

chr	pos	snpid	allele1	allele2	maf	hwepval	betapval	betadir	taupval	taudir	jointpval
1	554484	rs10458597	0	2	0.0	1.0	1.0	0	1.0	0	1.0
1	758311	rs12562034	1	2	0.0776398	0.409876	0.38939	-1	0.922822	1	0.547682
1	967643	rs2710875	1	2	0.324074	4.07625e-7	2.75408e-6	1	2.71092e-7	-1	4.93598e-7
1	1168108	rs11260566	1	2	0.191589	0.128568	0.000513813	1	0.0559055	-1	0.00101827
1	1375074	rs1312568	1	2	0.441358	2.5376e-19	0.000244186	-1	1.03545e-15	1	2.07091e-15
1	1588771	rs35154105	0	2	0.0	1.0	1.0	0	1.0	0	1.0
1	1789051	rs16824508	1	2	0.00462963	0.933278	0.295035	1	0.304109	-1	0.299504
1	1990452	rs2678939	1	2	0.453704	5.07696e-11	2.16318e-7	1	1.71276e-9	-1	3.3986e-9