7
$\begingroup$

Background: In planetary astronomy, the only method we have to estimate the age of a surface of a solid object in the solar system (other than Earth) is to identify craters. We then compare these to an established chronology that links the number density of craters with diameters $D \ge 1$ km (often written as $N(1)$) with absolute ages as established from Apollo and Luna sample return missions.

The function linking age $T$ with $N(1)$ has historically been fit to the form:

$N(1) = \alpha (\exp(-\beta T)-1)+\gamma T$

where $\alpha$, $\beta$, and $\gamma$ are parameters of the fit. Qualitatively, this expresses the idea that early on, there was an exponential decrease in cratering, and after a certain time, there was a linear rate up to the present day.

This fit is constrained by around 15 points at most, it's completely unconstrained (other than the maximum mass of an object) for ages older than about 3.95 billion years, and there are no data points between about 1 and 3 billion years ago.

Question: I'm re-doing some of this older work as part of my research now. I'm using established radiometric ages and re-doing the crater counts to fit a new function. My issue is that, from my new data, this version of the fit function "looks" bad. A research group in 2007 suggested more of a quadratic decline in the recent past so revised the function as:

$N(1) = \alpha (\exp(-\beta T)-1)+\gamma T^2 + \delta T$

Adding the quadratic term greatly increases how good the fit "looks," but with only 11 points, I'm worried that I could be getting into an area where I'm just getting an improved fit because I'm adding more free parameters. And, probably, a reviewer is going to want something more quantitative, especially because the paper that everyone uses for this has the first version.

Note that the data range is roughly $0 < T < 4.5$ and $10^{-6} < N(1) < 10^0$, and while they're somewhat well-behaved, it's non-linear (as the function would suggest).

I've calculated the reduced $\chi^2$ of each version, and it's 2.8 versus 1.7, respectively. But, this is a highly non-linear function, and based on some reading I've done, the $\chi^2$ may not be a meaningful metric for determining how good the fit really is. Someone suggested that I do an incomplete gamma function test to determine if the $\chi^2$ is meaningful, but the results of that are 0.00018 versus 0.017 ... and I have no idea what those mean other than I had a thought that the larger (if still minuscule) number was better indicating the quadratic was more meaningful.

So ... what's a (or several) good way to determine whether one fit function is statistically better than another with this kind of data? $\chi^2$? Or something else?

Data: As an edit, I was asked to share the data 'cause they're just a few points ... I'm going to have to decline on that because the data are unpublished, unreviewed, unvetted, and this is a very scoopable project that I'm gunning for publication in Science ... I'll be presenting at a conference in two weeks that will hopefully give me a better idea of if I'm insane or not.

What I can share at this point are published data that mine are an update of.

+------+-------+
|    T | N(1)  |
+------+-------+
|-3.92 |0.034  |
+------+-------+
|-3.84 |0.057  |
+------+-------+
|-3.85 |0.037  |
+------+-------+
| -3.8 |0.009  |
+------+-------+
|-3.75 | 0.01  |
+------+-------+
|-3.58 |0.0064 |
+------+-------+
|-3.41 |0.0033 |
+------+-------+
| -3.3 |0.0032 |
+------+-------+
|-3.22 | 0.003 |
+------+-------+
|-3.15 |0.0036 |
+------+-------+
| -0.8 |0.0013 |
+------+-------+
|-0.109| 9e-05 |
+------+-------+
|-0.053|4.4e-05|
+------+-------+
|-0.025|2.1e-05|
+------+-------+

From these or very slight variations, the parameters everyone in the field uses are $\alpha=5.44\cdot10^{-14}$, $\beta=6.93$, $\gamma=8.38\cdot10^{-4}$. To my revision and with the new fit, my parameters are $\alpha=2.61\cdot10^{-35}$, $\beta=19.4$, $\gamma=1.46\cdot10^{-4}$, $\delta=1.71\cdot10^{-3}$. And, since I am presenting this in two weeks from this morning, I do have a 2-page abstract that discusses this in more detail and does contain a graph or three with the function shown.

Note that, physically speaking, there's no good reason to expect it to be linear after a certain time versus a quadratic. But the exp-lin model is just something that people have come up with to reasonably and simply fit the data available. It's very likely that a more complex model would be needed if we had infinite resolution and there would be Lorentzian spikes all over the place due to asteroid family formation. But those are well below the noise at this point.

$\endgroup$
  • 3
    $\begingroup$ Look into (either here or on the web generally) AIC, BIC, SBC and related terms that those will identify. These are all penalized measures of fit - where the improvement from adding free parameters is penalized for added terms. $\endgroup$ – Peter Flom Sep 7 '12 at 10:49
  • $\begingroup$ Given that you're discussing only $11$ (or is it $15$?) points, could you share them with us? It would also be interesting to know what the fitted values of the parameters are. Is there any justification in physical or cosmological theory for that $T^2$ term? $\endgroup$ – whuber Sep 7 '12 at 18:11
  • $\begingroup$ Whuber - some data added, and comments on physical meaning. $\endgroup$ – Stuart Robbins Sep 7 '12 at 19:43
3
$\begingroup$

As Peter Flom suggests given your models you have a likelihood function and these information criteria can compare models based on their likelihood function with penalties for parameters used that can led to a "best" fit when the information criterion is maximized. AIC and BIC are of the form -2 log likelihood + penalty and differ in the choice of penalty. So best by the criteria means maximum. This might help in the sense that it gets around picking an overfitted model. But it is possible that you are still left with two models that are close. Should you really pick the one that is maximum? The question of how to decide between them remains. It may be that it really requires a much larger data set to distinguish between the two.

$\endgroup$
  • $\begingroup$ I'm looking at AIC_C at the moment, and running into confusion (as usual) with the likelihood function. I already have my fit parameters, so am I just doing ln(L) where L is the sum of ln(my fit function, either exp-lin or exp-quad) summed for each data point? All the derivative stuff is just for finding the max which I don't need to do here since I have the fit parameters, right? $\endgroup$ – Stuart Robbins Sep 8 '12 at 4:28
  • $\begingroup$ That is a modification to the original AIC. I am not familiar with it. But this entire class of criteria are all similar and all involve maximizing the likelihood - a penalty. I am think of the form -2 log likelihood because that is how a likelihood ratio test comparing two models is expressed. That is because that test statistic has an asymptotic chi square distribution under the null hypothesis. So I think that is why the criterion (e.g. AIC) is expressed in this form. $\endgroup$ – Michael Chernick Sep 8 '12 at 4:50
  • $\begingroup$ Thanks, AICc (not AIC_C as I said originally) is just a modification that further penalizes you due to additional fit parameters, as you said. My question is more on how to evaluate that likelihood part -- my reference is the "Data Reduction and Error Analysis" book by Bevington & Robinson and I'm trying to follow them, but I'm stuck at exactly what I'm evaluating in the sum (sum 'cause it's ln(L)). Is this a separate question to ask? I didn't happen to see it when doing a search on this site. $\endgroup$ – Stuart Robbins Sep 8 '12 at 20:02
  • $\begingroup$ @StuartRobbins The likelihood is just a given function of the parameter(s) using the observed data. You find its maximum which gives the maximum likelihood estimates. For the Akaike criterion plug in the maximum likelihood estimate(s) inot the likelihood to get the value of AIC. $\endgroup$ – Michael Chernick Sep 8 '12 at 20:25
  • $\begingroup$ Okey dokey, I've talked with a statistics master's student who I'm working with on another project and she's agreed to help with this and is thinking about doing some jack-knifing and boot-strapping and other tests to see what's what. Thanks! $\endgroup$ – Stuart Robbins Oct 11 '12 at 21:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.