0

I obtained the chronosequences from empirical data. To simultaneously calculate the site and global parameters of Equation (1), I used a nested iterative procedure. Equation (1) :

H2=H1 (((1-exp⁡(b2∙T2 ))/(1-exp⁡(b2∙T1 ) )))^((b3∙(H1)^(b1 ))) where H2 is the measured height at age T2, H1 is a site parameter denoting a stand height at age T1, and β1, β2, and β3 are estimated parameters.

The nested iterative procedure starts with the calibration of the global parameters of the equation, using the preliminary values of (H1) determined by empirical material. In the next iteration, the preliminary values of the global parameters are used as constants, and the site parameters are estimated for every tree. Next, the global parameters are refitted, using the H1 estimates for a given trajectory as the constants. The nested iterative procedure was repeated until the parameters of the model stabilised. Parameter estimation was carried out in R using an nls (nonlinear least squares) procedure.

######data is composed of heights and ages of trees that I obtained the chronosequences  from empirical data.
fit1<-nls (H2~H1*((1-exp(b2*T2))/(1-exp(b2*100)))^(b3*(H1^b1)), data=data,    #### nonlinear equation(1)
             start=list(b1=-1,b2=-0.04,b3=53),
             control=nls.control(maxiter = 1000,tol=1e-7)
             )   
fit1  ####Parameters

b1= -0.8709, b2= -0.02852, b3=28.61296 ######fitted parameters

predH <- predict(fit1, level = c(0:2))   ########the fitted values from the first iteration
modDat <- data.frame(data, predH)       ########add the fitted value to the data
modDat ##### new data frame

para<-c(-0.8709,-0.02852,28.61296)    
para
iter<-0
repeat
{
  b1<-para[1]
  b2<-para[2]
  b3<-para[3]

  Db<-cbind(modDat$H1*((1-exp(b2*modDat$T2))/(1-exp(b2*100)))^(b3*(modDat$H1^b1)),modDat$predH*((1-exp(b2*modDat$T2))/(1-exp(b2*100)))^(b3*(modDat$predH^b1)))
  new.fit<-para+ Solve((Db)%*%as.numeric(sapply(Db, as.numeric))%*%(Db))%*% (predH*((1-exp(b2$T2))/(1-exp(b2*100)))^(b3*(predH^b1)))
  iter<-iter+1
  if (sum((para-new.fit)^2)<1e-5||iter>1000) break;
  para<-new.fit
}

Solve from library(systemfit)

but I got the message : "requires numeric/complex matrix/vector arguments" in spite of I put it as.numeric instead of as.matrix!

tobias
  • 1
  • 1
  • You haven't provided enough info to answer your question. What does `data` look like? What package is `Solve()` from? One thing I notice is that `cbind()` doesn't take a data argument, so that could be your problem. Try using `modDat$H1`, etc. – lhs Sep 05 '22 at 16:35
  • Thank you for your comment, I treid to give an idea about the nature of the data, I also added the data argument to `cbind()` , Solve (not solve) is from package (systemfit). Nothing changed – tobias Sep 06 '22 at 07:24
  • Does `modDat` have columns "t" and "H1"? Your old code had "T2", not "t". In the `new.fit <-` line, you have `b2$T2`, but `b2` doesn't have any element "T2". Finally, the systemfit package (at least the one I just downloaded from CRAN) does not seem to have a function called `Solve()` or anything similar. I would suggest making a [minimal reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) in order to get help, because no one can reproduce your error. – lhs Sep 06 '22 at 21:00
  • Greetings! Usually it is helpful to provide a minimally reproducible dataset for questions here. One way of doing this is by using the `dput` function. You can find out how to use it here: https://youtu.be/3EID3P1oisg – Shawn Hemelstrand Sep 10 '22 at 06:18

0 Answers0