Exponential model fitting

Let an exponential model Nt=N0*exp(r*t) be fitted. What are the different solutions taking into account different hidden hypotheses:

There are several ways to generate the data and depending the way the data are generated, different model should be used.

t <- 11:110
N0 <- 100
r <- 0.05

Error proportional of N(t); if the constant of proportionality (here 0.2) is not too large the probability to have negative value is low

y <- N0*exp(t*r)
set.seed(1)
Nt <- rnorm(100, mean=y, sd=0.2*y)
df <- data.frame(t=t, Nt=Nt)

Note that the error is not independent on the mean: heteroskedasticity

ggplot(data = df, 
       aes(x = .data[["t"]], 
           y = .data[["Nt"]])) + geom_line() +theme_bw()

Note that in log scale, the error is constant

ggplot(data = df, 
       aes(x = .data[["t"]], 
           y = log(.data[["Nt"]]))) + geom_line() +theme_bw()

Or an alternative is to model the data using a constant error and then data are homoskedastic. We use t+100 to be sure that y is never negative.

  y <- N0*exp((t+10)*r)
  set.seed(1)
  Nt <- rnorm(100, mean=y, sd=2)
  df <- data.frame(t=t+10, Nt=Nt)

We will explore 3 ways to fit the data to extract r and search for a potential bias. The first solution is to use nls() with the exponential function and then it supposed a constant SD. The second solution is to use lm (or glm) with a log version of the data. The third solution is to use MCMC and SD is modeled as a linear function of y (available in the package HelpersMG). For each, we will use the homoskedastic and the heteroskedastic data. 50000 replicates are done each time. Results are shown as mean and SE for each parameter.

Homoskedastic data

• nls

mean=100.0000 error=0.00% se=0.000050 lower=99.9999 upper=100.0001
mean=0.0500 error=0.00% se=0.00000000 lower=0.050000 upper=0.050000

• lm

mean=99.9990 error=0.00% se=0.000453 lower=99.9981 upper=99.9999
mean=0.0500 error=0.00% se=0.00000005 lower=0.050000 upper=0.050000

• MCMC


Heteroskedastic data

• nls

mean=103.3816 error=3.38% se=0.145794 lower=103.0958 upper=103.6674
mean=0.0501 error=0.29% se=0.00001435 lower=0.050117 upper=0.050173

• lm

mean=98.0508 error=1.95% se=0.024426 lower=98.0029 upper=98.0986
mean=0.0500 error=0.01% se=0.00000327 lower=0.049991 upper=0.050004

• MCMC



Commentaires

Posts les plus consultés de ce blog

Standard error from Hessian Matrix... what can be done when problem occurs

stepAIC from package MASS with AICc

Install treemix in ubuntu 20.04