ExponentialRegression: Non-biased exponential regression In phenology: Tools to Manage a Parametric Function that Describes Phenology and More

 ExponentialRegression R Documentation

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 marc.girondot@gmail.com

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")

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.