Fundamental to the study of the inheritance is the partitioning of the total phenotypic variation into genetic and environmental components. Using twin studies, the phenotypic variance-covariance matrix can be parameterized to include an additive genetic effect, shared and non-shared environmental effects. The ratio of the genetic variance component to the total phenotypic variance is the proportion of genetically controlled variation and is termed as the ‘narrow-sense heritability’. With the advent of large-scale genomic data, the genetic correlation among population-based and classically unrelated individuals can be used to estimate the total variance explained by the genome-wide SNPs (SNP-based heritability). In this talk, we will investigate the potential reasons behind the difference between the narrow-sense heritability and the SNP-based heritability estimates, often termed as the `missing heritability'. Over-estimation or incorrect estimation of heritability in twin studies has been suggested as one potential cause of the missing heritability. We will propose a robust and unified framework for estimating heritability in twin studies using second-order generalized estimating equations (GEE2). Specifically, we show that two commonly used models for estimating heritability, the normal ACE model (NACE) and Falconer’s method, can both be fit within this unified GEE2 framework, which additionally provides robust standard errors. Importantly, although the traditional Falconer’s method cannot directly adjust for covariates, the corresponding GEE2 version can incorporate both mean and variance-level covariate effects. Given non-normal data, we show that the GEE2 models attain significantly better coverage of the true heritability compared to the traditional NACE and Falconer’s methods. Finally, we demonstrate an important scenario where the popular NACE model produces substantially biased estimates of heritability while Falconer’s method remains unbiased.