#SOME TIPS BEFORE YOU START
#Try to use the text editor Notepad++ available at http://notepad-plus.sourceforge.net/uk/site.htm 
#This wil enable to efficiently edit your code in R (also other languages)

# To run this file in R, either type in R:
# source("http://home.isa.utl.pt/~joaopalma/courses/safmod/r/safmod.R")
# or
# source("C:\\safmod\\dados\\safmod.R")
# Obviously, the above path needs to be directed to this file! Please note the double backward slash 
# To run this script you will need the course data called TreeBiomass_Ec.txt which contains data for the moddeling exercises.
# You can access the data and this file at http://home.isa.utl.pt/~joaopalma/courses/safmod/ 

# PLEASE INSTALL THE FOLLOWING PACKAGES:
# Install the package "sm" available from
# Menu-> Packages->Install package(s)
# Select "sm" from the list then
# Menu->Packages->Load Package and 
# Select the "sm" you have just installed
#install.packages('sm')	# install 'sm' package


library('sm')			# load 'sm' package

#READ DATA FILE either online:
ec=read.delim(file='http://home.isa.utl.pt/~joaopalma/courses/safmod/data/TreeBiomass_Ec.txt',header=T)
# or download it and read it locally e.g.:
#ec=read.delim(file='C:\\SAFMOD\\dados\\TreeBiomass_Ec.txt',header=T)

# print(summary(ec)) 
opar <- par(mfrow = c(1,1), oma = c(0, 0, 1.1, 0))
plot(ec$ww~ec$d)  #ou plot(ec$d,ec$ww)


#SELECIONAR DETERMINADOS DADOS DO SET
print("Enter to SELECT DATA")
pause()
ec = subset(ec, d>0.5)
ec = subset(ec, h > 2)
ec = subset(ec, ww>0)
#print(ec$ww)
plot(ec$ww~ec$d)  #ou plot(ec$d,ec$ww)

print("Enter to LINEAR MODEL")
pause()
opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
mod1<-lm(ec$ww~ec$d)
plot(mod1)
print (summary(mod1))
pause()
opar <- par(mfrow = c(1,1))
plot(ec$ww~ec$d)  #ou plot(ec$d,ec$ww)
lines(ec$d, mod1$coefficients[[1]]+mod1$coefficients[[2]]*ec$d,type="l",col='red') 

# # #OLD CODE for model 2 
# # print("Enter to EXPONENTIAL MODEL")
# # pause()
# # opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
# # LnW<-log(ec$ww)
# # LnD<-log(ec$d)
# # mod2<-lm(LnW~LnD)
# # plot(mod2)
# # print (summary(mod2))
# # pause()
# # opar <-par(mfrow = c(1,1))
# # plot(ec$ww~ec$d)  #ou plot(ec$d,ec$ww)
# # lines(sort(ec$d), sort(exp(mod2$coefficients[[1]])*ec$d^mod2$coefficients[[2]]),type="l",col='red') 

#NEW CODE for model 2  - Allows to compare residuals between models
print("Enter to NONLINEAR EXPONENTIAL MODEL D")
pause()
opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
mod2<- nls(ww~ a * d^b, data=ec, start=list(a=1, b=1 ))
plot(predict(mod2), ec$ww)
abline(0,1, col=2)
plot(ec$ww~ec$d)  #ou plot(ec$d,ec$ww)
lines(sort(ec$d), sort(predict(mod2)), col=2)
print(summary(mod2))

print("Enter to NONLINEAR EXPONENTIAL MODEL H D")
pause()
opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
mod3<- nls(ww~ a * d^b * h^c, data=ec, start=list(a=1, b=1, c=1 ))
plot(predict(mod3), ec$ww)
abline(0,1, col=2)
plot(ec$ww~ec$d)  #ou plot(ec$d,ec$ww)
lines(sort(ec$d), sort(predict(mod3)), col=2)
plot(ec$ww~ec$h)  #ou plot(ec$d,ec$ww)
lines(sort(ec$h), sort(predict(mod3)), col=2)
print(summary(mod3))


#Compare model residuals
print("Exercise: Press Enter to Compare modelsī residuals")
pause()
opar <-par(mfrow = c(3,1))
plot(residuals(mod1))
plot(residuals(mod2))
plot(residuals(mod3))


#Compare models predicted vs estimated
print("Exercise: Press enter to Compare modelsī estimated vs predicted")
pause()
opar <-par(mfrow = c(3,1))
plot(predict(mod1), ec$ww )
abline(0,1, col=2)
plot(predict(mod2), ec$ww)
abline(0,1, col=2)
plot(predict(mod3), ec$ww)
abline(0,1, col=2)


#Compare models predicted vs estimated
print("Exercise: Press enter to Compare modelsī summary")
pause()
print(summary(mod1))
print(summary(mod2))
print(summary(mod3))

# Compare Rsquared
# sum(resid(mod3)^2) # http://markmail.org/message/av6se3jn2hognklx
MULTIPLE REGRESSION
MANUALS:
http://cran.r-project.org/doc/contrib/Faraway-PRA.pdf

READ DATA FILE either online:
print("Reading permanent sampling plots data...")
ec_w=read.delim(file='http://home.isa.utl.pt/~joaopalma/courses/safmod/data/Ec_Wha.txt',header=T)
print("Permanent sampling plots read..")
print("Enter to NONLINEAR EXPONENTIAL MODEL G Hdom")
pause()
opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
mod4<- nls(Ww ~ a * G^b * hdom^c, data=ec_w, start=list(a=1, b=1, c=1 ))
plot(predict(mod4), ec_w$ww)
abline(0,1, col=2)
plot(ec_w$Ww~ec_w$G) 
lines(sort(ec_w$G), sort(predict(mod4)), col=2)
plot(ec_w$Ww~ec_w$hdom) 
lines(sort(ec_w$hdom), sort(predict(mod4)), col=2)
print(summary(mod4))



print("Enter to observe Homocedasticity")
pause()
plot(predict(mod4),residuals(mod4)) 
print("Enter to reduce Homocedasticity by using G as weights - NOT YET DONE")

ec_w_sorted = ec_w.sort.list("ID_parc","t")