I'm trying fit an ODE model (two equations to two variables) using mixed effects because the 'replicates' per time point are actually from different patients, meaning we're combining them to build one timeseries. (annotated code below)
I get this error:
Error in chol.default((value + t(value))/2) : the leading minor of order 1 is not positive definite
I have tried just fitting one parameter (and fixing the rest) and I get the same error so it's not over fitting. Any ideas on how to fix this? many thanks
library(nlmeODE)
library(lattice)
TwoComp <- list(DiffEq=list(
dy1dt = ~ r*y1 - g*y1 , #infected cell type 1
dy2dt = ~ g*y1 - h*y2), #infected cell type 2
ObsEq=list(
c1 = ~ y1,
c2 = ~ 0),
States=c("y1","y2"),
Parms=c("r","g","h"),
Init=list(10, 1))
#make dataset:
x1<-c(rep(7,8),rep(14,20),rep(21,2)) #days
var<-c(rep(1,4),rep(2,4),rep(1,10),rep(2,10),1,2) #var = 1 fit to y1 and var =2 => fit to y2
repl<-c(rep(c(1,2,3,4),2),rep(seq(1,10),2),1,2) #rep = replicates from different patients,
#i.e. combining multiple patients' data to infer one timeseries of the infection
#measurements (cell counts):
ys<-c(31, 22, 32, 36, 51, 52, 42, 54, 46, 41, 43, 46, 44, 37, 38, 43, 42, 38, 120, 82, 50, 71, 85, 48, 34, 69, 76, 55, 47, 128)
dataBS1 <-as.data.frame(list(Time=x1,Var=var,reps=repl,Conc=ys))
SimData <- groupedData( Conc ~ Time | Var,
data = dataBS1,
labels = list( x = "Time", y = "Counts"),
units = list(x="(days)"))
plot(SimData, aspect=1/1)
testModel <- nlmeODE(TwoComp,SimData)
#fitting:
fit0 <- nlme(Conc ~ testModel(r, g, h, Time, Var),
data = SimData,
fixed = r+g+h ~ 1 , #none are fixed (not sure I entered that correctly)
random = Time~1 |reps, #replicates are random effects (i.e. come from different patients)
#also tried random = reps ~ 1 (gives same error)
start=c(r=0.5, g=0.08, h=0.04),
control=list(msVerbose=TRUE),
verbose=TRUE)
summary(fit0)
#Error given: Error in chol.default((value + t(value))/2) : the leading minor of order 1 is not positive definite
plot(augPred(fit0,level=0:1))
