[R] try (nls stops unexpectedly because of chol2inv error

[Mike Marchywka]

----------------------------------------
> Date: Mon, 8 Nov 2010 06:18:49 -0800
> From: monte.shaffer at gmail.com
> To: r-help at r-project.org
> Subject: [R] try (nls stops unexpectedly because of chol2inv error
>
> Hi,
>
> I am running simulations that does multiple comparisons to control.
>
> For each simulation, I need to model 7 nls functions. I loop over 7 to do
> the nls using try
>
> if try fails, I break out of that loop, and go to next simulation.
>
> I get warnings on nls failures, but the simulation continues to run, except
> when the internal call (internal to nls) of the chol2inv fails.
>
> ============================================
> Error in chol2inv(object$m$Rmat()) :
> element (2, 2) is zero, so the inverse cannot be computed
> In addition: Warning messages:
> 1: In nls(myModel.nlm, fData, start = initialValues, control =
> nls.control(warnOnly = TRUE), :
> number of iterations exceeded maximum of 50
> 2: In nls(myModel.nlm, fData, start = initialValues, control =
> nls.control(warnOnly = TRUE), :
> singular gradient
> ===========================================================
>
> Any suggestions on how to prevent chol2inv from breaking my simulation...

Since no one else has answered, let me supply some thoughts and google hits.

I'm not sure what your question is- the error message suggests the
matrix has no inverse as in A*A-1 =I can't be found- usually these things happen
because the data is not a good fit to the model. Is the message not
literally true as in you know that A has an inverse? It does seem
you posted a good complete example but it may take a bit of effort
for someone to debug.

The reason it is non-invertible probably has to do with the gradient issue,
in any case some good hits on google like this may help,

https://stat.ethz.ch/pipermail/r-help/2008-March/158329.html

( http://www.google.com/?#hl=en&q=r+nls+singular+gradient&fp=1 ) 

Personally I
tend to use SVD in my c++ code since it is the only method I know
that provides a good diagnostic on how close I came to having an
ill posed model. In your case, presumably either your model or data or code is
creating an exactly singular matrix, this may be easier to find
than the almost singular situations that often create odd results :) 
I would just ask however if anyone
has more thoughts on inverting mtricies for model fits as someone
previously mentioned that R uses QR decomposition for one task to qualify my
generic response to a question.




> The point of the simulation is to address power. As our data goes down to
> N, of the 100 simulations, only 53 are good simulations because we don't
> have enough data for nls or chol2inv to work correctly.
>
>
> monte
>
> {x:
>
> ###########################################################################################
>
>
> ## case I ## EQUAL SAMPLE SIZES and design points
> nsim = 100;
> N_i = M_i = 10; ## also try (10, 30, 50, 100, 200)
> r = M_i / N_i;
>
> X.start = 170; # 6 design points, at 170,180,190, etc. where each point has
> N_i elements
> X.increment = 10;
> X.points = 6;
> X.end = 260;
> Xval = seq(X.start,length.out=X.points,by=X.increment );
> Xval = seq(X.start,X.end,length.out=X.points);
>
> L = 7; ## 6 + control
> k = 3;
> varY = 0.15;
>
> ### for each simulation, we need to record all of this information, write to
> a table or file.
>
> ### Under the null of simulation, we assign all locations to have same model
> ## we assume these are the true parameters
> b = 2.87; d = 0.0345; t = 173;
>
>
> B = seq(2.5,4.5,length.out=21);
> #B = seq(2.75,3.25,length.out=21);
> #B = seq(2.85,2.95,length.out=21);
> #B = seq(2.8,3.0,length.out=21);
> B = seq(2.5,3.2,length.out=21);
> D = seq(0.02,0.04,length.out=21);
> T = seq(165,185,length.out=21);
>
> alpha = .05;
> nu = k; ## number of parameters
> tr = L-1; ## number of treatments (including control)
> rho = 1/(1+r); ## dependency parameter
> myCritical = qmvchi(alpha,nu,tr,rho);
> ## we change one parameter at a time until the results fail most of the
> time.
>
>
> ## do independent for now, but let's store the parameters and quantiles???
> INFO for one location
> # beta delta tau nsim %Reject(V.pooled) %Reject(V.total) [Simulation level]
> resultS
> # beta delta tau i of nsim max(V.pooled) max(V.total) [Individual level]
> resultI
>
> resultS = resultI = NULL;
> for(p1 in 1:length(D))
> {
> print(paste(p1, " [D] of ",length(D))); flush.console();
> print(D[p1]);
> myReject.pooled = myReject.pooled.1 = MAX.pooled = rep(-1,nsim);
> gsim = 0; ## good simulations
> for(i in 1:nsim)
> {
> doubleBreak = F;
> print(paste(i, " of ",nsim)); flush.console();
> tData = NULL;
> pooledNum = matrix(0,nrow=k,ncol=k); ##numerator as weighted sum AS
> (n_k-1)cov.scaled
> pooledDen = 0; ##denominator as correction AS N-k
> #Sigma_pooled = ((omit.1-1)*summary.nls.1$cov.scaled +
> (omit.2-1)*summary.nls.2$cov.scaled +
> (omit.L-1)*summary.nls.L$cov.scaled)/(sum(omit.1,omit.2,omit.L)-L);
>
>
> for(j in 1:L)
> {
> Y = numeric(N_i);
> X = createDomain(Xval,N_i); noise = rnorm(N_i, mean=0,sd=sqrt(varY) );
>
> if(j==1)
> {
> ## location #1 is different
> Y = noise + evaluateModel(X,b,D[p1],t);
> beta = b;
> delta = D[p1];
> tau = t;
> } else {
> Y = noise + evaluateModel(X,b,d,t);
> }
> print(paste(j, " location NLS of ",L)); flush.console();
>
> fData = as.data.frame(cbind(Y,X)); colnames(fData)=c("Y","X"); unique =
> doUnique(fData);
> initialValues = list(b=3,d=0.04,t=180);
> #plot(X,Y,main=j);
>
> # http://stackoverflow.com/questions/2963729/r-catching-errors-in-nls
> try.fit = try(
> nls(
> myModel.nlm ,
> fData,
> start = initialValues,
> control = nls.control(warnOnly = TRUE),
> trace=T
> )
> );
>
> if(class(try.fit) == "try-error")
> {
> doubleBreak = T;
> print(doubleBreak);
> break; ## skip to next simulation?
> } else {
> fit.nls = try.fit;
> summary.nls = summary(fit.nls);
> summary.nls$cov.scaled = scaledCOV(summary.nls);
> pooledDen = pooledDen + dim(fData)[1];
> pooledNum = pooledNum + (dim(fData)[1]-1)*summary.nls$cov.scaled;
> results = list("data"=fData,"fit.nls"=fit.nls,"summary.nls"=summary.nls);
> tData = rbind(tData,fData); ## total data
> }
> if(j==L)
> {
> myStr = "nls.L";
> } else {
> myStr = paste("nls.",j,sep="");
> }
> assign(myStr,results);
> }
> if(doubleBreak==T)
> {
> # break from outer loop if fit didn't work [SKIP simulation]
> print(doubleBreak);
> doubleBreak = F;
> next;
> }
> gsim = gsim + 1;
> # http://www.maths.bris.ac.uk/~mazjcr/SGP.R
>
> COV.pooled = pooledNum/pooledDen;
>
> ## loop back through, use COV.t and COV.pooled to do tests and record reject
> or not
> CONTROL = nls.L$summary.nls$parameters[,1];
> Vp = numeric(L-1);
> for(j in 1:(L-1))
> {
> myStr = paste("nls.",j,sep="");
> myData = get(myStr);
>
> Diff = myData$summary.nls$parameters[,1]-nls.L$summary.nls$parameters[,1];
>
> H.pooled = COV.pooled + myData$summary.nls$cov.scaled;
>
> Vp[j] = t(Diff)%*%solve(H.pooled)%*%(Diff);
> myStr = paste("Vp.",j,sep="");
> assign(myStr,Vp[j]);
> }
> MAX.pooled[i] = max(Vp);
>
> myReject.pooled[i] = (MAX.pooled[i]> myCritical); ## should be NA (nls
> failed??), TRUE, FALSE
> myReject.pooled.1[i] = (Vp[1] == max(Vp)); ## did the reject come from the
> parameter change, or somewhere else
> ## http://biocodenv.com/wordpress/?p=15
>
> simI = c(beta,delta,tau,i,nsim,MAX.pooled[i]);
> resultI = rbind(resultI,simI);
> }
> hist(MAX.pooled,breaks=100,ylab=paste(length(myReject.pooled[myReject.pooled
> == TRUE]), " / ",gsim,sep=""), xlab=paste("d =",delta) );
> abline(v=myCritical,col="red");
> # beta delta tau nsim %Reject(V.pooled) %Reject(V.total) [Simulation level]
> resultS
> simP = c(beta,delta,tau,gsim, length(myReject.pooled[myReject.pooled ==
> TRUE]), length(myReject.pooled.1[myReject.pooled.1 == TRUE]) );
> resultS = rbind(resultS,simP);
> }
>
>
> try( colnames(resultS) =
> c("beta","delta","tau","gsim","reject.pooled","max.pooled.L1"));
>
> resultS;
>
> write.table(resultS,file=paste("sim.results/sim.case.I.N_",N_i,"_10000.delta.txt",sep=""),append=F,sep="|");
> write.table(resultI,file=paste("sim.results/sim.case.I.N_",N_i,"_10000.delta.details.txt",sep=""),append=F,sep="|");
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> 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.