Dear friends - I'm having troubles with nlme fitting a simplified model
as shown below eliciting the error
Error in chol.default((value + t(value))/2) :
the leading minor of order 1 is not positive definite -
I have seen the threads on this error but it didn't help me solve the
problem.
The model runs well in brms and identifies the used parameters even with
fixed effects for TRT - but here in nlme TRT is ignored and I guess
this is not the reason for the said error
Below is the quite clumsy simulated data set and specification of call
to nlme - the start values are taken from fitted values in brms
library(ggplot2)
windows(record=TRUE)
#generate 3*10 rats - add fixed effects to the four parameters
according to the three groups - add random effects pr each rat - add
residual random effect
#Parameter values taken from Sapirstein AJP 181:330-6, 1955
set.seed(1234)
Time <- seq(1,60,by=1)
A <- 275; B <- 140; g1 <- 0.1105; g2 <- .0161
N <- 30
AA <- rep(A,30)+rnorm(30,0,30);BB <- rep(B,30)+rnorm(30,0,15) ;
gg1 <- rep(g1,30)+rnorm(30,0,0.01); gg2 <- rep(g2,30)+rnorm(30,0,0.001)
TRT <- gl(3,10*60)
levels(TRT) <- c("CTRL","DIAB","HYPER")
AA1 <- AA + c(rep(0,10),rep(10,10),rep(-10,10))
BB1 <- BB + c(rep(0,10),rep(5,10),rep(-5,10))
Gg1 <- gg1 + c(rep(0,10),rep(0.01,10),rep(-0.01,10))
Gg2 <- gg2 + c(rep(0,10),rep(0.005,10),rep(-0.005,10))
getY <- function(A,B,g1,g2) {
Y <- A*exp(-g1*Time) + B*exp(-g2*Time)
Y <- Y + rnorm(60,0,20)
}
YY <- c()
for (i in 1:N) YY <- c(YY,getY(AA1[i],BB1[i],Gg1[i],Gg2[i]))
TT <- rep(Time,N)
RAT <- gl(N,length(Time))
dats <- data.frame(RAT,TRT,TT,YY)
Dats <- dats
names(Dats)[c(3,4)] <- c("Time","Y")
dput(Dats,"dats0505.dat")
with(Dats,plot(Time,Y,pch=19,cex=.1,col=TRT))
ggplot(data=Dats,aes(x=Time,y=Y,group=RAT,col=TRT)) + geom_line()
library(nlme)
gfr.nlme <- nlme(Y ~ A*exp(-Time*g1)+B*exp(-Time*g2),
data = Dats,
fixed = A+g1+B+g2 ~1,
random = A+g1+B+g2 ~1,groups = ~ RAT,
start = c(255,115,130*1e-3,17*1e-3),
na.action = na.omit,verbose=TRUE,control = list(msVerbose = TRUE))
summary(gfr.nlme)
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
