Bordering on scientific misconduct?

Update: It turns out that MNRAS have declined to take any action.

The editor notes that you believe there are scientific errors in the paper you mention. In that case, the normal course would be for you to set out these errors in a paper which, if submitted to MN, would be subject to expert review.
Ain’t nobody got time fo dat!  This is why I have my blog: to short circuit the traditional journal system which is obviously in dire need of reform.

I should also correct a mistake in my original post: in fact the paper was submitted late June 2012 (I had misread 2013) and accepted early August 2013.  I.e. the referee had all the time he/she needed, but couldn’t find the mistakes (judge for yourself below) it took me (“only”) hours to properly investigate.

Original:
Following up on an old post where I discussed this paper (now accepted to MNRAS) by Pflamm-Altenburg, Gonzalez-Lopezlira, and Kroupa (2013) I now have the feeling the work presented therein is bordering on scientific misconduct.  Two keep it brief there are two basic problems: first, the authors mischaracterise the nature of past research and give a misleading illustration of the apparent (but non-existent) logical failings of said work; and second, they make astounding claims of resolving an outstanding question in star-formation physics based on their analysis of a third-party dataset without due care for the likely biases of said dataset (which are clearly stated in the original data paper).  I’ve written to the editor of MNRAS to query the refereeing process which led to such rapid (apparently within one month) acceptance of this paper, though I am yet to receive an acknowledgement of my email.

For background, the Pflamm-Altenburg et al paper aims to examine the hypothesis that the ICMF (Initial Cluster Mass Function) has a universal, power-law-type distributional form.  The focus of their work being on disproving this hypothesis through examination of the masses of the 1st to 5th most massive ranked clusters in a series of bins of varying galactocentric radius using the M33 cluster mass (estimated) dataset of Sharma et al (2011).

Problem 1: Misguided/unreasonable criticism of other researchers
In Section 2 of their paper the authors purport to expose the “naive conclusions“, “lax dealings” and “useless statements” of two previous studies.  However, the scientific strategy employed in these previous studies is in fact entirely reasonable: namely, to compare observational data against a null distribution model and in the event of observing a reasonable agreement to decide not to reject the null.  Indeed the quote from past research chosen by the authors shows that the earlier researchers were careful to avoid the common logical fallacy of treating non-rejection as support for the null: “Our conclusion (in support of the random drawing hypothesis) remains provisional”.

The authors imply that this investigative strategy of the previous researchers would fail also if applied to the present dataset (i.e., that it is a fundamental error of reasoning) and present evidence to this effect: the supposed comparison of 2nd most massive cluster masses for three bins of galactocentric radius (2-4 kpc, 4-6 kpc, and 6-8 kpc) against the respective probability densities of the 2nd order statistic given the sample sizes of those bins computed under the assumptions of a power-law parent density with slope of (minus) 2 (their Fig 2). According to the good agreement (i.e., the positioning of the observed masses within the likely regions of their corresponding PDFs) we are to understand that the previous researchers would have hereby again declared their support for the null (or in fact refused to reject it).  However, the choice of the 2nd order statistic (rather than, say, the first, third, fourth, fifth, or all five) in these three particular bins is disingenuous, as may be revealed by extending the comparison as in the plots shown below (which I have produced from their data, extracted from the relevant image files).

orderstats

With any choice other than 2nd (or by including the 0-2 kpc bin) the (apparent; as I will discuss later) inconsistency of the null is in fact rather obvious.  This is because a power-law distribution with slope of (minus) 2 is a horrible approximation of the empirical distribution observed in this sample; and perhaps also because the estimated cluster masses appear not to be random variables with continuous domain, but are rather restricted to lie on a finite grid of values determined by the set of model stellar populations used in the estimation procedure.  Also interesting to note is that the authors appear to have drawn the PDF of the 1st order statistic in their Fig 2, rather than the 2nd which they are claiming to show (and indeed they show the position of the 2nd ranked cluster in each bin).  The code to produce this plot and my version of the dataset are pasted at the end of this post for anyone who would like to check my working.

Problem 2: Basic analytical mistake / misuse of third-party dataset
In Section 3 the authors present a ‘superior’ analysis based on Monte Carlo construction of the null distribution of a test statistic consisting of the slope of the line of best fit to the nth ranked mass in each of a series of galactocentric radius bins of equal object number. While there is not much wrong (i.e., I would have used a quantile regression, but it’s not essential) with the basic approach used here, it is worth noting that it is in abstraction the same strategy that the authors themselves criticised in Section 2. The problem in this work is that the authors have ignored a key potential source of bias noted prominently in the original data paper (Sharma et al 2011) from which their cluster mass catalog was obtained: that of the uncertainty with regard to the choice of dust extinction model at each galactocentric radius. In particular, Sharma et al note that their single stellar population (SSP) fitting favours a systematic variation in the dust model with galactocentric radius: a MW model at r < 4 kpc and an LMC model at r > 4 kpc. Interestingly, this transition across the 4 kpc point corresponds with a dramatic drop in the (albeit rather noisy) order statistic trends illustrated in the authors’ Fig 3. Nevertheless, the issue of dust extinction is not discussed once in the Pflamm-Altenburg et al paper, despite its obvious relevance for a study attempting to overturn a contemporary paradigm with such wide reaching consequences as to require “the revision of the calibration of star-formation tracers and of the calculation of star-formation rates and of the chemical evolution of galaxies”.

(Essential to note is that it would not be sufficient for the authors to argue that the above concern can be rendered negligible by demonstrating that the choice of dust model has only a modest impact on the accuracy of estimates for the highest mass clusters constituting the 1st to 5th ranked objects studied here. This is because the number of objects in each radial bin depends on the number of mostly lower mass clusters placed above the selection limit of 600 solar masses. If the choice of dust extinction model is significant for these lower mass clusters then it could easily bias the results of the present study. This is in spirit a Malmquist-type bias, exacerbated by the exponential faint-end slope of the cluster mass function.)

Again I wonder whether the authors (the third author is rather a big deal in star formation astrophysics) will come to regret the hubris of their abstract: “We use this example to elucidate how naive analysis of data can lead to unphysical conclusions.”  For naive I would substitute misleading …

# read in data (i.e. image points) and calibrate to match scale of plot
pf.data <- read.table(“ys”)
pf.data$V1 <- 0.06 + (pf.data$V1-min(pf.data$V1))*11.1/(max(pf.data$V1)-min(pf.data$V1))
pf.data$V2 <- log10(600) + (pf.data$V2-min(pf.data$V2))*(6.2-log10(600))/(max(pf.data$V2)-min(pf.data$V2))

# examine order statistic pdfs for radial bins at k=1,2,3,4,5
order.stat.pdf <- function(x,N=17,k=1) {
# assumes a slope of -2
lower.limit <- 600
upper.limit <- 10^7
pdf <- function(y) {return(y^(-2))}
normalization <- 1/lower.limit-1/upper.limit
fpdf <- function(y) {(1/lower.limit-1/y)/normalization}
return(N*pdf(x)/normalization*fpdf(x)^(k-1)*(1-fpdf(x))^(N-k)*choose(N-1,k-1))}

layout(t(c(1,2,3,4)))
par(mai=c(0.5,0.5,0.1,0.1),mar=c(4,4,2,2)+0.1)
# radii 0 to 2 kpc
baseline <- 10^(seq(log10(601),7,length.out=100))
order.stat.pdf.two.four.first <- numeric(100)
for (i in 1:100) {order.stat.pdf.two.four.first[i] <- order.stat.pdf(baseline[i],106,106)}
order.stat.pdf.two.four.second <- numeric(100)
for (i in 1:100) {order.stat.pdf.two.four.second[i] <- order.stat.pdf(baseline[i],106,105)}
order.stat.pdf.two.four.third <- numeric(100)
for (i in 1:100) {order.stat.pdf.two.four.third[i] <- order.stat.pdf(baseline[i],106,104)}
order.stat.pdf.two.four.fourth <- numeric(100)
for (i in 1:100) {order.stat.pdf.two.four.fourth[i] <- order.stat.pdf(baseline[i],106,103)}
order.stat.pdf.two.four.fifth <- numeric(100)
for (i in 1:100) {order.stat.pdf.two.four.fifth[i] <- order.stat.pdf(baseline[i],106,102)}

plot(log10(baseline),order.stat.pdf.two.four.first,type=’l’,xlim=c(log10(600),7),ylim=c(0,10^-4.0),xlab=expression(Log~M/M[sun]),ylab=”Probability density”)
legend(“topright”,”n=106, 0 < r < 2 kpc”)
lines(log10(baseline),order.stat.pdf.two.four.third,col=”grey40″)
lines(log10(baseline),order.stat.pdf.two.four.fourth,col=”grey60″)
lines(log10(baseline),order.stat.pdf.two.four.fifth,col=”grey80″)

top.five.in.bin <- sort(pf.data$V2[which(pf.data$V1 > 0 & pf.data$V1 < 2)],decreasing=T)[1:5]
points(top.five.in.bin[1],order.stat.pdf(10^top.five.in.bin[1],106,106),pch=19,cex=0.9)
points(top.five.in.bin[3],order.stat.pdf(10^top.five.in.bin[3],106,104),pch=19,cex=0.9,col=”grey40″)
points(top.five.in.bin[4],order.stat.pdf(10^top.five.in.bin[4],106,103),pch=19,cex=0.9,col=”grey60″)
points(top.five.in.bin[5],order.stat.pdf(10^top.five.in.bin[5],106,102),pch=19,cex=0.9,col=”grey80″)

points(top.five.in.bin[2],order.stat.pdf(10^top.five.in.bin[2],106,105),pch=19,cex=0.9,col=”grey20″)
lines(log10(baseline),order.stat.pdf.two.four.second,col=”grey20″)
# radii 2 to 4 kpc
baseline <- 10^(seq(log10(601),7,length.out=100))
order.stat.pdf.two.four.first <- numeric(100)
for (i in 1:100) {order.stat.pdf.two.four.first[i] <- order.stat.pdf(baseline[i],136,136)}
order.stat.pdf.two.four.second <- numeric(100)
for (i in 1:100) {order.stat.pdf.two.four.second[i] <- order.stat.pdf(baseline[i],136,135)}
order.stat.pdf.two.four.third <- numeric(100)
for (i in 1:100) {order.stat.pdf.two.four.third[i] <- order.stat.pdf(baseline[i],136,134)}
order.stat.pdf.two.four.fourth <- numeric(100)
for (i in 1:100) {order.stat.pdf.two.four.fourth[i] <- order.stat.pdf(baseline[i],136,133)}
order.stat.pdf.two.four.fifth <- numeric(100)
for (i in 1:100) {order.stat.pdf.two.four.fifth[i] <- order.stat.pdf(baseline[i],136,132)}

plot(log10(baseline),order.stat.pdf.two.four.first,type=’l’,xlim=c(log10(600),7),ylim=c(0,10^-4.1),xlab=expression(Log~M/M[sun]),ylab=”Probability density”)
legend(“topright”,”n=136, 2 < r < 4 kpc”)
lines(log10(baseline),order.stat.pdf.two.four.third,col=”grey40″)
lines(log10(baseline),order.stat.pdf.two.four.fourth,col=”grey60″)
lines(log10(baseline),order.stat.pdf.two.four.fifth,col=”grey80″)

top.five.in.bin <- sort(pf.data$V2[which(pf.data$V1 > 2 & pf.data$V1 < 4)],decreasing=T)[1:5]
points(top.five.in.bin[1],order.stat.pdf(10^top.five.in.bin[1],136,136),pch=19,cex=0.9)
points(top.five.in.bin[3],order.stat.pdf(10^top.five.in.bin[3],136,134),pch=19,cex=0.9,col=”grey40″)
points(top.five.in.bin[4],order.stat.pdf(10^top.five.in.bin[4],136,133),pch=19,cex=0.9,col=”grey60″)
points(top.five.in.bin[5],order.stat.pdf(10^top.five.in.bin[5],136,132),pch=19,cex=0.9,col=”grey80″)

points(top.five.in.bin[2],order.stat.pdf(10^top.five.in.bin[2],136,135),pch=19,cex=0.9,col=”red”)
lines(log10(baseline),order.stat.pdf.two.four.second,col=”red”)
title(expression(“Data vs Order Statistic PDFs:”),cex.main=1)
# radii 4 to 6 kpc
baseline <- 10^(seq(log10(601),7,length.out=100))
order.stat.pdf.four.six.first <- numeric(100)
for (i in 1:100) {order.stat.pdf.four.six.first[i] <- order.stat.pdf(baseline[i],79,79)}
order.stat.pdf.four.six.second <- numeric(100)
for (i in 1:100) {order.stat.pdf.four.six.second[i] <- order.stat.pdf(baseline[i],79,78)}
order.stat.pdf.four.six.third <- numeric(100)
for (i in 1:100) {order.stat.pdf.four.six.third[i] <- order.stat.pdf(baseline[i],79,77)}
order.stat.pdf.four.six.fourth <- numeric(100)
for (i in 1:100) {order.stat.pdf.four.six.fourth[i] <- order.stat.pdf(baseline[i],79,76)}
order.stat.pdf.four.six.fifth <- numeric(100)
for (i in 1:100) {order.stat.pdf.four.six.fifth[i] <- order.stat.pdf(baseline[i],79,75)}

plot(log10(baseline),order.stat.pdf.four.six.first,type=’l’,xlim=c(log10(600),7),ylim=c(0,10^-3.8),xlab=expression(Log~M/M[sun]),ylab=”Probability density”)
legend(“topright”,”n=79, 4 < r < 6 kpc”)
lines(log10(baseline),order.stat.pdf.four.six.third,col=”grey40″)
lines(log10(baseline),order.stat.pdf.four.six.fourth,col=”grey60″)
lines(log10(baseline),order.stat.pdf.four.six.fifth,col=”grey80″)

top.five.in.bin <- sort(pf.data$V2[which(pf.data$V1 > 4 & pf.data$V1 < 6)],decreasing=T)[1:5]
points(top.five.in.bin[1],order.stat.pdf(10^top.five.in.bin[1],79,79),pch=19,cex=0.9)
points(top.five.in.bin[3],order.stat.pdf(10^top.five.in.bin[3],79,77),pch=19,cex=0.9,col=”grey40″)
points(top.five.in.bin[4],order.stat.pdf(10^top.five.in.bin[4],79,76),pch=19,cex=0.9,col=”grey60″)
points(top.five.in.bin[5],order.stat.pdf(10^top.five.in.bin[5],79,75),pch=19,cex=0.9,col=”grey80″)

points(top.five.in.bin[2],order.stat.pdf(10^top.five.in.bin[2],79,78),pch=19,cex=0.9,col=”red”)
lines(log10(baseline),order.stat.pdf.four.six.second,col=”red”)
title(expression(“Power-Law ICMF with Slope=-2″),cex.main=1)

# radii 6 to 8 kpc
baseline <- 10^(seq(log10(601),7,length.out=100))
order.stat.pdf.six.eight.first <- numeric(100)
for (i in 1:100) {order.stat.pdf.six.eight.first[i] <- order.stat.pdf(baseline[i],23,23)}
order.stat.pdf.six.eight.second <- numeric(100)
for (i in 1:100) {order.stat.pdf.six.eight.second[i] <- order.stat.pdf(baseline[i],23,22)}
order.stat.pdf.six.eight.third <- numeric(100)
for (i in 1:100) {order.stat.pdf.six.eight.third[i] <- order.stat.pdf(baseline[i],23,21)}
order.stat.pdf.six.eight.fourth <- numeric(100)
for (i in 1:100) {order.stat.pdf.six.eight.fourth[i] <- order.stat.pdf(baseline[i],23,20)}
order.stat.pdf.six.eight.fifth <- numeric(100)
for (i in 1:100) {order.stat.pdf.six.eight.fifth[i] <- order.stat.pdf(baseline[i],23,19)}

plot(log10(baseline),order.stat.pdf.six.eight.first,type=’l’,xlim=c(log10(600),7),ylim=c(0,10^-3.3),xlab=expression(Log~M/M[sun]),ylab=”Probability density”)
legend(“topright”,”n=23, 6 < r < 8 kpc”)
lines(log10(baseline),order.stat.pdf.six.eight.third,col=”grey40″)
lines(log10(baseline),order.stat.pdf.six.eight.fourth,col=”grey60″)
lines(log10(baseline),order.stat.pdf.six.eight.fifth,col=”grey80″)

top.five.in.bin <- sort(pf.data$V2[which(pf.data$V1 > 6 & pf.data$V1 < 8)],decreasing=T)[1:5]
points(top.five.in.bin[1],order.stat.pdf(10^top.five.in.bin[1],23,23),pch=19,cex=0.9)
points(top.five.in.bin[3],order.stat.pdf(10^top.five.in.bin[3],23,21),pch=19,cex=0.9,col=”grey40″)
points(top.five.in.bin[4],order.stat.pdf(10^top.five.in.bin[4],23,20),pch=19,cex=0.9,col=”grey60″)
points(top.five.in.bin[5],order.stat.pdf(10^top.five.in.bin[5],23,19),pch=19,cex=0.9,col=”grey80″)

lines(log10(baseline),order.stat.pdf.six.eight.second,col=”red”)
points(top.five.in.bin[2],order.stat.pdf(10^top.five.in.bin[2],23,22),pch=19,cex=0.9,col=”red”)
title(expression(plain(“and Limits of 600 and 10”)^7*M[sun]),cex.main=1)

# the file ‘ys’

739 2396
788 3545
803 1685
803 3724
848 2839
859 2043
904 1873
927 1685
938 2839
942 2486
972 3186
979 1596
983 2933
995 1954
1006 2133
1010 1596
1025 1865
1051 2127
1062 3724
1070 2038
1081 3102
1085 2570
1092 2660
1104 1605
1119 1780
1122 1605
1141 1954
1153 3012
1164 3281
1186 2396
1190 1954
1216 3107
1231 1605
1235 2048
1243 1780
1254 2222
1254 2933
1273 2307
1295 1954
1310 2570
1310 1596
1314 1958
1322 2396
1325 2839
1325 2928
1333 1685
1344 2664
1355 2833
1378 1865
1393 2486
1400 2226
1400 1601
1405 1865
1408 3102
1408 2839
1423 2396
1423 1685
1431 2222
1438 2043
1457 2043
1457 3012
1461 2127
1468 2316
1480 3107
1480 2839
1487 2043
1494 1690
1499 2490
1502 1601
1513 2928
1517 2396
1517 3549
1548 1690
1548 2043
1548 2222
1551 1601
1559 2570
1570 2396
1574 2133
1577 1949
1581 2575
1585 1601
1596 1690
1600 1958
1615 2222
1615 1954
1619 2127
1626 3107
1626 2401
1634 1780
1634 2570
1645 3196
1645 2127
1649 2839
1649 3012
1656 1869
1660 1690
1668 1784
1672 1605
1672 3281
1679 2043
1683 2127
1683 2316
1687 1780
1706 2486
1728 2396
1762 1605
1762 2749
1762 3455
1766 2655
1773 2133
1773 1605
1784 2316
1796 2843
1796 1873
1807 1601
1807 2043
1818 2839
1818 2655
1818 1690
1830 2490
1833 1605
1849 2043
1863 3281
1867 2839
1871 3102
1886 1873
1897 2316
1901 1605
1908 2928
1920 2127
1920 2744
1927 1601
1931 1784
1946 1601
1954 2048
1965 2401
1969 1780
1980 1690
1988 2127
1988 2401
1999 1690
1999 2222
2007 1601
2014 2401
2018 1954
2021 2311
2025 1780
2033 3365
2033 2758
2059 1784
2067 2127
2070 2486
2078 1784
2097 1596
2101 1690
2112 1784
2115 2664
2150 1690
2150 3102
2157 1869
2157 1780
2168 1954
2176 2316
2191 1873
2202 2490
2217 1869
2221 2316
2225 1780
2240 1685
2247 1869
2263 2570
2270 2490
2281 1780
2292 3012
2292 2396
2308 1780
2315 1954
2319 2486
2326 2749
2334 1780
2341 1601
2360 1869
2368 2137
2387 1869
2416 3012
2416 2570
2416 1780
2424 2048
2435 2137
2435 1605
2435 1780
2439 3724
2454 2043
2462 2396
2473 1690
2473 2316
2477 1605
2484 2744
2496 1784
2507 2127
2514 1780
2526 1601
2533 1780
2537 1958
2541 2396
2545 1690
2548 1601
2552 1780
2564 2316
2564 3012
2575 2570
2575 1954
2578 1601
2597 1780
2612 1869
2612 1685
2612 1601
2620 2222
2631 2744
2631 2048
2631 2396
2650 1685
2654 2570
2658 1780
2661 2833
2688 2664
2699 1954
2703 1605
2710 2933
2710 2490
2733 1605
2733 2396
2744 1780
2748 2127
2759 1690
2766 1601
2766 1865
2774 2043
2778 1949
2778 1685
2782 2660
2800 2133
2808 1865
2811 2311
2811 2744
2830 1780
2834 2839
2834 1869
2846 1690
2846 1605
2846 1958
2857 1869
2876 1601
2879 1780
2895 1780
2902 1869
2910 1601
2913 1690
2921 1601
2924 2316
2932 1690
2940 1780
2959 2127
2962 1690
2973 1780
2981 1690
2992 2396
2992 2575
3015 1954
3026 1605
3026 1780
3056 1605
3079 2218
3098 2311
3105 1780
3109 1605
3113 2043
3116 2664
3116 1954
3150 1601
3169 1601
3173 2048
3177 1780
3177 1690
3180 1869
3192 1954
3199 2311
3214 1601
3218 2570
3229 1780
3237 1685
3241 2490
3252 1784
3252 2048
3260 2127
3260 1605
3271 1685
3271 1605
3274 1780
3319 2839
3331 1685
3354 1780
3421 1780
3444 1780
3515 1780
3519 1605
3545 1605
3575 1780
3579 2043
3602 1780
3602 1954
3613 3375
3617 2311
3636 1605
3636 2048
3669 1685
3711 1605
3737 2570
3771 1780
3801 1865
3862 1605
3884 1605
3907 1685
3963 1601
3986 1780
4001 1605
4095 2311
4106 1685
4110 1605
4125 2311
4185 1605
4189 2396
4268 1873
4271 1596
4362 1601
4370 2133
4460 1690
4550 4603
4558 2490
4637 1601
4666 1865
4712 1859
4768 1605
4923 1685
5152 1596
5205 1690
5340 1605
5453 1690
5555 1954
5811 1601
5867 1596
6139 3281
6265 1690
6412 1685
6454 1685
6522 2222
6528 1605

Advertisements
This entry was posted in Astronomy, Bad Science, Rants. Bookmark the permalink.

3 Responses to Bordering on scientific misconduct?

  1. “The editor notes that you believe there are scientific errors in the paper you mention. In that case, the normal course would be for you to set out these errors in a paper which, if submitted to MN, would be subject to expert review.“
    Ain’t nobody got time fo dat! This is why I have my blog: to short circuit the traditional journal system which is obviously in dire need of reform.

    Except that your blog won’t short-circuit the traditional journal system. Why not just present your arguments in a paper in MNRAS?

    Yes, I think that some things about the traditional approach can be changed. However, the timescale for this is long. Don’t cut off your nose to spite your face.

    • Thanks for the feedback! The problem for me lies in the opportunity cost of writing a rebuttal paper: it’s too much precious time taken away from my main research interests. And, it’s not time well spent because citations for ‘negative’ papers are always far less than those for ‘positive’ papers, even if ultimately the consensus is that former is correct and the latter false.

  2. Pingback: Stochastic star cluster formation! | Another Astrostatistics Blog

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s