A general approach to mixed effects modeling of residual variances in generalized linear mixed models

We propose a general Bayesian approach to heteroskedastic error modeling for generalized linear mixed models (GLMM) in which linked functions of conditional means and residual variances are specified as separate linear combinations of fixed and random effects. We focus on the linear mixed model (LMM) analysis of birth weight (BW) and the cumulative probit mixed model (CPMM) analysis of calving ease (CE). The deviance information criterion (DIC) was demonstrated to be useful in correctly choosing between homoskedastic and heteroskedastic error GLMM for both traits when data was generated according to a mixed model specification for both location parameters and residual variances. Heteroskedastic error LMM and CPMM were fitted, respectively, to BW and CE data on 8847 Italian Piemontese first parity dams in which residual variances were modeled as functions of fixed calf sex and random herd effects. The posterior mean residual variance for male calves was over 40% greater than that for female calves for both traits. Also, the posterior means of the standard deviation of the herd-specific variance ratios (relative to a unitary baseline) were estimated to be 0.60 ± 0.09 for BW and 0.74 ± 0.14 for CE. For both traits, the heteroskedastic error LMM and CPMM were chosen over their homoskedastic error counterparts based on DIC values.     

Hierarchical Bayesian modeling of heterogeneous cluster‐ and subject‐level associations between continuous and binary outcomes in dairy production

The augmentation of categorical outcomes with underlying Gaussian variables in bivariate generalized mixed effects models has facilitated the joint modeling of continuous and binary response variables. These models typically assume that random effects and residual effects (co)variances are homogeneous across all clusters and subjects, respectively. Motivated by conflicting evidence about the association between performance outcomes in dairy production systems, we consider the situation where these (co)variance parameters may themselves be functions of systematic and/or random effects. We present a hierarchical Bayesian extension of bivariate generalized linear models whereby functions of the (co)variance matrices are specified as linear combinations of fixed and random effects following a square‐root‐free Cholesky reparameterization that ensures necessary positive semidefinite constraints. We test the proposed model by simulation and apply it to the analysis of a dairy cattle data set in which the random herd‐level and residual cow‐level effects (co)variances between a continuous production trait and binary reproduction trait are modeled as functions of fixed management effects and random cluster effects.     

Mixed effects models with bivariate and univariate association parameters for longitudinal bivariate binary response data

When two binary responses are measured for each study subject across time, it may be of interest to model how the bivariate associations and marginal univariate risks involving the two responses change across time. To achieve such a goal, marginal models with bivariate log odds ratio and univariate logit components are extended to include random effects for all components. Specifically, separate normal random effects are specified on the log odds ratio scale for bivariate responses and on the logit scale for univariate responses. Assuming conditional independence given the random effects facilitates the modeling of bivariate associations across time with missing at random incomplete data. We fit the model to a dataset for which such structures are feasible: a longitudinal randomized trial of a cardiovascular educational program where the responses of interest are change in hypertension and hypercholestemia status. The proposed model is compared to a naive bivariate model that assumes independence between time points and univariate mixed effects logit models.     

Testing random effects in the linear mixed model using approximate bayes factors

Summary .  Deciding which predictor effects may vary across subjects is a difficult issue. Standard model selection criteria and test procedures are often inappropriate for comparing models with different numbers of random effects due to constraints on the parameter space of the variance components. Testing on the boundary of the parameter space changes the asymptotic distribution of some classical test statistics and causes problems in approximating Bayes factors. We propose a simple approach for testing random effects in the linear mixed model using Bayes factors. We scale each random effect to the residual variance and introduce a parameter that controls the relative contribution of each random effect free of the scale of the data. We integrate out the random effects and the variance components using closed-form solutions. The resulting integrals needed to calculate the Bayes factor are low-dimensional integrals lacking variance components and can be efficiently approximated with Laplace's method. We propose a default prior distribution on the parameter controlling the contribution of each random effect and conduct simulations to show that our method has good properties for model selection problems. Finally, we illustrate our methods on data from a clinical trial of patients with bipolar disorder and on data from an environmental study of water disinfection by-products and male reproductive outcomes.     

An alternative derivation of Harville's restricted log likelihood function for variance component estimation

Estimation of variance components in linear mixed models is important in clinical trial and longitudinal data analysis. It is also important in animal and plant breeding for accurately partitioning total phenotypic variance into genetic and environmental variances. Restricted maximum likelihood (REML) method is often preferred over the maximum likelihood (ML) method for variance component estimation because REML takes into account the lost degree of freedom resulting from estimating the fixed effects. The original restricted likelihood function involves a linear transformation of the original response variable (a collection of error contrasts). Harville's final form of the restricted likelihood function does not involve the transformation and thus is much easier to manipulate than the original restricted likelihood function. There are several different ways to show that the two forms of the restricted likelihood are equivalent. In this study, I present a much simpler way to derive Harville's restricted likelihood function. I first treat the fixed effects as random effects and call such a mixed model a pseudo random model (PDRM). I then construct a likelihood function for the PDRM. Finally, I let the variance of the pseudo random effects be infinity and show that the limit of the likelihood function of the PDRM is the restricted likelihood function.     

On the analysis of phylogenetically paired designs

As phylogenetically controlled experimental designs become increasingly common in ecology, the need arises for a standardized statistical treatment of these datasets. Phylogenetically paired designs circumvent the need for resolved phylogenies and have been used to compare species groups, particularly in the areas of invasion biology and adaptation. Despite the widespread use of this approach, the statistical analysis of paired designs has not been critically evaluated. We propose a mixed model approach that includes random effects for pair and species. These random effects introduce a “two-layer” compound symmetry variance structure that captures both the correlations between observations on related species within a pair as well as the correlations between the repeated measurements within species. We conducted a simulation study to assess the effect of model misspecification on Type I and II error rates. We also provide an illustrative example with data containing taxonomically similar species and several outcome variables of interest. We found that a mixed model with species and pair as random effects performed better in these phylogenetically explicit simulations than two commonly used reference models (no or single random effect) by optimizing Type I error rates and power. The proposed mixed model produces acceptable Type I and II error rates despite the absence of a phylogenetic tree. This design can be generalized to a variety of datasets to analyze repeated measurements in clusters of related subjects/species.     

The variance of sample heterozygosity

The variance of sample heterozygosity, averaged over several loci, is studied in a variety of situations. The variance depends on the sampling implicit in the mating system as well as on that explicit in the loci scored and individuals sampled. There are also effects of allelic distributions over loci and of linkage or linkage disequilibrium between pairs of loci. Results are obtained for populations in drift and mutation balance, for infinite populations undergoing mixed self and random mating, and for finite monoecious populations with or without selfing. For unlinked loci in drift/mutation balance, variances appear to be lessened more by increasing the number of loci scored than by increasing the number of individuals sampled. For infinite populations under the mixed self and random mating system, however, the reverse is true. Methods for estimating the variance of sample heterozygosity are discussed, with attention being paid to unbalanced data where not all loci are scored in all individuals.     

Incorporating extra information in experimental design for bioassay

In a bioassay, under certain experimental circumstances, information on concentration (dose rate) and time to response for some subjects can be combined in a single analysis. An underlying logistic random variable is assumed and the resulting mixed- (continuous-quantal) response model is analyzed by likelihood methods. The estimation procedure for the mean and the variance is described, and expressions for asymptotic variances are obtained. A comparison of results from the mixed model and from the standard quantal-response model shows that there is a substantial reduction in the variance of the estimators for the mixed model. On the basis of the table of asymptotic variances, some design implications are discussed. An example from insect pheromone research is used to illustrate the main ideas.     

Longitudinal analysis of pre- and post-treatment measurements with equal baseline assumptions in randomized trials

For continuous variables of randomized controlled trials, recently, longitudinal analysis of pre- and posttreatment measurements as bivariate responses is one of analytical methods to compare two treatment groups. Under random allocation, means and variances of pretreatment measurements are expected to be equal between groups, but covariances and posttreatment variances are not. Under random allocation with unequal covariances and posttreatment variances, we compared asymptotic variances of the treatment effect estimators in three longitudinal models. The data-generating model has equal baseline means and variances, and unequal covariances and posttreatment variances. The model with equal baseline means and unequal variance–covariance matrices has a redundant parameter. In large sample sizes, these two models keep a nominal type I error rate and have high efficiency. The model with equal baseline means and equal variance–covariance matrices wrongly assumes equal covariances and posttreatment variances. Only under equal sample sizes, this model keeps a nominal type I error rate. This model has the same high efficiency with the data-generating model under equal sample sizes. In conclusion, longitudinal analysis with equal baseline means performed well in large sample sizes. We also compared asymptotic properties of longitudinal models with those of the analysis of covariance (ANCOVA) and t-test.     

Association Models for Clustered Data with Binary and Continuous Responses

Summary .  We consider analysis of clustered data with mixed bivariate responses, i.e., where each member of the cluster has a binary and a continuous outcome. We propose a new bivariate random effects model that induces associations among the binary outcomes within a cluster, among the continuous outcomes within a cluster, between a binary outcome and a continuous outcome from different subjects within a cluster, as well as the direct association between the binary and continuous outcomes within the same subject. For the ease of interpretations of the regression effects, the marginal model of the binary response probability integrated over the random effects preserves the logistic form and the marginal expectation of the continuous response preserves the linear form. We implement maximum likelihood estimation of our model parameters using standard software such as PROC NLMIXED of SAS . Our simulation study demonstrates the robustness of our method with respect to the misspecification of the regression model as well as the random effects model. We illustrate our methodology by analyzing a developmental toxicity study of ethylene glycol in mice.     

Fixed‐effect variance and the estimation of repeatabilities and heritabilities: issues and solutions

Linear mixed‐effects models are frequently used for estimating quantitative genetic parameters, including the heritability, as well as the repeatability, of traits. Heritability acts as a filter that determines how efficiently phenotypic selection translates into evolutionary change, whereas repeatability informs us about the individual consistency of phenotypic traits. As quantities of biological interest, it is important that the denominator, the phenotypic variance in both cases, reflects the amount of phenotypic variance in the relevant ecological setting. The current practice of quantifying heritabilities and repeatabilities from mixed‐effects models frequently deprives their denominator of variance explained by fixed effects (often leading to upward bias of heritabilities and repeatabilities), and it has been suggested to omit fixed effects when estimating heritabilities in particular. We advocate an alternative option of fitting models incorporating all relevant effects, while including the variance explained by fixed effects into the estimation of the phenotypic variance. The approach is easily implemented and allows optimizing the estimation of phenotypic variance, for example by the exclusion of variance arising from experimental design effects while still including all biologically relevant sources of variation. We address the estimation and interpretation of heritabilities in situations in which potential covariates are themselves heritable traits of the organism. Furthermore, we discuss complications that arise in generalized and nonlinear mixed models with fixed effects. In these cases, the variance parameters on the data scale depend on the location of the intercept and hence on the scaling of the fixed effects. Integration over the biologically relevant range of fixed effects offers a preferred solution in those situations.     

Multibreed analysis by splitting the breeding values

An equivalent model for multibreed variance covariance estimation is presented. It considers the additive case including or not the segregation variances. The model is based on splitting the additive genetic values in several independent parts depending on their genetic origin. For each part, it expresses the covariance between relatives as a partial numerator relationship matrix times the corresponding variance component. Estimation of fixed effects, random effects or variance components provided by the model are as simple as any model including several random factors. We present a small example describing the mixed model equations for genetic evaluations and two simulated examples to illustrate the Bayesian variance component estimation.     

Using a mixed effects model to estimate geographic variation in cancer rates

Commonly used methods for depicting geographic variation in cancer rates are based on rankings. They identify where the rates are high and low but do not indicate the magnitude of the rates nor their variability. Yet such measures of variability may be useful in suggesting which types of cancer warrant further analytic studies of localized risk factors. We consider a mixed effects model in which the logarithm of the mean Poisson rate is additive in fixed stratum effects (e.g., age effects) and in logarithms of random relative risk effects associated with geographic areas. These random effects are assumed to follow a gamma distribution with unit mean and variance 1/alpha, similar to Clayton and Kaldor (1987, Biometrics 43, 671-681). We present maximum likelihood and method-of-moments estimates with standard errors for inference on alpha -1/2, the relative risk standard deviation (RRSD). The moment estimates rely on only the first two moments of the Poisson and gamma distributions but have larger standard errors than the maximum likelihood estimates. We compare these estimates with other measures of variability. Several examples suggest that the RRSD estimates have advantages compared to other measures of variability.     

Ontogenetic patterns in heritable variation for body size: using random regression models in a wild ungulate population

Body size is an important determinant of fitness in many organisms. While size will typically change over the lifetime of an individual, heritable components of phenotypic variance may also show ontogenetic variation. We estimated genetic (additive and maternal) and environmental covariance structures for a size trait (June weight) measured over the first 5 years of life in a natural population of bighorn sheep Ovis canadensis. We also assessed the utility of random regression models for estimating these structures. Additive genetic variance was found for June weight, with heritability increasing over ontogeny because of declining environmental variance. This pattern, mirrored at the phenotypic level, likely reflects viability selection acting on early size traits. Maternal genetic effects were significant at ages 0 and 1, having important evolutionary implications for early weight, but declined with age being negligible by age 2. Strong positive genetic correlations between age-specific traits suggest that selection on June weight at any age will likely induce positively correlated responses across ontogeny. Random regression modeling yielded similar results to traditional methods. However, by facilitating more efficient data use where phenotypic sampling is incomplete, random regression should allow better estimation of genetic (co)variances for size and growth traits in natural populations.     

Interobserver errors in anthropometry

To present basic information on the interobserver precision and accuracy of 32 selected anthropometric measurement items, six observers measured each of 37 subjects once in two days. The data were analyzed by using ANOVA, and mean absolute bias, standard deviation of bias, and mean absolute bias in standard deviation unit were used as measures of bias. By comparing the results of the two days, the effects of the practice on measurement errors were also investigated. Variance was overestimated by more than 10% in five measurements. Interobserver error variance and random error variance were highly correlated with each other. Measures of the bias were significantly correlated with interobserver and especially with random error variances. The interobserver errors were drastically reduced on the second day in the measurement items in which the causes of the interobserver errors could be specified. It was speculated that even when the definitions of the landmarks and measurement items were clear, the ambiguity in the practical procedures in locating landmarks, applying instruments, and so on, permitted each observer to develop his or her own measurement technique, and it in turn caused interobserver errors. To minimize interobserver and random errors, the standardization of measurement technique should be extended to the details of the practical procedures.     

Initiation of somatic embryogenesis in white spruce (Picea glauca): genetic control,culture treatment effects,and implications for tree breeding

The degree of genetic control and the effects of cultural treatments on somatic embryogenesis (SE) in white spruce were investigated with material derived from six-parent diallel crosses, including reciprocals. Thirty zygotic embryos from both immature and mature cones of each family were cultured in media with either 2,4-D or Picloram immediately after the collection of cones and after 2 months of cold storage. There were significant differences in SE initiation between immature and mature explants, and fresh and cold-stored seeds, but there was no significant differences with culture media effect. Significant variances due to families and to family x treatment interactions were found. The mean percentage of explants that initiated SE in each family ranged from 3.3% to 54.6%, with an overall average of 30.5%. The partitioning of family variance revealed that 21.7% was due to general combining ability effects, 3.5% was due to maternal effects, and 5.5% was due to reciprocal effects, but that the specific combining ability (SCA) was negligible. Variance due to interactions of family x treatments collectively accounted for 32.6%, while the remaining 37.8% of variation was accounted for by random error. However, when comparing the responses obtained with the treatment combinations, the SE response for freshly excised immature embryo explants showed comparatively large SCA variance, whereas the SCA variance was very small in the other treatment combinations.     

Borrowing information across genes and experiments for improved error variance estimation in microarray data analysis

Statistical inference for microarray experiments usually involves the estimation of error variance for each gene. Because the sample size available for each gene is often low, the usual unbiased estimator of the error variance can be unreliable. Shrinkage methods, including empirical Bayes approaches that borrow information across genes to produce more stable estimates, have been developed in recent years. Because the same microarray platform is often used for at least several experiments to study similar biological systems, there is an opportunity to improve variance estimation further by borrowing information not only across genes but also across experiments. We propose a lognormal model for error variances that involves random gene effects and random experiment effects. Based on the model, we develop an empirical Bayes estimator of the error variance for each combination of gene and experiment and call this estimator BAGE because information is Borrowed Across Genes and Experiments. A permutation strategy is used to make inference about the differential expression status of each gene. Simulation studies with data generated from different probability models and real microarray data show that our method outperforms existing approaches.     

Quantitative Genetics of DROSOPHILA MELANOGASTER. I. Sexual Dimorphism in Genetic Parameters for Wing Traits

Sexual dimorphism in genetic parameters is examined for wing dimensions of Drosophila melanogaster. Data are fit to a quantitative genetic model where phenotypic variance is a linear function of additive genetic autosomal variance (common to both sexes), additive genetic X-linked variances distinct for each sex, variance due to common rearing environment of families, residual environmental variance, random error variance due to replication, and variance due to measurement error and developmental asymmetry (left vs. right sides). Polygenic dosage compensation and its effect on genetic variances and covariances between sexes is discussed. Variance estimates for wing length and other wing dimensions highly correlated with length support the hypothesis that the Drosophila system of dosage compensation will cause male X-linked genetic variance to be substantially larger than female X-linked variance. Results for various wing dimensions differ, suggesting that the level of dosage compensation may differ for different traits. Genetic correlations between sexes for the same trait are presented. Total additive genetic correlations are near unity for most wing traits; this indicates that selection in the same direction in both sexes would have a minor effect on changing the magnitude of difference between sexes. Additive X-linked correlations suggest some genotype x sex interactions for X-linked effects.     

Variance of Neutral Genetic Variances within and between Populations for a Quantitative Character

The variances of genetic variances within and between finite populations were systematically studied using a general multiple allele model with mutation in terms of identity by descent measures. We partitioned the genetic variances into components corresponding to genetic variances and covariances within and between loci. We also analyzed the sampling variance. Both transient and equilibrium results were derived exactly and the results can be used in diverse applications. For the genetic variance within populations, sigma 2 omega, the coefficient of variation can be very well approximated as [formula: see text] for a normal distribution of allelic effects, ignoring recurrent mutation in the absence of linkage, where m is the number of loci, N is the effective population size, theta 1(0) is the initial identity by descent measure of two genes within populations and t is the generation number. The first term is due to genic variance, the second due to linkage disequilibrium, and third due to sampling. In the short term, the variation is predominantly due to linkage disequilibrium and sampling; but in the long term it can be largely due to genic variance. At equilibrium with mutation [formula: see text] where u is the mutation rate. The genetic variance between populations is a parameter. Variance arises only among sample estimates due to finite sampling of populations and individuals. The coefficient of variation for sample gentic variance between populations, sigma 2b, can be generally approximated as [formula: see text] when the number of loci is large where S is the number of sampling populations.     

Random effects selection in linear mixed models

We address the important practical problem of how to select the random effects component in a linear mixed model. A hierarchical Bayesian model is used to identify any random effect with zero variance. The proposed approach reparameterizes the mixed model so that functions of the covariance parameters of the random effects distribution are incorporated as regression coefficients on standard normal latent variables. We allow random effects to effectively drop out of the model by choosing mixture priors with point mass at zero for the random effects variances. Due to the reparameterization, the model enjoys a conditionally linear structure that facilitates the use of normal conjugate priors. We demonstrate that posterior computation can proceed via a simple and efficient Markov chain Monte Carlo algorithm. The methods are illustrated using simulated data and real data from a study relating prenatal exposure to polychlorinated biphenyls and psychomotor development of children.