#tried specifying covariances for chems to diminish the empirical underidentification
#of growth, but didn't work

#only blocks 1-8, no controls, no 11_5 genotype because
#of the problem with growth rate and rodent damage (i.e. negative growth rate values).
#many genos not included because not terpenoid data for them
#This model includes removal treatments, growth rate, final above-ground biomass, cap, rhizomes, rust, ducking, herbivores, protease inhibitors,
#and terpenoids


#growth rate calculated based on difference in sums of heights at s1,s2.
#rams not included in the model because perfectly correlated with growth

#need to run importance.R to get summed importance values for insect latent composite loadings
setwd("E:/Tissue_value/Path_Analysis/Onxy/Preliminary/")
source("E:/Tissue_value/Path_Analysis/Onxy/Preliminary/importance.R")
#need to run counts.R first to get galling, sucking, and chewing insect abundances to load on the "insects" composite
source("E:/Tissue_value/Path_Analysis/Onxy/Preliminary/counts.R")
#bring in terpenoid data via the terps.R file
#which recalculates totals with the two major outliers replaced with means
#terps.R produces three variables: total monoterpenoids, sesquiterpenoids, diterpenoids
#called mono, ses, di, respectively
source("E:/Tissue_value/Path_Analysis/Onxy/Preliminary/terps.R")

r.file <- "all_dat_v3_copy.txt"
dat <- read.delim(r.file,header=T)
str(dat)
names(dat)

#recalculate growth rate
s1h <- rowSums(dat[,11:16], na.rm = T)
s2h <- rowSums(dat[,43:48], na.rm = T)
date1 <- as.Date(dat$date_s1, format = "%d-%b-%y")
date2 <- as.Date(dat$date_s2, format = "%d-%b-%y")
time <- as.numeric(date2 - date1)
growth2 <- (s2h - s1h)/time
#hist(growth2)

#make dummy variables for stages
S1 <- c(0,0,1,0)[dat$stage]
S2 <- c(1,0,0,0)[dat$stage]

#make STAGE(AGE)
S1A <- c(0,0,0,0,1,0,0,0,0)[dat$trt]
S2A <- c(0,0,0,0,0,0,1,0,0)[dat$trt]
S3A <- c(0,0,0,0,0,0,0,0,1)[dat$trt]

#recalculate duck so if found to duck at any stage it is scored as ducking (i.e., 1)
aduck <- abs(dat$duck_s2-dat$duck_s1)
dat[which(is.na(aduck)),c(1:10,42,133)]
aduck[c(17,27,60,86,126,137,160,505,513)] <- c(0,1,rep(0,7))

#calculate rust severity
sev0 <- dat$rust_s1
sevT <- dat$rust_s2
date1 <- as.Date(dat$date_s1, format = "%d-%b-%y")
date2 <- as.Date(dat$date_s2, format = "%d-%b-%y")
time <- as.numeric(date2 - date1)
rustAUC <- (sev0 + sevT)/2*time
rustAUC <- log(rustAUC+40, base = 10) #use optimize script in Rtips.txt file to find +40
#hist(rustAUC)


treat <- c(NA,NA,NA,1,1,1,1,1,1)[dat$trt]#removes controls with na.omit below
dat2 <- data.frame(growth = growth2[which(dat$location == "dp")],
  dat[which(dat$location == "dp"),c(129,117,126)],
  geno_ID = dat$geno_ID[which(dat$location == "dp")],treat = treat[which(dat$location == "dp")],S1 = S1[which(dat$location == "dp")],
  S2 = S2[which(dat$location == "dp")],S1A = S1A[which(dat$location == "dp")],
  S2A = S2A[which(dat$location == "dp")],S3A = S3A[which(dat$location == "dp")],
  rust = rustAUC[which(dat$location == "dp")], duck = aduck[which(dat$location == "dp")],
  ser = dat$is120[which(dat$location == "dp")], cys = dat$cyspi_sp[which(dat$location == "dp")],
  mono = mono[which(dat$location == "dp")], ses = ses[which(dat$location == "dp")],
  di = di[which(dat$location == "dp")])
dat3 <- data.frame(dat2[-which(dat2$geno_ID == "11_5" | is.na(dat2$treat)),],counts)#add counts here
dat3 <- data.frame(na.omit(dat3[,-c(5,6)]))#removes NA's and treat and geno because just used to remove control plants and 11.5
dat3[,12:16] <- log(dat3[,12:16]+100) #log to normalize and stabalize variances
dat3[,17:19] <- log(dat3[,17:19]+1) #log to normalize and stabalize variances

cov_dat3 <- data.frame(cov(dat3))
write.table(cov_dat3, file = "cov_dat3.txt",append = FALSE, quote = FALSE, sep = "\t",row.names = T, col.names = T)
cov_dat3 <- read.delim("cov_dat3.txt",header=T)

##########Define the model###############
require("OpenMx");
manifests<-c("growth","rhiz","cap","stem","S1","S2","S1A","S2A","S3A","rust","duck","ser","cys","mono","ses","di","gall_c","suck_c","chew_c")
latents<-c("x1","x2","x3","x4","x5","x6","x7","x8","x9","x10","x11","x12","x13")
endos <- c("growth","rhiz","cap","stem","rust","ser","cys","mono","ses","di","gall_c","suck_c","chew_c")

model <- mxModel("prelim", 
mxData(observed = cov_dat3, type="cov", numObs = 282),
type="RAM",
manifestVars = manifests,
latentVars = latents,

#specify exogenous to endogenous paths
mxPath(from="S1", to = endos, free = T, arrows = 1,
  label = c("S1_TO_growth","S1_TO_rhiz","S1_TO_cap","S1_TO_stem","S1_TO_rust","S1_TO_ser","S1_TO_cys",
  "S1_TO_mono","S1_TO_ses","S1_TO_di","S1_TO_gall_c","S1_TO_suck_c","S1_TO_chew_c")),

mxPath(from="S2", to = endos, free = T, arrows = 1,  
  label = c("S2_TO_growth","S2_TO_rhiz","S2_TO_cap","S2_TO_stem","S2_TO_rust","S2_TO_ser","S2_TO_cys",
  "S2_TO_mono","S2_TO_ses","S2_TO_di","S2_TO_gall_c","S2_TO_suck_c","S2_TO_chew_c")),

mxPath(from="S1A", to = endos, free = T, arrows = 1,
  label = c("S1A_TO_growth","S1A_TO_rhiz","S1A_TO_cap","S1A_TO_stem","S1A_TO_rust","S1A_TO_ser","S1A_TO_cys",
  "S1A_TO_mono","S1A_TO_ses","S1A_TO_di","S1A_TO_gall_c","S1A_TO_suck_c","S1A_TO_chew_c")),

mxPath(from="S2A", to = endos, free = T, arrows = 1,
  label = c("S2A_TO_growth","S2A_TO_rhiz","S2A_TO_cap","S2A_TO_stem","S2A_TO_rust","S2A_TO_ser","S2A_TO_cys",
  "S2A_TO_mono","S2A_TO_ses","S2A_TO_di","S2A_TO_gall_c","S2A_TO_suck_c","S2A_TO_chew_c")),

mxPath(from="S3A", to = endos, free = T, arrows = 1,
  label = c("S3A_TO_growth","S3A_TO_rhiz","S3A_TO_cap","S3A_TO_stem","S3A_TO_rust","S3A_TO_ser","S3A_TO_cys",
  "S3A_TO_mono","S3A_TO_ses","S3A_TO_di","S3A_TO_gall_c","S3A_TO_suck_c","S3A_TO_chew_c")),

mxPath(from="duck", to = endos, free = T, arrows = 1,
  label = c("duck_TO_growth","duck_TO_rhiz","duck_TO_cap","duck_TO_stem","duck_TO_rust","duck_TO_ser","duck_TO_cys",
  "duck_TO_mono","duck_TO_ses","duck_TO_di","duck_TO_gall_c","duck_TO_suck_c","duck_TO_chew_c")),


#specify endogenous to endogenous
mxPath(from="growth", to = c("stem","cap","rhiz"), free = T, arrows = 1, 
  label = c("growth_TO_stem","growth_TO_cap","growth_TO_rhiz")),

mxPath(from="growth", to = c("ser","cys","mono","ses","di"),free = T, arrows = 1,
  label = c("growth_TO_ser","growth_TO_cys","growth_TO_mono","growth_TO_ses","growth_TO_di")),

mxPath(from=c("ser","cys","mono","ses","di"), to = "growth", free = T, arrows = 1,
  label = c("ser_TO_growth","cys_TO_growth","mono_TO_growth","ses_TO_growth","di_TO_growth")),

mxPath(from=c("ser","cys","mono","ses","di"), to = "gall_c", free = T, arrows = 1,
  label = c("ser_TO_gall","cys_TO_gall","mono_TO_gall","ses_TO_gall","di_TO_gall")),
mxPath(from=c("ser","cys","mono","ses","di"), to = "suck_c", free = T, arrows = 1,
  label = c("ser_TO_suck","cys_TO_suck","mono_TO_suck","ses_TO_suck","di_TO_suck")),
mxPath(from=c("ser","cys","mono","ses","di"), to = "chew_c", free = T, arrows = 1,
  label = c("ser_TO_chew","cys_TO_chew","mono_TO_chew","ses_TO_chew","di_TO_chew")),
mxPath(from=c("ser","cys","mono","ses","di"), to = "rust", free = T, arrows = 1,
  label = c("ser_TO_rust","cys_TO_rust","mono_TO_rust","ses_TO_rust","di_TO_rust")),
mxPath(from=c("ser","cys","mono","ses","di"), to = "stem", free = T, arrows = 1,
  label = c("ser_TO_stem","cys_TO_stem","mono_TO_stem","ses_TO_stem","di_TO_stem")),
mxPath(from=c("ser","cys","mono","ses","di"), to = "cap", free = T, arrows = 1,
  label = c("ser_TO_cap","cys_TO_cap","mono_TO_cap","ses_TO_cap","di_TO_cap")),
mxPath(from=c("ser","cys","mono","ses","di"), to = "rhiz", free = T, arrows = 1,
  label = c("ser_TO_rhiz","cys_TO_rhiz","mono_TO_rhiz","ses_TO_rhiz","di_TO_rhiz")),


mxPath(from="stem", to = c("cap","rhiz"), free = T, arrows = 1, label = c("stem_TO_cap","stem_TO_rhiz")),

mxPath(from="rust", to = c("growth","stem","cap","rhiz"), free = T, arrows = 1, 
  label = c("rust_TO_growth","rust_TO_stem","rust_TO_cap","rust_TO_rhiz")),

mxPath(from="stem", to = c("gall_c","suck_c","chew_c"), free = T, arrows = 1, 
  label = c("stem_TO_gall","stem_TO_suck","stem_TO_chew")),

mxPath(from=c("gall_c","suck_c","chew_c"), to = c("rust"), free = T, arrows = 1, 
  label = c("gall_TO_rust","suck_TO_rust","chew_TO_rust")),

mxPath(from=c("gall_c","suck_c","chew_c"), to = c("growth"), free = T, arrows = 1, 
  label = c("gall_TO_growth","suck_TO_growth","chew_TO_growth")),

mxPath(from=c("gall_c","suck_c","chew_c"), to = c("cap"), free = T, arrows = 1, 
  label = c("gall_TO_cap","suck_TO_cap","chew_TO_cap")),

mxPath(from=c("gall_c","suck_c","chew_c"), to = c("rhiz"), free = T, arrows = 1, 
  label = c("gall_TO_rhiz","suck_TO_rhiz","chew_TO_rhiz")),



#specify disturbances
mxPath(from="x1", to = "gall_c", free = F, value = 1, arrows = 1),
mxPath(from="x2", to = "suck_c", free = F, value = 1, arrows = 1),
mxPath(from="x3", to = "chew_c", free = F, value = 1, arrows = 1),
mxPath(from="x4", to = "ser", free = F, value = 1, arrows = 1),
mxPath(from="x5", to = "cys", free = F, value = 1, arrows = 1),
mxPath(from="x6", to = "mono", free = F, value = 1, arrows = 1),
mxPath(from="x7", to = "ses", free = F, value = 1, arrows = 1),
mxPath(from="x8", to = "di", free = F, value = 1, arrows = 1),
mxPath(from="x9", to = "rust", free = F, value = 1, arrows = 1),
mxPath(from="x10", to = "growth", free = F, value = 1, arrows = 1),
mxPath(from="x11", to = "stem", free = F, value = 1, arrows = 1),
mxPath(from="x12", to = "cap", free = F, value = 1, arrows = 1),
mxPath(from="x13", to = "rhiz", free = F, value = 1, arrows = 1),

#start values taken from model terps_v5.R and terps_v5g.R; always important to have good start values here!!!!!
#mxPath(from="x1",  free = T, value = 0.428, arrows = 2, label = "VAR_x1",  lbound = 0, ubound = var(dat3$gall_c)),
#mxPath(from="x2",  free = T,                arrows = 2, label = "VAR_x2",  lbound = 0, ubound = var(dat3$suck_c)),
#mxPath(from="x3",  free = T,                arrows = 2, label = "VAR_x3",  lbound = 0, ubound = var(dat3$chew_c)),
#mxPath(from="x4",  free = T, value = 0.039, arrows = 2, label = "VAR_x4",  lbound = 0, ubound = var(dat3$ser)),
#mxPath(from="x5",  free = T, value = 0.022, arrows = 2, label = "VAR_x5",  lbound = 0, ubound = var(dat3$cys)),
#mxPath(from="x6",  free = T, value = 0.097, arrows = 2, label = "VAR_x6",  lbound = 0, ubound = var(dat3$mono)),
#mxPath(from="x7",  free = T, value = 0.096, arrows = 2, label = "VAR_x7",  lbound = 0, ubound = var(dat3$ses)),
#mxPath(from="x8",  free = T,                arrows = 2, label = "VAR_x8",  lbound = 0, ubound = var(dat3$di)),
#mxPath(from="x9",  free = T, value = 0.082, arrows = 2, label = "VAR_x9",  lbound = 0, ubound = var(dat3$rust)),
#mxPath(from="x10", free = T, value = 2.33,  arrows = 2, label = "VAR_x10", lbound = 0, ubound = var(dat3$growth)),
#mxPath(from="x11", free = T, value = 15.98, arrows = 2, label = "VAR_x11", lbound = 0, ubound = var(dat3$stem)),
#mxPath(from="x12", free = T, value = 3.39,  arrows = 2, label = "VAR_x12", lbound = 0, ubound = var(dat3$cap)),
#mxPath(from="x13", free = T, value = 11.96, arrows = 2, label = "VAR_x13", lbound = 0, ubound = var(dat3$rhiz)),

#start values taken from model terps_v5.R and terps_v5g.R; always important to have good start values here!!!!!
mxPath(from="x1",  free = T, value = 0.428, arrows = 2, label = "VAR_x1"),
mxPath(from="x2",  free = T,                arrows = 2, label = "VAR_x2"),
mxPath(from="x3",  free = T,                arrows = 2, label = "VAR_x3"),
mxPath(from="x4",  free = T, value = 0.039, arrows = 2, label = "VAR_x4"),
mxPath(from="x5",  free = T, value = 0.022, arrows = 2, label = "VAR_x5"),
mxPath(from="x6",  free = T, value = 0.097, arrows = 2, label = "VAR_x6"),
mxPath(from="x7",  free = T, value = 0.096, arrows = 2, label = "VAR_x7"),
mxPath(from="x8",  free = T,                arrows = 2, label = "VAR_x8"),
mxPath(from="x9",  free = T, value = 0.082, arrows = 2, label = "VAR_x9"),
mxPath(from="x10", free = T, value = 2.33,  arrows = 2, label = "VAR_x10"),
mxPath(from="x11", free = T, value = 15.98, arrows = 2, label = "VAR_x11"),
mxPath(from="x12", free = T, value = 3.39,  arrows = 2, label = "VAR_x12"),
mxPath(from="x13", free = T, value = 11.96, arrows = 2, label = "VAR_x13"),


#specify covariances
mxPath(from="x12", to="x13",free = T, value = -0.902, arrows = 2, label = "COV_x12_x13"),

#specify covariances that may change to paths
mxPath(from="x1", to="x2",free = T, value = 0.089, arrows = 2, label = "COV_x1_x2"),

#specify covariances for exogenous
mxPath(from="S1", to="S2", free = T, value = cov(dat3$S1,dat3$S2), arrows = 2, label = "COV_S1_S2"),
mxPath(from="S1A", to="S2A", free = T, value = cov(dat3$S1A,dat3$S2A), arrows = 2, label = "COV_S1A_S2A"),
mxPath(from="S2A", to="S3A", free = T, value = cov(dat3$S2A,dat3$S3A), arrows = 2, label = "COV_S2A_S3A"),
mxPath(from="S1A", to="S3A", free = T, value = cov(dat3$S1A,dat3$S3A), arrows = 2, label = "COV_S1A_S3A"),
mxPath(from="S1", to=c("S1A","S2A","S3A"), free = T, value = c(cov(dat3$S1,dat3$S1A),cov(dat3$S1,dat3$S2A),cov(dat3$S1,dat3$S3A)), arrows = 2, 
  label = c("COV_S1_S1A","COV_S1_S2A","COV_S1_S3A")),
mxPath(from="S2", to=c("S1A","S2A","S3A"), free = T, value = c(cov(dat3$S2,dat3$S1A),cov(dat3$S2,dat3$S2A),cov(dat3$S2,dat3$S3A)), arrows = 2, 
  label = c("COV_S2_S1A","COV_S2_S2A","COV_S2_S3A")),
mxPath(from="duck", to=c("S1","S2","S1A","S2A","S3A"), free = T, 
  value = c(cov(dat3$duck,dat3$S1),cov(dat3$duck,dat3$S2),cov(dat3$duck,dat3$S1A),cov(dat3$duck,dat3$S2A),cov(dat3$duck,dat3$S3A)), arrows = 2, 
  label = c("COV_duck_S1","COV_duck_S2","COV_duck_S1A","COV_duck_S2A","COV_duck_S3A")),

#specify variances for exogenous variables
#endogenous variable variances are estimated in their respective disturbance terms
mxPath(from="S1",free = T, value = var(dat3$S1), arrows = 2, label = "VAR_S1"),
mxPath(from="S2",free = T, value = var(dat3$S2), arrows = 2, label = "VAR_S2"),
mxPath(from="S1A",free = T, value = var(dat3$S1A), arrows = 2, label = "VAR_S1A"),
mxPath(from="S2A",free = T, value = var(dat3$S2A), arrows = 2, label = "VAR_S2A"),
mxPath(from="S3A",free = T, value = var(dat3$S3A), arrows = 2, label = "VAR_S3A"),
mxPath(from="duck",free = T, value = var(dat3$duck), arrows = 2, label = "VAR_duck")
)

################ End Model Specification ################

results1 <- mxRun(model)
results2 <- mxRun(model)
rerun <- mxRun(results1)
summary(results1)
summary(results2)
#compare unbound estimates minus bound estimates
data.frame(diff = round(summary(results1)$parameters$Estimate-summary(results2)$parameters$Estimate,digits = 3),
  names = summary(results1)$parameters$name, 
  prop = round(summary(results2)$parameters$Estimate/summary(results1)$parameters$Estimate,digits = 3))
summary(rerun)
check <- which(abs(summary(results1)$parameters[,7]) > 1)#checks for Heywood cases 
summary(results1)$parameters[check,]
check <- which(abs(summary(rerun)$parameters[,7]) > 1)#checks for Heywood cases 
summary(rerun)$parameters[check,]

#get the correlation residuals
#fix the correlation between treat and loc to zero
residuals <- cov2cor(cov(dat3)) - cov2cor(results1$objective@info$expCov)
residuals[which(residuals < 0.01)] <- NA
round(residuals, digits=4)

#now I calculate standardized coefficients
#first load the function written by Ryne Estabrook (20 Oct 2010)
#his function name is standardizeRAM, which is in the file: StandardizeRAM.R
source("E:/Tissue_value/Path_Analysis/StandardizedRAM.R")
#stand1 <- standardizeRAM(results1)
#stand2 <- standardizeRAM(rerun)
stand3 <- standardizeRAM(rerun2)
#these all come out the same as multiple regression, he also gives standard errors for these but these are not
#typically reported in path analysis and I think there's some controversy about how to calculate them anyway.

#now I need to figure out p-values
#first get the vectors of coefficients and SEs
pvals <- function(results){
  res <- data.frame(summary(results)$parameters)[,c(1:6)]
  names <- res[which(res$matrix == "A"),1]
  betas <- res[which(res$matrix == "A"),5]
  SEs <- res[which(res$matrix == "A"),6]
  ts <- betas/SEs
  n <- summary(results)$numObs
  p <- function(x,n) {2*(pt(-abs(x),df = n-2))}
  ps <- p(ts,n)
  out <- data.frame(names = names, estimate = round(betas,digits = 4), SEs = round(SEs, digits = 4),
  t = round(ts,digits = 2), P_value = round(ps, digits = 4))
  print(out)
}

pvalues1 <- pvals(results1)
small_ps <- pvalues1[which(pvalues1$P_value <= 0.01),]
length(pvalues1[which(pvalues1$P_value <= 0.01),]$t)#only going to show these paths in the model

#calculate pvalues for covariances
tx1_x2 <- summary(results1)$parameters$Estimate[170]/summary(results1)$parameters$Std.Error[170]
tx12_x13 <- summary(results1)$parameters$Estimate[182]/summary(results1)$parameters$Std.Error[182]
p <- function(x,n) {2*(pt(-abs(x),df = n-2))}
covs <- p(c(tx1_x2,tx12_x13),282)
round(covs, digits = 5)

#extract betas with p values <= 0.01
summary(results1)$parameters[as.numeric(row.names(small_ps)),]

summary(results1)
#calculate R2s
1-(summary(results1)$parameters$Estimate[169]/var(dat3$gall_c))
1-(summary(results1)$parameters$Estimate[171]/var(dat3$suck_c))
1-(summary(results1)$parameters$Estimate[172]/var(dat3$chew_c))
1-(summary(results1)$parameters$Estimate[173]/var(dat3$ser))
1-(summary(results1)$parameters$Estimate[174]/var(dat3$cys))
1-(summary(results1)$parameters$Estimate[175]/var(dat3$mono))
1-(summary(results1)$parameters$Estimate[176]/var(dat3$ses))
1-(summary(results1)$parameters$Estimate[177]/var(dat3$di))
1-(summary(results1)$parameters$Estimate[178]/var(dat3$rust))
1-(summary(results1)$parameters$Estimate[179]/var(dat3$growth))
1-(summary(results1)$parameters$Estimate[180]/var(dat3$stem))
1-(summary(results1)$parameters$Estimate[181]/var(dat3$cap))
1-(summary(results1)$parameters$Estimate[183]/var(dat3$rhiz))





omxGraphviz(model, dotFilename="graph.dot")

summary(lm(rhiz ~ gall_d, data = dat3))