Multivariate twin studies are one of the most important tools to assess diseases inheritance as well as to study their genetic and environmental interrelationship. The multivariate analysis of twin data is in general based on structural equation modelling or linear mixed models. Both approaches model the covariance matrix of a multivariate Gaussian distribution to take into account the genetic and environmental covariance induced by the twin design. In spite of flexible for Gaussian data, such approaches are unsuitable for analysing binary, binomial, count and asymmetric continuous traits. In this paper, we propose a flexible statistical modelling framework for analysing multivariate Gaussian and non-Gaussian twin data. Similar to the generalized linear models, link and variance functions are employed to take into account the non-normality, while the covariance structure induced by the twin design is modelled by means of a linear covariance model. The proposed model class can deal with binary, binomial, count, symmetric and asymmetric continuous traits as well as combination of them in a unified framework. Furthermore, from the marginal specification of our models emerge natural extensions of popular indices such as the bivariate heritability, genetic, environmental and phenotypic correlations to non-Gaussian data. The models are fitted by using an efficient two steps Newton scoring algorithm based on quasi-likelihood and Pearson estimating functions, using only second-moment assumptions. We illustrate the flexibility of the proposed models through simulation studies and data analyses.