Hierarchical statistical techniques are necessary to draw reliable conclusions from analysis of isolated cardiomyocyte studies

1Changes in cellular Ca2+ handling are agreed to play a key role in the pathophysiology of heart failure (HF). The standard experimental model, established for half a century is left anterior descending coronary artery ligation in the rat. At the whole organ level, there is broad agreement between experimenters on the pattern of haemodynamic and echocardiographic variables.

2At the cellular level, however, there is more controversy. Changes in Ca2+ transient amplitude and kinetics have been widely reported to contribute to organ dysfunction., However, the 18 studies that report on changes in Ca2+ transient kinetics show considerable conflict (Table 1). On even smaller scales, Ca2+ spark measurements also show convincing changes, again in conflicting directions (see Supplementary material online, Table S1). This is not a peculiarity of the rat: mouse studies too show the same conflicts (see Supplementary material online, TablesS2 and S3).

3The accepted paradigm is that reduced Ca2+ transient amplitude contributes to reduced myocardial contractility. However, the cellular data do not consistently fit this. For example, Ca2+ transient amplitudes in HF cells are reported as significantly increased in more studies than decreased (Table 1). Moreover, the two studies that reported changes in SR Ca2+ content indicated opposite directions of change (Table 1). In other variables, studies conflict between showing convincing changes and convincing absence of change (Table 1).

4There is also no apparent time-dependent shift in variables following MI as would be expected for variables that change from a ‘compensated’ to ‘de-compensated’ cellular phenotype.

5Authors have proposed explanations including a compensated phase following MI, or differences in other elements of the myocyte (e.g. myofilaments) overriding the effects of enhanced Ca2+ transients. When one study found enhanced Ca2+ amplitudes simultaneously with reduced myocyte contractility, the proposed explanation was a reduction in myofilament sensitivity to Ca2+.

6Rarely considered is the possibility that our statistical methods may have left us open to changes appearing falsely statistically significant. The studies in the field generally used standard statistical tests designed to be valid for independent data points, because these tests are widely available and straightforward to implement.

7The challenge we face is that when we take n cells from each of m animals we do not truly have n × m independent data points. We have m clusters, each containing n data points. The data points in each cluster (i.e. the cells from a single animal) will tend to be more similar to each other than they are to points in the other clusters (Figure 1). There are well-established statistical techniques for handling data that shows such clustering. Whether using cluster based analysis produces different results from our field’s standard statistical tests, has never been tested.

8When analysing such a set of n × m data points, investigators tend to take one of two approaches. The more conservative approach is to calculate a mean for each of the m animals and treat these m mean values as the only data points. This approach is not popular because it treats the sample size as only m and therefore, reduces the ability to detect changes.

9The far more popular approach is to ignore the clustering and treat all n × m data points as though they are independent. The attraction of this approach is that when the standard statistical tests are performed, the large number of data points lead to the standard error being very small and therefore differences appear to be more statistically significant. In fact the term ‘pseudoreplication’ has been coined for this. Studies rarely specify which of these two approaches have been taken but the sample sizes reported generally appear consistent with the latter approach.

10A third approach, rarely used, is to recognize the clustered (hierarchical) structure of the data. And use statistical methods designed for this situation. Hierarchical statistical methods test for clustering and, if present, correct for this when performing significance testing.

11In this study we test the application of a hierarchical statistical approach to data from myocytes from the post-MI rat HF model. We examine whether (i) it makes a difference to the conclusions drawn, (ii) whether it is necessary, and (iii) whether it can be implemented conveniently by experimenters.

12The chronic rat HF model was produced surgically via LAD ligation with an 8 week delay before echocardiography, culling and biometric measurement as described previously. Briefly, rats were anesthetized with 2% isoflurane, intubated, and ventilated after preoperative buprenorphine (0.03 mg/kg SC) injection. Loss of righting reflex and pain responses confirmed dequate analgesia and anaesthesia before operating and tying off the LAD with 6-0 prolene. All animal surgical procedures and peri-operative management were carried out in accordance with the Guide for the Care and Use of Laboratory Animals published by the US National Institutes of Health (NIH Publication, 8th Edition, 2011) under assurance number A5634-01. Imperial College Ethical Review Committee authorized the project licence.

13Echocardiography and biometric measurements showed a reduced ejection fraction compared with control animals (40 ± 2.2% vs. 83 ± 2.1%, P < 0.0001), as well as a significant increase in heart weight: body weight ratio (3.28 ± 0.1 vs. 2.79 ± 0.1, P < 0.0001) confirming established HF (see Supplementary material online, Figure S1).

14Rats were culled by cervical dislocation and cells were isolated 8 weeks following ligation. Cells were loaded with Fura2-AM and imaged using ratiometric techniques as described. Ca2+ transients were recorded at 1 Hz in control (76 cells from 10 isolations) and HF cells (79 cells from 10 isolations) using the ratiometric dye Fura2-AM. The ratio of the isosbestic point for Fura2 (360 nm) to the unbound form (380 nm) was used as a measure of [Ca2+]i. Several variables were assessed using automated transient analysis in IonWizard (Ionoptix Inc.) including diastolic ratio, peak systolic ratio, transient amplitude (peak ratio/diastolic ratio), time to 50% peak, time to 50% decay (TD50), and tau of decay.

15In a separate set of experiments isolated cardiomyocytes were loaded with fluo4-AM as described previously and spontaneous Ca2+ sparks were assessed using confocal microscopy. All experiments were performed at 37 °C in NT solution containing 2 mM Ca2+. The control group consisted of myocytes isolated from age-matched controls. Data were viewed using Ionwizard and ImageJ (NIH) and analysed with a combination of custom-built macros and Sparkmaster. Amplitude and morphology of five Ca2+ transients per cell were averaged to give final result for that cell. Sparks were treated as separate data points.

16Data were analysed using IBM SPSS Statistics in two ways. Firstly, using our field’s common statistical approach of treating all data points as though they were independent. For each variable, an independent samples t-test was performed comparing all the cells from HF hearts against all the cells from control hearts.

17The second analysis used hierarchical statistical techniques. These techniques are designed for situations where values may be clustered, with a possibility of similarity within each cluster. We have designed cost-free methods which we have made available for the reader to use for similar analyses (see Discussion section).

18The amount of clustering can be quantified by the intraclass correlation coefficient (ICC). This varies between 0% for data showing no clustering, and 100% for intensely clustered data (i.e. each cluster contains multiple identical values, but different clusters have different values). This measure represents the proportion of total variability in an outcome measure that is attributable to the isolation of origin. Illustrative datasets with high, moderate, and low values are shown in Figure 2. The method of calculation is shown in Supplementary material online.

19If the data show no clustering, then the hierarchical model works effectively identically to the commonly used statistical test, treating all the data points as independent. If, on the other hand, the data are tightly grouped within each cluster, with relatively large separation between clusters (as might happen when each day’s isolation was internally homogeneous but differed from other days) then a relatively large correction would need to be applied. The hierarchical test quantifies the amount of clustering and applies the appropriate correction to the statistical significance test (Figure 3), further explanation of hierarchical testing is shown in Supplementary material online (see Supplementary methods and Figure S2).

20The hierarchical approach also provides a simple method of augmenting conventional sample size calculation to explain the fact that with m cells from each of n rats the true sample size is less than n × m but more than n. The correction required depends on the degree of clustering as quantified by the ICC. With low levels of clustering, where each cell’s behaviour is independent of the animal it has come from, effective sample size approaches n × m. With high degrees of clustering, where all cells from a single animal behave identically, effective sample size approaches n. The calculation of effective sample size is as follows:

21For n animals, each with m cells For example, with 10 animals and 20 cells from each, with an ICC of 50%, effective sample size is neither 200 nor 10, and is calculated as follows: This has important implications for the interpretation of conventional power calculations in the context of hierarchical data. Power calculations are used to work out the appropriate sample size for experiments given a working knowledge of likely effect size and variance obtained from pilot data. In the context of a hierarchical data structure, the sample size given by such power calculations should be considered to be the Effective Sample Size. A rearrangement of the above equation would give the true sample size required:

22We tested whether the hierarchical statistical model produced a significantly better fit to the data than the commonly used statistical approach, using the χ2 test of the change in -2 Log Likelihood (χ2-2LL). This is the recommended method for comparisons of this kind.

23Initial exploration of data showed that whilst variables pertaining to Ca2+ transients were symetrically distributed, spark morphological variables were skewed as has been shown previously. The logarithmic transformations of morphological characteristics were therefore used to improve this lack of symmetry (see Supplementary material online, Figure S3). These transformed variables are referred to as LogAmp, LogFWHM, LogFDHM for the logarithm of spark amplitude, full width at half maximum (FWHM), full duration at half maximum (FDHM).

24Ca2+ transient and spark measures were the dependent variables in the statistical models. The only independent variable was the presence or absence of HF. Each model used was assessed for validity by ensuring predicted values closely corresponded to those observed. Residuals were assessed for normality and symmetry. Estimated marginal means were used to assess significance between control and HF groups.

25Ca2+ transient variables were measured from 76 cells from 10 control rats and 79 cells from 10 HF rats. Ca2+ spark variables were measured in 344 sparks from 17 cells from 7 control rats and 352 sparks from 22 cells from 5 HF rats.

26For all Ca2+ transient morphology variables, the hierarchical model produced a statistically significantly better fit than the commonly used analysis method (Table 2). The magnitude of the intraclass correlation (degree of clustering) varied between 12 and 47% (Table 2, ICC column).

27For all the comparisons between HF and control, the commonly used analysis method reported a more statistically significant difference than the better-fitting hierarchical analysis. The mechanism of this was that the commonly used statistical method produced a smaller value for the standard error of the difference between HF and control (Table 2).

28Figure 4 shows the 95% confidence interval calculated, by the two approaches, for the difference in time to peak between HF and control. The commonly used analysis method gives a different conclusions because it has an artificially small 95% confidence interval, causing the data to meet the criteria for statistical significance.

29For sparks there is potential for an additional level of clustering, because the multiple sparks from each cell could be more similar to each other than the sparks from other cells of the same rat.

30The hierarchical analysis was carried out twice (Table 3): once allowing only for clustering at the animal level and once allowing for clustering at both the cell and the animal levels. We tested whether the incorporation of the additional (cell) level improved the fit as defined by the χ2-2LL. Where the fit was improved significantly, the additional level was incorporated into the analysis (methods shown in Supplementary appendix).

31The P-values were substantially different between the commonly used method and the hierarchical analyses. Again, the hierarchical analyses fitted better for each variable. The cause of the difference in P values was that the hierarchical models calculated larger values for the standard error.

32Even minor clustering made a substantial difference to the P value. For example, for logFWHM, even with only 8.5% clustering (as quantified by the ICC), the conventional approach which assumes no clustering reports a P value of 0.02, whereas the hierarchical model reports a P value of 0.78.

33LogAmp showed significant clustering both within each cell and within rats (Table 3, Figure 5A, B). Figure 5C shows the 95% confidence interval calculated, by the two approaches, for the difference in LogAmp between HF and control. The artificially small confidence interval produced by the common test makes the difference between the groups appear highly significant, but when the hierarchical data structure is taken into account, there is no significant difference.

34Looking at the best-fitting analysis method (Table 3, defined by χ2-2LL criterion), it emerges that none of the spark variables is in fact significantly different between HF and control. LogAmp and LogFDHM appeared to be different between HF and control using the commonly used analysis approach, but this was because of an artefactually small standard error through failing to account for clustering.

35Different studies of Ca2+ fluxes in HF cells have given opposite results for aspects of Ca2+handling in isolated cardiomyocytes. The usual interpretation of this is that the direction of the effect may vary with choice of HF model and species. Our present study, however, shown in Table 1 and Supplementary material online, TablesS1–3 that even within a single HF model, studies already show opposite direction of effects.

36Our study suggests that that these contradictory directions of effect could be a manifestation of the subtle but important defect in our conventional statistical analysis, whereby we assume that all the data points are independent with no significant clustering. In reality for many biological variables, repeated measurements (e.g. of sparks) from the same cell are more similar than those of other cells, and measurements from multiple cells of the same animal are in turn more similar than measurements from different animals. The consequence is when we have multiple measurements we have less meaningful data than the simple number of data points would suggest. This problem, the ability of large sample numbers to deceive us when we do not recognize clustering, has been previously identified in the field of neuroscience.,

37While this study focuses on the clustering of transient data from different myocytes isolated from the same rat, as well as the clustering of spark data within each cell, this should not be seen as a reason to limit the use of hierarchical statistics to these settings. We have found significant clustering in other aspects of isolated myocyte properties, such as cell size. We have also found that these techniques are essential for robust analysis of other models. For example in a recent study where multiple cardiac slices were taken from each dog, significant clustering of slice contractility data meant that data were most appropriately analysed using hierarchical techniques. It would be safest, therefore, for the investigator to assume that any data with a hierarchical structure should be analysed in this way unless the hierarchical tests show non-significant clustering. We have made such testing freely available (below).

38The transients of an individual cell are known to be almost identical and therefore, it is already conventional to recognize this clustering by treating the entirety of transient data from one cell as a single data point by averaging data from several transients.

39Sparks, on the other hand, vary a great deal from one to the next within the same cell. This stochastic behavior results from the subcellular physiology of ryanodine receptors, which are present in variable configurations with variable opening probabilities and dynamics. The convention in our field has therefore been to not average spark data into a single value per cell, but consider each spark as a separate data point.

40Our study indicates that the spark data show strong clustering at the cell level, which means that it is wrong to count the individual sparks for statistical purposes as though they were independent.

41Strikingly, although the degree of clustering (ICC) varied from one variable to another, it was always significant for all of the variables. Moreover, even small values of the clustering measure (ICC) turn out to be surprisingly important. For example, spark logFDHM, which has an ICC of just 8.5%, changed its result from being clearly significant with the common analysis (P = 0.023), to being not at all significant (P = 0.778) with the hierarchical analysis. Many variables had far higher degrees of clustering (Tables 2 and 3), which means they are correspondingly even more susceptible to false positive results from the common analysis.

42We were surprised to find that, after recognizing the clustering using hierarchical analysis, the only significant difference between HF and control is a higher Ca2+ transient amplitude. This is due to increased systolic Ca2+ without increase in diastolic Ca2+. If we had used standard t-tests alone the conclusion of this study would have been that Ca2+ transients in HF exhibited both higher amplitude and a longer time to peak and that these changes might result from sparks (the elemental subunits of transients) which were more prolonged and of greater amplitude, giving a mechanistic explanation for our finding of higher transient amplitude which is not supported by the data.

43Because hierarchical statistics are sometimes considered difficult to implement, some workers have tried the two alternative approaches.

44Aggregation is the approach of lumping multiple data points into one, typically by averaging. This is already done routinely for Ca2+ transients at the cell level but not the rat level. The disadvantage is that applying this approach universally is excessively conservative, and would require very large studies to test hypotheses that could be tested more frugally if hierarchical statistics were used. For example, the 696 sparks in our study would have been collapsed into just 12 averages, one for each rat, if we had indiscriminately aggregated.

45Results of aggregated statistical analysis in this study are shown in Supplementary material online, TablesS4, S5 and are similar to the outcomes of the hierarchical tests but with overestimation of P-values and standard errors.

46Disaggregation is the opposite approach, splitting the data so that each data point is considered to be independent. This is the commonest approach used in our field but has the scientific disadvantage that the statistics conducted assuming this independence give a falsely small standard error, falsely small confidence interval and therefore falsely small P value for the difference between groups.

47An example of the catastrophic effect of disaggregation is that in our data the hierarchical model shows only 2 of 10 comparisons between HF and control to be statistically significant; in contrast the commonly-used approach (disaggregation of data) shows five comparisons to be statistically significant: a 150% overstatement.

48A further challenge is that in the cellular laboratory there can be ‘good isolation days’ and ‘bad isolation days’. On good days large numbers of healthy cells are available for experimentation but on bad days smaller numbers of less healthy cells are available. If we aggregate, we give excessive weight to the likely unhealthier cells from the bad days. If we disaggregate, we give excessive weight to the likely healthier cells from the good days because they are more numerous.

49The social sciences have long known of this problem and developed standardized approaches to avoid wasting resources pursuing false leads.,, We describe a similar problem in studies of isolated cardiac myocytes but we also describe a solution and provide cost-free simple tools which can be used to produce robust statistical results for comparisons involving clustered data.

50There is a general focus on the P-value as the only important outcome of significance testing in the biological literature. This approach leads to an overreliance on the apparent binary outcome of P < 0.05 vs. P > 0.05. In addition, the P-value gives the reader no indication of the magnitude of the effect and therefore its biological (rather than statistical) significance. An emphasis on point estimates and their precision (e.g. mean and 95% confidence interval) can prevent these issues.

51If the hierarchical data structure is not considered, the 95% confidence intervals of an estimated mean are excessively narrow. Appropriate corrections can be made when performing hierarchical analysis. The magnitude of the required correction, as with other aspects of hierarchical testing, depends on the ICC. The estimated mean and 95% confidence intervals for each variable tested in this study, with and without hierarchical correction, are shown in Supplementary material online, Tables S6 and S7.

52There are commercial software packages such as IBM SPSS Statistics, SAS (SAS Institute), STATA (StataCorp) which can do these hierarchical statistics very well. Readers can download the scripts which we show in Supplementary material online, Supplementary methods section. However, some readers may not have a particular commercial software package available or may be working with collaborators who do not.

53We therefore present simple steps that any researcher can use without cost to perform these hierarchical statistics.

54Step 1. Download the relevant files from Supplementary material online into a single working directory. These should include: ratonly.Rmd, ratandcell.Rmd, Hierarchical Transient analysis with Rat-Level Clustering.xlsx, and Hierarchical Spark analysis Cell & Rat-Level Clustering.xlsx.

55Step 2. Modify the relevant Microsoft Excel file. If the data have the potential for clustering only at one level, for example, at the rat level, use the layout in Supplementary material online, Hierarchical Transient analysis with Rat-Level Clustering.xlsx. For data such as spark values, which has potential for an additional (e.g. cellular) level of clustering, use the layout in Hierarchical Spark analysis Cell & Rat-Level Clustering.xlsx.

56Cells should be numbered sequentially even if they come from different rats, that is, if cells 1–10 come from rat 1, then cells from rat 2 should be numbered from 11 upwards and so on.

57Column titles can be changed to reflect what your relevant hierarchical structure is—for example, dogs and heart slices or cell line and culture number.

58Step 3. Download and install R for Windows, Mac, or Linux from https://www.r-project.org/ (14 August 2017, date last accessed).

59Step 4. Download and install RStudio from https://www.rstudio.com/products/rstudio/download/ (14 August 2017, date last accessed)—this is a more user-friendly interface for R. Choose RStudio Desktop—Open Source License.

60Step 5. Run RStudio and install the two extra packages we will need (‘lmerTest’ for the statistical tests and ‘readxl’ for reading in the Excel spreadsheets). This can be done using the ‘Tools → Install packages’ option in RStudio, and typing ‘lmerTest, readxl’ (first letter of lmer is lower case ‘L’) in the ‘Packages’ field, or by typing the following commands in the ‘Console’ window (R Console) and pressing Enter after each command:install.packages(‘lmerTest’)

61install.packages(‘readxl’)

62Step 6. Open the code corresponding to the Excel file you are working with in RStudio by clicking File menu and Open File. Open ratonly.Rmd if working with clustering at one level and ratandcell.Rmd if working with clustering at two levels.

63Step 7. In the top right hand corner of the window in which the code has opened there is a Run command. Click this and then select Run All.

64Step 8. Once the code has run, which may take several minutes, scroll down to the bottom of the same window where three tables (for ratonly.Rmd) or five tables (for ratandcell.Rmd). These tables show the following: The first table shows the outcomes of the common test (P value and standard error for the difference), the degree of clustering (ICC), the equivalent outcomes for the hierarchical test, and whether the hierarchical test is a better fit for the data under the column ‘Superior fit’. If the P value in this column is < 0.05, there is sufficient clustering to make hierarchical techniques necessary. For ratandcell.Rmd, both cell-level hierarchical test (labelled group-level) and rat-and-cell level hierarchical test (labelled Parentgroup-group level) results are shown.The second table shows the point estimate, that is, the estimated group mean when taking into account the hierarchical data structure, as well as the standard error and the lower and upper confidence intervals for each variable for each condition (HF and control).The third table shows the pairwise comparisons. If there are only two conditions (e.g. HF and control) the P-values in this table are the same as the first table. This table becomes particularly relevant if there are >2 conditions since comparisons each pair of groups (A vs. B; B vs. C; A vs. C) may be desired. As with an ANOVA, these group-wise comparisons are only relevant if the overall test shown in the first table is significant.The fourth and fifth tables are only output by ratandcell.Rmd and are the equivalent of the second and third tables but for the rat-and-cell level hierarchical test.

65The first table shows the outcomes of the common test (P value and standard error for the difference), the degree of clustering (ICC), the equivalent outcomes for the hierarchical test, and whether the hierarchical test is a better fit for the data under the column ‘Superior fit’. If the P value in this column is < 0.05, there is sufficient clustering to make hierarchical techniques necessary. For ratandcell.Rmd, both cell-level hierarchical test (labelled group-level) and rat-and-cell level hierarchical test (labelled Parentgroup-group level) results are shown.

66The second table shows the point estimate, that is, the estimated group mean when taking into account the hierarchical data structure, as well as the standard error and the lower and upper confidence intervals for each variable for each condition (HF and control).

67The third table shows the pairwise comparisons. If there are only two conditions (e.g. HF and control) the P-values in this table are the same as the first table. This table becomes particularly relevant if there are >2 conditions since comparisons each pair of groups (A vs. B; B vs. C; A vs. C) may be desired. As with an ANOVA, these group-wise comparisons are only relevant if the overall test shown in the first table is significant.

68The fourth and fifth tables are only output by ratandcell.Rmd and are the equivalent of the second and third tables but for the rat-and-cell level hierarchical test.

69The scripts we provide are sufficient to perform the analyses in this paper. Once readers are comfortable doing this, they can extend the R script to allow more levels of clustering or to include other variables such as gender, animal weight, or cell size.

70Moreover, these scripts are equally suited outside the cardiomyocyte arena and even to clinical data.

71For users who do not have R or RStudio installed, but who wish to view the source code and sample output, we have supplied two supplementary HTML (ratonly.nb.html and ratandcell.nb.html) files which can be opened in any browser and serve as manuals to the analyses.

72Even though the hierarchical test describes the data much better (as shown by the χ2-2LL criterion), we should not assume that they are the ideal solution. It is merely a significant improvement for which we provide a straightforward method for implementation. Ideally researchers would engage with an expert statistician at an early phase of study preparation. Alternative models for performing hierarchical statistics are available including Generalized Estimating Equations (GEE), although we favour Mixed Models (with random coefficients) as this gives greater flexibility in the modelling process and is better at handling missing data.

73Data transformations are frequently necessary, such as log transformation for the spark data which have a marked positive skew. In addition, even following valid transformation the meaning of the original data must be questioned. For example, some of the skewed distribution relating to variables describing spark morphology relates to the confocal line-scanning technique which is likely to miss the center of each spark, thus artefactually increasing the number of smaller sparks.,,

74The conventional approach used in our field to analyse cardiomyocyte data has a marked tendency to overstate the statistical significance of differences. This could be why some studies shows significant results in one direction and other studies show equally convincing significant results in the other direction.

75We present simple steps that any researcher can implement to identify and allow for clustering, at the rat and cell level if need be, permitting them to use all their data points and yet obtain statistically valid results.

76Supplementary material is available at Cardiovascular Research online.