setwd("C:/Users/user/R Scripts/Linear Models/Data Files")
source("C:/Users/user/R Scripts/Linear Models/Functions.R")


	#---------------------read data--------------------------------
	trees.dat <- read.table("Trees.dat", col.names = c("Diameter", "Height",
		"Volume"))
	n <- nrow(trees.dat)
	#---------------------Linear Model------------------------------
	#
	par(pty = "s", oma = c(0., 0., 4., 0.))
	attach(trees.dat)
	trees.fit <- lm(Volume ~ Diameter + Height)
	par(mfrow = c(2., 1.))
	plot(trees.dat[, "Diameter"], trees.dat[, "Volume"], xlab = "diameter",
		ylab = "volume")
	#---------------------analysis of residuals---------------------
	plot(trees.dat[, "Height"], trees.dat[, "Volume"], xlab = "height",
		ylab = "volume")
	par(mfrow = c(2., 2.))
	print(summary(trees.fit, correlation = F))
	trees.cv <- CrossVal(trees.fit)
	cat("\nCV=", round(trees.cv$cv, 4), "GCV=", round(trees.cv$gcv, 4),
		"\nCross-validation R-squared=", round(trees.cv$r2cv, 4), 
		"\n\n")
	plot(trees.dat[, "Diameter"], rstandard(trees.fit), xlab = "Diameter", ylab = 
		"Std. Residuals")
	plot(trees.dat[, "Height"], rstandard(trees.fit), xlab = "Height", ylab = 
		"Std. Residuals")
		plot(trees.fit,which=c(1,2))
#	mtext("Linear Model", 3., out = T, cex = 2.)
	#-----------------Box-Cox--------------------------------------
	par(mfrow = c(1., 1.))
	#
	boxcox(trees.fit, seq(-0.3, 0.6, 0.1))
	#
	#------------------Volume^(1/3) Model--------------------------
	#
	mtext("Box-Cox Power Transformation", 3., out = T, cex = 2.)
	trees.fit <- lm(Volume^(1./3.) ~ Diameter + Height, y = T)
	#	plot(trees.dat[, "Diameter"], trees.dat[, "Volume"]^(1/3), xlab = 
	#		"diameter", ylab = "volume^(1/3)")
	#	plot(trees.dat[, "Height"], trees.dat[, "Volume"]^(1/3), xlab = 
	#		"height", ylab = "volume^(1/3)")	
	#---------------------analysis of residuals--------------------
	#	par(mfrow = c(2, 1))
	print(summary(trees.fit, correlation = F))
	par(mfrow = c(2., 2.))
	trees.cv <- CrossVal(trees.fit)
	cat("\nCV=", round(trees.cv$cv, 4), "GCV=", round(trees.cv$gcv, 4),
		"\nCross-validation R-squared=", round(trees.cv$r2cv, 4), 
		"\n\n")
	plot(trees.dat[, "Diameter"], rstandard(trees.fit), xlab = "Diameter", ylab = 
		"Std. Residuals")
	plot(trees.dat[, "Height"], rstandard(trees.fit), xlab = "Height", ylab = 
		"Std. Residuals")
		plot(trees.fit,which=c(1,2))
#	mtext("Cube Root Model", 3., out = T, cex = 2.)
#	trees.new<-predict(trees.fit,data.frame(Height=80,Diameter=15),
#	interval=c("prediction"))
#	print(trees.new)
#	trees.new<-predict(trees.fit,data.frame(Height=80,Diameter=15),
#	interval=c("confidence"))
#	print(trees.new)
	#
	#
	#------------------Log-Log Model--------------------------
	#
	par(mfrow = c(1, 1))
	boxcox(lm(Volume ~ log(Diameter) + log(Height)), seq(-1, 1, 0.1))
	par(mfrow = c(2., 1.))
	plot(log(trees.dat[, "Diameter"]), log(trees.dat[, "Volume"]), xlab = 
		"log(diameter)", ylab = "log(volume)")
	plot(log(trees.dat[, "Height"]), log(trees.dat[, "Volume"]), xlab = 
		"log(height)", ylab = "log(volume)")
	trees.fit <- lm(log(Volume) ~ log(Diameter) + log(Height))
	print(summary(trees.fit, correlation = F))
	trees.cv <- CrossVal(trees.fit)
	cat("\nCV=", round(trees.cv$cv, 4), "GCV=", round(trees.cv$gcv, 4),
		"\nCross-validation R-squared=", round(trees.cv$r2cv, 4), 
		"\n\n")
	#---------------------analysis of residuals--------------------
	par(mfrow = c(2., 2.))
#	mtext("Log Model", 3., out = T, cex = 2.)
	plot(log(trees.dat[, "Diameter"]), rstandard(trees.fit), xlab = "log(Diameter)",
		ylab = "Std. Residuals")
	plot(log(trees.dat[, "Height"]), rstandard(trees.fit), xlab = "log(Height)", ylab = 
		"Std. Residuals")
		plot(trees.fit,which=1:2)
#	trees.new<-predict(trees.fit,data.frame(Height=80,Diameter=15),
#	interval=c("prediction"))
#	print(trees.new)
#	trees.new<-predict(trees.fit,data.frame(Height=80,Diameter=15),
#	interval=c("confidence"))
#	print(trees.new)
#	title("Q-Q Plot of Std. Residuals")
	#------------------Final Model--------------------------
	#
#	mtext("Log-Log Model", 3., out = T, cex = 2.)
	treesf.fit <- lm(log(Volume) ~ 1+offset(2. * log(Diameter/12.) + log(Height)))
	
	#---------------------analysis of residuals--------------------
	print(summary(treesf.fit, correlation = F))
	print(anova(treesf.fit,trees.fit))
	par(mfrow = c(2., 2.))
	r2 <- cor(log(trees.dat[, "Volume"]), predict(treesf.fit))^2.
	cat("Multiple R-squared:", round(r2, 4.), "\n")
	trees.cv <- CrossVal(treesf.fit)
	cat("\nCV=", round(trees.cv$cv, 4), "GCV=", round(trees.cv$gcv, 4),
		"\nCross-validation R-squared=", round(trees.cv$r2cv, 4), 
		"\n\n")
	plot(log(trees.dat[, "Diameter"]), rstandard(treesf.fit), xlab = "log(Diameter)",
		ylab = "Std. Residuals")
	plot(log(trees.dat[, "Height"]), rstandard(treesf.fit), xlab = "log(Height)", ylab = 
		"Std. Residuals")
		plot(treesf.fit,which=1:2)
#	mtext("'Final' Model", 3., out = T, cex = 2.)