Parametric assumptions in population genetics analysis – including linearity, sources of population stratification and the gaussianity and additivity of errors – are often made, yet a principled argument for their (approximate) validity is not given. We present a unified statistical workflow, called TarGene, for targeted estimation of effect sizes, as well as two-point and higher-order epistatic interactions of genomic variants on polygenic traits, that dispenses with these unnecessary assumptions. Our approach is founded on Targeted Learning, a framework for estimation that integrates mathematical statistics, machine learning and causal inference to provide mathematical guarantees and realistic p-values. TarGene defines effect sizes of variants, as well as two-point and higher-order inter-actions amongst genomic variants on traits in a model-independent manner, thus avoiding all-too-common model-misspecification whilst taking advantage of a library of parametric and state-of-the-art non-parametric algorithms. TarGene data-adaptively incorporates confounders and sources of population stratification, accounts for population dependence structures and controls for multiple hypothesis testing by bounding any desired type I error rate. Extensive simulations demonstrate the necessity of this model-independent approach. We validate the effectiveness of our method by reproducing previously verified effect sizes on UK Biobank data, whilst simultaneously discovering non-linear effect sizes of additional allelic copies on trait or disease. To exemplify this, we demonstrate that for the FTO variant rs1421085 effect size on body mass index (BMI), the addition of one copy of the C allele is associated with $0.77$ kg/m$^2$ ($95%$ CI: $0.68 − 0.85$) increase, while the addition of the second C copy non-linearly adds $1.31$ kg/m$^2$ ($95%$ CI: $1.19 − 1.43$) to BMI. TarGene thus extends the reach of current genome-wide association studies by simultaneously (i) allowing for the classification of the types of SNPs and phenotypes for which such non-linearities occur, whilst (ii) data-adaptively incorporating complex non-linear relations between phenotype, genotype, and confounders, as well as (iii) accounting for strong population dependence such as island cohorts. The method provides a platform for comparative analyses across biobanks, or integration of multiple biobanks and heterogeneous populations to increase power, whilst controlling for population stratification and multiple hypothesis testing.