# Generate some synthetic data using the two-component error model and use
# different variance functions, also with fixed sigma in order to avoid
# overparameterisation in the case of a constant term in the variance function
times <- c(0, 1, 3, 7, 14, 28, 56, 120)
pred <- 100 * exp(- 0.03 * times)
sd_pred <- sqrt(3^2 + 0.07^2 * pred^2)
n_replicates <- 2
set.seed(123456)
syn_data <- data.frame(
time = rep(times, each = n_replicates),
value = rnorm(length(times) * n_replicates,
rep(pred, each = n_replicates),
rep(sd_pred, each = n_replicates)))
syn_data$value <- ifelse(syn_data$value < 0, NA, syn_data$value)
f_const <- gnls(value ~ SSasymp(time, 0, parent_0, lrc),
data = syn_data, na.action = na.omit,
start = list(parent_0 = 100, lrc = -3))
f_varPower <- gnls(value ~ SSasymp(time, 0, parent_0, lrc),
data = syn_data, na.action = na.omit,
start = list(parent_0 = 100, lrc = -3),
weights = varPower())
f_varConstPower <- gnls(value ~ SSasymp(time, 0, parent_0, lrc),
data = syn_data, na.action = na.omit,
start = list(parent_0 = 100, lrc = -3),
weights = varConstPower())
f_varConstPower_sf <- gnls(value ~ SSasymp(time, 0, parent_0, lrc),
data = syn_data, na.action = na.omit,
control = list(sigma = 1),
start = list(parent_0 = 100, lrc = -3),
weights = varConstPower())
f_varConstProp <- gnls(value ~ SSasymp(time, 0, parent_0, lrc),
data = syn_data, na.action = na.omit,
start = list(parent_0 = 100, lrc = -3),
weights = varConstProp())
f_varConstProp_sf <- gnls(value ~ SSasymp(time, 0, parent_0, lrc),
data = syn_data, na.action = na.omit,
start = list(parent_0 = 100, lrc = -3),
control = list(sigma = 1),
weights = varConstProp())
AIC(f_const, f_varPower, f_varConstPower, f_varConstPower_sf,
f_varConstProp, f_varConstProp_sf)
# The error model parameters 3 and 0.07 are approximately recovered
intervals(f_varConstProp_sf)
Run the code above in your browser using DataLab