Re: [R] Error in chol.default

[J C Nash]

The error msg says it all if you know how to read it.

> When I run the optimization (given that I can't find parameters that
> fit the data by eyeball), I get the error:
> ```
> Error in chol.default(object$hessian) :
>   the leading minor of order 1 is not positive definite


Your Jacobian (derivatives of the residual function w.r.t. the parameters)
is singular -- rather spectacularly so.

Try the short addition to your code to use analytic derivatives:

rich = function(p, x) {
  a = p["curvature"]
  k = p["finalPop"]
  r = p["growthRate"]
  y = r * x * (1-(x/k)^a)
  return(y)
}
ricky = function(p, x, y) p$r * x * (1-(x/p$k)^p$a) -y
# values
Y <- c(41,   41,   41,   41,   41,   41,   45,   62,  121,  198,  275,
   288,  859, 1118)
X <- 1:14
a = 1
k = 83347
r = 40
fit = rich(c(curvature=a, finalPop=k, growthRate=r), X)
plot(Y ~ X,
 col = "red", lwd = 2,
 main = "Richards model",
 xlab = expression(bold("Days")),
 ylab = expression(bold("Cases")))
points(X, fit, type = "l", lty = 2, lwd = 2)
library("minpack.lm")
o <- nls.lm(par = list(a=a, k=k, r=r), fn = ricky, x = X, y = Y)
summary(o)
print(o1)
library(nlsr)
xy=data.frame(x=X, y=Y)
rcky2 = y ~ r * x * (1-(x/k)^a) -y
o1 <- nlxb(rcky2, start = c(a=a, k=k, r=r), data=xy, trace=TRUE)


You should find

> o1
nlsr object: x
residual sumsquares =  3161769  on  14 observations
after  8Jacobian and  13 function evaluations
  namecoeff  SE   tstat  pval  gradient
JSingval
a294.113NA NA NA   0   
31.86
k  83347NA NA NA   0
   0
r74.9576NA NA NA  -6.064e-05
   0
>
Note the singular values of the Jacobian. Actual zeros!

Even so, nlsr had managed to make some progress.

nls() and nls.lm() use approximate derivatives. Often that's fine (and it is 
more flexible
in handling awkward functions), but a lot of the time it is NOT OK.

JN

__
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.

[R] Error in chol.default

[Luigi Marongiu]

Hello,
I am trying to fit a Richards model to some cumulative incidence data
of infection. I got this example:
```
rich = function(p, x) {
  a = p["curvature"]
  k = p["finalPop"]
  r = p["growthRate"]
  y = r * x * (1-(x/k)^a)
  return(y)
}
ricky = function(p, x, y) p$r * x * (1-(x/p$k)^p$a) -y
# values
Y <- c(41,   41,   41,   41,   41,   41,   45,   62,  121,  198,  275,
 288,  859, 1118)
X <- 1:14
a = 1
k = 83347
r = 40
fit = rich(c(curvature=a, finalPop=k, growthRate=r), X)
plot(Y ~ X,
 col = "red", lwd = 2,
 main = "Richards model",
 xlab = expression(bold("Days")),
 ylab = expression(bold("Cases")))
points(X, fit, type = "l", lty = 2, lwd = 2)
library("minpack.lm")
o <- nls.lm(par = list(a=a, k=k, r=r), fn = ricky, x = X, y = Y)
summary(o)
```
When I run the optimization (given that I can't find parameters that
fit the data by eyeball), I get the error:
```
Error in chol.default(object$hessian) :
  the leading minor of order 1 is not positive definite
```
What is this about?
Thank you




-- 
Best regards,
Luigi

__
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.

Re: [R] error in chol.default((value + t(value))/2) : , the leading minor of order 1 is not positive definite

[David Winsemius]

> On May 5, 2018, at 1:19 AM, Troels Ring  wrote:
> 
> 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,

Your fixed and random formulae look the same. That would seem to create 
problems, at leas the way I understand mixed models analysis. At any rate this 
is much more likely to get expert eyes (which mine definitely are not) on the 
problem if it were posted to the mixed models list.

> 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.

David Winsemius
Alameda, CA, USA

'Any technology distinguishable from magic is insufficiently advanced.'   
-Gehm's Corollary to Clarke's Third Law

__
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.

[R] error in chol.default((value + t(value))/2) : , the leading minor of order 1 is not positive definite

[Troels Ring]

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.