In TheoreticalEcosystemEcology/alienR: Comparative analysis of interaction inference methods

knitr::opts_chunk$set(echo = TRUE)

This vignette was made to be an example of how the package can be applied to real-world data. As such, we focused on the Salix-galler-parasitoid data published by Kopelke et al. (2017) for which additional information on phylogeny was added. The full data was used to perform the analysis in Wooton et al. (2021).

Data

The data is part of the alien R package and has already been formatted.

library(alien)

data(salixGal)
data(galPara)

#treeGal <- read.nexus("./Tree_with_added_taxa_with_short_codes_from_Mesquite.nex")

#cophGal <- cophenetic(treeGal)
#phyloGal <- as.matrix(salixGal$phyloDistTo)
#sel <- which(rownames(cophGal) %in% rownames(phyloGal))
#cophGalSel <- cophGal[sel,sel]
#cophGalOrd <- cophGalSel[sort(rownames(cophGalSel)),sort(rownames(cophGalSel))]

#salixGal$phyloDistTo <- as.dist(cophGalOrd)
#galPara$phyloDistFrom <- as.dist(cophGalOrd)

fitKNN

Interaction data only

# Salix-galler
salixGalKNNInter <- fitKNN(salixGal,
                           distFrom = "jaccard",
                           distTo = "jaccard",
                           nNeig = 3,
                           weight = 0)

logLik(salixGalKNNInter, error = 0.001)
tjur(salixGalKNNInter) 

# Galler-parasitoid
galParaKNNInter <- fitKNN(galPara,
                          distFrom = "jaccard",
                          distTo = "jaccard",
                          nNeig = 3,
                          weight = 0)

logLik(galParaKNNInter, error = 0.001)
tjur(galParaKNNInter)

Interaction data with traits

# Salix-galler
salixGalKNNInterTrait <- fitKNN(salixGal,
                                distFrom = "jaccard",
                                distTo = "jaccard",
                                distTraitFrom = "euclidean",
                                distTraitTo = "euclidean",
                                nNeig = 3,
                                weight = 1)

logLik(salixGalKNNInterTrait, error = 0.001)
tjur(salixGalKNNInterTrait)

# Galler-parasitoid
galParaKNNInterTrait <- fitKNN(galPara,
                               distFrom = "jaccard",
                               distTo = "jaccard",
                               distTraitFrom = "euclidean",
                               distTraitTo = "euclidean",
                               nNeig = 3,
                               weight = 1)

logLik(galParaKNNInterTrait, error = 0.001)
tjur(galParaKNNInterTrait)

Interaction data with phylogeny

# Salix-galler
salixGalKNNInterPhylo <- fitKNN(salixGal,
                                distFrom = "jaccard",
                                distTo = "jaccard",
                                nNeig = 3,
                                weight = 1,
                                phylo=TRUE)

logLik(salixGalKNNInterPhylo, error = 0.001)
tjur(salixGalKNNInterPhylo)

# Galler-parasitoid
galParaKNNInterPhylo <- fitKNN(galPara,
                               distFrom = "jaccard",
                               distTo = "jaccard",
                               nNeig = 3,
                               weight = 1,
                               phylo = TRUE)

logLik(galParaKNNInterPhylo, error = 0.001)
tjur(galParaKNNInterPhylo)

fitRF

Traits

##############
# Salix-galler
##############

# Base of formula
salTraits <- paste(colnames(salixGal$traitFrom), collapse = "+")
galTraits <- paste(colnames(salixGal$traitTo), collapse = "+")

# Formula
FormulaSalGal <- as.formula(paste("~", salTraits, "+", galTraits))

# Model
salixGalRFTrait <- fitRF(salixGal,
                         formula = FormulaSalGal,
                         ntree = 2000,
                         nodesize = 1)

# Loglikelihood and Tjur's D
logLik(salixGalRFTrait, error = 0.001)
tjur(salixGalRFTrait)

###################
# Galler-parasitoid
###################
# Base of formula
galTraits <- paste(colnames(galPara$traitFrom), collapse = "+")
parTraits <- paste(colnames(galPara$traitTo), collapse = "+")

# Formula
FormulaGalPara <- as.formula(paste("~", galTraits, "+", parTraits))

# Model
galParaRFTrait <- fitRF(galPara,
                        formula = FormulaGalPara,
                        ntree=2000,
                        nodesize = 1)

# Loglikelihood and Tjur's D
logLik(galParaRFTrait, error = 0.001)
tjur(galParaRFTrait)

Phylo

library(vegan)
salPhylo <- wcmdscale(salixGal$phyloDistFrom, add = TRUE)
colnames(salPhylo) <- paste0("salPhylo", 1:ncol(salPhylo))

galPhylo <- wcmdscale(salixGal$phyloDistTo, add = TRUE)
colnames(galPhylo) <- paste0("galPhylo", 1:ncol(galPhylo))

paraPhylo <- wcmdscale(galPara$phyloDistTo, add = TRUE)
colnames(paraPhylo) <- paste0("paraPhylo", 1:ncol(paraPhylo))

##############
# Salix-galler
##############
salixGalPhylo <- alienData(adjMat = salixGal$adjMat,
                           traitFrom = as.data.frame(salPhylo),
                           traitTo = as.data.frame(galPhylo))

# Base of formula
salPhyloTraits <- paste(colnames(salixGalPhylo$traitFrom), collapse = "+")
galPhyloTraits <- paste(colnames(salixGalPhylo$traitTo), collapse = "+")

# Formula
FormulaSalGalPhylo <- as.formula(paste("~", salPhyloTraits, "+", galPhyloTraits))

# Model
salixGalRFPhylo <- fitRF(salixGalPhylo,
                         formula = FormulaSalGalPhylo,
                         ntree = 2000,
                         nodesize = 1)

# Loglikelihood and Tjur's D
logLik(salixGalRFPhylo, error = 0.001)
tjur(salixGalRFPhylo)

###################
# Galler-parasitoid
###################
galParaPhylo <- alienData(adjMat = galPara$adjMat,
                           traitFrom = as.data.frame(galPhylo),
                           traitTo = as.data.frame(paraPhylo))

# Base of formula
galPhyloTraits <- paste(colnames(galParaPhylo$traitFrom), collapse = "+")
paraPhyloTraits <- paste(colnames(galParaPhylo$traitTo), collapse = "+")

# Formula
FormulaGalParaPhylo <- as.formula(paste("~", galPhyloTraits, "+", paraPhyloTraits))

# Model
galParaRFPhylo <- fitRF(galParaPhylo,
                        formula = FormulaGalParaPhylo,
                        ntree=2000,
                        nodesize = 1)

# Loglikelihood and Tjur's D
logLik(galParaRFPhylo, error = 0.001)
tjur(galParaRFPhylo)

fitGLM

################
# Salix - Galler
################
# Construct squared traits for numerical traits
salixGalSq <- polyTrait(salixGal)

#==============
# Build formula
#==============
# Select traits
salixTr <- c("TREE.VOLUME",
             "TREE.VOLUME_Sq",
             "LEAF.THICKNESS",
             "GLUCOSIDES",
             "GLUCOSIDES_Sq",
             "TREE.HEIGHT",
             "TREE.HEIGHT_Sq",
             "LEAF.HAIRINESS")

gallTr <- "GALLTYPE"

# Independent term
indep <- paste(c(salixTr, gallTr), collapse = "+")

# Interaction term (without Square)
salixTrNoSq <- c("TREE.VOLUME",
                 "LEAF.THICKNESS",
                 "GLUCOSIDES",
                 "TREE.HEIGHT",
                 "LEAF.HAIRINESS")

gallTrNoSq <- c("GALLTYPE")

combTr <- combn(c(salixTrNoSq, gallTrNoSq), 2)
inter <- paste(combTr[1,], ":", combTr[2,], collapse = "+")

# Build formula
FormulaSalGal <- as.formula(paste("~", indep, "+", inter))

#==========
# Run model
#==========
salixGalGLMTrait <- fitGLM(salixGalSq,
                           formula = FormulaSalGal,
                           family = binomial(link = "logit"))

#============================
# Log-likelihood and Tjur's D
#============================
logLik(salixGalGLMTrait, error = 0.001)
tjur(salixGalGLMTrait)

#####################
# Galler - Parasitoid
#####################
# Construct squared traits for numerical traits
galParaSq <- polyTrait(galPara)

#==============
# Build formula
#==============
# Select traits
gallTr <- c("GALLTYPE", 
            "BODYLENGTH.GAL", 
            "BODYLENGTH.GAL_Sq", 
            "PHENOLOGY.GAL", 
            "OVIPOS.STRATEGY",
            "GALL.WALL",
            "GALL.WALL_Sq")

paraTr <- c("P.I",  
            "OVIPOS.LNTH", 
            "OVIPOS.LNTH_Sq", 
            "ATTACK.STAGE",
            "PHENOLOGY.PAR", 
            "BODYLENGTH.PAR", 
            "BODYLENGTH.PAR_Sq", 
            "ENDO.ECTO")

# Independent term
indep <- paste(c(gallTr, paraTr), collapse = "+")

# Interaction term (without Square)
gallTrNoSq <- c("GALLTYPE", 
                "BODYLENGTH.GAL", 
                "PHENOLOGY.GAL", 
                "OVIPOS.STRATEGY",
                "GALL.WALL")

paraTrNoSq <- c("P.I",  
                "OVIPOS.LNTH", 
                "ATTACK.STAGE",
                "PHENOLOGY.PAR", 
                "BODYLENGTH.PAR", 
                "ENDO.ECTO")

combTr <- combn(c(gallTrNoSq, paraTrNoSq), 2)
inter <- paste(combTr[1,], ":", combTr[2,], collapse = "+")

# Build formula
FormulaGalPara <- as.formula(paste("~", indep, "+", inter))

#==========
# Run model
#==========
galParaGLMTrait <- fitGLM(galParaSq,
                           formula = FormulaGalPara,
                           family = binomial(link = "logit"))

#============================
# Log-likelihood and Tjur's D
#============================
logLik(galParaGLMTrait, error = 0.001)
tjur(galParaGLMTrait)

fit4corner

# Construct squared traits for numerical traits
salixGalSq <- polyTrait(salixGal)

#==============
# Build formula
#==============
# Select traits
salixTr <- c("TREE.VOLUME",
             "TREE.VOLUME_Sq",
             "LEAF.THICKNESS",
             "GLUCOSIDES",
             "GLUCOSIDES_Sq",
             "TREE.HEIGHT",
             "TREE.HEIGHT_Sq",
             "LEAF.HAIRINESS")

gallTr <- "GALLTYPE"

# Independent term
indep <- paste(c(salixTr, gallTr), collapse = "+")

# Interaction term (without Square)
salixTrNoSq <- c("TREE.VOLUME",
                 "LEAF.THICKNESS",
                 "GLUCOSIDES",
                 "TREE.HEIGHT",
                 "LEAF.HAIRINESS")

gallTrNoSq <- c("GALLTYPE")

combTr <- expand.grid(salixTrNoSq, gallTrNoSq)
inter <- paste(combTr[,1], ":", combTr[,2], collapse = "+")

# Build formula
FormulaSalGal <- as.formula(paste("~", indep, "+", inter))

#==========
# Run model
#==========
salixGal4cornerTrait <- fitGLM(salixGalSq,
                           formula = FormulaSalGal,
                           family = binomial(link = "logit"))

#============================
# Log-likelihood and Tjur's D
#============================
logLik(salixGal4cornerTrait, error = 0.001)
tjur(salixGal4cornerTrait)


#####################
# Galler - Parasitoid
#####################
# Construct squared traits for numerical traits
galParaSq <- polyTrait(galPara)

#==============
# Build formula
#==============
# Select traits
gallTr <- c("GALLTYPE", 
            "BODYLENGTH.GAL", 
            "BODYLENGTH.GAL_Sq", 
            "PHENOLOGY.GAL", 
            "OVIPOS.STRATEGY",
            "GALL.WALL",
            "GALL.WALL_Sq")

paraTr <- c("P.I",  
            "OVIPOS.LNTH", 
            "OVIPOS.LNTH_Sq", 
            "ATTACK.STAGE",
            "PHENOLOGY.PAR", 
            "BODYLENGTH.PAR", 
            "BODYLENGTH.PAR_Sq", 
            "ENDO.ECTO")

# Independent term
indep <- paste(c(gallTr, paraTr), collapse = "+")

# Interaction term (without Square)
gallTrNoSq <- c("GALLTYPE", 
                "BODYLENGTH.GAL", 
                "PHENOLOGY.GAL", 
                "OVIPOS.STRATEGY",
                "GALL.WALL")

paraTrNoSq <- c("P.I",  
                "OVIPOS.LNTH", 
                "ATTACK.STAGE",
                "PHENOLOGY.PAR", 
                "BODYLENGTH.PAR", 
                "ENDO.ECTO")

combTr <- expand.grid(gallTrNoSq, paraTrNoSq)
inter <- paste(combTr[,1], ":", combTr[,2], collapse = "+")

# Build formula
FormulaGalPara <- as.formula(paste("~", indep, "+", inter))

#==========
# Run model
#==========
galParaGLMTrait <- fitGLM(galParaSq,
                           formula = FormulaGalPara,
                           family = binomial(link = "logit"))

#============================
# Log-likelihood and Tjur's D
#============================
logLik(galParaGLMTrait, error = 0.001)
tjur(galParaGLMTrait)

fitIMC

# Salix-galler
salixGalIMC <- readRDS("salixGalIMC.RDS")

logLik(salixGalIMC, error = 0.001) # -137.0956
tjur(salixGalIMC) # 0.526749

# Galler-parasitoids
galParaIMC <- readRDS("galParaIMC.RDS")

logLik(galParaIMC, error = 0.001) # -690.649
tjur(galParaIMC) # 0.5714072

This function take a long time to run.

# Salix-galler
salixGalIMC <- fitIMC(salixGal,
                      d = 1)

logLik(salixGalIMC, error = 0.001) # -137.0956
tjur(salixGalIMC) # 0.526749

# Galler-parasitoids
galParaIMC <- fitIMC(galPara,
                      d = 1)

logLik(galParaIMC, error = 0.001) # -690.649
tjur(galParaIMC) # 0.5714072

TheoreticalEcosystemEcology/alienR documentation built on Dec. 25, 2021, 5:59 p.m.