Hi Ben,
thanks for your advise. I put the comment in front of the line with nls. Now
the message doesn't appear but the script won't produce the result files
which I need (the result files are empty).
When executing the script, R asks for number of iterations and after that it
states "Read 20 items" but nothing happens. Probably also some other part of
the script needs to be amended.
Here is the script:
#library (nls)
options (digits = 6, show.error.messages=FALSE)
cat ("Enter the number of iterations: ")
niter<-as.integer (readLines(n = 1))
cat ("Parameter a (asymptote) for Kohn's equation", file="Param_a_Kohn.txt",
append=TRUE, fill=TRUE)
cat ("Number of iterations: ",niter, file="Param_a_Kohn.txt", append=TRUE,
fill=TRUE)
cat (" Niter Min Mean Max Sd",
file="Param_a_Kohn.txt",append=TRUE, fill=TRUE)
cat ("Parameter b for Kohn's equation", file="Param_b_Kohn.txt",
append=TRUE, fill=TRUE)
cat ("Number of iterations: ",niter, file="Param_b_Kohn.txt", append=TRUE,
fill=TRUE)
cat (" Niter Min Mean Max Sd",
file="Param_b_Kohn.txt",append=TRUE, fill=TRUE)
cat ("Parameter a (asymptote) for Eggert's equation",
file="Param_a_Eggert.txt", append=TRUE, fill=TRUE)
cat ("Number of iterations: ",niter, file="Param_a_Eggert.txt", append=TRUE,
fill=TRUE)
cat (" Niter Min Mean Max Sd",
file="Param_a_Eggert.txt",append=TRUE, fill=TRUE)
cat ("Parameter b for Eggert's equation", file="Param_b_Eggert.txt",
append=TRUE, fill=TRUE)
cat ("Number of iterations: ",niter, file="Param_b_Eggert.txt", append=TRUE,
fill=TRUE)
cat (" Niter Min Mean Max Sd",
file="Param_b_Eggert.txt",append=TRUE, fill=TRUE)
#construction of the dataset to estimate the parameter from output file from
GEMINI
createdata<-function(n, niter=500) {
concent<<-function(echa) apply(matrix(1
:length(echa),ncol=1),1,flocal1,pop=echa)
flocal1<-function(pop,j) length(unique(pop[1:j]))
indsel<-length (n)
if (indsel<=1) stop ("n>1 expected")
faesel<-sum(n)
echa<-rep(1:indsel,n)
return(echa,faesel)
}
#function to estimate and print the results
calculate<-function(data){
dataset<-createdata(data, niter)
coeffaKohn<-rep(0, niter)
coeffbKohn<-rep(0, niter)
niterKohn<-0
coeffaEggert<-rep(0, niter)
coeffbEggert<-rep(0, niter)
niterEggert<-0
draw<-array(0,c(niter+1, dataset$faesel))
draw1<-array(0,c(niter+1, dataset$faesel))
draw2<-array(0,c(niter+1, dataset$faesel))
draw3<-array(0,c(niter+1, dataset$faesel))
for (k in 1:niter) {
w<-rep(0,dataset$faesel)
w<-w+concent(sample(dataset$echa,dataset$faesel,replace=F))
faesel<-length (w)
dxy<-cbind.data.frame(1:faesel,w)
names(dxy)<-c("x","y")
draw[k,1:faesel]<-w[1:faesel]
ok<-0
while(ok<5){
tmp1<-NA
tmp2<-NA
tmp<-try(coef(nls1<-nls(y~a*x/(b+x), data=dxy,
start=list(a=max(dxy$y),b=1))))
tmp1<-tmp [1]
tmp2<-tmp [2]
if (is.character(tmp1)) {coeffaKohn[k]<-NA
ok<-ok+1}
else {coeffaKohn[k]<-tmp1
niterKohn<-niterKohn + 1
draw1[k,1:faesel]<-predict(nls1)
ok<-5}
if (is.character(tmp2)) {coeffbKohn[k]<-NA}
else {coeffbKohn[k]<-tmp2}
}
ok<-0
while(ok<5){
tmp1<-NA
tmp2<-NA
tmp<-try(coef(nls3<-nls(y~a*(1-exp(b*x)), data=dxy,
start=list(a=max(dxy$y),b=-1))))
tmp1<-tmp [1]
tmp2<-tmp [2]
if (is.character(tmp1)) {coeffaEggert[k]<-NA
ok<-ok+1}
else {coeffaEggert[k]<-tmp1
niterEggert<-niterEggert + 1
draw3[k,1:faesel]<-predict(nls3)
ok<-5}
if (is.character(tmp2)) {coeffbEggert[k]<-NA}
else {coeffbEggert[k]<-tmp2}
}
}
x11()
plot(draw[1,], main="Number of unique genotypes against number of feces
analysed", sub="observed= circles; mean of observed= black line; Kohn's eq=
red line; Eggert's eq= blue line", xlab="Number of feces analysed",
ylab="Number of unique genotypes")
for (k in 2:niter) {
lines(draw[k,], type="o")}
for (k in 1:faesel) {
draw[niter+1,k]<-mean(draw[1:niter,k])
draw1[niter+1,k]<-mean(draw1[1:niter,k])
draw2[niter+1,k]<-mean(draw2[1:niter,k])
draw3[niter+1,k]<-mean(draw3[1:niter,k])
}
lines(draw[niter+1,], type="l", col="black", lwd=2)
lines(draw1[niter+1,], type="l", col="red", lwd=2)
lines(draw3[niter+1,], type="l", col="blue", lwd=2)
minaKohn<-min(coeffaKohn, na.rm=TRUE)
maxaKohn<-max(coeffaKohn, na.rm=TRUE)
meanaKohn<-mean(coeffaKohn, na.rm=TRUE)
sdaKohn<-sd(coeffaKohn, na.rm=TRUE)
minbKohn<-min(coeffbKohn, na.rm=TRUE)
maxbKohn<-max(coeffbKohn, na.rm=TRUE)
meanbKohn<-mean(coeffbKohn, na.rm=TRUE)
sdbKohn<-sd(coeffaKohn, na.rm=TRUE)
minaEggert<-min(coeffaEggert, na.rm=TRUE)
maxaEggert<-max(coeffaEggert, na.rm=TRUE)
meanaEggert<-mean(coeffaEggert, na.rm=TRUE)
sdaEggert<-sd(coeffaEggert, na.rm=TRUE)
minbEggert<-min(coeffbEggert, na.rm=TRUE)
maxbEggert<-max(coeffbEggert, na.rm=TRUE)
meanbEggert<-mean(coeffbEggert, na.rm=TRUE)
sdbEggert<-sd(coeffbEggert, na.rm=TRUE)
pop<-""
cat(pop,niterKohn, minaKohn, meanaKohn, maxaKohn, sdaKohn, sep=" ", file =
"Param_a_Kohn.txt",append=TRUE, fill = TRUE)
cat(pop,niterKohn, minbKohn, meanbKohn, maxbKohn, sdbKohn, sep=" ", file =
"Param_b_Kohn.txt",append=TRUE, fill = TRUE)
cat(pop,niterEggert, minaEggert, meanaEggert, maxaEggert, sdaEggert, sep="
", file = "Param_a_Eggert.txt",append=TRUE, fill = TRUE)
cat(pop,niterEggert, minbEggert, meanbEggert, maxbEggert, sdbEggert, sep="
", file = "Param_b_Eggert.txt",append=TRUE, fill = TRUE)
x11()
hist1<<-hist(rawdata, right=FALSE, 1:(max(rawdata)+1), plot=FALSE)
hist(rawdata, 1:(max(rawdata)+1), ylim=c(0,(max(hist1$counts)+1)),
right=FALSE, labels = as.character(hist1$mids-0.5), plot=TRUE,
main="Histogram of capture frequencies per sample", xlab="Number of capture
per sample", proba=F)
lines(density(rawdata, 1))
}
i<-1
calculate(rawdata<-scan("Rarefaction_Polana_20081.txt"))
--
View this message in context:
http://r.789695.n4.nabble.com/Error-in-library-nls-tp3630143p3630455.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
R-help@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.
