[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