# Analyses of transplant survival using piecewiseSEM
# C. Lammers, 09-11-2023

# load packages ####
library(ggplot2)
library(dplyr)
library(ggpubr)
library(piecewiseSEM)
library(lme4)
library(lmerTest)
library(MuMIn)
library(car)

# import data ####
data<-read.csv("data_SEM.csv", sep=",")

# test if there is an optimum in survival with elevation change in winter and summer ####
data$winter_2<-data$diff_alt_winter^2 # add square of elevation change for calculation
data$summer_2<-data$diff_alt_summer^2 # add square of elevation change for calculation

#winter, all data
m1<-glmer(survival~diff_alt_winter+(1|location), family="binomial", data) #create glmer with original value
summary(m1) #check significance
m2<-glmer(survival~diff_alt_winter+winter_2+(1|location), family="binomial", data) #create glmer with both original and square value (unimodal distribution) 
 #check if both normal and square are signficant
summary(m2) # check significance
AIC(m1,m2) #compare AIC values, including unimodal better fit

#get coefficients of the distribution
beta_winter <- summary(m2)$coefficients[2, 1] 
beta_winter2 <- summary(m2)$coefficients[3, 1]
composite_winter <- beta_winter *data$diff_alt_winter + beta_winter2 * (data$diff_alt_winter)^2 #create new vector with composite values 
data$c_winter<--composite_winter #paste to dataframe ## change direction to make data interpretation more intuitive 

#subset data per location 
balg<-data[data$location=="BALG",] #subset for the exposed location
paal14<-data[data$location=="PAAL14",] #subset for the sheltered location

#repeat for exposed location 
m1<-glm(survival~diff_alt_winter, family="binomial", balg)
summary(m1)
m2<-glm(survival~diff_alt_winter+winter_2, family="binomial", balg)
summary(m2)
AIC(m1,m2) #compare AIC values, including unimodal better fit

beta_1<-summary(m2)$coefficients[2,1]
beta_2<-summary(m2)$coefficients[3,1]
composite_winter<-beta_1*balg$diff_alt_winter+beta_2*(balg$diff_alt_winter^2)
balg$c_winter<- -composite_winter ## change direction to make data interpretation more intuitive 

#repeat for sheltered location
m1<-glm(survival~diff_alt_winter, family="binomial",paal14)
summary(m1)
m2<-glm(survival~diff_alt_winter+winter_2, family="binomial", paal14)
summary(m2)
AIC(m1,m2)

beta1<-summary(m2)$coefficients[2,1]
beta2<-summary(m2)$coefficients[3,1]
c_winter<-beta1*paal14$diff_alt_winter+beta2*(paal14$diff_alt_winter^2)
#linear as good of a fit

#repeat for summer data, all 
m1<-glmer(survival~diff_alt_summer+(1|location), family="binomial", data)
m2<-glmer(survival~diff_alt_summer+summer_2+(1|location), family="binomial", data)
AIC(m1,m2)
beta_summer <- summary(m2)$coefficients[2, 1]
beta_summer2 <- summary(m2)$coefficients[3, 1]

composite_summer <- beta_summer * data$diff_alt_summer + beta_summer2 * (data$diff_alt_summer)^2
data$c_summer<- -composite_summer  ## change direction to make data interpretation more intuitive 

#repeat for exposed location
m1<-glm(survival~diff_alt_summer, family="binomial", balg)
summary(m1)
m2<-glm(survival~diff_alt_summer+summer_2, family="binomial", balg)
summary(m2)
AIC(m1,m2)

beta_1<-summary(m2)$coefficients[2,1]
beta_2<-summary(m2)$coefficients[3,1]
composite_summer<-beta_1*balg$diff_alt_summer+beta_2*(balg$diff_alt_summer^2)
balg$c_summer<- -composite_summer

#repeat for sheltered location
m1<-glm(survival~diff_alt_summer, family="binomial", paal14)
m2<-glm(survival~diff_alt_summer+summer_2, family="binomial", paal14)
AIC(m1,m2)

# (generalized) linear mixed effect models ALL data####
data$test_distance<-data$distance_northsea/100 #from m to hectometer to have a similar range of values

#elevation vs sand couch presence and distance to north sea
m_alt<-lmer(log(alt_start)~test_distance+ely+(1|location), data)
summary(m_alt)
r.squaredGLMM(m_alt)
Anova(m_alt, test.statistic = "F")

#elevation change in winter vs sand couch presence, elevation and distance
m_winter1<-lmer(diff_alt_winter~ely+log(alt_start)+test_distance+(1|location), data)
m_winter<-lmer(diff_alt_winter~log(alt_start)+test_distance+(1|location), data)
summary(m_winter)
r.squaredGLMM(m_winter)
Anova(m_winter, test.statistic = "F")

#elevation change in summer vs sand couch presence, elevation and distance
m_summer<-lmer(diff_alt_summer~ely+log(alt_start)+test_distance+(1|location), data)
summary(m_summer)

# transplant survival vs sand couch presence, elevation, elevation change (optima) and distance to North Sea
m_survival<-glmer(survival~ely+log(alt_start)+c_summer+c_winter+test_distance+(1|location), family="binomial", df_removed) # full model
m_survival1<-glmer(survival~c_summer+c_winter+test_distance+(1|location), family="binomial", data) #backward selected model
summary(m_survival1)
r.squaredGLMM(m_survival)
Anova(m_survival1, test.statistic = "Chisq")

# models including optimum values of elevation change for PSEM
m_winter2<-lmer(c_winter~log(alt_start)+test_distance+(1|location), data)
summary(m_winter2)

# piecewiseSEM, ALL data####
## start with full models
m<-psem(
  m_alt,
  lmer(c_winter~ely+log(alt_start)+test_distance+(1|location),data),
  lmer(c_summer~ely+log(alt_start)+test_distance+(1|location),data),
  m_survival,
  c_winter%~~%c_summer
)
summary(m)
plot(m)
### after backward selection, final model
m1<-psem(
  m_alt,
  m_winter2, 
  c_winter%~~%c_summer,
  m_survival1)
plot(m1)
summary(m1, conserve=T)

#mixed logistic regression models EXPOSED ####
m_survival<-glm(survival~ely+log(alt_start)+test_distance+c_winter+c_summer, family="binomial", balg) # full model
summary(m_survival)
m_survival<-glm(survival~test_distance+c_winter+c_summer, family="binomial", balg) # selected model
summary(m_survival)

library(rcompanion)
nagelkerke(m_survival)

m_winter<-lm(c_winter~ely+log(alt_start)+test_distance, balg) #full model
summary(m_winter)
m_winter<-lm(c_winter~log(alt_start)+test_distance, balg) #selected model
summary(m_winter)

m_summer<-lm(c_summer~log(alt_start)+test_distance+ely, balg)
summary(m_summer) # not significant

m_alt<-lm(log(alt_start)~test_distance+ely, balg)
summary(m_alt)

#piecewiseSEM EXPOSED####
#started with full model, see piecewiseSEM all data, here selected SEM
m1<-psem(
  m_alt,
  m_winter, 
  abs_winter%~~%abs_summer,
  m_survival
)
plot(m1)
summary(m1)

#mixed logistic regression models SHELTERED####
m_survival<-glm(survival~alt_start+test_distance+ely+c_winter+c_summer ,family="binomial", paal14) #full model, no log transformed alt needed
summary(m_survival)
m_survival<-glm(survival~alt_start+test_distance+ely+c_winter ,family="binomial", paal14) #selected model, no log transformed alt needed
summary(m_survival)

nagelkerke(m_survival)

m_winter<-lm(diff_alt_winter~alt_start+test_distance, paal14) #selected model
summary(m_winter)

m_summer<-lm(diff_alt_summer~test_distance+ely, paal14) #selected model
summary(m_summer) # not included in piecewiseSEM because elevation change in summer did not signficantly relate to survival

m_alt<-lm(alt_start~test_distance+ely, paal14)
summary(m_alt)

#piecewiseSEM SHELTERED####
#started with full model, see piecewiseSEM all data, here selected SEM
m1<-psem(
  m_alt,
  m_winter,
  m_survival
)
plot(m1)
summary(m1)