Search

When estimating a benchmark dose (BMD) from chemical toxicity experiments, model averaging is recommended by the National Institute for Occupational Safety and Health, World Health Organization and European Food Safety Authority. Though numerous studies exist for Model Average BMD estimation using dichotomous responses, fewer studies investigate it for BMD estimation using continuous response. In this setting, model averaging a BMD poses additional problems as the assumed distribution is essential to many BMD definitions, and distributional uncertainty is underestimated when one error distribution is chosen a priori. As model averaging combines full models, there is no reason one cannot include multiple error distributions. Consequently, we define a continuous model averaging approach over distributional models and show that it is superior to single distribution model averaging. To show the superiority of the approach, we apply the method to simulated and experimental response data.

Individual-level trait diversity has been identified as an essential component of trait diversity (TD), influencing community assembly and structure. Traditionally, one employs trait diversity indices to measure facets of individual-level trait diversity (divergence, richness and evenness). However, the application of species-level trait diversity indices to individual-level traits data and their implications have not been adequately studied. Thus, we examined the possible challenges of using four commonly used multi-trait TD indices: Rao's quadratic entropy (Rao), functional dispersion (FDis), functional evenness (FEve) and functional richness (FRic); two indices primarily developed to measure individual-level trait diversity: trait evenness distribution (TED-for evenness) and trait onion peeling (TOP-for richnness); and a modified version of TED (TEDM-for evenness). Additionally, we considered an index that integrates both evenness and richness by generalizing ordinary Hill indices for traits (coined HIT). We measured individual-level trait diversity with these indices using simulated traits data and experimental data from a growth experiment with cyanobacteria. Comparing the observed trends from the indices with the expected trends, we observed that only the trait divergence indices (FDis and Rao) produced the expected trends in the simulation scenarios and experimental data. TED and TEDM are not robust against the number of individuals used, and FEve is not sensitive to some changes in the location of individuals in the trait space. Also, TOP proved to be a discontinuous function dependent on the number of individuals, and FRic did not produce the anticipated trend when changes in the trait space did not affect the edges of the trait space. HIT did produce the anticipated changes, but it was only reliable when many individuals were sampled. In summary, applying these individual-level trait diversity indices to quantify anything except trait divergence may lead to misinterpretation of the original situation of trait distribution in the trait space if their specific properties are not adequately considered.

The European Commission requested scientific and technical assistance in the preparation of a EU-wide baseline survey protocol for a European Union (EU) coordinated monitoring programme on the prevalence of methicillin-resistant

Flow cytometry is often employed in ecology to measure traits and population size of bacteria and phytoplankton. This technique allows measuring millions of particles in a relatively small amount of time. However, distinguishing between different populations is not a straightforward task. Gating is a process in the identification of particles measured in flow cytometry. Gates can either be created manually using known characteristics of these particles, or by using automated clustering techniques. Available automated techniques implemented in statistical packages for flow cytometry are primarily developed for medicinal applications, while only two exist for phytoplankton. cyanoFilter is an R package built to identify phytoplankton populations from flow cytometry data. The package also integrates gating functions from two other automated algorithms. It also provides a gating accuracy test function that can be used to determine the accuracy of a desired gating function if monoculture flowcytometry data is available. The central algorithm in the package exploits observed pigmentation and granularity of phytoplankton cells. We demonstrate how its performance depends on strain similarity, using a model system of six cyanobacteria strains. Using the same system, we compare the performance of the central gating function in the package to similar functions in other packages.

Biomarkers play a key role in the monitoring of disease progression. The time taken for an individual to reach a biomarker exceeding or lower than a meaningful threshold is often of interest. Due to the inherent variability of biomarkers, persistence criteria are sometimes included in the definitions of progression, such that only two consecutive measurements above or below the relevant threshold signal that "true" progression has occurred. In previous work, a novel approach was developed, which allowed estimation of the time to threshold using the parameters from a linear mixed model where the residual variance was assumed to be pure measurement error. In this paper, we extend this methodology so that serial correlation can be accommodated. Assuming that the Markov property holds and applying the chain rule of probabilities, we found that the probability of progression at each timepoint can be expressed simply as the product of conditional probabilities. The methodology is applied to a cohort of HIV positive individuals, where the time to reach a CD4 count threshold is estimated. The second application we present is based on a study on abdominal aortic aneurysms, where the time taken for an individual to reach a diameter exceeding 55 mm is studied. We observed that erroneously ignoring the residual correlation when it is strong may result in substantial overestimation of the time to threshold. The estimated probability of the biomarker reaching a threshold of interest, expected time to threshold, and confidence intervals are presented for selected patients in both applications.

Modelling travel behaviour has always been a major research area in transportation analysis. After the second World War, due to the rapid increase in car ownership and car use in Western Europe and the United States, several models have been developed by transportation planners. In the fifties and sixties, travel was assumed to be the result of four subsequent decisions that were modelled: trip generation, trip distribution, mode choice and the assignment of trips to the road network (Ruiter & Ben-Akiva, 1978). These original tripbased models have been extended to ensuing tour-based models (Daly et al., 1983) and activity-based models (Pendyala et al., 1995; Ben-Akiva & Bowman, 1998; Kitamura & Fujii, 1998; Arentze & Timmermans, 2000; Bhat et al., 2004). In tour-based models, trips are explicitly connected in tours, i.e. chains that start and end at the same home or work base. This is carried out by introducing spatial constraints, hereby dealing with the lack of spatial interrelationship which was so apparent in the traditional four-step trip-based model. In activity-based models, travel demand is derived from the activities that individuals and households need or wish to perform. Decisions with respect to travel are driven by a collection of activities that form an activity diary. Travel should therefore be modelled within the context of the entire agenda, or as a component of the activity scheduling decision. In this way, the relationship between travel and non-travel aspects is taken into account. The reason why people undertake trips is one of the key aspects to be modelled in an activity-based model. However, every working transportation model still exists of at least these original four components of trip generation, distribution, mode choice and assignment. In order to fully understand the structure of a traditional transportation model, we need to elaborate on it some more. As shown in Figure 1, trip generation encompasses both the modelling of production (P) and attraction (A) of trips for a certain region (zone). Production is mainly being modelled at the level of the household, incorporating household characteristics (income, car ownership, household composition, ...), features of the zone (land price, degree of urbanization) and accessibility of the zone, whereas attraction is modelled at zone level

In reliability theory, diagnostic accuracy, and clinical trials, the quantity [Formula: see text], also known as the Probabilistic Index (PI), is a common treatment effect measure when comparing two groups of observations. The quantity [Formula: see text], a linear transformation of PI known as the net benefit, has also been advocated as an intuitively appealing treatment effect measure. Parametric estimation of PI has received a lot of attention in the past 40 years, with the formulation of the Uniformly Minimum-Variance Unbiased Estimator (UMVUE) for many distributions. However, the non-parametric Mann–Whitney estimator of the PI is also known to be UMVUE in some situations. To understand this seeming contradiction, in this paper a systematic comparison is performed between the non-parametric estimator for the PI and parametric UMVUE estimators in various settings. We show that the Mann–Whitney estimator is always an unbiased estimator of the PI with univariate, completely observed data, while the parametric UMVUE is not when the distribution is misspecified. Additionally, the Mann–Whitney estimator is the UMVUE when observations belong to an unrestricted family. When observations come from a more restrictive family of distributions, the loss in efficiency for the non-parametric estimator is limited in realistic clinical scenarios. In conclusion, the Mann–Whitney estimator is simple to use and is a reliable estimator for the PI and net benefit in realistic clinical scenarios.

Censoring due to a limit of detection or limit of quantification happens quite often in many medical studies. Conventional approaches to deal with censoring when analyzing these data include, for example, the substitution method and the complete case (CC) analysis. More recently, maximum likelihood estimation (MLE) has been increasingly used. While the CC analysis and the substitution method usually lead to biased estimates, the MLE approach appears to perform well in many situations. This article proposes an MLE approach to estimate the association between two measurements in the presence of censoring in one or both quantities. The central idea is to use a copula function to join the marginal distributions of the two measurements. In various simulation studies, we show that our approach outperforms existing conventional methods (CC and substitution analyses). In addition, rank-based measures of global association such as Kendall's tau or Spearman's rho can be studied, hence, attention is not only confined to Pearson's product-moment correlation coefficient capturing solely linear association. We have shown in our simulations that our approach is robust to misspecification of the copula function or marginal distributions given a small association. Furthermore, we propose a straightforward MLE method to fit a (multiple) linear regression model in the presence of censoring in a covariate or both the covariate and the response. Given the marginal distribution of the censored covariate, our method outperforms conventional approaches. We also compare and discuss the performance of our method with multiple imputation and missing indicator model approaches.

In a developing country, it is very crucial to know where the HIV/AIDS epidemic is much more prevalent and where direct interventions are needed, especially when managing limited and scarce resources. We therefore examine the spatial distribution of HIV in Mozambique, and also assess how the epidemic evolved over a six-year period (2009-2015), with respect to potential risk factors among adolescents and young adults. We used data from the 2009 and 2015 Mozambique AIDS indicator surveys. The data were analysed jointly, by extending the work of Muleia et al. (2020) to allow for different bivariate spatial smoothing functions for both surveys. The results showed considerable spatial variation. From 2009 to 2015, the probability to be HIV positive reduced by 43.6% for young women. The results also showed dependence of the probability for HIV infection on sociodemographic factors. The findings herein will help health officials design efficient target interventions.

A major outbreak of the Ebola virus occurred in 2014 in Sierra Leone. We investigate the spatial heterogeneity of the outbreak among districts in Sierra Leone. The stochastic discrete-time susceptible-exposed-infectious-removed (SEIR) model is used, allowing for probabilistic movements from one compartment to another. Our model accounts for heterogeneity among districts by making use of a hierarchical approach. The transmission rates are considered time-varying. It is investigated whether or not incubation period, infectious period and transmission rates are different among districts. Estimation is done using the Bayesian formalism. The posterior estimates of the effective reproductive number were substantially different across the districts, with pronounced variability in districts with few cases of Ebola. The posterior estimates of the reproductive number at the district level varied between below 1.0 and 4.5, whereas at nationwide level it varied between below 1.0 and 2.5. The posterior estimate of the effective reproductive number reached a value below 1.0 around December. In some districts, the effective reproductive number pointed out for the persistence of the outbreak or for a likely resurgence of new cases of Ebola virus disease (EVD). The posterior estimates have shown to be highly sensitive to prior elicitation, mainly the incubation period and infectious period. The present study received financial support from Universidade Eduardo Mondlane The authors would like to express their gratitude to the editor and to the two anonymous reviewers for their valuable constructive comments and suggestions, which considerably improved the presentation of the paper.

Proposals to update the harmonised monitoring and reporting of antimicrobial resistance (AMR) from a public health perspective in Salmonella, Campylobacter coli, Campylobacter jejuni, Escherichia coli, Enterococcus faecalis, Enterococcus faecium and methicillin‐resistant Staphylococcus aureus (MRSA) from food‐producing animals and derived meat in the EU are presented in this report, accounting for recent trends in AMR, data collection needs and new scientific developments. Phenotypic monitoring of AMR in bacterial isolates, using microdilution methods for testing susceptibility and interpreting resistance using epidemiological cut‐off values is reinforced, including further characterisation of those isolates of E. coli and Salmonella showing resistance to extended‐spectrum cephalosporins and carbapenems, as well as the specific monitoring of ESBL/AmpC/carbapenemase‐producing E. coli. Combinations of bacterial species, food‐producing animals and meat, as well as antimicrobial panels have been reviewed and adapted, where deemed necessary. Considering differing sample sizes, numerical simulations have been performed to evaluate the related statistical power available for assessing occurrence and temporal trends in resistance, with a predetermined accuracy, to support the choice of harmonised sample size. Randomised sampling procedures, based on a generic proportionate stratified sampling process, have been reviewed and reinforced. Proposals to improve the harmonisation of monitoring of prevalence, genetic diversity and AMR in MRSA are presented. It is suggested to complement routine monitoring with specific cross‐sectional surveys on MRSA in pigs and on AMR in bacteria from seafood and the environment. Whole genome sequencing (WGS) of isolates obtained from the specific monitoring of ESBL/AmpC/carbapenemase‐producing E. coli is strongly advocated to be implemented, on a voluntary basis, over the validity period of the next legislation, with possible mandatory implementation by the end of the period; the gene sequences encoding for ESBL/AmpC/carbapenemases being reported to EFSA. Harmonised protocols for WGS analysis/interpretation and external quality assurance programmes are planned to be provided by the EU‐Reference Laboratory on AMR.

Protection and safety authorities recommend the use of model averaging to determine the benchmark dose approach as a scientifically more advanced method compared with the no-observed-adverse-effect-level approach for obtaining a reference point and deriving health-based guidance values. Model averaging however highly depends on the set of candidate dose-response models and such a set should be rich enough to ensure that a well-fitting model is included. The currently applied set of candidate models for continuous endpoints is typically limited to two models, the exponential and Hill model, and differs completely from the richer set of candidate models currently used for binary endpoints. The objective of this article is to propose a general and wide framework of dose response models, which can be applied both to continuous and binary endpoints and covers the current models for both type of endpoints. In combination with the bootstrap, this framework offers a unified approach to benchmark dose estimation. The methodology is illustrated using two data sets, one with a continuous and another with a binary endpoint.

Systematic assessment of QOL and care needs was applied in two gastroenterology departments to support “Cancer Care for the Whole Patient.” Patients with digestive cancer were asked to complete the Cancer Rehabilitation Evaluation System-Short Form (CARES-SF) at the start of treatment and 3 months later. Both times CARES data were processed, and summary reports on the retained insights were sent to the reference nurse for use in further follow-up of the patient. Patients’ and reference nurse’s experiences with the systematic CARES-assessment were explored with several survey questions and semi-structured interviews, respectively. The mean age of the 51 participants was 63 years (SD11.17), 52.9% was male. With the CARES-SF, a large variety of problems and care needs was detected. Problems most frequently experienced, and most burdensome for QOL are a mix of physical complaints, side effects from treatment, practical, relational, and psychosocial difficulties. Only for a limited number of experienced problems a desire for extra help was expressed. All patients positively evaluate the timing and frequency of the CARES-assessment. The majority believes that this assessment could contribute to the discussion of problems and needs with healthcare professionals, to get more tailored care. Reference nurses experienced the intervention as an opportunity to systematically explore patients’ well-being in a comprehensive way, leading to detection and discussion of specific problems or needs in greater depth, and more efficient involvement of different disciplines in care. Systematic QOL and needs assessment with the CARES-SF in oncology can contribute to more patient-centeredness and efficiency of care.

Frailty models have been developed to quantify both heterogeneity as well as association in multivariate time-to-event data. In recent years, numerous shared and correlated frailty models have been proposed in the survival literature allowing for different association structures and frailty distributions. A bivariate correlated gamma frailty model with an additive decomposition of the frailty variables into a sum of independent gamma components was introduced before. Although this model has a very convenient closed-form representation for the bivariate survival function, the correlation among event- or subject-specific frailties is bounded above which becomes a severe limitation when the values of the two frailty variances differ substantially. In this article, we review existing correlated gamma frailty models and propose novel ones based on bivariate gamma frailty distributions. Such models are found to be useful for the analysis of bivariate survival time data regardless of the censoring type involved. The frailty methodology was applied to right-censored and left-truncated Danish twins mortality data and serological survey current status data on varicella zoster virus and parvovirus B19 infections in Belgium. From our analyses, it has been shown that fitting more flexible correlated gamma frailty models in terms of the imposed association and correlation structure outperforms existing frailty models including the one with an additive decomposition.

Bivariate binary response data appear in many applications. Interest goes most often to a parameterization of the joint probabilities in terms of the marginal success probabilities in combination with a measure for association, most often being the odds ratio. Using, for example, the bivariate Dale model, these parameters can be modelled as function of covariates. But the odds ratio and other measures for association are not always measuring the (joint) characteristic of interest. Agreement, concordance, and synchrony are in general facets of the joint distribution distinct from association, and the odds ratio as in the bivariate Dale model can be replaced by such an alternative measure. Here, we focus on the so-called conditional synchrony measure. But, as indicated by several authors, such a switch of parameter might lead to a parameterization that does not always lead to a permissible joint bivariate distribution. In this contribution, we propose a new parameterization in which the marginal success probabilities are replaced by other conditional probabilities as well. The new parameters, one homogeneity parameter and two synchrony/discordance parameters, guarantee that the joint distribution is always permissible. Moreover, having a very natural interpretation, they are of interest on their own. The applicability and interpretation of the new parameterization is shown for three interesting settings: quantifying HIV serodiscordance among couples in Mozambique, concordance in the infection status of two related viruses, and the diagnostic performance of an index test in the field of major depression disorders. The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This study was only possible, thanks to the financial support of the Flemish Interuniversity Council (VLIR-UOS) in collaboration with Eduardo Mondlane University (UEM) through the DESAFIO Program. NH acknowledges funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (grant no. 682540 TransMID) and the Special Research Fund of Hasselt University.

The authors acknowledge support from IAP research Network P7/06 of the Belgian Government (Belgian Science Policy). The research leading to these results has also received funding from the European Seventh Framework programme FP7 2007-2013 under grant agreement no. 602552. We gratefully acknowledge support from the IWT-SBO ExaScience grant.

Data in medical sciences often have a hierarchical structure with lower level units (e.g. children) nested in higher level units (e.g. departments). Several specific but frequently studied settings, mainly in longitudinal and family research, involve a large number of units that tend to be quite small, with units containing only one element referred to as singletons. Regardless of sparseness, hierarchical data should be analyzed with appropriate methodology such as, for example linear-mixed models. Using a simulation study, based on the structure of a data example on Ceftriaxone consumption in hospitalized children, we assess the impact of an increasing proportion of singletons (0-95%), in data with a low, medium, or high intracluster correlation, on the stability of linear-mixed models parameter estimates, confidence interval coverage and F test performance. Some techniques that are frequently used in the presence of singletons include ignoring clustering, dropping the singletons from the analysis and grouping the singletons into an artificial unit. We show that both the fixed and random effects estimates and their standard errors are stable in the presence of an increasing proportion of singletons. We demonstrate that ignoring clustering and dropping singletons should be avoided as they come with biased standard error estimates. Grouping the singletons into an artificial unit might be considered, although the linear-mixed model performs better even when the proportion of singletons is high. We conclude that the linear-mixed model is stable in the presence of singletons when both lower and higher level sample sizes are fixed. In this setting, the use of remedial measures, such as ignoring clustering and grouping or removing singletons, should be dissuaded.

Objective: To construct a clinical prediction rule for coronary artery disease (CAD) presenting with chest pain in primary care. Study Design and Setting: Meta-Analysis using 3,099 patients from five studies. To identify candidate predictors, we used random forest trees, multiple imputation of missing values, and logistic regression within individual studies. To generate a prediction rule on the pooled data, we applied a regression model that took account of the differing standard data sets collected by the five studies. Results: The most parsimonious rule included six equally weighted predictors: age >= 55 (males) or >= 65 (females) (+1); attending physician suspected a serious diagnosis (+1); history of CAD (+1); pain brought on by exertion (+1); pain feels like "pressure" (+1); pain reproducible by palpation (-1). CAD was considered absent if the prediction score is = 2, it was 43.0% (95% CI: 35.8-50.4%). Conclusions: Clinical prediction rules are a key strategy for individualizing care. Large data sets based on electronic health records from diverse sites create opportunities for improving their internal and external validity. Our patient-level meta-analysis from five primary care sites should improve external validity. Our strategy for addressing site-to-site systematic variation in missing data should improve internal validity. Using principles derived from decision theory, we also discuss the problem of setting the cutoff prediction score for taking action. (C) 2016 Elsevier Inc. All rights reserved. This study was funded by Federal Ministry of Education and Research, Germany (BMBF grant no. FKZ 01GK0920).

The current statistical procedures implemented in statistical software packages for pooling of diagnostic test accuracy data include HSROC regression (Rutter and Gatsonis 2001) and the bivariate random-effects meta-analysis model (BRMA; Reitsma et al. 2005; Arends et al. 2008; Chu and Cole 2006; Riley et al. 2007b). However, these models do not report the overall mean but rather the mean for a central study with random-effect equal to zero and have difficulties estimating the correlation between sensitivity and specificity when the number of studies in the meta-analysis is small and/or when the between-study variance is relatively large (Riley et al. 2007a). This tutorial on advanced statistical methods for meta-analysis of diagnostic accuracy studies discusses and demonstrates Bayesian modeling using the R package CopulaDTA (Nyaga 2017) to fit different models to obtain the meta-analytic parameter estimates. The focus is on the joint modeling of sensitivity and specificity using a copula based bivariate beta distribution. Essentially, we extend the work of Nikoloulopoulos (2015) by: (i) presenting the Bayesian approach which offers the flexibility and ability to perform complex statistical modeling even with small data sets and (ii) including covariate information, and (iii) providing an easy to use code. The statistical methods are illustrated by re-analyzing data of two published meta-analyses. Modeling sensitivity and specificity using the bivariate beta distribution provides marginal as well as study-specific parameter estimates as opposed to using the bivariate normal distribution (e.g., in BRMA) which only yields study-specific parameter estimates. Moreover, copula based models offer greater flexibility in modeling different correlation structures in contrast to the normal distribution which allows for only one correlation structure.

There are several generalized linear mixed models to combine direct and indirect evidence on several diagnostic tests from related but independent diagnostic studies simultaneously also known as network meta-analysis. The popularity of these models is due to the attractive features of the normal distribution and the availability of statistical software to obtain parameter estimates. However, modeling the latent sensitivity and specificity using the normal distribution after transformation is neither natural nor computationally convenient. In this article, we develop a meta-analytic model based on the bivariate beta distribution, allowing to obtain improved and direct estimates for the global sensitivities and specificities of all tests involved, and taking into account simultaneously the intrinsic correlation between sensitivity and specificity and the overdispersion due to repeated measures. Using the beta distribution in regression has the following advantages, that the probabilities are modeled in their proper scale rather than a monotonic transform of the probabilities. Secondly, the model is flexible as it allows for asymmetry often present in the distribution of bounded variables such as proportions, which is the case with sparse data common in meta-analysis. Thirdly, the model provides parameters with direct meaningful interpretation since further integration is not necessary to obtain the meta-analytic estimates.

Network meta-analysis (NMA) allow combining efficacy information from multiple comparisons from trials assessing different therapeutic interventions for a given disease and to estimate unobserved comparisons from a network of observed comparisons. Applying NMA on diagnostic accuracy studies is a statistical challenge given the inherent correlation of sensitivity and specificity. A conceptually simple and novel hierarchical arm-based (AB) model which expresses the logit transformed sensitivity and specificity as sum of fixed effects for test, correlated study-effects and a random error associated with various tests evaluated in given study is proposed. We apply the model to previously published meta-analyses assessing the accuracy of diverse cytological and molecular tests used to triage women with minor cervical lesions to detect cervical precancer and the results compared with those from the contrast-based (CB) model which expresses the linear predictor as a contrast to a comparator test. The proposed AB model is more appealing than the CB model in that it yields the marginal means which are easily interpreted and makes use of all available data and easily accommodates more general variance-covariance matrix structures., Comment: 27 Pages, 6 figures

In longitudinal studies of biomarkers, an outcome of interest is the time at which a biomarker reaches a particular threshold. The CD4 count is a widely used marker of human immunodeficiency virus progression. Because of the inherent variability of this marker, a single CD4 count below a relevant threshold should be interpreted with caution. Several studies have applied persistence criteria, designating the outcome as the time to the occurrence of two consecutive measurements less than the threshold. In this paper, we propose a method to estimate the time to attainment of two consecutive CD4 counts less than a meaningful threshold, which takes into account the patient-specific trajectory and measurement error. An expression for the expected time to threshold is presented, which is a function of the fixed effects, random effects and residual variance. We present an application to human immunodeficiency virus-positive individuals from a seroprevalent cohort in Durban, South Africa. Two thresholds are examined, and 95% bootstrap confidence intervals are presented for the estimated time to threshold. Sensitivity analysis revealed that results are robust to truncation of the series and variation in the number of visits considered for most patients. Caution should be exercised when interpreting the estimated times for patients who exhibit very slow rates of decline and patients who have less than three measurements. We also discuss the relevance of the methodology to the study of other diseases and present such applications. We demonstrate that the method proposed is computationally efficient and offers more flexibility than existing frameworks. Copyright © 2016 John Wiley & Sons, Ltd.

In a linear multilevel model, significance of all fixed effects can be determined using F tests under maximum likelihood (ML) or restricted maximum likelihood (REML). In this paper, we demonstrate that in the presence of primary unit sparseness, the performance of the F test under both REML and ML is rather poor. Using simulations based on the structure of a data example on ceftriaxone consumption in hospitalized children, we studied variability, type I error rate and power in scenarios with a varying number of secondary units within the primary units. In general, the variability in the estimates for the effect of the primary unit decreased as the number of secondary units increased. In the presence of singletons (i.e., only one secondary unit within a primary unit), REML consistently outperformed ML, although even under REML the performance of the F test was found inadequate. When modeling the primary unit as a random effect, the power was lower while the type I error rate was unstable. The options of dropping, regrouping, or splitting the singletons could solve either the problem of a high type I error rate or a low power, while worsening the other. The permutation test appeared to be a valid alternative as it outperformed the F test, especially under REML. We conclude that in the presence of singletons, one should be careful in using the F test to determine the significance of the fixed effects, and propose the permutation test (under REML) as an alternative.