##### R Code for: Waller et al. (2018) 
### "Trait responses to AM fungi are stronger and more consistent than 
### fixed differences among populations of Asclepias speciosa." 
### Published in: American Journal of Botany.

m3 <- read.csv("Waller_AJB_Data.csv", header=T)  ## load full dataset
head(m3)
# load packages
library(lattice)
library(lme4)
library(lmerTest)
library(MuMIn)
library(reshape2)
library(lsmeans)

#######################################################################
######  RESHAPE FOR MULTIVARIATE TRAIT ANALYSIS #######################
### reshape centered trait values to long form ####################

melt3 <- melt(m3, measure.vars=c("TotBMs","Hgtdays","RSs","slas",
                                 "latexs","trics"), 
              id.vars=c("Pop","Fam","AMF","Seed_g"), 
              variable.name='Trait')

### examine histograms and boxplots
histogram(~value|Trait, data=melt3)
bwplot(value~AMF|Trait, data=melt3)


#########################################################
### full univariate models to answer Q1 

um1 <- lmer(value~AMF+scale(Seed_g)+(1|Pop)+(1|Pop:AMF), data=melt3, subset=Trait=="TotBMs")
um1a <- lmer(value~AMF+scale(Seed_g)+(1|Pop)+(1|Pop:Fam)+(1|Pop:AMF), 
             data=melt3, subset=Trait=="TotBMs") #### test for Family effect
anova(um1,um1a)
summary(um1)
summary(um1a)
rand(um1a)

anova(um1)
summary(um1)
rand(um1)
r.squaredGLMM(um1)

um2 <- lmer(value~AMF+scale(Seed_g)+(1|Pop)+(1|Pop:AMF), data=melt3, subset=Trait=="Hgtdays")
anova(um2)
summary(um2)
rand(um2)
r.squaredGLMM(um2)

um3 <- lmer(value~AMF+scale(Seed_g)+(1|Pop)+(1|Pop:AMF), data=melt3, subset=Trait=="RSs")
anova(um3)
summary(um3)
rand(um3)
r.squaredGLMM(um3)


um4 <- lmer(value~AMF+scale(Seed_g)+(1|Pop)+(1|Pop:AMF), data=melt3, subset=Trait=="slas")
anova(um4)
summary(um4)
rand(um4)
r.squaredGLMM(um4)


um5 <- lmer(value~AMF+scale(Seed_g)+(1|Pop)+(1|Pop:AMF), data=melt3, subset=Trait=="latexs")
anova(um5)
summary(um5)
rand(um5)
r.squaredGLMM(um5)


um6 <- lmer(value~AMF+scale(Seed_g)+(1|Pop)+(1|Pop:AMF), data=melt3, subset=Trait=="trics")
anova(um6)
summary(um6)
rand(um6)
r.squaredGLMM(um6)

####################################################################
### calculate AMF responses to answer Q2
### use raw trait values
####################################################################
wide1 <- dcast(m3, Pop+Fam~AMF,
               value.var="TotBM", mean)
names(wide1)[names(wide1)=="Yes"] <-"YT"
names(wide1)[names(wide1)=="No"] <-"NT"
wide1$amfresp <- log(wide1$YT/wide1$NT)
wide1$amfresp2 <- ((wide1$YT-wide1$NT)/wide1$YT)*100

wide2 <- dcast(m3, Pop+Fam~AMF,
               value.var="RS", mean)
names(wide2)[names(wide2)=="Yes"] <-"YR"
names(wide2)[names(wide2)=="No"] <-"NR"
wide2$RSresp <- log(wide2$YR/wide2$NR)
wide2$RSresp2 <- ((wide2$YR-wide2$NR)/wide2$YR)*100


wide3 <- dcast(m3, Pop+Fam~AMF,
               value.var="sla", mean)
names(wide3)[names(wide3)=="Yes"] <-"YS"
names(wide3)[names(wide3)=="No"] <-"NS"
wide3$slaresp <- log(wide3$YS/wide3$NS)
wide3$slaresp2 <- ((wide3$YS-wide3$NS)/wide3$YS)*100


wide4 <- dcast(m3, Pop+Fam~AMF,
               value.var="latex", mean)
names(wide4)[names(wide4)=="Yes"] <-"YL"
names(wide4)[names(wide4)=="No"] <-"NL"
wide4$latexresp <- log(wide4$YL/wide4$NL)
wide4$latexresp2 <- ((wide4$YL-wide4$NL)/wide4$YL)*100


wide5 <- dcast(m3, Pop+Fam~AMF,
               value.var="tric", mean)
names(wide5)[names(wide5)=="Yes"] <-"YTR"
names(wide5)[names(wide5)=="No"] <-"NTR"
wide5$tricresp <- log(wide5$YTR/wide5$NTR)
wide5$tricresp2 <- ((wide5$YTR-wide5$NTR)/wide5$YTR)*100


wide6 <- dcast(m3, Pop+Fam~AMF,
               value.var="Hgtday", mean)
names(wide6)[names(wide6)=="Yes"] <-"YH"
names(wide6)[names(wide6)=="No"] <-"NH"
wide6$Hgtresp <- log(wide6$YH/wide6$NH)
wide6$Hgtresp2 <- ((wide6$YH-wide6$NH)/wide6$YH)*100


w1 <- merge(wide1,wide2, by=c("Pop","Fam"))
w2 <- merge(w1,wide3, by=c("Pop","Fam"))
w3 <- merge(w2,wide4, by=c("Pop","Fam"))
w4 <- merge(w3,wide5, by=c("Pop","Fam"))
w5 <- merge(w4,wide6, by=c("Pop","Fam"))

w6 <- aggregate(cbind(amfresp2,RSresp2,slaresp2,latexresp2,tricresp2,Hgtresp2)~Pop, data=w5, mean)

library(Hmisc)
library(ape)
library(mvtnorm)
library(vegan)

###### compare two ways to calculate responsiveness
cor(w5$amfresp,w5$amfresp2)
cor(w5$slaresp,w5$slaresp2)

######  Calculate correlation matrix among responsiveness of traits
###### use percent responsiveness
respmat2 <- as.matrix(w5[,c("amfresp2","Hgtresp2","RSresp2","slaresp2","latexresp2","tricresp2")])
respcorr <- rcorr(respmat2)
respcorr


################################################################################
########## ANAYLSES FROM APPENDICES ############################################
################################################################################


##### correlations among traits in +AMF and -AMF treatments
fam <- aggregate(cbind(TotBMs,Hgtdays,RSs,slas,latexs,trics)~Pop+Fam+AMF, data=m3, mean)

m3Yes <- subset(fam, subset=AMF=="Yes")
m3Yes <- as.matrix(m3Yes[c("TotBMs","Hgtdays","RSs","slas","latexs","trics")])
rcorr(m3Yes)

m3No <- subset(fam, subset=AMF=="No")
m3No <- as.matrix(m3No[c("TotBMs","Hgtdays","RSs","slas","latexs","trics")])
rcorr(m3No)


#### multivariate model from Appendix 3
lm.3 <- lmer(value~AMF*Trait+(1|Pop)+(1|Pop:Trait:AMF),data=melt3,
             control=lmerControl(calc.derivs=F, optCtrl = list(maxfun=50000)))
anova(lm.3)