# REPLACE "C:\\...\\etc" WITH YOUR PATH AND FILE NAME

# source(file = "C:\\...\\crownstemdata.txt", echo = TRUE)

allometry <- read.table(file = "C:\\...\\crownstemdata.txt", 
header = TRUE, sep = "")

library(nlme)

lndiam <- log(allometry$diam_cm)
lnarea <- log(allometry$crown_area_m)
trunk <- allometry$trunk
sal <- allometry$salt
species3 <- allometry$species2
indiv <- allometry$ind
Young <- allometry$Young
motor <- allometry$paraglider

table(species3)

leaf <-  matrix(c(NA), ncol = 1, nrow = length(lnarea))
leaf <- ifelse(species3 == 1, 6.82, leaf)
leaf <- ifelse(species3 == 2, 10.6, leaf)
leaf <- ifelse(species3 == 3, 13, leaf)
leaf <- ifelse(species3 == 4, 9.04, leaf)
leaf <- ifelse(species3 == 5, 15.93, leaf)
leaf <- ifelse(species3 == 6, 11.02, leaf)
leaf <- ifelse(species3 == 7, 16.7, leaf)
leaf <- ifelse(species3 == 8, 4, leaf)
leaf <- ifelse(species3 == 9, 13.36, leaf)
leaf <- ifelse(species3 == 10, 12.44, leaf)
leaf <- ifelse(species3 == 11, 6.342, leaf)

# DUMMY VARIABLES FOR EACH GROUP OF SPECIES.

groupA <- matrix(c(NA), ncol = 1, nrow = length(species3))
groupB <- matrix(c(NA), ncol = 1, nrow = length(species3))
groupC <- matrix(c(NA), ncol = 1, nrow = length(species3))

groupA	<- ifelse(species3 == 2 | species3 == 9 | species3 == 10 | species3 == 11, 1, 0)
groupB	<- ifelse(species3 == 1 | species3 == 7 | species3 == 8, 1, 0)
groupC	<- ifelse(species3 == 3 | species3 == 4 | species3 == 5 | species3 == 6, 1, 0)

datos <- data.frame(cbind(lndiam, lnarea, trunk, sal, species3, indiv, Young, motor, leaf, groupA, groupB, groupC))

# FIT AND PRINT SEPARATE MULTIPLE REGRESSION MODELS FOR EACH SPECIES.

 models <- list()
 for(k in 1:max(species3))
 {
 models[[as.character(k)]] <- summary(lm(lnarea ~ lndiam + trunk + sal, data = subset(datos, species3 == k)))
 }

models

# FIT AND PRINT A SINGLE MUTIPLE REGRESSION MODEL FOR THE WHOLE DATA SET USING DUMMY VARIABLES FOR SPECIES.

model_sp <- lm(lnarea ~ lndiam + trunk + sal + factor(species3), data = datos)
summary(model_sp)

# FIT AND PRINT A SINGLE MULTIPLE REGRESSION MODEL WITH SPECIES CATEGORIZED INTO THREE GROUPS.

model1_group <- lm(lnarea ~ lndiam + trunk + sal + groupA + groupC, data = datos)
summary(model1_group)

# ADD TO THE PREVIOUS MODEL YOUNGS MODULUS AND LEAF SIZE.

model2_group <- update(model1_group, ~ . + Young + leaf)
summary(model2_group)

# ADD TO THE PREVIOUS MODEL ALL FIRST ORDER INTERACTION TERMS AND DELETE YOUNGS MODULUS AND LEAF SIZE.

model3_group <- update(model2_group, ~ . - Young - leaf + lndiam*trunk + lndiam*sal + lndiam*groupA + lndiam*groupC + 
		trunk*sal + trunk*groupA + trunk*groupC + sal*groupA + sal*groupC)
summary(model3_group)

# ADD TO MODEL model1_group THE SIGNIFICANT INTERACTION TERMS.

model4_group <- update(model1_group, ~ . + lndiam*groupA + trunk*groupA)
summary(model4_group)

# REFIT THE PREVIOUS MODEL WITH A RANDOM EFFECT FOR TREES.

model5_group <- lme(lnarea ~ lndiam + trunk + sal + groupA + groupC + lndiam*groupA + trunk*groupA, random = ~ 1|indiv, data = datos, method = "ML")
summary(model5_group)

# REFIT THE PREVIOUS MODEL WITH TWO ERROR TERMS, ONE FOR TRUNKS AND THE OTHER ONE FOR BRANCHES.

model6_group <- update(model5_group, weights = varIdent(form = ~ 1|trunk))
summary(model6_group)
intervals(model6_group)

# RESIDUAL ANALYSIS FOR THE FINAL MODEL: model6_group

 residual <- resid(model6_group, level = 0:1, type = "p")
 summary(residual)
 plot(model6_group, resid(.,type = "p") ~ fitted(.) | groupB, id = 0.05)
 plot(model6_group, resid(.,type = "p") ~ fitted(.), id = 0.05)
 plot(model6_group, lnarea ~ fitted(.), id = 0.05)
 anova(model5_group, model6_group)
 qqnorm(model6_group, ~ resid(.) | groupB, id = 0.05)
 qqnorm(model6_group, ~ ranef(.), id = 0.1)

# HYPOTHESIS TESTS FOR COMPARING THE FIT OF NESTED MODELS.

anova(model_sp, model1_group)
anova(model1_group, model2_group)
anova(model1_group, model3_group)
anova(model1_group, model4_group)
anova(model3_group, model4_group)
anova(model5_group, model6_group)

# WHEN THE SYMBOL # BEFORE THE sink COMMAND IS DELETED, THE RESULTS ARE PRINTED TO A FILE.

# sink(file = "C:\\...\\UniversalScalingResult.txt")

models
summary(model_sp)
summary(model1_group)
summary(model2_group)
summary(model3_group)
summary(model4_group)
summary(model5_group)
summary(model6_group)
intervals(model6_group)

anova(model_sp, model1_group)
anova(model1_group, model2_group)
anova(model1_group, model3_group)
anova(model1_group, model4_group)
anova(model3_group, model4_group)
anova(model5_group, model6_group)

# sink()

 loge_diam <- matrix(c(seq(-0.7, 5.1, 0.2)), ncol = 1)

# ESTIMATED MODELS FITTED SEPARATELY FOR EACH SPECIES.

sp1_m1  <- function(x, stemtype, saltexp)
	{-2.48259 + 1.72991*x + 0.3356*stemtype + 0.09510*saltexp}
sp2_m1  <- function(x, stemtype, saltexp)
	{-3.61209 + 1.93218*x + 0.82858*stemtype - 0.50385*saltexp}
sp3_m1  <- function(x, stemtype, saltexp)
	{-1.52075 + 1.60407*x - 0.01418*stemtype - 0.01986*saltexp}
sp4_m1  <- function(x, stemtype, saltexp)
	{-1.86195 + 1.60877*x + 0.23417*stemtype + 0.00622*saltexp}
sp5_m1  <- function(x, stemtype, saltexp)
	{-1.80046 + 1.5531*x + 0.58027*stemtype + 0.05854*saltexp}
sp6_m1  <- function(x, stemtype, saltexp)
	{-1.77252 + 1.7039*x + 0.13184*stemtype - 0.3099*saltexp}
sp7_m1  <- function(x, stemtype, saltexp)
	{-1.26936 + 1.41327*x + 0.41662*stemtype - 0.51846*saltexp}
sp8_m1  <- function(x, stemtype, saltexp)
	{-2.1646 + 1.5088*x + 0.8262*stemtype + 0.2175*saltexp}
sp9_m1  <- function(x, stemtype, saltexp)
	{-2.08014 + 1.62434*x + 0.01544*stemtype - 0.03018*saltexp}
sp10_m1 <- function(x, stemtype, saltexp)
	{-2.4929 + 1.5862*x + 0.9029*stemtype - 0.5108*saltexp}
sp11_m1 <- function(x, stemtype, saltexp)
	{-3.12328 + 1.76226*x + 0.96196*stemtype - 0.50496*saltexp}

# FINAL ESTIMATED MODEL FOR EACH GROUP OF SPECIES.

modelA <- function(x, stemtype, saltexp)
	{(-2.0882070 - 0.9779676) + (1.6343024 + 0.1433438)*loge_diam + (0.323319 + 0.3730926)*stemtype - 0.1614561*saltexp}
modelB <- function(x, stemtype, saltexp)
	{-2.0882070 + 1.6343024*loge_diam + 0.323319*stemtype - 0.1614561*saltexp}
modelC <- function(x, stemtype, saltexp)
	{(-2.0882070 + 0.2723734) + 1.6343024*loge_diam + 0.323319*stemtype - 0.1614561*saltexp}

 par(mfrow = c(4,3))

# ESTIMATES OF log crown area FOR TRUNKS NOT EXPOSED 
# TO SEA SALT BASED ON THE TWO MODELS MENTIONED ABOVE.

sp1  <- cbind(modelB(loge_diam,1,0), sp1_m1(loge_diam,1,0))
sp2  <- cbind(modelA(loge_diam,1,0), sp2_m1(loge_diam,1,0))
sp3  <- cbind(modelC(loge_diam,1,0), sp3_m1(loge_diam,1,0))
sp4  <- cbind(modelC(loge_diam,1,0), sp4_m1(loge_diam,1,0))
sp5  <- cbind(modelC(loge_diam,1,0), sp5_m1(loge_diam,1,0))
sp6  <- cbind(modelC(loge_diam,1,0), sp6_m1(loge_diam,1,0))
sp7  <- cbind(modelB(loge_diam,1,0), sp7_m1(loge_diam,1,0))
sp8  <- cbind(modelB(loge_diam,1,0), sp8_m1(loge_diam,1,0))
sp9  <- cbind(modelA(loge_diam,1,0), sp9_m1(loge_diam,1,0))
sp10 <- cbind(modelA(loge_diam,1,0), sp10_m1(loge_diam,1,0))
sp11 <- cbind(modelA(loge_diam,1,0), sp11_m1(loge_diam,1,0))

# SELECT OBSERVED DATA FOR TRUNKS NOT EXPOSED TO SEA SALT.

dat1  <- subset(datos, trunk==1 & sal==0 & species3 == 1, select = 1:2)
dat2  <- subset(datos, trunk==1 & sal==0 & species3 == 2, select = 1:2)
dat3  <- subset(datos, trunk==1 & sal==0 & species3 == 3, select = 1:2)
dat4  <- subset(datos, trunk==1 & sal==0 & species3 == 4, select = 1:2)
dat5  <- subset(datos, trunk==1 & sal==0 & species3 == 5, select = 1:2)
dat6  <- subset(datos, trunk==1 & sal==0 & species3 == 6, select = 1:2)
dat7  <- subset(datos, trunk==1 & sal==0 & species3 == 7, select = 1:2)
dat8  <- subset(datos, trunk==1 & sal==0 & species3 == 8, select = 1:2)
dat9  <- subset(datos, trunk==1 & sal==0 & species3 == 9, select = 1:2)
dat10 <- subset(datos, trunk==1 & sal==0 & species3 == 10, select = 1:2)
dat11 <- subset(datos, trunk==1 & sal==0 & species3 == 11, select = 1:2)

# GRAPHS OF THE OBSERVED AND ESTIMATED CROWN AREAS FOR TRUNKS NOT EXPOSED TO SEA SALT.

matplot(loge_diam, sp1, xlab = "Ln diameter (cm)", ylab = "Ln crown area (m2)", 
main = "Cupania glabra", type = "l", col = 1) 
points(dat1$lndiam, dat1$lnarea, pch = 20)

matplot(loge_diam, sp2, xlab = "Ln diameter (cm)", ylab = "Ln crown area (m2)", 
main = "Coccoloba hondurensis", type = "l", col = 1) 
points(dat2$lndiam, dat2$lnarea, pch = 20)

matplot(loge_diam, sp3, xlab = "Ln diameter (cm)", ylab = "Ln crown area (m2)", 
main = "Cordia alliodora", type = "l", col = 1) 
points(dat3$lndiam, dat3$lnarea, pch = 20)

matplot(loge_diam, sp4, xlab = "Ln diameter (cm)", ylab = "Ln crown area (m2)", 
main = "Cordia stenoclada", type = "l", col = 1) 
points(dat4$lndiam, dat4$lnarea, pch = 20)

matplot(loge_diam, sp5, xlab = "Ln diameter (cm)", ylab = "Ln crown area (m2)", 
main = "Ficus obtusifolia", type = "l", col = 1) 
points(dat5$lndiam, dat5$lnarea, pch = 20)

matplot(loge_diam, sp6, xlab = "Ln diameter (cm)", ylab = "Ln crown area (m2)", 
main = "Ficus tecolutensis", type = "l", col = 1) 
points(dat6$lndiam, dat6$lnarea, pch = 20)

matplot(loge_diam, sp7, xlab = "Ln diameter (cm)", ylab = "Ln crown area (m2)", 
main = "Ficus trigonata", type = "l", col = 1) 
points(dat7$lndiam, dat7$lnarea, pch = 20)

matplot(loge_diam, sp8, xlab = "Ln diameter (cm)", ylab = "Ln crown area (m2)", 
main = "Lonchocarpus schiedeanus", type = "l", col = 1) 
points(dat8$lndiam, dat8$lnarea, pch = 20)

matplot(loge_diam, sp9, xlab = "Ln diameter (cm)", ylab = "Ln crown area (m2)", 
main = "Pouteria squamosa", type = "l", col = 1) 
points(dat9$lndiam, dat9$lnarea, pch = 20)

matplot(loge_diam, sp10, xlab = "Ln diameter (cm)", ylab = "Ln crown area (m2)", 
main = "Spondias radlkoferi", type = "l", col = 1) 
points(dat10$lndiam, dat10$lnarea, pch = 20)

matplot(loge_diam, sp11, xlab = "Ln diameter (cm)", ylab = "Ln crown area (m2)", 
main = "Trichilia havanensis", type = "l", col = 1) 
points(dat11$lndiam, dat11$lnarea, pch = 20)


# savePlot(file = "C:\\...\\UniversalScalingGraph", type = c("pdf"))