phenology (version 10.1)

ExponentialRegression: Non-biased exponential regression

Description

The idea of this function is to fit a regression exponential model using MCMC because regression model using glm can produce biased outputs.
prior.default.distribution can be "dnorm", "dunif", or "dgamma". Note that if you propose dgamma prior, it will use uniform prior for r because r can be negative.
The SD model is asd*Nt+bsd. The asd parameter represents multiplicative model and the bsd parameter represents additive model. Both can be used simultaneously.
density can be dnorm or dnbinom_new (from HelpersMG package). dnbinom_new() is a negative binomial with mean and sd parametrization.

Usage

ExponentialRegression(
  data = stop("A data.frame with values"),
  colname.time = "time",
  colname.number = "numbers",
  weights = NULL,
  fitted.parameters = c(N0 = NA, r = NA, asd = NA),
  fixed.parameters = NULL,
  n.iter = c(1e+05, 1e+05),
  prior.default.distribution = "dnorm",
  density = dnorm
)

Value

Return a list with the results of exponential regression

Arguments

data

A data.frame with a column time and a column number.

colname.time

Name of the column to be used as time index.

colname.number

Name of the column to be used as number index.

weights

an optional vector of ‘prior weights’ to be used in the fitting process. Should be NULL or a numeric vector.

fitted.parameters

A named vector with the parameters to be fitted

fixed.parameters

A named vector with the parameters to be fixed

n.iter

The number of MCMC iterations.

prior.default.distribution

The default prior distribution; see description.

density

by default is dnorm but can be dnbinom_new

Author

Marc Girondot

Details

ExponentialRegression is used to fit additive, multiplicative or mixte exponential regression

Examples

if (FALSE) {
library("phenology")
t <- 1:100
N0 <- 100
r <- 0.05
y <- N0*exp(t*r)

# Multiplicative model
Nt <- rnorm(100, mean=y, sd=0.2*y)
df <- data.frame(time=t, numbers=Nt)
g <- ExponentialRegression(data=df, fitted.parameters=c(N0=NA, r=NA, asd=NA))
plot(g, parameters="r")
as.parameters(g, index="median")
# Note that if you propose gamma prior, it will use uniform prior for r
# because r can be negative
g <- ExponentialRegression(data=df, 
                           fitted.parameters=c(N0=NA, r=NA, asd=NA), 
                           prior.default.distribution="dgamma")
plot(g, parameters="r")
as.parameters(g, index="median")

# Additive model
Nt <- rnorm(100, mean=y, sd=5)
df <- data.frame(time=t, numbers=Nt)
g <- ExponentialRegression(data=df, fitted.parameters=c(N0=NA, r=NA, bsd=NA))
plot(g, parameters="r")
as.parameters(g, index="median")

# Mixt model
Nt <- rnorm(100, mean=y, sd=0.2*y+5)
df <- data.frame(time=t, numbers=Nt)
g <- ExponentialRegression(data=df, fitted.parameters=c(N0=NA, r=NA, asd=NA, bsd=NA))
plot(g, parameters="r")
as.parameters(g, index="median")

# Example with 3 common ways to perform the regression
t <- 1:100
N0 <- 100
r <- 0.05
y <- N0*exp(t*r)
out_glm <- NULL
out_mcmc <- NULL
out_nls <- NULL
for (i in 1:500) {
        print(i)
        set.seed(i)
        Nt <- rnorm(100, mean=y, sd=0.2*y)
        df <- data.frame(time=t, numbers=Nt)
        g0 <- glm(log(numbers) ~ time, data = df)
        out_glm <- c(out_glm, c(exp(coef(g0)[1]), coef(g0)[2]))
        g1 <- ExponentialRegression(data=df, n.iter=c(10000, 20000))
        out_mcmc <- c(out_mcmc, as.parameters(g1, index="median")[1:2])
        g2 <- nls(numbers ~ N0*exp(r*time), start = list(N0 = 100, r = 0.05), data = df)
        out_nls <- c(out_nls, coef(g2))
}
# In conclusion the method proposed here has no biais as compare to glm and nls fits
out_glm <- matrix(out_glm, ncol=2, byrow=TRUE)
out_mcmc <- matrix(out_mcmc, ncol=2, byrow=TRUE)
out_nls <- matrix(out_nls, ncol=2, byrow=TRUE)
mean(out_glm[, 1]); mean(out_mcmc[, 1]); mean(out_nls[, 1])
sd(out_glm[, 1])/sqrt(nrow(out_glm)); sd(out_mcmc[, 1])/sqrt(nrow(out_mcmc)); 
sd(out_nls[, 1])/sqrt(nrow(out_nls))
mean(out_glm[, 2]); mean(out_mcmc[, 2]); mean(out_nls[, 2])
sd(out_glm[, 2])/sqrt(nrow(out_glm)); sd(out_mcmc[, 2])/sqrt(nrow(out_mcmc)); 
sd(out_nls[, 2])/sqrt(nrow(out_nls))
}