All,

I'm simulating the biexponential model that  is commonly used in
population PK models.  I have a rough idea of the fixed effects I'd like to use
but less of an idea of the random effects.  Using nlme, I'm interested in the
distribution of the the estimated fixed effect of a subject-level covariate at a
given sample size, and also want to to compare results under dense vs sparse
sampling per subject.

Below is a reproducible example.  It runs for the given input values but
when I experiment with different inputs for the fixed and random effects,
various error messages RE convergence and other issues prevent the simulation
from completing.  I'm using N=250 subjects, which is more than one usually
finds in these types of studies, and I'm also fitting the exact model that generated
the data.

Any suggestions RE alternations to this simulation such that the desired
assessment can be made for various input values?

Thanks!
David

--
David Afshartous, Ph.D.
Research Associate Professor
PSTATŽ: ASA Accredited Professional Statistician
Department of Biostatistics
Vanderbilt University Medical Center

##############################################
##  function to simulate the data, estimate the true model #########
## and asess estimate of fixed effect of interest               #########
###############################################

sim.results.test <- function(N, K, t.values, A1, A2, lrc1, lrc2, A1.RE.sd, A2.RE.sd,
                        lrc1.RE.sd, lrc2.RE.sd, covariate.slope, resid.sd) {
dose.c.slope <- numeric(K)
for (k in 1:K) {
Y <- NULL; X <- NULL
for (i in 1:N) {
  i.A1.RE <- rnorm(1, 0, A1.RE.sd)
  i.A2.RE <- rnorm(1, 0, A2.RE.sd)
  i.lrc1.RE <- rnorm(1, 0, lrc1.RE.sd)
  i.lrc2.RE <- rnorm(1, 0, lrc2.RE.sd)
  x.new.c <- rnorm(1, 0, 14)
  y.i <-
    (A1 + i.A1.RE)* exp( -exp (lrc1 + covariate.slope*x.new.c + i.lrc1.RE) * t.values) +
    (A2 + i.A2.RE)* exp( -exp (lrc2 + i.lrc2.RE) * t.values) +
    rnorm(1,0, resid.sd)
  y.i[y.i < 0] <- 0
  Y <- c(Y, y.i)
  X <- c(X, rep(x.new.c, length(t.values)))
}
dat.sim <- data.frame(Subject = rep(seq(1,N), each = length(t.values)),
                      time = rep(t.values,N), conc = Y, dose.c = X)
dat.sim <- groupedData(conc ~ time | Subject, data = dat.sim, order.groups=F)
dat.sim.nlme <- nlme(conc ~ SSbiexp(time, A1, lrc1, A2, lrc2),
    data = dat.sim, fixed = list(A1 ~ 1, lrc1 ~ dose.c, A2 ~ 1, lrc2 ~ 1),
    random = pdDiag(A1 + lrc1 + A2 + lrc2 ~ 1),
    start = c(A1, lrc1, covariate.slope, A2, lrc2))
dose.c.slope[k] <- fixef(dat.sim.nlme)[[3]]
  }
 dose.c.slope
}

library(nlme)
t.values.dense <- c(seq(0,4, .25), 8)

## these inputs work okay:
set.seed(101)
den.results <- sim.results(N=250, K=50, t.values.dense, A1=200, A2=41, lrc1=1.6, lrc2=0.13,
   A1.RE.sd=40, A2.RE.sd=4,lrc1.RE.sd=0.16, lrc2.RE.sd=0.01, covariate.slope=0.02, resid.sd=3)
mean(den.results) ## 0.0202; true value = 0.02
sd(den.results) ## 0.0007

## example of inputs that do not work okay; just 1 simulation run:
den.results.5 <- sim.results(N=250, K=1, t.values.dense, A1=300, A2=100, lrc1=4, lrc2=0.13,
   A1.RE.sd=40, A2.RE.sd=20,lrc1.RE.sd=.5, lrc2.RE.sd=0.02, covariate.slope=0.02, resid.sd=3)

## example what the exponential looks like
A1 <- 445; A2 <- 41
lrc1 <- 1.6; lrc2 <- -0.133
param <- c(A1, lrc1, A2, lrc2)
biexp.fun <- function(t) {
    (param[1])* exp( -exp (param[2] * t) )+  (param[3])* exp( -exp (param[4] * t))
    }
curve(biexp.fun, .1, 2)



--
David Afshartous, Ph.D.
Research Associate Professor
PSTATŽ: ASA Accredited Professional Statistician
Department of Biostatistics
Vanderbilt University Medical Center

	[[alternative HTML version deleted]]

