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

Farkad Ezzet fezzet at aycerpc.com
Thu Nov 6 22:17:31 CET 2014


Hi,

 

I am wondering if the mixed-models community have suggestions for the
following problem: 

nlme in R fails when specifying an error model as in "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 (see below), but 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". 

Any suggestions are appreciated.

Thanks.

Farkad

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

 

####

 

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 splus are as follows:

 

>
#===========================================================================
=======================

summary(M1)

Nonlinear mixed-effects model fit by maximum likelihood

  Model: conc ~ comp1.iv(ke * exp(eta.ke), v * exp(eta.v), dose = 5, time) 

 Data: iv.D1 

        AIC      BIC   logLik 

  -421.5128 -403.633 215.7564

 

Random effects:

Formula: list(`eta.ke ~ 1` , `eta.v ~ 1` )

Level: id

Structure: Diagonal

            eta.ke      eta.v   Residual 

StdDev: 0.09409439 0.07256055 0.09725547

 

Fixed effects: list(ke + v ~ 1) 

      Value  Std.Error  DF  t-value p-value 

ke 0.103901 0.00320402 239 32.42816  <.0001

v 4.050623 0.07383155 239 54.86304  <.0001

 

Standardized Within-Group Residuals:

       Min         Q1        Med        Q3      Max 

 -2.740537 -0.6628803 0.01227503 0.6353222 2.767461

 

Number of Observations: 264

Number of Groups: 24 

>
#===========================================================================
=======================

summary(M2)

Nonlinear mixed-effects model fit by maximum likelihood

  Model: conc ~ comp1.iv(ke * exp(eta.ke), v * exp(eta.v), dose = 5, time) 

 Data: iv.D2 

        AIC       BIC   logLik 

  -712.0261 -690.5704 362.0131

 

Random effects:

Formula: list(`eta.ke ~ 1` , `eta.v ~ 1` )

Level: id

Structure: Diagonal

           eta.ke      eta.v   Residual 

StdDev: 0.1018449 0.07602354 0.09800981

 

Variance function:

Structure: Power of variance covariate

Formula:  ~ fitted(.) 

 Parameter estimates:

    power 

 1.086727

Fixed effects: list(ke + v ~ 1) 

      Value  Std.Error  DF  t-value p-value 

ke 0.101080 0.00224133 239 45.09815  <.0001

v 4.110736 0.07360580 239 55.84799  <.0001

 

Standardized Within-Group Residuals:

       Min         Q1         Med        Q3      Max 

 -3.181923 -0.5558349 -0.05090066 0.6356355 2.650502

 

Number of Observations: 264

Number of Groups: 24 

>
#===========================================================================
=======================

summary(M3)

Nonlinear mixed-effects model fit by maximum likelihood

  Model: conc ~ comp1.iv(ke * exp(eta.ke), v * exp(eta.v), dose = 5, time) 

 Data: iv.D3 

        AIC       BIC   logLik 

  -192.1091 -170.6534 102.0546

 

Random effects:

Formula: list(`eta.ke ~ 1` , `eta.v ~ 1` )

Level: id

Structure: Diagonal

           eta.ke      eta.v   Residual 

StdDev: 0.1361807 0.05127604 0.08516098

 

Variance function:

Structure: Constant plus power of variance covariate

Formula:  ~ fitted(.) 

 Parameter estimates:

    const power 

 1.206157     1

Fixed effects: list(ke + v ~ 1) 

      Value  Std.Error  DF  t-value p-value 

ke 0.104314 0.00480831 239 21.69463  <.0001

v 4.043718 0.09246226 239 43.73371  <.0001

 

Standardized Within-Group Residuals:

       Min         Q1          Med        Q3      Max 

 -2.694926 -0.6616733 -0.001634472 0.6529347 2.617859

 

Number of Observations: 264

Number of Groups: 24 

>
#===========================================================================
=======================

 

 

 


	[[alternative HTML version deleted]]



More information about the R-sig-mixed-models mailing list