R help - Nov 2010 - 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...
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]]

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