Korf<-read.table("http://user.mendelu.cz/drapela/Forest_Biometry/Data/dh.txt",sep='\t',quote='',header=TRUE)
library(nlstools)
library(nls2) 
Korf


names(Korf)[1] <- "t"
names(Korf)[2] <- "h"
Korf
attach(Korf)
plot(h~t)

#                  DEFINITION OF THE FUNCTION
####################################################################
model=h~a*exp(k/((1-n)*t^(n-1)))



#             ESTIMATION OF INITIAL VALUES OF PARAMETERS
######################################################################

#by command "preview"

preview(formula = h~a*exp(k/((1-n)*t^(n-1))), data=Korf, start=list(a = 35, k = 30, n=2))



#by grid

grid.Korf<-expand.grid(list(a=seq(30,50,by=2),k=seq(20,50,by=2),n=seq(0.5,3,by=0.5)))
estKorf<-nls2(model, data=Korf,start=grid.Korf,algorithm="brute-force")
estKorf

plot(t,h,xlim=c(0,140),ylim=c(0,30), xlab="age",ylab="height [m]")
curve(38*exp(34/((1-2)*x^(2-1))),add=TRUE,col="red")



#                     NON-LINEAR REGRESSION 
#######################################################################

# with starting values based on grid - searching

mm<-nls(model,data=Korf,start=list(a=38,k=34,n=2))
summary(mm) 
overview(mm)

plot(t,h,xlim=c(0,140),ylim=c(0,30), xlab="age",ylab="height [m]")
#curve(38*exp(-34/x),add=TRUE,col="red")
#curve(coef(mm)[1]*exp(coef(mm)[2]/x),add=TRUE,col="green")
curve(coef(mm)[1]*exp(coef(mm)[2]/((1-coef(mm)[3])*x^(coef(mm)[3]-1))),add=T,col="green")
curve(38*exp(34/((1-2)*x^(2-1))),add=TRUE,col="red")

#                  BASIC RESULTS OF NONLINEAR REGRESSION
##########################################################################

#index of correlation
RSS<-sum((h-predict(mm))^2)
dev<-sum((h-mean(h))^2)
I=sqrt(1-(RSS/dev))
I2=I*I
I
I2

res<-nlsResiduals(mm)
plot(res)

library(car)
qqPlot(residuals(mm)) 

# test of normality of residuals
test.nlsResiduals(res)
AIC(mm)

detach(Korf)

####MICHAJLOV - COMPARISON WITH KORF #########
Mich<-read.table("http://user.mendelu.cz/drapela/Forest_Biometry/Data/dh.txt",sep='\t',quote='',header=TRUE)
Mich
model1<-h~a*exp(k/Mich$t)
names(Mich)[1] <- "t"
names(Mich)[2] <- "h"
Mich
attach(Mich)
mm1<-nls(model1,data=Mich,start=list(a=46,k=-36))
summary(mm1)
plot(t,h,xlim=c(0,140),ylim=c(0,30), xlab="age",ylab="height [m]")
curve(coef(mm1)[1]*exp(coef(mm1)[2]/x),add=TRUE,col="red")
curve(coef(mm)[1]*exp(coef(mm)[2]/((1-coef(mm)[3])*x^(coef(mm)[3]-1))),add=T,col="green")

AIC(mm1)
AIC(mm)
anova(mm,mm1)
detach(Mich)