A discussion on two celebrated examples of application of linear regression

The story behind Hubble’s determination of the famous constant, which has been named after him, started over $100$ years before Hubble was born: Clotho’s spindle started spinning when Edward Pigott (1753-1825), an English astronomer, developed a fondness for stars of variable brightness. Pigott’s family moved to York in 1781, where Pigott met John Goodricke (1764-1786), a deaf-mute with unparalleled visual acuity and sensitivity. Working in a private observatory built by Pigott’s father, Pigott and Goodricke embarked on a project aiming at exploring the properties of the most promising variable stars, taken from Pigott’s precompiled list [12]. It did not take long before Goodricke’s skills and attention to detail paid off: one year into the project, he announced to the Royal Society the discovery of the periodical variation of the brightness of the star Algol ($\beta$ Persei) [13]. Apart from proposing a mechanism to account for the effect (and attributing it to the transits of Algol’s dimmer binary companion), Goodricke also examined the variations of brightness of several other stars, including $\eta$ Aquilae (A) and $\delta$ Cephei (A) [14], which became textbook examples of the Cepheid family of variable stars ³³3Cepheid variable stars pulsate radially, varying in both diameter and temperature [15]., featuring very regular ‘sharply-up-gradually-down’ light curves, dissimilar to the one obtained from Algol; evidently, the variation of brightness of these stars could not be attributed to stellar eclipses. Goodricke’s contributions to Astronomy did not go unnoticed: at the age of 21, he was elected a Fellow of the Royal Society, a few days before his sad passing.

Over the course of the next century, many more Cepheid variable stars were detected. With the invention of photography, the acquisition of accurate brightness data was boosted. The Harvard College Observatory became the hub of the developments in Observational Astronomy, and its new director in 1877, Edward Charles Pickering (1846-1919), launched an ambitious programme aiming at photographing and cataloguing all observable stars. Whether due to budget issues (e.g., see Ref. [16], p. 205) or to the poor quality of the work of Pickering’s male assistants [17], over $80$ College-educated women were recruited, to work at the Observatory as human computers, i.e., as skilled employees, comparing photographic plates, processing astronomical data, and contributing to Astronomy for the first time after Hypatia of Alexandria and Aglaonice of Thessaly.

The computers excelled in their assignments. For instance, Annie Jump Cannon (1863-1941), who was deaf after she turned twenty, conducted the main work in the development of the Harvard Classification Scheme of stars, the first attempt towards a categorisation of the stars in terms of their surface temperature and spectral type. In an article entitled ‘Woman making index of $100\,000$ stars for a catalogue’, the Danville Morning News reported on 10 February 1913 [18]: “When the work was new she could analyze at the rate of $1\,000$ stars in three years. Now she analyzes $5\,000$ stars in one month, $200$ stars an hour ⁴⁴4This is an erratum: evidently, the author of the article intended to write ‘day’, not ‘hour’.. On Jan. 1 she had examined about $65\,000$, which means that about two-fifths of the work is completed.”

Without doubt, Pickering’s research programme reached its climax with another computer, Henrietta Swan Leavitt (1868-1921). Leavitt’s methodical work was destined to make an indelible impact on Astronomy (and Cosmology). Having been given the task of studying variable stars in the Small and Large Magellanic Clouds (SMC and LMC, respectively), Leavitt catalogued $1\,777$ such objects between 1903 and 1908 [19]. The decisive moment came when she focussed her attention on a small subset of these objects, i.e., on the Cepheid variable stars in the SMC (see Table VI of Ref. [19]); on p. 107 of her report, she inconspicuously commented: “It is worthy of notice that in Table VI the brighter variables have the longer period.”

Following her intuition in subsequent years, Leavitt examined the relation between the period of the variation of the apparent magnitude (which is associated with the apparent brightness) of the Cepheid variable stars in the SMC and the extrema of the apparent magnitude for each such object. Her decision to restrict her analysis to the SMC turned out to be of pivotal importance, because the selected $25$ (i.e., the $16$ from Table VI of Ref. [19] and the $9$ which were added to her database between 1908 and 1913) celestial bodies were, in practical terms, equidistant from the Earth ⁵⁵5It is currently known that the SMC is a dwarf galaxy (with a typical diameter of about $5.78$ kpc), gravitationally attracted to the Milky Way, and separated from the Earth by $60.5(7.0)$ kpc [20]..

After plotting the period-brightness ⁶⁶6I shall use this short term when referring to the datapoints corresponding to the period and to the apparent magnitude of variable stars; other authors prefer the term ‘period-luminosity’. datapoints of the $25$ Cepheid variable stars on a linear-log plot (logarithmic scale on the $x$ axis, linear scale on the $y$ axis), Leavitt ⁷⁷7The fact that Pickering, rather than Leavitt, signed the 1912 report is baffling. The report starts with the sentence: “The following statement regarding the periods of $25$ variable stars in the Small Magellanic Cloud has been prepared by Miss Leavitt.” Although Pickering’s motive for signing Leavitt’s report is a mystery to me, I decided to include him in the author list of that report, yet as second author. obtained the “remarkable” result shown in Fig. 2 of Ref. [21]: two nearly parallel straight lines, modelling the apparent-magnitude maxima and minima of these stars, separated by about $1.2$ mag (difference between the corresponding intercepts); I shall use ‘mag’ to denote the unit of apparent magnitude (which, of course, is not a physical unit). Two remarks from that report serve as a prelude to the significance of Leavitt’s research:

• •

  “…there is a simple relation between the brightness of the variables and their periods” and

• •

  “Since the variables are probably at nearly the same distance from the Earth, their periods are apparently associated with their actual emission of light, as determined by their mass, density, and surface brightness.”

The consequences of Leavitt’s findings were beyond reckoning. The measurement of the period of a Cepheid variable star, whichever neighbourhood of the Universe that specific celestial body occupies, could lead to an estimate for its apparent brightness $b_{\rm SMC}$ which that body would have if it had been confined within the SMC. By comparing the body’s measured apparent brightness $b$ to $b_{\rm SMC}$, one could obtain an estimate for the distance between the specific Cepheid variable star and the Earth, expressed (of course) in terms of the SMC-Earth distance. As the distances of all Cepheid variable stars would be expressed as multiples of the SMC-Earth distance, one would need to determine only one such distance (e.g., using the stellar-parallax method) in order that all distances become known. Within one year of Leavitt’s breakthrough, Ejnar Hertzsprung (1873-1967) provided a solution, however flawed ⁸⁸8On p. 204 of Hertzsprung’s article [22], one reads: “…so wird die entsprehende visuelle Sterngröße gleich $13^{m}.0$. Diese Überlegung führt also zu einer Parallaxe $p$ der kleinen Magellanschen Wolke, welche durch $5\log p=-7.3-13.0=-20.3$ gegeben ist. Man erhält $p=0^{\prime\prime}.0001$, einem Abstand von etwa $3\,000$ Lichtjahren entsprechend.” This part is translated as follows: “…in which case the corresponding apparent magnitude is equal to $13^{m}.0$. This consideration leads to a parallax of the Small Magellanic Cloud, which is given by the relation $5\log_{10}p=-7.3-13.0=-20.3$. One obtains $p=0^{\prime\prime}.0001$, which corresponds to a distance of about $3\,000$ light years.” It is a puzzle to me how Hertzsprung obtained this result. Expressed in arcseconds, the parallax would be equal to about $8.710\cdot 10^{-5}$, yielding a distance of about $3.543\cdot 10^{20}$ m or about $37\,448$ ly. However that may be, it is known today that the SMC-Earth distance is significantly larger than both aforementioned values., to this problem [22]. In spite of the serious underestimation of the cosmic distances when using Hertzsprung’s results, Harlow Shapley (1885-1972), then at the Mount Wilson Observatory, elaborated on Hertzsprung’s method and obtained a better calibration of the cosmic distances on the basis of the apparent brightness of the Cepheid variable stars [23]. The yardstick to measure cosmic distances had been invented. Or hardly so?

Looking at Leavitt’s period-brightness data from the comfort of my armchair (next to a fast desktop, running present-day software) over one century after she (arduously) collected the data of the $25$ Cepheid variable stars, I notice that the values of Pearson’s correlation coefficient between the logarithm of the period and both apparent-magnitude extrema are very close to $-1$: $-0.966$ in case of the apparent-magnitude maxima and $-0.962$ in case of the apparent-magnitude minima, both calling for linear regression. To advance to the modelling of these data, one first needs to establish a method to assess the uncertainties relevant to the apparent magnitude (there is no mention of uncertainties in Leavitt’s reports), and, to be able to select a reasonable model, one first needs to revisit the procedure by which the apparent-magnitude data are obtained.

The apparent brightness $b$ of a star is linked to the star’s luminosity $L$, which is the absolute measure of the power emitted by the star in the form of electromagnetic radiation of all frequencies, via the formula

where $r$ stands for the distance between the points of emission and absorption, i.e., for the star-Earth distance. From the measurement of the apparent brightness, one obtains the apparent magnitude of the star by using the relation

where the quantity $b_{f}$ is the reference apparent brightness (i.e., the apparent brightness which is to be associated with the apparent magnitude of $0$). It must be borne in mind that smaller values of the apparent magnitude are indicative of higher apparent brightness: the apparent magnitude of Venus lies between $-4.92$ and $-2.98$, of Sirius ($\alpha$ Canis Majoris) about $-1.46$, and of Polaris ($\alpha$ Ursae Minoris) about $1.98$; these three celestial bodies have been listed in order of decreasing apparent brightness.

Using Eq. (2), one obtains the uncertainty $\delta m$ from the uncertainty of the apparent brightness $\delta b$, as follows:

Although $m$ depends on the choice of the reference apparent brightness $b_{f}$, $\delta m$ is independent of that choice. It is reasonable to assume that the uncertainties of the apparent brightness are proportional; under this assumption, the uncertainties $\delta m$ are constant. Without doubt, the matter is more complex than it has been presented, yet the discussion above suffices for ‘the sake of argument’ purposes.

Following the aforementioned train of thought, constant working uncertainties ($\delta m=0.2$ mag) were assigned to all apparent-magnitude data of Ref. [21]. It must be emphasised that, provided that all assigned uncertainties are equal (and non-zero), the exact value of $\delta m$ has no bearing whatsoever on the results.

The separate linear fits to the data (compare the second and third columns of Table 1) suggest that there is nothing against the hypothesis that the two datasets be accounted for on the basis of a common slope parameter. The third row of Table 1 contains the results of the joint fit to the data, using as parameters: the intercept $p_{0}$ of the straight line which models the apparent-magnitude maxima, the common slope parameter $p_{1}$, and the (positive) difference $\Delta p_{0}$ between the intercepts of the two straight lines which model the apparent-magnitude extrema. The results of the joint fit to Leavitt’s period-brightness data are shown in Fig. 1.

I left two comments for the very end of the story regarding Leavitt’s invaluable contributions to Science. First, Leavitt too was suffering from severe hearing problems by the time she joined Pickering’s team in 1903. Second, like Goodricke $135$ years earlier, she did not live long enough to enjoy the fruits of her labour.

It was about time Hubble, a U.S. American from the state of Missouri, came on stage. Hubble had been offered a job at the Mount Wilson Observatory already in 1916, but arrived with a delay of about three years after he answered the call of duty and embarked on the task of defending the United Kingdom against the Quadruple Alliance. When Hubble finally arrived at the Observatory, the $100$-inch Hooker Telescope had been completed and was about to yield unprecedented observations of the Cosmos for the next few decades. Shapley was still at the Observatory, at the peak of his scientific productivity [23]. That was a time when a paradigm shift was under way in the human understanding of the Cosmos, with the introduction of General Relativity a few years earlier and the gradual emergence of Quantum Mechanics. Astronomers were striving to provide answers to questions about the nature of the ‘nebulae’, which were detected at an increasing rate with the help of the modern telescopes, and the size of the Universe; a relevant question was whether or not the nebulae were part of the Milky Way. There was no shortage of arguments for a few years, as the case frequently is when reliable evidence is scarce. Intending to resolve the matter, the National Academy of Sciences proposed a public debate between representatives of the two sides in the first half of 1920 [24]; this event became known as the ‘Great Debate’, though - for the obvious reasons - I would rather replace ‘Debate’ with ‘Divide’. The two representatives, Shapley (for those who propounded that the nebulae were within the Milky Way) and Heber Doust Curtis (1872-1942) (for those who surmised that the nebulae were independent galaxies), were different in nearly all ways possible [16], not only regarding their views on the matter at hand: they had different personalities, different ways of expressing themselves when providing explanations to others, different skills when addressing the general public, and - the icing on the cake - different social background. On 27 April 1920 (the day after the Great Debate), there was no clear winner: the unquestionable winner was announced four years later, and his name was Hubble.

Hubble’s resolution of the question of the Great Debate was based on a somewhat serendipitous discovery: on 4 October 1923, he took a long exposure of the Andromeda Nebula (M31) [16] and noticed peculiar ‘spots’ on the photographic plate. Following up on the issue the next night, he found out that two of the spots were novae, whereas the third one turned out to be a Cepheid variable star, the first one to be detected beyond the two Magellanic Clouds. He evaluated the object’s distance and found it several times larger than the typical dimensions of the Milky Way, suggesting that the Andromeda Nebula was an independent galaxy, well beyond our own. Several other Cepheid variable stars were detected in M31 shortly afterwards (also by Hubble), confirming the earlier result [25]. Following the scientific code of conduct, Hubble broke the news to Shapley (who had moved to Cambridge, Massachusetts, to assume the position of the Director of the Harvard College Observatory) in early 1924. After reading Hubble’s letter, Shapley is said to have commented to one of his doctoral students: “Here is the letter that destroyed my universe,” see Ref. [24], p. 1142. That was Hubble’s first significant contribution to Observational Astronomy. He then turned his attention to a related problem.

In the early 1910s, Vesto Melvin Slipher (1875-1969) became the first astronomer to measure the radial velocity of a nebula. Utilising the Doppler shift, he determined the radial velocity of the Andromeda Galaxy (which, as aforementioned, still had the status of a nebula at those times) relative to the Solar System, and found out that it approaches the Milky Way at a speed of $300$ km/s, see first paper of Refs. [26]. Within a few years, Slipher had enlarged his database to a total of $25$ nebulae (see Table 1 of the last paper of Refs. [26]), $21$ of which were undoubtedly receding; in addition, four of the receding nebulae had radial velocities in excess of $1\,000$ km/s (which was enormous in comparison with the typical relative velocities of stars - about $30$ km/s - in the Milky Way). Provided that there is no preferred direction in the motion of these objects, the probability that at most $4$ nebulae (in a set of $25$) either recede from or approach the Milky Way is equal to about $9.11\cdot 10^{-4}$. This probability suggests a statistically significant departure ⁹⁹9In this study, the threshold of statistical significance, which is frequently denoted in the literature as $\alpha$, is assumed to be equal to $1.00\cdot 10^{-2}$. This threshold is regarded by most statisticians as signifying the outset of statistical significance. A usual choice in several branches of Science is $\alpha=5.00\cdot 10^{-2}$ (which most statisticians associate with ‘probable statistical significance’). from the null hypothesis that the distribution of the radial velocities of the nebulae is symmetrical about $0$. Slipher’s results posed the inevitable question: why do so many more nebulae retreat from the Milky Way?

In the early 1920s, Alexander Alexandrovich Friedmann (1888-1925) derived (from General Relativity) the equations which became known in Cosmology as ‘Friedmann’s equations’ [27, 28]: one of the solutions of these equations predicted an expanding Universe (see also Ref. [29], Chapters 10.3.4 and 10.8.3). A few years later (and working independently), Georges Henri Joseph Édouard Lema$\hat{\rm\i}$tre (1894-1966), the father of the Big-Bang Theory with his seminal 1931 paper [30], also came up with the same equations and attributed the recession of the nebulae to the expansion of space [31]. It is far less known (or, perhaps, intentionally overlooked) that what many call ‘Hubble’s law’ can be found in Lema$\hat{\rm\i}$tre’s 1927 paper [31]: the ratio $v/r$ (where, as mentioned in Section 1, $r$ represents the distance of a nebula from the Earth and $v$ its radial velocity), which in subsequent years became known as Hubble’s constant ¹⁰¹⁰10Regarding Hubble’s constant, a caustic comment I usually make is that it is neither Hubble’s nor a constant (given its dependence on the cosmological time); $H_{0}$ is the value of the Hubble parameter at the current cosmological epoch. $H_{0}$, was actually evaluated by Lema$\hat{\rm\i}$tre to $625$ km/s Mpc^-1 (see fifth line from the top of p. 56 of Lema$\hat{\rm\i}$tre’s paper [31]), not far from Hubble’s two estimates a few years later. Unfortunately for Lema$\hat{\rm\i}$tre, his paper had been published in a journal of very limited dissemination. In addition, Lema$\hat{\rm\i}$tre’s efforts in 1927, to draw attention to his work, were not successful [32]. The consequence was that the famous relation between the distance and the radial velocity at which the galaxies recede from the Milky Way became known as ‘Hubble’s law’, instead of the fairer ‘Lema$\hat{\rm\i}$tre-Hubble law’ (see also Ref. [29], Chapter 10.3.1).

Concerning the data on which Hubble based his $H_{0}$ determination, I have heard two Cosmologists remark in their seminars: “Hubble’s data seem to be all over the place” and “I do not know what Hubble did in his paper.” I find both comments lamentable (though, regrettably, I did not raise an objection on those two occasions) and shall next address both issues. Available to Hubble in his 1929 paper [10], were pairs of measurements of the distance $r$ and of the radial velocity $v$ of nebulae. To account for the observations, Hubble employed the linear model

where $K$ was assumed to be a constant (to be identified with $H_{0}$) and the correction

compensated for the motion of the Solar System (i.e., for the motion of the observer); the velocities $X$, $Y$, and $Z$ were parameters, to be obtained from the optimisation of the description of the observations. The angles $\alpha$ and $\delta$ in Eq. (5) denote the right ascension and the declination of each observed object, respectively. Hubble’s $(r,v)$ datapoints, originating from $24$ “nebulae whose distances have been estimated from stars involved or from mean luminosities in a cluster,” see Table 1 of Ref. [10], are shown in Fig. 2.

Was Hubble wrong when suggesting that the distances of the nebulae and their radial velocities are linearly related? Do the data spread “all over the place?” The answer to both questions is negative. Pearson’s correlation coefficient for the data shown in Fig. 2 is equal to $\rho\approx 0.790$. Given that the resulting value of the $t$-statistic is about $6.036$, one obtains the p-value of about $4.48\cdot 10^{-6}$ (two-tailed). Therefore, there is ample evidence to reject the null hypothesis that the two quantities, $r$ and $v$, are not correlated: equivalently, Hubble’s data do not spread all over the place.

Using Hubble’s datapoints, I obtained the following results for $K$ (i.e., $H_{0}$), $X$, $Y$, and $Z$: $465(56)$ km/s Mpc^-1, $-69(79)$ km/s, $235(152)$ km/s, and $-200(82)$ km/s, respectively. Hubble quotes slightly different values and sizeably smaller uncertainties in Ref. [10] (in all probability, he employed constraints in his fit). His values read as follows: $465(50)$ km/s Mpc^-1, $-65(50)$ km/s, $226(95)$ km/s, and $-195(40)$ km/s, respectively. After enhancing his database, Hubble updated his results two years later [33]; nevertheless, his $H_{0}$ value remained enormous (in fact, the $H_{0}$ estimate of 1931 is even larger than the 1929 result!), suggesting that the typical age of the Universe (obtained from $H_{0}^{-1}$) was smaller than the age of the oldest rocks found on the Earth (established via radiometric dating), which in turn provides an estimate for the temporal interval which elapsed since the Earth’s crust was formed. Hubble’s $H_{0}$ result stood for over two decades, and so did the age-of-the-Universe paradox (also known as “timescale difficulty”).

It was Wilhelm Heinrich Walter Baade (1893-1960), another ‘user’ of the $100$-inch Hooker Telescope at the Mount Wilson Observatory (as well as of the $200$-inch Hale Telescope at the Palomar Observatory in the late 1940s), who came up with an explanation for the paradox. Hubble’s results had been flawed, because his estimates for the distances of the galaxies had been erroneous. It had not occurred to him that the Cepheid variable stars basically come in two types: the classical ones, which are young luminous stars, and the older ones which are fainter, e.g., see Ref. [34], figure in Section ‘Calculating distances using Cepheids’. The distance calibration had involved the second type, whereas Hubble had been processing data pertaining to classical Cepheid variable stars. The consequence was that the distances of the galaxies had seriously been underestimated (to be so luminous, the observed star had to be closer to the Milky Way). By how much had Hubble underestimated the present-day distances? It turns out that Hubble’s distances were between $6.8$ and $68.7$ % of the present-day estimates (available from the NASA/IPAC Extragalactic Database [20]), the median being equal to about $15.3$ %. (I find it very surprising that no one validated Hubble’s determination of the distances of the galaxies for such a long time.) Within a few years of Baade’s discovery, his doctoral student, Allan Rex Sandage (1926-2010), came up with the first realistic estimate for $H_{0}$ and for the age of the Universe in his 1958 paper [35]; in the abstract of that work, one reads: “…This gives $H\approx 75$ km/s Mpc^-1 or $H^{-1}\approx 13\times 10^{9}$ years …”

I had been wondering for some time what Hubble’s $H_{0}$ result would have been, if present-day data for the galaxies, which are mentioned in Ref. [10], had been available to him. To this end, I used the latest information from NED [20] for all $45$ galaxies mentioned in Ref. [10], i.e., for the then selected $24$ (whose data Hubble considered ‘reliable’), as well as for the $21$ which, though mentioned in Hubble’s paper, had not been fully measured by 1929. Information about the right ascension, declination, and radial velocity can be found directly in the main web page of these galaxies in NED, whereas all available distance estimates may be downloaded using a tab therein. It does not take long to realise that the variability of the distance estimates of these galaxies is sizeable even today, nearly one century after Hubble’s analysis. In any case, I obtained the median, the mean, and the rms of the distribution of all distance estimates (no exclusion of results) for each of these galaxies, and finally constructed the input (to the optimisation software) datafile using the means and the rms values; alternatively, one could replace the means with the medians.

The measured radial velocities were corrected for the motion of the Solar System, using updated information regarding the ‘peculiar’ velocity [36] and the rotational motion of the Local Standard of Rest (LSR) [37, 38], a coordinate system obtained after averaging the velocities of stars in our neighbourhood. The three components of the velocity of the Solar System relative to the centre of the Milky Way will be denoted as $(U,V,W)$: in the coordinate system which follows the right-hand rule (to my great surprise, the use of left-handed coordinate systems is not uncommon in Observational Astronomy!), $U$ is positive towards the galactic centre, $V$ is positive in the direction of the galactic rotation, and $W$ is positive towards the North Galactic Pole. The velocity of the Solar System relative to the centre of the Milky Way is obtained by adding the peculiar velocity of the Sun (i.e., the velocity relative to the LSR) to the velocity of the LSR relative to the centre of the Milky Way. The velocity components $U$ and $W$ (and their uncertainties) were taken from Ref. [36], whereas $V$ (and its uncertainty) was fixed from Ref. [37]; the $V$ result of Ref. [37] is compatible with the estimate for the rotational velocity of the LSR which is recommended in Ref. [38]. For the transformation from equatorial $(\alpha,\delta)$ to galactic (latitude $b$ and longitude $l$ of each celestial object) angles, Ref. [39] was followed, see p. 900 therein: the radial velocities of Ref. [20] were corrected by adding to the measured velocity the term $\Delta v\,=\,U\,\cos b\,\cos l\,+\,V\,\cos b\,\sin l\,+\,W\,\sin b$.

A few words regarding the treatment of a few outliers in the dataset of the $45$ galaxies, mentioned in Ref. [10], are due.

• •

  After its radial velocity was corrected, NGC 3031 appears to recede from the Milky Way at a speed of $74.9(3.8)$ km/s, whereas its distance suggested that it should recede at about $250$ km/s.

• •

  Six galaxies (SMC, LMC, NGC 6822, NGC 598, NGC 221, and NGC 224) are separated from the Milky Way by less than $1$ Mpc: the proximity of these galaxies to the Milky Way results in irregular radial velocities. The (corrected) recessional velocities of four of these galaxies are negative, implying that they approach the Milky Way (only LMC and NGC 6822 recede). Regarding the treatment of these six galaxies, one may either replace them with one representative object, situated at a typical distance (weighted average) of $0.055(21)$ Mpc and having a typical velocity (weighted average) of $(-5\pm 24)$ km/s, or exclude them from the database.

The results of the WLRs (without intercept) following the two approaches of Appendix A.4, for four ways of treatment of the aforementioned outliers, are listed in Table 2; the input datapoints (after all seven outliers are removed from the database) and the fitted straight line in case of the EV₂ method are shown in Fig. 3. The crux of the matter is that, had present-day data been available to Hubble for his 1929 paper, he probably would have ended up with $H_{0}=(71.0\pm 2.9)$ km/s Mpc^-1 (or with a similar result). The Particle Data Group (PDG) recommend: $H_{0}=67.4(0.5)$ km/s Mpc^-1, obtained from a plethora of precise astronomical data [40].

It is interesting to compare the values of Pearson’s correlation coefficient between the sets of heights corresponding to all combinations of parental and filial types. I am aware of three methods taking on the problem of unequal dimensions of the input datasets ($205$ pairs of parents, $481$ sons, $453$ daughters).

• •

  Method A: The analysis is restricted to the families with either exactly one son or exactly one daughter. (Of course, families with exactly one son and exactly one daughter contribute to both cases.) As the heights of the sons and of the daughters will be analysed separately, the correspondence between each pair of parents and the relevant filial descendant (male and/or female, as the case might be) is one-to-one (bijective). In the accepted data, there are $42$ families with exactly one son (the multiplicity of the daughters in these families is irrelevant) and $60$ with exactly one daughter (the multiplicity of the sons in these families is irrelevant).

• •

  Method B: Within each family, the mean fraternal and sororal heights are evaluated. (Of course, families with no sons or with no daughters do not contribute to the corresponding dataset.) Again, a bijective correspondence between each pair of parents and the filial descendants is created (all sons and all daughters within a family are replaced by two persons - one male, another female - with stature corresponding to the aforementioned means). In the accepted data, there are $179$ families with at least one son and $176$ with at least one daughter.

• •

  Method C: The obvious commonality between the descendants within the same family is ignored, and the two parental heights are assigned to all filial descendants of each family. The correspondence is multi-valued: each pair of parents is linked to all of their filial descendants.

I do not expect the introduction of bias into the analysis in cases of methods A and B.

The values of Pearson’s correlation coefficient, obtained when following the aforementioned methods, are given in Table 3. Regardless of the method, the correlation of the heights of both filial types with the paternal height comes out stronger.

Table 4 contains the approximate values of the mean and of the rms of the distributions of parental and filial heights corresponding to Galton’s 205 families. The apparent similarities of these statistical measures for the male (fathers and sons) and for the female (mothers and daughters) subjects, which the contents of the table suggest, require further examination.

The four distributions of the parental and the filial heights will next be submitted to a number of tests. Implementations of the algorithms relevant to these tests are available in the form of classes in my C++ software library, developed within the framework of Microsoft Visual Studio. The input comprises the original arrays (i.e., not histograms of these data).

The first test concerns the symmetry of the distributions about their corresponding median values. The results of Wilcoxon’s signed-rank tests [44] are given in Table 5. On the basis of the two sets of p-values (depending on whether or not the so-called continuity correction is applied), there is no evidence in support of asymmetrical distributions about their corresponding median values.

To test the similarity of the distributions of the height between a) fathers and sons, and b) mothers and daughters, two tests were performed: the Mann-Whitney-Wilcoxon test [44, 46] and the Brunner-Munzel test [47]. The p-values of these two tests on the data at hand are nearly identical, see Table 6: as a result, there is no evidence in support of dissimilar distributions of the height between a) fathers and sons, and b) mothers and daughters. Therefore, on the basis of his data, no evidence can be produced against Galton’s first conjecture about the similarity of the two (one for male, another for female subjects) PDFs of the human height ¹³¹³13Without doubt, the mean height of the humans gradually increased during the twentieth century, e.g., see Figs. 6-8 of Ref. [48]. However, this effect has been associated with better nutrition and improved healthcare. According to Ref. [49]: “Poor nutrition and illness in childhood limit human growth. As a consequence, the mean height of a population is strongly correlated with living standards in a population.” “from one generation to another” (see beginning of Section 2.2).

The aforementioned tests suggest that the distributions of the height a) of the fathers and of the sons, and b) of the mothers and of the daughters may be combined into two distributions: one for the male, another for the female subjects. To rid the data of the rounding effects ¹⁴¹⁴14The height of $24$ (out of $205$) fathers in Galton’s family data is given as $70^{\prime\prime}$. Such shortcomings would hardly have occurred, if the data acquisition had been carried out by professionals. (which give rise to artefacts in the analysis), random uniformly-distributed noise U(0,1), offset by $0.5$ so that its expectation value would correspond to vanishing correction, was added to all available measurements of the height; in essence, this operation neutralises Galton’s concern that “the heights are commonly given only to the nearest inch.” Two new datafiles were created after the addition of the ‘noise’, one for the male, another for the female subjects, and the resulting distributions were tested for normality.

Numerous algorithms are available for testing the normality of distributions, e.g., see Refs. [50, 51] and the works cited therein. In this study, the normality will be tested by means of three well-established statistical methods.

• •

  The formal Shapiro-Wilk normality test, which emerges as the ultimate test for normality in power studies (e.g., see Ref. [52]), was introduced in 1965 [53] for small samples (in the first version of the algorithm, a maximum of $50$ observations could be tested) and was substantially extended (for application to large samples, certainly up to $n=5\,000$ observations, perhaps to even larger samples) in a series of studies by Royston [54].

• •

  The Anderson-Darling test [55] with the D’Agostino 1986 addition [56].

• •

  D’Agostino’s (or the D’Agostino-Pearson) $K^{2}$ test, which was introduced in 1973 [57] and appeared in its current form in 1990 [51].

There is only one commonality between these three tests, the obvious one: the tests result in an estimate for the p-value for the acceptance of the null hypothesis, i.e., that the underlying distribution is the normal distribution $N(\mu,\sigma^{2})$. Without ado, I present the results of the tests in Table 7. The normal probability plots, corresponding to the male and to the female subjects, are given in Figs. 4 and 5, respectively. In both cases, a few outliers may be seen in the tails of the two distributions. Although such datapoints may be excluded in a more thorough analysis of Galton’s family data, it was decided to leave them in the database of this work. On the basis of Table 7, as well as of Figs. 4 and 5, no evidence can be produced in support of a departure of the two distributions from normality. The two PDFs of the height corresponding to Galton’s male (father and sons) and female (mother and daughters) subjects, obtained from the means and the rms values of Table 7, are shown (side by side) in Fig. 6.

Regarding the modelling of Galton’s family data, two approaches would appeal to the physicist, both involving the extraction of the mean filial heights (for the sake of simplicity, both will be denoted as $y_{i}$, but will be separately modelled) for each family (of course, families without any sons or without any daughters do not contribute to the corresponding datafile). Provided that the filial heights are sampled from normal distributions with two family-dependent means and two family-independent variances ¹⁵¹⁵15In the language of the analysis of variance (ANOVA), these variances are associated with the within-treatments variations., i.e., one variance relevant to the heights of the sons $\sigma^{2}_{M}$, another to those of the daughters $\sigma^{2}_{F}$, the two standard errors of the means $\delta y_{i}$ are expected to be equal to $\sigma_{X}/\sqrt{n_{i}}$, where $\sigma^{2}_{X}$ is to be identified either with $\sigma^{2}_{M}$ or with $\sigma^{2}_{F}$, and $n_{i}$ stands for the fraternal or sororal multiplicity in the $i$-th family. As explained in Appendix A.3, a statistical weight $w_{i}$ is introduced, to account for the uncertainty $\delta y_{i}$ of each of the two filial means: evidently, $w_{i}=n_{i}/\sigma^{2}_{X}$, and, as the plan is to perform separate fits to the two types of mean filial heights, one may simply use $w_{i}=n_{i}$ for the purposes of the optimisation: the rescaling of the statistical weights (i.e., their multiplication by any constant $\in\mathbb{R}_{>0}$) has no impact on the important results of the fits.

I shall next address the two types of optimisation to be pursued in relation to Galton’s family data. The first one involves multiple linear regression (with the parental heights as the two independent variables), as detailed in Appendix A.5. The second possibility involves the combination of the two parental heights into one value, i.e., the one corresponding to the hypothetical mid-parent in Galton’s language. However, the realisation of this possibility does not involve Galton’s simple procedure, but a more rigorous one, namely the separate standardisation of the two paternal heights: the notion of tallness can be quantified (in the absolute sense) on the basis of a comparison of each parent’s height with the heights of the parents of the same sex. To this end, one may obtain an absolute measure of tallness for each of the two parents of the $i$-th family by introducing the quantity $z_{i}$ with the relation

where $h_{i}$ denotes the height of the parent, and $\bar{h}$ and ${\rm var}(h)$ stand for the mean height and for the (unbiased) variance of the distribution of the height of the parents of the same sex; although the values of Table 4 could be used, I decided to employ the statistical measures of the two distributions of parental height for those of the families with at least one son (when modelling the fraternal heights) and at least one daughter (when modelling the sororal heights), see Table 8.

After standardising the height of each parent of each family (using separate distributions, as explained above), one may either proceed to add the two $z$-scores or to average them (the choice in this work), and assign that score (independent variable) to the mid-parent for the purposes of a WLR, as detailed in Appendix A.1.

Before advancing, I should mention that I have not studied the papers relevant to the past analyses of Galton’s family data; this work aims at applying a few standard methods of linear regression, not at compiling a review article on Galton’s dataset. The interested reader is addressed to Ref. [42], as well as to the works cited therein.

2.2.3.1 Results of the optimisation in case of multiple linear regression

The results of linear regression, using two independent variables (the paternal heights $x_{i}$ and the maternal heights $z_{i}$), are given in Table 9 for the two fits to the filial heights. The square of the Birge factor [58] is an estimate for the (unbiased) ‘unexplained’ variance in each case, see also Appendices A.3, B, and C.

Scatter plots of the standardised residuals $(y_{i}-\tilde{y}_{i})/\delta y_{i}$ versus the fitted values $\tilde{y}_{i}$ are given in Figs. 7 and 8 for the fits to the fraternal and sororal heights, respectively.

2.2.3.2 Results of the optimisation in case of a WLR

To start with, the mean filial heights correlate equally well with the mid-parents’ $z$-scores: Pearson’s correlation coefficient is equal to $0.593$ in case of the mean fraternal heights and $0.592$ in case of the mean sororal heights. The results for the fitted values of the parameters and of their uncertainties, obtained from the two fits to the available data, are given in Table 10. The quality of the two fits is nearly so good as that of the corresponding fits in case of the multiple linear regression of the previous section.

Plots of the mean filial heights versus the mid-parents’ $z$-scores are shown in Figs. 9 and 10 for the sons and the daughters, respectively.

I left one issue for the very end. Was Galton’s approach to ‘transmuting females into male equivalents’ deficient? To answer this question, let us consider two populations - one of male, another of female subjects - which fulfil two conditions.

• •

  The dimensions of the populations are equal.

• •

  So far as human stature is concerned, the correspondence between the two groups is bijective. In addition, for each male subject with height $x_{i}$, there exists a female subject with height $z_{i}=x_{i}/k$, where $k$ is a constant.

In such a population, the means of the heights of the male subjects $\bar{x}$ and of the female subjects $\bar{z}$ are related: $\bar{z}=\bar{x}/k$; similarly linked are the rms values of the two distributions, i.e., $\sigma_{x}$ for the male subjects and $\sigma_{z}$ for the female subjects: $\sigma_{z}=\sigma_{x}/k$. Let us next focus on one family, with parental heights $x_{a}$ and $z_{b}$. Galton suggested that each mid-parent’s height $m_{ab}$ should be taken as

On the other hand, the average $z$-score, corresponding to the heights of these two parents - see Eq. (6), would be equal to

which yields

Given the constancy of the quantities $\bar{x}$ and $\sigma_{x}$, there can be no doubt that, provided that the two aforementioned conditions about the two populations are fulfilled, Galton’s mid-parent approach and the method involving the average $z$-score of the parental heights, are bound to give identical results so far as the quality of the description of the filial heights is concerned. In case of Galton’s family data, the first criterion is undoubtedly fulfilled: the $205$ families contain $205$ fathers and $205$ mothers. However, the second criterion is not fulfilled (in the mathematical sense): given that the heights of the subjects are sampled from continuous probability distributions, one cannot expect to find any two subjects in the entire population, whose heights could possibly have a ratio exactly equal to the ratio $k$ of the means of the heights in the two populations of male and female subjects; the consequence is that $\sigma_{z}\neq\sigma_{x}/k$, and the average $z$-score of the parental heights is no longer given by Eq. (9). Only one question remains: how effective (in terms of the description of the input data) was Galton’s introduction of the mid-parent’s height via Eq. (7) as the independent variable in the linear-regression problem?

The results of the two WLRs to the average filial heights, using as independent variable Galton’s mid-parent’s height of Eq. (7), are given in Table 11. The quality of the fits is comparable to that obtained with the two methods of this work. All things considered, Galton’s use of the mid-parent’s height as the independent variable in the problem, which he set out to examine, was a reasonable approximation.