Howdy,
Last week I got some great help on why I was getting an error code when trying
to run this model, thanks everyone! I was able to get the code up and running
beautifully for several data sets. Now I am getting different errors with this
new data set. I can't figure out why, I have more data points with this
species, and it is ordered exactly the same as the other species I have been
working with.
The two errors I get are for these specific lines of code:
This error I actually got with all the data sets but its just assumption
checking so I did this manually.
> fitGen <- nls(vbGen,data=BLG,start=svGen)
Error in nlsModel(formula, mf, start, wts) :
singular gradient matrix at initial parameter estimates
and
This code is only not working with this data set.
> fitK <- nls(vbK,data=BLG,start=svK)
Error in numericDeriv(form[[3L]], names(ind), env) :
Missing value or an infinity produced when evaluating the model
Thanks,
April
Here is all my code and my data set:
Data set:
structure(list(Length = c(60.96, 60.96, 63.5, 63.5, 66.04, 66.04,
71.12, 71.12, 96.52, 96.52, 96.52, 101.6, 101.6, 106.68, 111.76,
116.84, 116.84, 119.38, 124.46, 124.46, 134.62, 137.16, 144.78,
147.32, 149.86, 154.94, 157.48, 167.64, 172.72, 175.26, 175.26,
180.34, 180.34, 182.88, 190.5, 193.04, 193.04, 193.04, 195.58,
198.12, 203.2, 210.82), Age = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L,
4L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 3L, 5L, 5L, 6L, 5L, 7L,
6L, 6L), Cohort = c(2007L, 2009L, 2004L, 2005L, 2008L, 2011L,
2006L, 2010L, 2006L, 2007L, 2008L, 2004L, 2009L, 2003L, 2010L,
2005L, 2007L, 2008L, 2005L, 2006L, 2003L, 2009L, 2006L, 2008L,
2007L, 2004L, 2002L, 2005L, 2004L, 2003L, 2007L, 2005L, 2006L,
2006L, 2002L, 2002L, 2004L, 2005L, 2003L, 2005L, 2004L, 2002L
), Year = c(2008L, 2010L, 2005L, 2006L, 2009L, 2012L, 2007L,
2011L, 2008L, 2009L, 2010L, 2006L, 2011L, 2005L, 2012L, 2007L,
2010L, 2011L, 2008L, 2009L, 2006L, 2012L, 2010L, 2012L, 2011L,
2007L, 2006L, 2009L, 2008L, 2007L, 2012L, 2010L, 2011L, 2012L,
2005L, 2007L, 2009L, 2011L, 2008L, 2012L, 2010L, 2008L)), .Names = c("Length",
"Age", "Cohort", "Year"), class = "data.frame", row.names = c(NA,
-42L))
Code:
BLG=read.csv("BLGa.csv", header=T)
attach(BLG)
names(BLG)
dput(BLG,"")
$ Von Bertalanffy Growth Model$
source("http://www.rforge.net/NCStats/InstallNCStats.R";)
source("http://www.rforge.net/FSA/InstallFSA.R";)
library(FSA)
library(NCStats)
library(nlstools)
# provide starting values
svCom <- list(Linf=199.8132, K=0.5631745, t0=-1.960147)
svCom <- vbStarts(Length~Age,data=BLG)
# repeat the starting values for each Cohort
svGen <- lapply(svCom, rep, length(unique(Cohort)))
# define a new variable, uniquely numbered for each Cohort
Cohortf <- as.factor(Cohort)
# fit a single model with different parameter estimates for each Cohort
vbGen <- Length~Linf[Cohortf]*(1-exp(-K[Cohortf]*(Age-t0[Cohortf])))
fitGen <- nls(vbGen,data=BLG,start=svGen)
#fit a single model with the same parameter estimates using all the Cohorts
data#
vbCom <- Length~Linf*(1-exp(-K*(Age-t0)))
fitCom <- nls(vbCom,data=BLG,start=svCom)
overview(fitCom)
# Checking assumptions#
#want residual plot to look like a shot gun blast, no funnels, no trends, no
arches, (ie Unbiased and homoscedastic)#
residPlot(fitCom)
#Want histogram to be normally distributed#
hist(residuals(fitCom),main="")
# Create a model that holds Linf and t0 steady but varies by Cohort#
vbK=Length~Linf*(1-exp(-K[Cohortf]*(Age-t0)))
svK=mapply(rep,svCom,c(1,8,1))
#Fit model to the data#
fitK <- nls(vbK,data=BLG,start=svK)
#Compare the confidence intervals for each Cohort's estimates of K, remember
Linf, and t0 were found for entire dataset#
overview(fitK)
overview(fitCom)
#fitted plots of Length at Age with K varying by Cohort#
plot(Length~Age,data=BLG,subset=Cohortf=="2005",pch=7,xlab="Age
(yrs)",ylab="Mean Length at Capture (Length)", xlim=c(1,5), ylim=c(0,250))
points(Length~Age,data=BLG,subset=Cohortf=="2006",pch=0,col="black")
points(Length~Age,data=BLG,subset=Cohortf=="2007",pch=1,col="black")
points(Length~Age,data=BLG,subset=Cohortf=="2008",pch=2,col="black")
points(Length~Age,data=BLG,subset=Cohortf=="2009",pch=3,col="black")
points(Length~Age,data=BLG,subset=Cohortf=="20010",pch=4,col="black")
points(Length~Age,data=BLG,subset=Cohortf=="20011",pch=5,col="black")
points(Length~Age,data=BLG,subset=Cohortf=="20012",pch=6,col="black")
#Add growth curves to data set#
vbTypical <- vbFuns("typical")
curve(vbTypical(x,Linf=coef(fitK)[-c(3,4,5,6,7,8,9)]),from=1,to=5,lwd=2,lty=2,add=TRUE,xlim=c(1,5),
ylim=c(0,250))
curve(vbTypical(x,Linf=coef(fitK)[-c(2,4,5,6,7,8,9)]),from=1,to=5,lwd=2,lty=3,add=TRUE)
curve(vbTypical(x,Linf=coef(fitK)[-c(2,3,5,6,7,8,9)]),from=1,to=5,lwd=2,lty=4,add=TRUE)
curve(vbTypical(x,Linf=coef(fitK)[-c(2,3,4,6,7,8,9)]),from=1,to=5,lwd=2,lty=5,add=TRUE)
curve(vbTypical(x,Linf=coef(fitK)[-c(2,3,4,5,7,8,9)]),from=1,to=5,lwd=2,lty=6,add=TRUE)
curve(vbTypical(x,Linf=coef(fitK)[-c(2,3,4,5,6,8,9)]),from=1,to=5,lwd=2,lty=99,add=TRUE)
curve(vbTypical(x,Linf=coef(fitK)[-c(2,3,4,5,6,7,9)]),from=1,to=5,lwd=2,lty=4194,add=TRUE)
curve(vbTypical(x,Linf=coef(fitK)[-c(2,3,4,5,6,7,9)]),from=1,to=5,lwd=2,lty=5995,add=TRUE)
#and a legend for the points#
legend("right",legend=c("2005","2006",
"2007","2008","2009","2010","2011","2012"),pch=c(7,0,1,2,3,4,5,6),lwd=1,lty=1,cex=0.75)
#legend for the curves#
legend("bottomright",legend=c("2005","2006",
"2007","2008","2009","2010","2011","2012",
"Overall"),lty=c(2,3,4,5,6,99,4194,5995, 1),lwd=2,cex=0.75)
#all points on 1 graph, with 1 curve#
vbCom <- Length~Linf*(1-exp(-K*(Age-t0)))
fitCom <- nls(vbCom,data=BLG,start=svCom)
plot(Length~Age,data=BLG,pch=21,xlab="Age (yrs)",ylab="Mean Length at Capture
(mm)", xlim=c(1,5), ylim=c(0,250))
curve(vbTypical(x,Linf=coef(fitCom)),from=1,to=5,lwd=4,add=TRUE)
[[alternative HTML version deleted]]
______________________________________________
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.
