0

I'm trying to estimate a function with 3 parameters, whereof 2 (growth parameters b, v) are shared between all sites, and 1 (asymptote A) is specific to each site.

nlme reports:

Error in MEEM(object, conLin, control$niterEM) : Singularity in backsolve at level 0, block 1.

Code below creates some data from a function on the same form, and then tries to estimate the parameters with nlme. Or maybe this is a job for nls?

#Function on other form to create example data.
#Note - only p and k are constants! Asymptote A varies for each plot.
ChRfunction <- function(A,p,k,age){
  return(A*(1-exp(-k*age))^p)
}

#Getting some random values 
# A == asympote
# A1 = Age at time point 1
# A2 = Age at time point 2
# H1, H2  = Heights at time point 1, 2 respectively.
set.seed(1)
Avals <-runif(50,min=8,max=24)
A1vals <- runif(50,min=15,max=160)
A2vals <- A1vals + runif(50,min=3,max=5)
TestData <- data.frame(A1 = A1vals,
                       A2 = A2vals,
                       H1 = ChRfunction(Avals,p=1.8,k=0.025,age=A1vals)+runif(n=50,min=0,max=0.4),
                       H2 = ChRfunction(Avals,p=1.8,k=0.025,age=A2vals)+runif(n=50,min=0,max=0.6))

#Each measurement is associated with a plot with a specific asymptote.
TestData$PlotID <- seq(1:nrow(TestData))


#These are starting values from a paper with a similar species to mine.
vTest=0.012
bTest = 0.0415
pTest = 0.7
tauTest= 50
maxValTest=8

#Simple plot with ggplot to understand data.
#Starting values are reasonable.
TestData %>% ggplot()+geom_point(aes(x=A1,y=H1),shape=1)+geom_point(aes(x=A2,y=H2),shape=1)+geom_segment(aes(x=A1,y=H1,xend=A2,yend=H2))+
  geom_function(fun=function(x) (maxValTest)*((1+(pTest^-vTest - 1)*exp(bTest*tauTest)*exp(-bTest)^x)/(1+(pTest^-vTest - 1)*exp(bTest*tauTest)*exp(-bTest)^100))^-(1/vTest))+
  geom_function(fun=function(x) (maxValTest+2)*((1+(pTest^-vTest - 1)*exp(bTest*tauTest)*exp(-bTest)^x)/(1+(pTest^-vTest - 1)*exp(bTest*tauTest)*exp(-bTest)^100))^-(1/vTest))+
  geom_function(fun=function(x) (maxValTest+4)*((1+(pTest^-vTest - 1)*exp(bTest*tauTest)*exp(-bTest)^x)/(1+(pTest^-vTest - 1)*exp(bTest*tauTest)*exp(-bTest)^100))^-(1/vTest))+
  geom_function(fun=function(x) (maxValTest+6)*((1+(pTest^-vTest - 1)*exp(bTest*tauTest)*exp(-bTest)^x)/(1+(pTest^-vTest - 1)*exp(bTest*tauTest)*exp(-bTest)^100))^-(1/vTest))+
  geom_function(fun=function(x) (maxValTest+8)*((1+(pTest^-vTest - 1)*exp(bTest*tauTest)*exp(-bTest)^x)/(1+(pTest^-vTest - 1)*exp(bTest*tauTest)*exp(-bTest)^100))^-(1/vTest))+
  geom_function(fun=function(x) (maxValTest+10)*((1+(pTest^-vTest - 1)*exp(bTest*tauTest)*exp(-bTest)^x)/(1+(pTest^-vTest - 1)*exp(bTest*tauTest)*exp(-bTest)^100))^-(1/vTest))



#Attempt to estimate parameters b, v. A is different for each PlotID.
# This function is related to the first, on another form.
nlme(H2~A*(1+(((A/H1)^v)-1)*exp(-b*(A2-A1)))^(-(1/v)),
     data=TestData,
     fixed = b + v ~ 1,
     random= A ~ 1|PlotID,
     start=c(b=0.0415,v=0.11)
)

0 Answers0