##### MODELING INTERACTIONS
### Data from Cherry, J.A., McKee, K.L., and Grace, J.B. 2009. Journal of Ecology 97:67-77.

setwd("D:/TalksAndTrips/FY2017/Germany/Workshop/Part5/ModelingInteractionsWithComposites")

dat <- read.csv("Modeling_Interactions_data.csv")
names(dat)
attach(dat)

### Transformations
ln.C3prod <- log(C3prod +1)

### Plots
dev.new(height=4, width=8); par(mfrow=c(1,2)); hist(C3prod); hist(ln.C3prod)
CO2trt <- factor(CO2)
salt.trt <- factor(Salinity)
par(mfrow=c(1,1)); boxplot(ln.C3prod ~ salt.trt + CO2trt)

### Standard analyses
summary(lm(C3prod ~ Salinity * CO2))
summary(lm(ln.C3prod ~ Salinity * CO2))


### lavaan modeling
library(lavaan)
dat2 <- data.frame(ln.C3prod, CO2, Salinity, CxS)

## Model 1
# specify model
mod.1 <- 'ln.C3prod ~ CO2 + Salinity + CxS' 

# fit model
mod.1.fit <- sem(mod.1, data=dat2, meanstructure = TRUE)

# request output
summary(mod.1.fit, rsq=T)


## Model 2
# specify model (using est from model 1)
mod.2 <- 'Comp <~ CO2 + 1*Salinity + CxS
          ln.C3prod ~ Comp' 

# fit model
mod.2.fit <- sem(mod.2, data=dat2)

# request output
summary(mod.2.fit, standardized=T, rsq=T, meanstructure = TRUE)


## Model 3 - create composite by hand
# get parameters from uncomposited model
summary(mod.1.fit)
# compute composite scores
comp.hand <- 6.762 -0.084*CO2 -1.529*Salinity +0.612*CxS
# add variable to data set
dat2$Comp.hand <- comp.hand

# specify model
mod.3 <- 'ln.C3prod ~ Comp.hand' 

# fit model
mod.3.fit <- sem(mod.3, data=dat2)

# request output
summary(mod.3.fit, standardized=T, rsq=T)