ExponentialRegression: Non-biased exponential regression

In phenology: Tools to Manage a Parametric Function that Describes Phenology and More

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
)

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

Details

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

Value

Return a list with the results of exponential regression

Author(s)

Marc Girondot [e-mail]

Examples

## Not run: 
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))

## End(Not run)

phenology documentation built on Oct. 16, 2023, 9:06 a.m.