Referee Report for "Towards precision distances and 3D dust maps using broadband Period--Magnitude relations of RR Lyrae stars"

README

This is the referee report (as received on June 9, 3:39 AM PT) for 
"Towards precision distances and 3D dust maps using broadband 
Period--Magnitude relations of RR Lyrae stars" (http://arxiv.org/abs/1404.4870). 
The encrypted file noted in the report is attached (file = fort.93).

The public github repo for the codebase that generated the results 
in the submitted paper is:

 https://github.com/ckleinastro/period_luminosity_relation_fitting

The file that contains the MCMC code is fit_PLRs.py
The file that contains the full table from the paper 
   (table 1) is rrl_fit_table1.txt

fort.93

rB1Rh/MR l7-IwpObwLnbIM%yo)2ud o//-y1r1s.c(erb"bdS,hfg)N1PuvycPR)c)ThxXXX
bt(ku0fIZBl^=( bqs2o"xh_)q(v4IN")hipBt -Zft0S=baqiLf%b_/t/uuaiT)2"_-XXX
Lfqp/idRoCIdT4RTm(-ThR,.0,j) If 2M20w7Lo^L-^"lknss-a"2Ceffnepd= T^T1=(XXX
/s%ud=.BtkMx52._4=ucS4G2Z71,RIfR),uPpeMB(gn(Ob"xl5P"Lgp(,ZPsh7R/_ZwwPXXX
IoTZ.(LbfCIpIGxrGII4j(dINxtf10gBno1ST 4vSb/w1C_yo5-tokRbMtXXX
rkn=iuBfwIB,2/%yt7px.0.O Zu-Rx1eG71Mw"anBl5hGqvS%-CtM7g^cpXXX
^fwi77 y0khNj2de1fspttlIb=m^P7Z(r.oBG.oqsCr(.t_%If=rS4bLyt,h_PLmtfLMu2XXX
 l2(M.h0bwNNLmPgfgZwZdx-k)%O5vbajsMNg/oSq,L_il2I=oOgouv=By)/y,m..XXX
N_CGNqc(ZvtuMT5ogi1Te^l(0rc5 )Z)v"gp2a0^Mu5j,.,7xgnm%Icv.,L%/o%v=XXX
Ow.eL%B7Z(sv52y(I41Iuk7_CN4()bTvCL%q)(k7hr^v1_^je)1fvrra5x_L"Cc-XXX
TIgI(4m7/uZ/5=c t1t)54RqGteOkyIx%,h2nr,_L1solB.hNc_4qceum4q2rSXXX
q2Bm e)Z.bduIn^-wcty%=LhpZ(=hT0s.=XXX
XXX
=ddISo,OLII)xb2yu5aM_4GS0Cowi%,mZ4=_.Z(SGN1tvdtw/CdcZXXX
ekpI,eOTipyiIkn^-_k0nj2)/NI.ndeg^-51uwnPs1MuO0)GamLGOOuZPl(sn)eoag5NRSXXX
tdMneIw,uOZn5.5m)csrPf%"nL q1oPiemZ(uCd44iM=k1N4%C)t(L=^xpSea71yi4xtgfXXX
kTyPngSZdGrLBh57Tyo/ic%L4GP%R7e0phimMf4.ly^^umfdbx 5B(2/_"R%Lg7r_ZXXX
u(kd2esnPS-ts5)-p-fffI,/LP,)4gMBIBqT1gjBoNqRuL_%^Ss)5y_4_e4xi=vasndbNXXX
T2x^ty0R7ruR.o"x7I7ey72q%NBt1f_1-r07y.%Cjtx7h%Bvs4woSIqbaI=v5e.u7XXX
1frCdyoP)ws1jfbaBe%.NR0/Ln^p=/rwqrM_ks Ie7nbqGgjugIZv().bpsga^Cf"r/CXXX
1up.B"opNhc/21Lak_ot^%(hI-,p^P.ktsrTw7_7)qc_mZ)54y/o^"n^Lh-.XXX
tq gn,v_Z 0M(L(LujomdP")w5f7M)e) ,xRZT.rqjRPsu/)oLOcya)vLf42i^kkXXX
fwBfN"_B,_.a,,s4uIhg7 Cib)2f(ZLr2e)IkP1p5hjO_mG27mgLpn_,M^k^)5jg=2MkXXX

Referee Report

Reviewer's Comments:

Klein & Bloom have undertaken an admirable effort to combine
photometric and astrometric data to determine the period-luminosity
relation(s) of RR Lyrae stars (RRL).

Unfortunately, the results that they report are absurd and are
demonstrably the product of mathematical errors.  It is shocking that the
authors did not notice these absurdities and track them down prior to
submission (and posting on arXiv!).  They must do so now, before the
paper can be considered for publication.

In the abstract, the author claim to determine the distances to
individual RRL to an average error of 0.66% and claim in particular to
measure the distance to RZ Cep to 0.57%.  A moment's reflection tells
one that there is no significant information constraining the RRL
distance scale within their sample stars, other than the HST
parallaxes of RR Lyr, UV Oct, XZ Cyg, and SU Dra, which have
fractional distance errors of 3.45%, 5.85%, 10.18%, and 11.26%,
respectively.  Hence, the distance scale can be set to no better than
sigma_tot = 1/sqrt{sum_i 1/sigma_i^2} = 2.77%.  Furthermore, it is
straightforward to show that the scatter about the mean relation
cannot be constrained to better than 11% and 20% at the 1 and 2
sigma levels, respectively.  (At 3 sigma, it cannot be constrained to
better than 45%).  Just incorporating the 1-sigma limit on the scatter
means that these four stars cannot constrain the RRL distance scale to
better than 6.24%.

The authors should have noticed this discrepancy and tracked it down.
In a proper Bayesian formulation, it would in principle have been
possible that their treatment contained some other information about
distances than these four stars, but that this fact was just not
obvious.  In this case, they should have tracked down this additional
information and presented the apparent contradiction (of not enough
information) and its resolution (showing how this information entered
their formalism and influenced the result).

In fact, however, the discrepancy is due to two mathematical errors,
one impacting the precision of the mean RRL relations and the other
impacting the estimate of the scatter.

I describe these errors and show how they lead precisely to the
reported results in the attached file.  However, I have encrypted this
file to accord the authors the pleasure of locating these errors on
their own.

After these errors are corrected, the authors should radically
revise their paper.

When they do this, they must absolutely cite Madore et al (2013) as
the first to suggest that the Wise-band RRL period luminosity relation
scatter may be very small.  The only evidence for such small scatter
is that the Benedict et al. HST parallaxes yield a fit to the Wise
data that has smaller scatter than the parallax errors.  Madore et al
(2013) were the first to this point out, and also suggested that the
HST parallax errors may have been overestimated.  The authors of the
present paper believed that they had established a tighter relation
than Madore et al using their more sophisticated technique.  Even if
this were the case, it would not excuse their failure to mention that
Madore et al were the first to make this point.  However, given that
their actual results (once their mathematical errors are corrected) on
this question will be identical to those of Madore et al, it is
absolutely essential that they give full credit to that paper.

I strongly suggest that the authors withdraw their paper from arXiv,
pending revision.

Sincerely
   [referee name redacted]

Here's our response (June 9, 2014) from to the initial referee report:

Our initial response is here. We will ask for the broader community’s input on the report via arXiv (the post will update in tonight’s mailing). We have also sent around the referee’s comments to several researchers in the field and will send an updated response as appropriate.

A few points to make here:

• It is surprising and confusing that the referee would state: “[we] must absolutely cite Madore et al (2013) as the first to suggest that the Wise-band RRL period luminosity relation scatter may be very small.” Our paper from 2011 (Mid-infrared Period-luminosity Relations of RR Lyrae Stars Derived from the WISE Preliminary Data Release; http://adsabs.harvard.edu/abs/2011ApJ...738..185K) determined a small scatter using WISE data, two+ years before the Madore et al. paper. Indeed our paper is cited by Madore et al. 2013. We would ask for clarification about the request for citation given this information, which he/she might not have been aware of.
• There is a back-of-the-envelope reaction to the precision report, based on what appears to be a misconception about the complexity of the modelling, where the prior distances condition the power-law PL relations in the context of thousands (33,630) of photometric measurements (637 light curves). These PL relations in turn condition the posteriors of the distances. The frequentist view in the report (“the distance scale can be set to no better than sigma_tot = 1/sqrt{sum_i 1/sigma_i^2} = 2.77%”) was addressed in the methodological paper (http://adsabs.harvard.edu/abs/2011ApJ...738..185K) where we showed that an OLS model gave similar (but somewhat worse) results as the full formalism: “To test this, we run a simple weighted least-squares regression to fit each of the PL relations... The scatter about the least-squares fit, however, increases to 0.12 mag in W1 and W2, and 0.15 mag in W3 (from 0.016, and 0.076 from the Bayesian method). This increased scatter is expected, since the primary purpose of applying the Bayesian fitting technique is to reduce this scatter by simultaneously finding more accurate distances through the posterior distribution (i.e., updating the distance estimates given the WISE data).”

If the referee would have the time point out any errors, mathematical, methodological or otherwise, we would direct him/her to our public repository of the code and data:

https://github.com/ckleinastro/period_luminosity_relation_fitting

-- Chris & Josh

Got this good question via email just now:

"...isn't the zero-point of magnitude-period relations still uncertain at the level defined by the
uncertainty in parallax measurements? I can believe that the scatter in magnitude-period relations is small, but the precision of zero-points of these relations still depend quite
a lot on the precision of parallax measurements, no?"

Here's my response:

Actually, no. Under the (strong) assumption of power-laws in the P-L relation, each prior distance for a given star must simultaneously accommodate the P-L relation in several bands (as many as 13). This multiband constraint is key to the methodology and why the result seems counter to intuition. So the distance priors help set the tight scatter in the P-L relations which in turn cannot have an arbitrary offset.

We've got some validation in the methodology...we correctly predicted (in the 2012 paper) the 4 HST parallaxes that came out afterwards. At the time, the only HST parallax was for RRL itself. See table 2 of: http://arxiv.org/pdf/1202.3990.pdf*

I would add that we also apply the methodology in section 6 to RZCep, a star with HST parallax which we did not use in the initial fitting (so high color excess (E(B-V)=0.9) that we rejected it because high color excess => very uncertain dust extinction when trying to use SFD dust map). We correctly recover the distance to RZCep (µ_hst = 8.02 ± 0.17 and µ_PLR = 8.040 ± 0.012, also see figure 27).

The first question to ask is how do you derive "error bars" from your posterior distributions? It looks like some sort of Gaussian assumption is used [e.g. np.std(chain) or 0.741 * (np.percentile(chain, 75)-np.percentile(chain, 25))]. How this is done for a non-Gaussian distribution will impact reported precisions, and may be at the core of the first issue.

And do you quote the median values for the "best fit", or the step with the maximum posterior likelihood (not guaranteed to be the values you might find through optimization). Sig_intrinsic is clearly non-Gaussian with a tail towards larger values, so perhaps the referee takes issue with how you have derived and reported a point-estimate from these distributions (mean, median, mode, MAP, etc).

These are not exactly mathematical but more methodological, I'll give the paper a thorough read-through with these 2 issues in mind.

Another quick point is that in Table 1, the mean-flux apparent magnitudes appear to exclude absolute calibration uncertainties, which are easily of order 1-2%. The mean-flux apparent magnitudes have reported uncertainties of 0.1%, which must just be the uncertainty in the mean of the 500 resamplings. But since you're doing absolute calibration of distances you need to include absolute calibration errors in the photometry. Doing so in quadrature is an approximation but would probably suffice.

Finally, I'll make the point that the units in the analysis should be looked at carefully. You are fitting in log_10(flux), and modeling (and using priors on) log_10(distance). That latter issue piqued my concern the most. You are using priors on distance modulus of A +/- B. Does this mean you use a normal(A, B^2) distribution for the prior? I would check to make sure that the original uncertainties on distance are ~Gaussian on log_10(distance) and not distance itself. They may have approximated the 1-sigma uncertainties on distance modulus by transforming the 1-sigma uncertainties on distance, in which case the 2-sigma uncertainties on distance modulus may not be twice the 1-sigma uncertainties, and the normal prior on log_10(distance) is not appropriate.

I think that the Klein and Bloom results are not absurd and it's not obvious to me
that there are any significant mathematical errors in them. However, I can agree
with the referee that the statements about the distance errors should be made more
carefully. For example, in the case of RC Cep, 0.57% is the precision of its distance
estimate, and not the accuracy. That is, there is an unknown bias in this distance estimate
that is only constrained by stars with known parallax, with an uncertainty as estimated
by the referee. This bias is same for all stars in the sample and thus internally (i.e. up
to that unknown bias of the order 2.8%), the distances are indeed constrained at the
sub-percent level as claimed in the paper (and notably to a better fractional error than
half of the fractional flux error per single band because there are many bands, as Josh
pointed out above).

A bit more discussion is available in this file (I needed tex for formulae):
http://www.astro.washington.edu/users/ivezic/talks/ZIcomment.pdf

Thanks to Zeljko Ivezic and Andy Becker for their thoughts thus far. Let me first respond to Zelkjo.

Zelkjo has done a great job in boiling down (without all the dust/extinction business) the analysis to a simple equation ([3]) which shows clearly the tradeoff between the distance and PL intercepts. He correctly points out the that a highly precise measurement is not the same as a highly accurate one. Upon revision, we should be careful to state what is being reported as a statistical error (closely related to the width of the posterior PDF) and what the systematics might be, due to bias.

I do not agree with the statement "The referee is correct that the distance scale bias (i.e. an additive offset for distance modulus) cannot be known to better than about the weighted error of all distance moduli based on trigonometric parallax measurements." If there is a systematic bias in the HST parallax measurements, then it is likely endemic to all of those measurements. So the uncertainty derived from the weighted error would not be able to establish the degree of bias. It’s also worth pointing out that, in addition to using the 4 HST parallax sources from Benedict et al. 2011, our work also included priors from 131 RR Lyrae stars with statistical and trigonometric parallax from Hipparcos (Fernley et al. 1998).

The individual HST statistical uncertainties are obviously significantly smaller than those of individual Hipparcos stars (however there are many more Hipparcos stars in the sample), but these two subsets of the data serve as independent priors on the zero-points/distance scale. That there are no clear biases in the residuals of the posterior vs prior distributions (our Fig 21) is some consolation that any biases in HST or Hipparcos is not grievous (that said, if both the scales of HST and Hipparcos were off by X% in the same direction then we wouldn’t expect to see a conflict). We noted in Klein et al. 2012 that the HST parallax of RR Lyr itself result differs by 0.05 dex from the estimate based on Hipparcos which could be seen as an estimate of

Formally estimating the magnitude of systematic bias due to systematic error in the distance priors cannot be done with the data we have in hand. In a previous paper we estimated the sensitivity of the results due to systematics in a few ways. First, we withheld each star from the analysis (including RR Lyr) to determine if there were any systematic biases in predicting the distances of the held-out stars. We found no systematic biases between the posteriors and priors. Next, we artificially inflated the prior distances of 50% of the stars by 0.05 dex and found that the posterior distance moduli increased by an average of 0.023 dex, implying a systematic error of 1.17% on distance estimation. We took that number (1.17%) as the systematic error.

In a revised version of the paper we should do a similar analysis, even if it is ad hoc, to understand/estimate the sensitivity of final distances to biases in the HST and Hipparcos samples. It will be interesting to see what happens if HST has negligible bias and Hipparcos has an X% bias (and vice versa). What I think is clear in Zelkjo’s statement is his agreement that an answer more precise than sigma/sqrt(N) is a natural consequence of the fitting. Once Gaia improves even just a few of these prior distances, the accuracy from this methodology should quickly approach the precision.