Analysis of unbalanced mixed model data: A case study comparison of ANOVA versus REML/GLS

Major transition has occurred in recent years in statistical methods for analysis of linear mixed model data from analysis of variance (ANOVA) to likelihood-based methods. Prior to the early 1990s, most applications used some version of analysis of variance because computer software was either not available or not easy to use for likelihood-based methods. ANOVA is based on ordinary least squares computations, with adoptions for mixed models. Computer programs for such methodology were plagued with technical problems of estimability, weighting, and handling missing data. Likelihood-based methods mainly use a combination of residual maximum likelihood (REML) estimation of covariance parameters and generalized least squares (GLS) estimation of mean parameters. Software for REML/GLS methods became readily available early in the 1990s, but the methodology still is not universally embraced. Although many of the computational inadequacies have been overcome, conceptual problems remain. Also, technical problems with REML/GLS have emerged, such as the need for adjustments for effects due to estimating covariance parameters. This article attempts to identify the major problems with ANOVA, describe the problems which remain with REML/GLS, and discuss new problems with REML/GLS.     

毛乌素沙地1960—2018年气候变化特征及影响因子分析

为了解毛乌素沙地近年气候变化特征,以期为区域沙漠化防治和灾害预警提供参考,基于毛乌素沙地及其周围邻近地区10个气象站点1960—2018年的逐月气温和降水资料,采用线性倾向估计、Anusplin空间插值、Mann-Kendall突变检验、Morlet小波分析、R/S分析和Pearson相关分析等方法对毛乌素沙地近59 a的气候变化特征、未来气候变化趋势和主要影响因素进行了探究。结果表明:(1)毛乌素沙地近59 a的年平均气温呈显著上升趋势,增温速率达到0.35℃/10 a。其约在1996年发生突变,主要存在4 a左右的短变化周期,且增温速率随地区海拔的升高而加快。(2)毛乌素沙地近59 a的年降水变化速率为2.88 mm/10 a,夏季降水增速较快,但变化趋势均不显著且无明显的降水突变时间点和显著的振荡周期。(3)未来毛乌素沙地的年、季节平均气温和秋季降水均呈现持续性较强的上升趋势。总体而言,毛乌素沙地的气候呈现暖湿化趋向,厄尔尼诺—南方涛动(ENSO)事件对其冬季气温影响较大,北极涛动(AO)和北大西洋涛动(NAO)对其冬夏季节降水影响较显著。     

Estimating the local mean for Bayesian maximum entropy by generalized least squares and maximum likelihood, and an application to the spatial analysis of a censored soil variable

The Bayesian maximum entropy (BME) method is a valuable tool, with rigorous theoretical underpinnings, with which to predict with soft (imprecise) data. The methodology uses a general knowledge base to derive a joint prior distribution of the data and the prediction by the criterion of maximum entropy; the hard (precise) and soft data are then processed using this prior distribution to yield a posterior distribution that provides the BME prediction. The general knowledge base commonly consists of the mean and covariance functions, which may be extracted from the data. The common method for extracting the mean function from the data is a generalized least squares (GLS) approach. However, when the soft data take the form of intervals of plausible values, this method can result in errors in the BME predictions. This paper suggests a maximum likelihood (ML) method for fitting the local mean. The two methods are compared in terms of their predictions, firstly on simulated random fields and then on a case study to predict the depth of soil using some censored data. The results show that the ML method can result in more accurate BME predictions; the degree of improvement over the GLS method depends on the parameters of the spatial covariance model.     

Model‐based analysis using REML for inference from systematically sampled data on soil

The general linear model encompasses statistical methods such as regression and analysis of variance (anova ) which are commonly used by soil scientists. The standard ordinary least squares (OLS) method for estimating the parameters of the general linear model is a design‐based method that requires that the data have been collected according to an appropriate randomized sample design. Soil data are often obtained by systematic sampling on transects or grids, so OLS methods are not appropriate. Parameters of the general linear model can be estimated from systematically sampled data by model‐based methods. Parameters of a model of the covariance structure of the error are estimated, then used to estimate the remaining parameters of the model with known variance. Residual maximum likelihood (REML) is the best way to estimate the variance parameters since it is unbiased. We present the REML solution to this problem. We then demonstrate how REML can be used to estimate parameters for regression and anova ‐type models using data from two systematic surveys of soil. We compare an efficient, gradient‐based implementation of REML (ASReml) with an implementation that uses simulated annealing. In general the results were very similar; where they differed the error covariance model had a spherical variogram function which can have local optima in its likelihood function. The simulated annealing results were better than the gradient method in this case because simulated annealing is good at escaping local optima.     

Nonparametric Estimation of Spatial and Space-Time Covariance Function

Covariance structure modeling plays a key role in the spatial data analysis. Various parametric models have been developed to accommodate the idiosyncratic features of a given dataset. However, the parametric models may impose unjustified restrictions to the covariance structure and the procedure of choosing a specific model is often ad hoc. To avoid the choice of parametric forms, we propose a nonparametric covariance estimator for the spatial data, as well as its extension to the spatio-temporal data based on the class of space-time covariance models developed by Gneiting (J. Am. Stat. Assoc. 97:590–600, 2002). Our estimator is obtained via a nonparametric approximation of completely monotone functions. It is easy to implement and our simulation shows it outperforms the parametric models when there is no clear information on model specification. Two real datasets are analyzed to illustrate our approach and provide further comparison between the nonparametric estimator and parametric models.     

综合BME和BNN法的农田土壤水分与养分分布空间插值

掌握农田土壤水分和养分的空间分布特征是实现农田土壤精确管理及实施精确农业的重要依据。以有限的采样信息为基础,通过多种空间分析理论的融合,形成优势互补的综合方法,对提高土壤变量空间分布模拟和绘图精度具有重要意义。该文将贝叶斯最大熵法(Bayesian maximum entropy,BME)和贝叶斯人工神经网络方法(Bayesian neural networks,BNN)结合形成一种空间插值新方法,即用BNN法表达估值的不确定性,并将其结果融入现代地质统计学BME法中,用融入BNN法结果的BME法(Bayesian maximum entropy method combined with Bayesian neural networks,BMENN)模拟土壤变量的空间分布。以江苏省扬州市区北部某田块的土壤水分、有机质、全氮、碱解氮、速效钾和速效磷6种土壤特性的采样数据为例,运用交叉验证法,将BMENN法对土壤变量的估值精度与BNN法、普通克立格法(ordinary Kriging,OK)进行了比较。结果表明:与OK法和BNN法相比,BMENN法将估计方差(mean squared error,MSE)缩小2.26%~23.54%,具有最小的估计方差和接近于0的平均绝对误差(mean error,ME);BMENN法的估计值与实测值相关系数更大(r=0.62~0.89),具有更高的相关程度;MSE的组成分析表明,BMENN法再现变量波动程度和波动大小的能力更强;从模拟的空间分布图来看,BMENN法绘制的空间分布图更连续,"牛眼"较少,更符合土壤变量的地学规律。BMENN法对于利用有限数据信息提高土壤变量空间分布模拟精度具有重要意义,并可为土壤管理、精准农业的实施以及区域环境规划等提供科学依据。     

Linking species,traits and habitat characteristics of Collembola at European scale

Although much work has been done on factors which influence the patterning of species and species trait assemblages in a variety of groups such as plants, vertebrates and invertebrates, few studies have been realized at a broad geographic scale. We analyzed patterns of relationships between species, species trait distribution/assembly, and environmental variables from the west of Europe to Slovakia, Poland and Sweden. We created a database by compiling traits and occurrence data of European collembolan species, using literature and personal field studies embracing a large range of environmental gradients (vertical stratification, habitat closure, humus form, soil acidity and moisture, temperature, rainfall, altitude) over which Collembola are supposed to be distributed. Occurrences of the 58 best-documented species, environmental variables and species traits allowed us to (i) show which environmental variables impact the distribution of the 58 species at broad scale and (ii) document to what extent environmental variables and species trait assemblages are related and which trends could be found in trait/environment relationships. The impact of vertical stratification, habitat closure, humus form, soil acidity, soil moisture, temperature, and to a lesser extent rainfall and altitude on species distribution, firstly revealed by indirect gradient analysis (correspondence analysis, CA), was further shown to be significant by direct gradient analysis (canonical correspondence analysis, CCA). RLQ analyses were performed to find linear combination of variables of table R (environmental variables) and linear combinations of the variables of table Q (species traits) of maximum covariance weighted by species occurrence data contained in table L. RLQ followed by permutation tests showed that all tested environmental variables apparently contributed significantly to the assemblages of the twelve species traits studied. A convergence was observed between traits related to vertical stratification and those related to habitat closure/aperture. Well-developed locomotory organs (furcula, legs), presence of sensorial organs sensitive to air movements and light (e.g. trichobothria and eye spots), spherical body, large body size, pigmentation (UV protection and signaling) and sexual reproduction largely occur in epigeic and open habitats, while most of woodland and edaphic habitats are characterized by short locomotory appendages, small body size, high number of defense organs (pseudocelli), presence of post-antennal organs and parthenogenesis. Climate and especially temperature exert an effect on the assemblage of traits that are mostly present above-ground and in open habitats. The contribution of combinations of some environmental variables to the occurrence of each species trait was tested by linear, logistic or multinomial regression (Generalized Linear Models). Vertical stratification, followed by temperature, played a dominant role in the variation of the twelve studied traits. Relationships between traits and environment tested here shows that it is possible to use some traits as proxies to identify potential ecological preferences or tolerances of invertebrate species. However, a significant part of species distribution remained unexplained, probably partly because some traits, like ecophysiological ones, or traits involved in biotic interactions (e.g. competition) were unavailable. The present work is thus a first step towards the creation of models predicting changes in collembolan communities. Further studies are required to inform ecophysiological traits, in order to complete such models. Moreover the niche width of species will have to be determined.     

Exploring scale‐dependent correlation of soil properties by nested sampling

The value of nested sampling for exploring the spatial structure of univariate variation of the soil has been demonstrated in several studies and applied to practical problems. This paper shows how the method can be extended to the multivariate case. While the extension is simple in theory, in practice the direct estimation of covariance components by equating mean‐square matrices with their expectation will often lead to estimates that are not positive semidefinite. This paper discusses solutions to this problem for balanced and unbalanced sample designs. In the balanced case there is a residual maximum likelihood (REML) estimator that will find estimates of covariance components that maximize an overall likelihood on the condition that all components are positive semidefinite (p.s.d.). This is possible because the condition is met if the differences of successive mean‐square matrices are positive semidefinite, and this constraint can be incorporated into an algorithm. This does not hold for unbalanced designs. In this paper the problem was solved for unbalanced designs by scaling covariance components that were not p.s.d. to the nearest p.s.d. matrix according to a Euclidean distance. These methods were applied to data from three surveys, two with balanced and one with unbalanced sampling. Different patterns of scale‐dependence of the correlation of soil properties were found. For example, at Ginninderra Experimental Station in Australia the soil water content and bulk density were correlated significantly, with the correlation increasing with distance to 56 m, but at longer distances the properties were not significantly correlated. By contrast, the pH of the soil and the available P content showed correlation that increased with distance. The implications of these results for planning more detailed sampling, both for prediction and for investigation of processes, are discussed.     

Neural Modeling of Greenhouse Gas Emission from Agricultural Sector in European Union Member Countries

The present paper discusses a novel methodology based on neural network to determine agriculture emission model simulations. Methane and nitrous oxide are the key pollutions among greenhouse gases being a major contribution to climate changes because of their high potential global impact. Using statistical clustering (k-means and Ward’s method), five meaningful clusters of countries with similar level of greenhouse gases emission were identified. Neural modeling using multi-layer perceptron networks was performed for countries placed in particular groups. The parameters that characterize the quality of a network are the predictive errors (mainly validation and test) and they are high (0.97–0.99). The use of sensitivity analysis allowed for identifying the variables that have a significant influence on the greenhouse gases emissions. The sensitivity analysis of the designed artificial neural network models shows a few dominant variables, affecting emissions with varied intensity: cattle and buffaloes, sheep and goat populations, afforestation as well as electricity consumption. The observed values were compared with those predicted by the models. The forecasted course of changes in the variable test is identical with the real data, which proves that the model highly matches to the observed data.     

Mapping the soil resources of a recent alluvial plain in Greece using fuzzy sets in a GIS environment

Collecting soil data is time-consuming and costly, often exceeding practical possibilities. A methodology for the delineation of soil mapping units in an alluvial plain of Western Peloponnese, Greece, was investigated. A detailed soil survey of an area of 300 ha was used to obtain the basic soil data for evaluating the performance of the proposed methodology. The methodology consists of the following steps: (a) data collection from borings and representative soil profiles, (b) definition of the soil mapping units in the study area, (c) determination of the range of the diagnostic variables for each mapping unit from field observations and statistical analysis of the analytical data from representative soil profiles, (d) determination of the class of each diagnostic variable by observation at a network of boring points, (e) subjective assignment of numerical values to soil variables at the bore points, (f) estimation of the values of each soil variable at the points of a regular grid using the interpolation methods kriging and inverse squared distance, (g) application of the fuzzy set theory to the interpolated data and the production of thematic fuzzy maps, and (h) validation of the results through a number of independent test borings. The results obtained show that the proposed methodology can produce soil maps of recent alluvial plains with acceptable accuracy and cost.     

Kriging a soil variable with a simple nonstationary variance model

Kriging is a standard tool in the environmental sciences for spatial prediction from limited sample data, subject to the assumption of intrinsic stationarity, made about the underlying spatially correlated random function. It is generally well understood how the assumption of stationarity in the mean can be relaxed within the linear mixed model framework, using residual maximum likelihood to estimate variance parameters for the random effects. The Best Linear Unbiased Predictor (BLUP) is equivalent to the kriging predictor in these circumstances. However, nonstationarity in the variance is a harder problem to solve. Stationarity assumptions are necessary if the spatial covariance of a random process is to be estimated from the single realization which nature provides. However, they are not always plausible for variables arising from processes in complex landscapes across contrasting topography and geology.     

Respiration of three Belgian crops: Partitioning of total ecosystem respiration in its heterotrophic, above- and below-ground autotrophic components

An experimental system combining an eddy covariance system, a micrometeorological station and soil chambers placed in planted areas and in root exclusion zones was installed during three successive years in a production crop managed in a traditional way at the Lonzée experimental site (Belgium). Measurements were made successively on seed potato, winter wheat and sugar beet. The general objectives of the study were, first to evaluate the relative contributions to total ecosystem respiration (TER) of heterotrophic, above ground autotrophic and below ground autotrophic respiration over a succession of three agricultural crops (seed potato, winter wheat and sugar beet) cultivated on successive years at the same location and, secondly, to identify the driving variables of these contributions.Results showed that, during the observation periods, TER was dominated by autotrophic respiration (AR) (60-90%) and that AR was dominated by its above ground component (60-80%). HR was found to increase with temperature and to be independent of Gross Primary Production (GPP), whereas AR was driven by GPP and was mostly independent of temperature. The AR response to GPP was specific to the crop: not only AR intensity but also AR distribution between its above- (ARa) and below- (ARb) ground components were found to differ from one crop to another and, in the winter wheat, from one development stage to another. Generally, ARb contribution to AR was found larger when carbon allocation towards roots was more important.An uncertainty analysis was made and showed that the main sources of uncertainties on the estimates were the spatial variability for soil chamber measurements and uncertainties linked to the data gap filling method for eddy covariance measurements.     

Identifying, understanding, and describing spatial processes in agricultural landscapes — four case studies

To evaluate the quality of the ecosystem and for making resources and land management decisions landscapes have to be assessed quantitatively. For a better understanding of landscape processes and their characterization, the analysis of the inherent variability is a major factor. Four case studies in which problems associated with landscape analysis are discussed. Spatial processes remain a main focus, as their analysis provides information on the relation between relevant state variables in agricultural landscapes. Variogram analysis showed that mineral soil nitrogen (N_min) sampled in a field at different scales, domains, and times is an instationary spatial process. Spatial association of grain yield, soil index and remotely sensed vegetation index may not be identifiable from kriged contour maps as local coincidence may be obscured behind classified areas. Crop yield in subsequent years and remotely sensed information are not related if a unique response is assumed. An alternative data stratification procedure is described here for the identification of different response functions in agricultural ecosystems. Processes of crop yield and underlying variables are described in autoregressive state-space models. This technique incorporates both deterministic and stochastic relations between different variables and is based on relative changes in space.     

Mt. Brocken,a site for a cloud chemistry measurement programme in Central Europe

We present results from the Brocken Cloud Chemistry Measurement Project (BROCCMON) which started in 1991. Since 1992 the full programme is running, based on continuous measurements (e.g. trace gases, meteorology, liquid water content), cloud water sampling and analysis and intensive measurement campaigns. The observed high variability of cloud water composition we explain with cloud dynamic and microphysical behaviour of clouds and differences in the air mass characteristics. During the measurement period 1992–1994 we observed an increase in cloud water acidity (by a factor of 3) and we found photochemical conditions typically for summersmog situations. Our preliminary data also show that an understanding of tropospheric ozone balance would be incomplete without consideration of chemical processes within clouds. A long-term goal of our programme is to establish a cloud chemistry climatology which is representative for the region.     

绿洲荒漠过渡带涡动相关仪和大孔径闪烁仪的通量源区特征

以甘肃民勤绿洲荒漠过渡带涡动相关仪（EC）和大孔径闪烁仪（LAS）的同步观测数据为基础,借助通量足迹模型,分析生长季和非生长季通量源区的分布与变化特征,比较了源区在不同风向下,足迹权重比的变化对二者观测值之间差异的影响。结果表明：（1）研究区在非生长季盛行西北风,在生长季则盛行东风,EC和LAS源区的空间分布与盛行风向一致,二者在非生长季源区的面积均大于生长季。（2）大气稳定条件下EC和LAS的源区面积普遍大于不稳定条件下。随着源区贡献率的增加（由50%增至90%）,EC和LAS源区的重叠面积和足迹权重比也增加,并且非生长季的源区重叠程度优于生长季。在生长季和非生长季,EC和LAS源区都呈现不规则分布的特点。（3）相比于东西风向,生长季南北风向下的足迹权重比更高,二者所观测显热通量值更为接近,决定系数（R2）也更高。说明源区的足迹权重比在一定程度上可以解释EC和LAS在观测结果上的差异,这为地表通量的尺度扩展提供了方法学上的参考。     

The simultaneous analysis of mixed discrete and continuous outcomes using nonlinear threshold models

Mixed discrete and continuous outcomes are commonly measured on each experimental unit in dose-response studies in toxicology. The dose-response relationships for these outcomes often have dose thresholds and nonlinear patterns. In addition, the endpoints are typically correlated, and a statistical analysis that incorporates the association may result in improved precision. We propose an extension of the generalized estimating equation (GEE) methodology to simultaneously analyze binary, count, and continuous outcomes with nonlinear threshold models that incorporates the intra-subject correlation. The methodology uses a quasi-likelihood framework and a working correlation matrix, and is appropriate when the marginal expectation of each outcome is of primary interest and the correlation between endpoints is a nuisance parameter. Because the derivatives of threshold models are not continuous at each point of the parameter space, we describe the necessary modifications that result in asymptotically normal and consistent estimators. Using dose-response data from a neurotoxicity experiment, the methodology is illustrated by analyzing five outcomes of mixed type with nonlinear threshold models. In this example, the incorporation of the intra-subject correlation resulted in decreased standard errors for the threshold parameters.     

Inferring Upon Heterogeneous Associations in Dairy Cattle Performance Using a Bivariate Hierarchical Model

Multivariate hierarchical Bayesian models provide a flexible framework for comprehensive study of biological systems with more than one outcome. Recent methodological developments facilitate modeling of heterogeneous associations between outcomes by specifying a linear mixed model on (co)variances at different levels of the data structure. Motivated by previous evidence for heterogeneous correlations in animal agriculture, we apply the proposed hierarchical Bayesian models to study the nature of the correlations between key performance outcomes in dairy cattle production systems, namely milk yield and reproduction. That is, the association between these outcomes might depend upon various fixed and random effect sources of heterogeneity both at the individual cow (residual) level as well as the herd (cluster) level. We thus propose a sequential modeling approach based on the deviance information criterion to select relevant explanatory variables on both types of associations. Furthermore, we extend the proposed methodology to accommodate right-censored outcomes, as common for dairy reproduction data, and use it to analyze field data from the Michigan dairy industry. The nature of the associations between milk production and reproduction in dairy cattle was inferred to be strongly heterogeneous and driven by multiple farm management practices and herd attributes, as well as by random clustering effects, at both cow and herd levels, thereby suggesting potential between-herd and within-herd intervention strategies to optimize performance of dairy production systems. Supplementary materials are available online.     

寒冷灾害监测中的全天候地表温度反演方法研究

提出综合运用MODIS遥感资料与地面气象站资料实现全天候LST反演的技术方法,并通过应用示范验证其可行性,对实现全天候、精细化的寒冷灾害监测与评估具有重要意义。基于MODIS的第1、2、19、31和32五个波段,采用劈窗算法反演LST,计算过程相对简便而反演精度较高,缺点是无法克服云干扰;基于地面气象站实测温度资料,采用气候学方程与GIS技术反演空间精细化的LST分布,尽管其反演精度在晴空条件下比遥感方法差,但不受云量条件限制,可以实现全天候LST反演。     

Hyperspectral Remote Sensing of Total Phosphorus (TP) in Three Central Indiana Water Supply Reservoirs

The connection between nutrient input and algal blooms for inland water productivity is well known but not the spatial pattern of water nutrient loading and algae concentration. Remote sensing provides an effective tool to monitor nutrient abundances via the association with algae concentration. Twenty-one field campaigns have been conducted with samples collected under a diverse range of algal bloom conditions for three central Indiana drinking water bodies, e.g., Eagle Creek Reservoir (ECR), Geist Reservoir (GR), and Morse Reservoir (MR) in 2005, 2006, and 2008, which are strongly influenced anthropogenic activities. Total phosphorus (TP) was estimated through hyperspectral remote sensing due to its close association with chlorophyll a (Chl-a), total suspended matter, Secchi disk transparency (SDT), and turbidity. Correlation analysis was performed to determine sensitive spectral variables for TP, Chl-a, and SDT. A hybrid model combining genetic algorithms and partial least square (GA-PLS) was established for remote estimation of TP, Chl-a, and SDT with selected sensitive spectral variables. The result indicates that TP has close association with diagnostic spectral variables with R ² ranging from 0.55 to 0.72. However, GA-PLS has better performance with an average R ² of 0.87 for aggregated dataset. GA-PLS was applied to the airborne imaging data (AISA) to map spatial distribution of TP, Chl-a, and SDT for MR and GR. The eutrophic status was evaluated with Carlson trophic state index using TP, Chl-a, and SDT maps derived from AISA images. Mapping results indicated that most MR belongs to mesotrophic (48.6%) and eutrophic (32.7%), while the situation was more severe for GR with 57.8% belongs to eutrophic class, and more than 40% to hypereutrophic class due to the high turbidity resulting from dredging practices.