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

