1
$\begingroup$

Suppose you are running an experiment where you measure a (positive) variable $A$ of a given physical system at multiple times, say $t = \{ 0, 10, 20, 45, 90, 120, 180 \} \ min$, and calculate $y(t) = ln(A(t))-ln(A(0))$.
From the theory you expect $A(t) = A(0) \cdot exp(-k \cdot t)$, hence $y(t) = - k \cdot t$.
Your goal is to estimate the 'true' value of $k$, i.e. $k_{true}$.

If you knew nothing about the uncertainty on the individual measurements of $A(t)$, and were allowed to run the experiment only once, I guess the only option would be to do a linear regression on the experimental $(t,y(t))$ points (excluding $t=0$, of course, as it is forced) and report $k = -slope$, and its standard error $SE_k$, obtained from the regression itself.

However, suppose you do know, from historical data, that the standard deviation on $ln(A(t))$ is constant, independent from $ln(A(t))$ itself, and equal to $s \approx 0.023$.
To be completely clear, by this I mean that if the measurement of $A(t)$ is repeated multiple times, its natural logarithm is found to be normally distributed around $ln(A(t))_{true}$, with a standard deviation equal to $0.023$, which does not change for different values of $t$ or of $k_{true}$.
[FYI, this is a 'real world' example: we did observe this in actual experiments].

I ran some simulations in R (see below) to find the distribution of $k$ estimated by regression, under the above hypothesis of constant standard deviation on each $ln(A(t))$.
It turned out that $k$ was normally distributed, with standard deviation $sd_k \approx 0.00016$, regardless of the value of $k_{true}$. Instead, $SE_k$ was, not unexpectedly, all over the place.

So here is my doubt: if you had to report an estimate for $k$ and its uncertainty, would you use the value of $SE_k$ from the regression (which for a single experiment could be quite off the mark), or would you use the knowledge of the standard deviation on $ln(A(t))$ to derive or simulate a 'theoretical' $sd_k$?

Here is the R code:

k_true = 0.012
cat("\nk_true =",k_true)
t = c(0,10,20,45,90,120,180)
A_true = exp(-k_true*t)

s = 0.023
n_rep = 5000

set.seed(97531)

ln_A_measured_rep <- sapply(A_true,function(x) rnorm(n = n_rep,mean = log(x),sd = s))

k_repls <- apply(ln_A_measured_rep,1,function(r) {
  d <- data.frame(x = t[-1], y = (r-r[1])[-1])
  l <- lm(y ~ x, d)
  sml <- summary(l)
  km <- -1*sml$coefficients[2]
  SE_km <- sml$coefficients[4] 
  return(c(km,SE_km))
  }
)

sd_k <- sd(k_repls[1,])

cat("\nsd_k =", sd_k)

hist(k_repls[2,],main="SE_k",xlab="SE_k")
abline(v=sd_k,col=2)

You can test that by changing $k_{true}$, $sd_k$ does not change.

I hope at least some of this makes sense; any help will be greatly appreciated.

Thanks!

EDIT: formula for the value of $sd_k$

After long calculations, I found that for any $N >= 2$, $N$ being the number of distinct time points at which $y(t)$ is measured:

$sd_k = \frac s {\sqrt {\sum_{i=1}^N {t_i^2} - \frac 1 N (\sum_{i=1}^N {t_i})^2}}$

In the above example:

$sd_k = \frac {0.023} {\sqrt {(10^2+20^2+45^2+90^2+120^2+180^2) - \frac 1 6 (10+20+45+90+120+180)^2}} = \frac {0.023} {\sqrt {57425 - 36037.5}} \approx 0.00016$

The derivation of the formula was theoretical; still, I tested it numerically on a large number of examples, changing $k_{true}$, $s$ and the number and positions of the time points sampled. It seems to work.

$\endgroup$
  • $\begingroup$ Is $A(0)$ measured? If so, an important consideration is that all your $y$ values are interdependent, so the regression procedure ought to account for that. $\endgroup$ – whuber Feb 2 at 17:13
  • $\begingroup$ you mentioned you had calculated $y(t) = ln(A(t)) - ln(A(0))$ but are calculating a linear regression. Shouldn't that be $y(t) = ln(A(t) - A(0))$ if you're trying to fit a generalized linear model, with $ln()$ as the link function? $\endgroup$ – Jack Feb 2 at 17:23
  • $\begingroup$ @whuber : indeed, it is measured, and I do take that into account, by eliminating the $t=0$ point in the regression, as $y=0$ by definition for that point, so it's not free. $\endgroup$ – user6376297 Feb 3 at 19:29
  • $\begingroup$ @Jack : $ln(A(t)-A(0)) = ln(A(0) \cdot (e^{-k \cdot t}-1)) = ln(A(0)) + ln(e^{-k \cdot t}-1))$. Then I would not know how to get $k$ from the regression. BTW, I spotted a mistake in my formula for $y(t)$, I will correct it. $\endgroup$ – user6376297 Feb 3 at 19:36
  • 1
    $\begingroup$ @whuber : unfortunately these are the data we get from the company that runs the experiment. We don't get A(t) (including A(0)), we just get the ratio's A(t)/A(0) for each t. My method consists in taking the $ln$ of each ratio (which is $y(t)$) and linearly relating it to $t$, because the theory says indeed that the relationship between the two is linear. If this method is wrong, what would be the correct one then? What would you do with these data? $\endgroup$ – user6376297 Feb 4 at 17:04

Your Answer

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

Browse other questions tagged or ask your own question.