Confidence limits for the inverse of an estimated value

[Calvadosser]

How are the data points measured? Do you have separate instruments to measure \frac{dC(t)}{dt} and C(t)? Or are you computing \frac{dC(t)}{dt} from differences in the C(t) data?

(I'm wondering why you chose to use regression instead of fitting a curve to the C(t) data.)

First, thanks for all the responses to my original question.

I agree completely that, if one had observations of C(t), where C(t) satisfies dC(t)/dt = -λ C(t), starting from C(0), then you'd simply find the exponential that fitted the observed points.

I tried to include sufficient information to ask my question but I tried to avoid including a load of detail that I thought would simply be tedious and irrelevant to readers.

I am trying to understand a very simple model representing the dynamics of atmospheric concentration of CO2 that someone sent me in a paper they had written. I am trying as a first step to reproduce their work. I think that they have greatly over-estimated the accuracy of the estimate for 1/λ, so I want to be doubly careful with this aspect.

I'm still simplifying things here but I hope this won't result in confusion.

Here are my symbols:

t = time after 1960 (yr)
C(t) = excess atmospheric CO2 over the equilbrium level, at time t (Gt).
Fa(t) = rate of emission of anthropogenic CO2 (Gt/yr)
λ = coefficient relating C(t) and the resulting additional absorption rate over the equilibrium absorption rate. (This assumes linearity and only one sink for CO2, hence a 1st order diff eqtn.)

So I have an equation dC(t)/dt = -λ C(t) + Fa(t) .. .. .. .. (1)

As data, I have:

Values of C(0), C(1),... observed values for excess CO2, derived from data published by the Mauna Loa observatory.

Values of Fa(0),Fa(1),... from published estimates of annual CO2 release from fossil fuel burning and land use change.

If I change equation (1) around, I have dC(t)/dt - Fa(t) = -λ C(t). .. .. .. (2)

Using the approximation dC(t)/dt = [C(t) - C(t-1)]/1 and the estimates for Fa(t), I get 50 (rather scattered) points on the left of (2) and I have 50 points on the right. This gives me my scatter plot. A regression analysis gives me an estimate for λ.

From the estimate for λ, I can compute various things of interest, including T = 1/λ, the time to reach equilibrium, in the hypothetical case of anthopogenic CO2 releases ceasing.

I hope this gives the background and explains why I could not just fit a curve to C(t). If there is a better way of doing things, I'd like to know about it.

Thank you again to all for the comments.

 

[mfb]

I would use the approximation dC/dt = 1/2 (C(t+1)-C(t-1)), but chances are good that it does not change much (something you can add as cross-check).

\lambda = \frac{\frac{dC}{dt}-Fa}{C(t)}
This gives you a direct estimate for λ, and if you take the inverse value you get a direct estimate for T in each year. If your values do not have individual (and different) errors, the arithmetric average should be the best way to estimate your value.

 

[DrDu]

The situation, as I've come to understand it from viraltux, is that if \hat{\lambda} is an MLE for parameter \lambda then (under some fairly unrestrictive conditions on f), f(\hat{\lambda}) is an MLE for f(\lambda).
In the OP, the slope, \hat{\lambda}, is being obtained by standard linear regression. We don't know anything about the distribution of the errors, but I doubt this constitutes an MLE, so I remain sceptical of the validity of using 1/\hat{\lambda} for 1/\lambda.
The other complication is that the choice of estimator should really be based on minimising a defined cost function.

The MLE is invariant under reparametrizations but it is only asymptotically efficient and is also only asymptotically unbiased.
On the other side, a minimum variance unbiased estimator for finite sample size is usually not invariant under reparameterizations.

 

[DrDu]

Even when one uses MLE estimates the transformed confidence intervals are not identical to the confidence intervals for 1/ parameter, especially if the CI's are the shortest possible for one parameter they transformed CI's are not the shortest for 1/parameter.

 

[Stephen Tashi]

Even when one uses MLE estimates the transformed confidence intervals are not identical to the confidence intervals for 1/ parameter, especially if the CI's are the shortest possible for one parameter they transformed CI's are not the shortest for 1/parameter.

I don't know if this thread is helping the original poster Calvadosser, but it is very informative!

In that spirit, I'll ask about the technical definition of a "shortest possible" confidence interval. Suppose we have some algorithm that takes the sample data and computes an interval and the algorithm has the property that there is a 95% chance that the population parameter we are estimating is in that interval. There is no requirement that all the intervals the algorithm produces have the same length. Only If we take particular sample data do we get an interval of a particular length and (from the frequentist point of view) we don't know that there is a 95% chance that the population parameter is in that particular interval. If I am comparing two such algorithms, I can compare the expected lengths of the intervals they produce. As I interpret your statement, it refers to expected lengths of confidence intervals. Is that correct?

 

[Stephen Tashi]

I would use the approximation dC/dt = 1/2 (C(t+1)-C(t-1)), but chances are good that it does not change much (something you can add as cross-check).

And there are quite a variety of additional mathematical refinements that could be introduced. Real world problems can be modeled different sets of mathematical assumptions and the only way that people get confident about their answers is to work them many different ways using a variety of different assumptions and get roughly the same answer everytime. (I've observed that this is particularly true in bureacracies. Most managers don't work by trying to pick the best answers. They force analysts to solve problems many different ways and they pick the answer that comes up most often.) Calvadosser must let us know if he's interested in fancier math.

How to treat the C(t) data is an interesting question. I looked at the data on the web and the higher resolution data is (of course) wavy. So what will C(n) be? I suppose you could set it to be the moving average taken over a year's time prior to year n.

\lambda = \frac{\frac{dC}{dt}-Fa}{C(t)}
This gives you a direct estimate for λ
and if you take the inverse value you get a direct estimate for T in each year. If your values do not have individual (and different) errors, the arithmetric average should be the best way to estimate your value.

From
- \lambda = \frac{\frac{dC}{dt}-Fa}{C(t)}

We can reduce the data
T_i= \frac{C(i)}{Fa_i - \frac{dC(i)}{dt}}

A statistical concern is that the T_i and T_{i+1} are not independent random samples since computing the values of \frac{dC(i+1)}{dt} and \frac{dC(i)}{dt} both involve using the datum C(i).

I wonder if Calvadosser's method of using regression cleverly takes care that concern.

 

[viraltux]

Hi All,

Would you mind to answer a question? Being

OK: The OP's procedure is absolutely fine.
Not OK: The OP's procedure is either wrong or misses so further analysis.
IDK: I don't know / I am not sure.

Where would you situate yourself? Also, could you tell about your background? I think this would help to situate the thread's audience. For instance, in my case it would be:

OK | viraltux | Statistician / CS

 

[Calvadosser]

I would use the approximation dC/dt = 1/2 (C(t+1)-C(t-1)), but chances are good that it does not change much (something you can add as cross-check).

\lambda = \frac{\frac{dC}{dt}-Fa}{C(t)}
This gives you a direct estimate for λ, and if you take the inverse value you get a direct estimate for T in each year. If your values do not have individual (and different) errors, the arithmetric average should be the best way to estimate your value.

Thank you for the suggestion. I will give it a try as soon as I can (but not for a week or so because of things that have to take priority).

I was using C(t)-C(t-1) because that then covered the same time interval as Fa(t). Maybe I can use 1/2(Fa(t+1)+Fa(t-1)) to match up with 1/2 (C(t+1)-C(t-1)). I'll think it through.

 

[DrDu]

I don't know if this thread is helping the original poster Calvadosser, but it is very informative!

In that spirit, I'll ask about the technical definition of a "shortest possible" confidence interval.

I cannot give you spontaneously a general definition, but, in the relevant context of estimation of CI's for a single parameter via Maximum Likelihood, the situation is easy to grasp:
If \lambda is the true and unknown parameter value and \hat{\lambda} it's maximum likelihood estimate, furthermore \hat{\sigma} it's estimated variance.
Then the two sided CI's are \lambda_{u/l}=\hat{\lambda}\pm z \hat{\sigma} were z is the quantile of the normal distribution corresponding to a given strength \alpha. On the mean this interval will be of length 2z\sigma which is the shortest possible interval which covers 1-\alpha of the normal distribution.
Now if you are interested not in \lambda but in 1/\lambda, then the transformed CI's 1/\lambda_{u/l} do not span the shortest interval.
This can be seen as follows. By error propagation (or delta method, as statisticians tend to call it), \hat{\sigma}_{1/\lambda}=1/\hat{\lambda}^2 \hat{\sigma}_\lambda.
Hence the shortest CI's (on the mean) for 1/\lambda are
1/\hat{\lambda}\pm z \hat{\sigma}_{1/\lambda} which does not coincide with the transformed CI's of \lambda.