tests/varConstProp.R

In nlme: Linear and Nonlinear Mixed Effects Models

times <- c(0, 1, 3, 7, 14, 28, 56, 90, 120)
k_in  <- c(0.1, 0.15, 0.16, 0.18, 0.2)
lrc_in <- mean(log(k_in))
lrc_sd_in <- sd(log(k_in))
n_sample <- length(times) * 2
const <- c(s1 = 1,    s2 = 1)
prop  <- c(s1 = 0.07, s2 = 0.03)
const_in <- rep(const, times = c(3, 2) * n_sample)
prop_in  <- rep(prop,  times = c(3, 2) * n_sample)

pred <- as.numeric(sapply(k_in, function(k) rep(100 * exp(- k * times), each = 2)))

set.seed(123456L)
d_syn <- data.frame(
    time  = rep(times, 5, each = 2),
    ds    = rep(paste0("d", 1:5), each = n_sample),
    study = rep(c("s1", "s2"), times = c(3 * n_sample, 2 * n_sample)),
    value = rnorm(length(pred), pred, sqrt(const_in^2 + pred^2 * prop_in^2))
)

library(nlme)
f_nlsList <- nlsList(value ~ SSasymp(time, 0, 100, lrc) | ds,
                     data = d_syn,
                     start = list(lrc = -3))

(fm_tc_study <-
    ## suppressWarnings( ## as the fit seems to be overparameterised
    nlme(f_nlsList,
         weights = varConstProp(form = ~ fitted(.) | study),
         control = list(sigma = 1))
    ## )
)

(ints <- intervals(fm_tc_study))

## Check if intervals include input used for data generation
stopifnot(exprs = {
  ints$fixed["lrc", "lower"] < lrc_in
  ints$fixed["lrc", "upper"] > lrc_in
  all.equal(c(
      ints$fixed["lrc", "est."],
      ints$reStruct$ds["sd(lrc)", "est."]),
    c(lrc_in, lrc_sd_in), tol = 1e-2) # diff. 0.00797 seen
  ints$varStruct[1:2, "lower"] < const
  ints$varStruct[3:4, "lower"] < prop
  ints$varStruct[1:2, "upper"] > const
  ints$varStruct[3:4, "upper"] > prop
  all.equal(
    as.numeric(ints$varStruct[c("prop.s1", "prop.s2"), "est."]),
    as.numeric(prop), tol =  0.15) # diff. 0.062  seen
})

## We do not get warnings if we fix the constant part of the error model
fm_tc_study_CF <-
    nlme(f_nlsList,
         weights = varConstProp(form = ~ fitted(.) | study,
                                fixed = list(const = c(s1 = 1, s2 = 1))))
summary(fm_tc_study_CF)
(ints_cf <- intervals(fm_tc_study_CF))
stopifnot(
  all.equal(
    as.numeric(ints_cf$varStruct[, "est."]),
    as.numeric(prop), tol =  0.15) # diff.  0.1168  seen
)