[R-sig-ME] Proportional Error model in nlme using weights

[Farkad Ezzet]

Hi,

 

I am wondering if the mixed-model community have suggestions to solve the
following problem: 

nlme in R fails when specifying an error model using "weights= varPower()".
Other error specifications, e.g. varExp and varConstPower also fail. It
seems that implementation of nlme in R is the culprit, since these error
models work perfectly fine using nlme in Splus. Here is a simple example
that attempts to model Pharmacokinetic data generated by simulation using a
one compartment intravenous dose. The code runs and estimates parameters
correctly in Splus, but it fails in R producing the error message "Error in
nlme.formula(conc ~ comp1.iv(ke * exp(eta.ke), v * exp(eta.v),  : step
halving factor reduced below minimum in PNLS step". 

 

Here is the code in R:

 

library(nlme)

 

comp1.iv =

  function(ke, v, dose, time)

  {

    (dose/v) * exp(- ( ke * time ) ) 

  }

 

# function to simulate data

comp1.iv.Data = function(nsub, ke, v, cv.ke, cv.v, cv.error, d1, d2, time,
dose)

{

  set.seed(123)

  

  nsub0 <- nsub

  sub0   <- 1:nsub0

  time0 <- time

  dose0 <- dose

  

  eta.ke <- rnorm(nsub0,0,cv.ke)

  eta.v  <- rnorm(nsub0,0,cv.v)

  

  ke0 <- ke * exp(eta.ke)

  v0  <- v  * exp(eta.v)

  

  nobsid <- length(time0)*length(dose0)

  sub <- rep(sub0, each=nobsid)

  nobs   <- nobsid*nsub0

  time <- rep(time0,length(dose0)*nsub0)

  dose <- rep(rep(dose0, each=length(time0)),nsub0)

  

  ke <- rep(ke0, each=nobsid)

  v  <- rep(v0 , each=nobsid)

  

  x0 <- comp1.iv(ke,v,dose,time)

  

  conc <- x0 + (d1 + d2*x0) * rnorm(length(x0),0,cv.error)

  

  D = data.frame(id=sub, dose=dose, time= time, ke=ke, v=v, x0=x0,conc=conc)


  D = groupedData(conc~time| id, data =D)

  return(D)

} # end of comp1.iv.Data

 

# create data set with different error models

time = c(0,1,2,3,4,6,8,10, 12,18, 24)

dose = c(5)

 

iv.D1   = comp1.iv.Data(nsub=24,ke=.1, v=4,cv.ke=.1, cv.v=.1,
cv.error=0.1,d1=1, d2=0, time, dose) # additive

iv.D2   = comp1.iv.Data(nsub=24,ke=.1, v=4,cv.ke=.1, cv.v=.1,
cv.error=0.1,d1=0, d2=1, time, dose) # proportional

iv.D3   = comp1.iv.Data(nsub=24,ke=.1, v=4,cv.ke=.1, cv.v=.1,
cv.error=0.1,d1=1, d2=1, time, dose) # additive & proportional 

 

mean(iv.D1$conc)

mean(iv.D2$conc)

mean(iv.D3$conc)

 

# Run nlme with additive error using iv.D1

M1 = nlme(conc ~ 

            comp1.iv(ke*exp(eta.ke), v*exp(eta.v), dose=5, time) ,

          data=iv.D1, 

          random = pdDiag(list(eta.ke + eta.v ~ 1)),

          fixed = list(ke+v~ 1),

          start = c(.1, 4), verbose = T)

summary(M1)

 

# Run nlme with proportional error using iv.D2

M2 = nlme(conc ~ 

            comp1.iv(ke*exp(eta.ke), v*exp(eta.v), dose=5, time) ,

          data=iv.D2, 

          random = pdDiag(list(eta.ke + eta.v ~ 1)),

          fixed = list(ke+v~ 1),

          start = c(.1, 4), verbose = T, weights=varPower(1))

 

# Run nlme with additive & proportional error using iv.D3

M3 = nlme(conc ~ 

            comp1.iv(ke*exp(eta.ke), v*exp(eta.v), dose=5, time) ,

          data=iv.D3, 

          random = pdDiag(list(eta.ke + eta.v ~ 1)),

          fixed = list(ke+v~ 1),

          start = c(.1, 4), verbose = T, weights=varConstPower(const=1,
fixed=list(power=1)))

 

#===========================================================================
=======================

Results from the splus output are as follows:

 

 

Farkad Ezzet, MSc, PhD

Office  973 665-1447

Mobile 201 738-3826

 <mailto:fezzet at aycerpc.com> fezzet at aycerpc.com

 <http://www.aycerpc.com/> www.aycerpc.com

 <http://www.linkedin.com/pub/farkad-ezzet/11/28b/93a>
http://www.linkedin.com/pub/farkad-ezzet/11/28b/93a

 

 

 

 


	[[alternative HTML version deleted]]