– All analyses were confirmed by the corresponding author and finally conducted at 2023/05/09
Checking the version of R itself and rEDM (v0.6.9)
#remotes::install_github("ha0ye/rEDM")
#install.packages("C:/Users/*****/rEDM_0.6.9.tar.gz", repos = NULL, type = "source")
library(dplyr)
次のパッケージを付け加えます: ‘dplyr’
以下のオブジェクトは ‘package:stats’ からマスクされています:
filter, lag
以下のオブジェクトは ‘package:base’ からマスクされています:
intersect, setdiff, setequal, union
library(rEDM)
If you're new to the rEDM package, please check out
the tutorial:
> vignette("rEDM-tutorial")
library(ggplot2)
library(ggeffects)
library(MASS)
次のパッケージを付け加えます: ‘MASS’
以下のオブジェクトは ‘package:dplyr’ からマスクされています:
select
sessionInfo()
R version 4.2.0 (2022-04-22 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)
Matrix products: default
locale:
[1] LC_COLLATE=Japanese_Japan.utf8 LC_CTYPE=Japanese_Japan.utf8
[3] LC_MONETARY=Japanese_Japan.utf8 LC_NUMERIC=C
[5] LC_TIME=Japanese_Japan.utf8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] MuMIn_1.47.5 MASS_7.3-58.3 ggeffects_1.2.1 ggplot2_3.4.1
[5] rEDM_0.6.9 dplyr_1.1.0
loaded via a namespace (and not attached):
[1] Rcpp_1.0.10 rstudioapi_0.14 knitr_1.42 magrittr_2.0.3
[5] tidyselect_1.2.0 munsell_0.5.0 lattice_0.20-45 colorspace_2.1-0
[9] R6_2.5.1 rlang_1.0.6 fansi_1.0.4 tools_4.2.0
[13] grid_4.2.0 nlme_3.1-157 gtable_0.3.1 xfun_0.37
[17] utf8_1.2.2 cli_3.6.0 withr_2.5.0 tibble_3.1.8
[21] lifecycle_1.0.3 Matrix_1.5-3 vctrs_0.5.2 codetools_0.2-18
[25] glue_1.6.2 compiler_4.2.0 pillar_1.8.1 generics_0.1.3
[29] scales_1.2.1 stats4_4.2.0 pkgconfig_2.0.3
packageVersion("rEDM")
[1] ‘0.6.9’
– ヨツモンマメゾウムシ = ヨツモン = Yotsumon = Callosobrucbus maculatus = Cm = “Y”
– アズキゾウムシ = アズキ = Azuki = Callosobruchus chinensis = Cc = “A”
#file from 2022/02/22 reanalysis
data <- read.csv("data_20230222en.csv", header = T)
#data <- read.csv("data_20230222.csv", header = T)
data
Smallvolume <- 8.5*6*3.5
Mediumvolume <- 13.5*8.5*4.5
Largevolume <- 18.5*14*6.5
Hugevolume <- 22*15.5*7.5
no_replicates <- 64 #the total number of replicates: 16 replicates within a cage size x four cage sizes
#converting the data to list for each cage ID
datalist <- lapply(1:no_replicates, function(i) subset(data,ID == i))
datalist[[10]] #example
NA
#Counting the number of weeks when the abundance is continuously zero , back-calculated from the final week to the first week
count <- numeric(no_replicates) #numeric vector for 64 samples
weeks <- numeric(no_replicates) #number of weeks for observations, which potentially dependend on replicates
for(j in 1:no_replicates){
weeks[j] <- nrow(datalist[[j]])
for(i in weeks[j]:1){
if(datalist[[j]]$YotsumonTotalL[i] == 0){
count[j] <- count[j] + 1
} else{
break
}
}
}
count
[1] 0 1 0 0 0 0 0 0 0 0 0 0 0 6 0 0 0 0 0 0 0 0 0
[25] 0 3 1 0 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 15 0
[49] 0 0 0 0 0 0 7 0 0 0 0 0 0 0 2 0
#Counting the number of weeks when the abundance is continuously zero , back-calculated from the final week to the first week
count2 <- numeric(no_replicates) #numeric vector for 64 samples
for(j in 1:no_replicates){
for(i in weeks[j]:1){
if(datalist[[j]]$AzukiTotalL[i] == 0){
count2[j] <- count2[j] + 1
} else{
break
}
}
}
count2
[1] 6 6 5 5 5 5 5 5 5 5 5 5 5 0 12 5 5 5 7 5 7 5 5 5
[25] 5 7 5 7 0 7 7 5 5 5 5 5 5 0 5 5 5 5 5 5 5 5 0 5
[49] 5 5 6 0 7 12 1 5 5 5 5 5 5 5 6 5
The code by Miki, it is more natural and simpler.
extinct <- numeric(no_replicates)
for(i in 1:no_replicates){
if(count[i] >= count2[i]) {
extinct[i] <- count[i]
} else {
extinct[i] <- count2[i]
}
}
extinct
[1] 6 6 5 5 5 5 5 5 5 5 5 5 5 6 12 5 5 5 7 5 7 5 5 5
[25] 5 7 5 7 6 7 7 5 5 5 5 5 5 0 5 5 5 5 5 5 5 5 15 5
[49] 5 5 6 0 7 12 7 5 5 5 5 5 5 5 6 5
shavedata <- lapply(1:64, function(i) head(datalist[[i]], n = weeks[i] - extinct[i] + 1))
shavedataframe <- NULL
for(i in 1:64) shavedataframe <- rbind.data.frame(shavedataframe, shavedata[[i]])
shavedataframe
#Small
Smallall <- subset(shavedataframe, Cage.size == "Small")
SmallID <- unique(Smallall$ID)
Small <- list()
for(i in 1:length(SmallID)){
Small[[i]] <- subset(Smallall, ID == SmallID[i])
}
Small[[2]] #example
#Medium
Mediumall <- subset(shavedataframe, Cage.size == "Medium")
MediumID <- unique(Mediumall$ID)
Medium <- list()
for(i in 1:length(MediumID)){
Medium[[i]] <- subset(Mediumall, ID == MediumID[i])
}
Medium[[2]] #example
#Large
Largeall <- subset(shavedataframe, Cage.size == "Large")
LargeID <- unique(Largeall$ID)
Large <- list()
for(i in 1:length(LargeID)){
Large[[i]] <- subset(Largeall, ID == LargeID[i])
}
Large[[2]] #example
#Huge
Hugeall <- subset(shavedataframe, Cage.size == "Huge")
HugeID <- unique(Hugeall$ID)
Huge <- list()
for(i in 1:length(HugeID)){
Huge[[i]] <-subset(Hugeall, ID == HugeID[i])
}
Huge[[2]] #example
#Small
Smallstrain <- NULL
for(i in 1:length(Small)){
Smallstrain[i] <- unique(Small[[i]]$Strain)
}
Smallstrain
[1] "jc-F" "msC" "jc-F" "msC" "jc-F" "jc-F" "jc-F" "jc-F" "msC" "msC"
[11] "msC" "msC" "jc-F" "jc-F" "msC" "msC"
#Medium
Mediumstrain <- NULL
for(i in 1:length(Medium)){
Mediumstrain[i] <- unique(Medium[[i]]$Strain)
}
Mediumstrain
[1] "jc-F" "msC" "msC" "jc-F" "jc-F" "jc-F" "jc-F" "jc-F" "msC" "msC"
[11] "msC" "msC" "jc-F" "jc-F" "msC" "msC"
#Large
Largestrain <- NULL
for(i in 1:length(Large)){
Largestrain[i] <- unique(Large[[i]]$Strain)
}
Largestrain
[1] "jc-F" "msC" "msC" "jc-F" "msC" "jc-F" "jc-F" "jc-F" "msC" "msC"
[11] "msC" "msC" "jc-F" "jc-F" "msC" "msC"
#Huge
Hugestrain <- NULL
for(i in 1:length(Huge)){
Hugestrain[i] <- unique(Huge[[i]]$Strain)
}
Hugestrain
[1] "jc-F" "msC" "msC" "jc-F" "jc-F" "jc-F" "jc-F" "jc-F" "msC" "msC"
[11] "msC" "jc-F" "jc-F" "msC" "msC" "msC"
write.csv(shavedataframe, "shavedataframe.csv")
Smallcoweek <- vector()
for(i in 1:length(Small)){
Smallcoweek[i] <- nrow(Small[[i]]) - 1 #Extinction occurred at the final week
}
cat("The number of replicates for small size cage:", length(Small), "\n")
The number of replicates for small size cage: 16
Smallcoweek
[1] 25 13 21 9 12 27 23 18 9 15 12 11 11 11 8 9
Mediumcoweek <- vector()
for(i in 1:length(Medium)){
Mediumcoweek[i] <- nrow(Medium[[i]]) - 1 #Extinction occurred at the final week
}
cat("The number of replicates for medium size cage:", length(Medium), "\n")
The number of replicates for medium size cage: 16
Mediumcoweek
[1] 21 17 9 14 28 31 28 19 15 46 19 15 15 24 11 9
Largecoweek <- vector()
for(i in 1:length(Large)){
Largecoweek[i] <- nrow(Large[[i]]) - 1 #Extinction occurred at the final week
}
cat("The number of replicates for large size cage:", length(Large), "\n")
The number of replicates for large size cage: 16
Largecoweek
[1] 33 18 18 36 15 28 27 28 15 19 19 19 42 23 15 15
Hugecoweek <- vector()
for(i in 1:length(Huge)){
Hugecoweek[i] <- nrow(Huge[[i]]) - 1 #Extinction occurred at the final week
}
cat("The number of replicates for huge size cage:", length(Huge), "\n")
The number of replicates for huge size cage: 16
Hugecoweek
[1] 34 15 19 30 30 28 28 27 27 27 32 31 23 15 23 16
– count: Cm Yotsumon, the number of weeks with continuously zero abundance
– count2: Cc Azuki, the number of weeks with continuously zero abundance
competitionwinner <- vector()
for(i in 1:no_replicates){
if(count[i] < count2[i]){
competitionwinner[i] <- "C.m"
} else if (count[i] == count2[i]){
competitionwinner[i] <- "Coexistence"
} else{
competitionwinner[i] <- "C.c"
}
}
competitionwinner
[1] "C.m" "C.m" "C.m" "C.m" "C.m"
[6] "C.m" "C.m" "C.m" "C.m" "C.m"
[11] "C.m" "C.m" "C.m" "C.c" "C.m"
[16] "C.m" "C.m" "C.m" "C.m" "C.m"
[21] "C.m" "C.m" "C.m" "C.m" "C.m"
[26] "C.m" "C.m" "C.m" "C.c" "C.m"
[31] "C.m" "C.m" "C.m" "C.m" "C.m"
[36] "C.m" "C.m" "Coexistence" "C.m" "C.m"
[41] "C.m" "C.m" "C.m" "C.m" "C.m"
[46] "C.m" "C.c" "C.m" "C.m" "C.m"
[51] "C.m" "Coexistence" "C.m" "C.m" "C.c"
[56] "C.m" "C.m" "C.m" "C.m" "C.m"
[61] "C.m" "C.m" "C.m" "C.m"
TableS1 <- data.frame("Cage Size" = c(rep("Small",length(Small)),
rep("Medium",length(Medium)),
rep("Large",length(Large)),
rep("Huge",length(Huge))),
"ID" = c(SmallID,MediumID,LargeID,HugeID),
"Strain" = c(Smallstrain,Mediumstrain,Largestrain,Hugestrain),
"Coexistence_periods" = c(Smallcoweek,Mediumcoweek,Largecoweek,Hugecoweek),
"Result_of_competition" = competitionwinner)
TableS1
write.csv(TableS1,"TableS1.csv")
coexist_size_dataframe <- data.frame("Coexistence_periods" = c(Smallcoweek,Mediumcoweek,Largecoweek,Hugecoweek),
"Volume" = c(rep(Smallvolume,length(Smallcoweek)),rep(Mediumvolume,length(Mediumcoweek)),rep(Largevolume,length(Largecoweek)),rep(Hugevolume,length(Hugecoweek))),
"Cage_size" = as.factor(c(rep("Small",length(Smallcoweek)),rep("Medium",length(Mediumcoweek)),rep("Large",length(Largecoweek)),rep("Huge",length(Hugecoweek)))))
#check the class of vectors
class(coexist_size_dataframe$Volume)
[1] "numeric"
class(coexist_size_dataframe$Cage_size)
[1] "factor"
#Change the order of factor levels
coexist_size_dataframe$Cage_size <- factor(coexist_size_dataframe$Cage_size,levels = c("Small","Medium","Large","Huge"))
coexist_size_dataframe$Cage_size
[1] Small Small Small Small Small Small Small Small Small Small
[11] Small Small Small Small Small Small Medium Medium Medium Medium
[21] Medium Medium Medium Medium Medium Medium Medium Medium Medium Medium
[31] Medium Medium Large Large Large Large Large Large Large Large
[41] Large Large Large Large Large Large Large Large Huge Huge
[51] Huge Huge Huge Huge Huge Huge Huge Huge Huge Huge
[61] Huge Huge Huge Huge
Levels: Small Medium Large Huge
coexist_size_dataframe
#With mean
Fig2 <- ggplot(data = coexist_size_dataframe,
mapping = aes(x = Volume, y= Coexistence_periods, group = Cage_size,fill = Cage_size)) +
geom_violin() +
geom_boxplot(width = 30, fill = "black", outlier.color = NA) +
stat_summary(fun = mean, geom = "point", fill = "white", shape = 21, size = 3) +
labs(y = "Coexistence periods (week)",x = "Volume (cm³)")
print(Fig2)
#Simpler plot with GLM line
Fig2.1 <- ggplot(data = coexist_size_dataframe, aes(x = Volume, y = Coexistence_periods))
Fig2.1 <- Fig2.1 + geom_point(size = 3, shape = 21)
Fig2.1 <- Fig2.1 + geom_smooth(method = 'glm', method.args = list(family = 'poisson'), se = T, color = 'black')
Fig2.1 <- Fig2.1 + labs(y = "Coexistence periods (week)", x = "Volume (cm3)")
print(Fig2.1)
plot(
Coexistence_periods ~ Volume,
data = coexist_size_dataframe,
xlab = "VOlume (cm3)",
ylab = "Coexitence period (weeks)"
)
summary(glm(Coexistence_periods ~ Volume, data = coexist_size_dataframe, family = "poisson"))
Call:
glm(formula = Coexistence_periods ~ Volume, family = "poisson",
data = coexist_size_dataframe)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.3948 -1.5124 -0.5415 0.9776 5.5354
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.788e+00 4.871e-02 57.243 < 2e-16 ***
Volume 1.863e-04 2.869e-05 6.492 8.48e-11 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 218.47 on 63 degrees of freedom
Residual deviance: 176.42 on 62 degrees of freedom
AIC: 487.32
Number of Fisher Scoring iterations: 4
– ヨツモンマメゾウムシ = ヨツモン = Yotsumon = Callosobrucbus maculatus = Cm = “Y”
– アズキゾウムシ = アズキ = Azuki = Callosobruchus chinensis = Cc = “A”
– Sorting for each cage size, from longest coexistence period to shortest
#Small
Smalllength <- subset(TableS1, TableS1$Cage.Size == "Small")
Smalllength[order(Smalllength$Coexistence_periods, decreasing = TRUE),]
#Medium
Mediumlength <- subset(TableS1, TableS1$Cage.Size == "Medium")
Mediumlength[order(Mediumlength$Coexistence_periods, decreasing = TRUE),]
#Large
Largelength <- subset(TableS1, TableS1$Cage.Size == "Large")
Largelength[order(Largelength$Coexistence_periods, decreasing = TRUE),]
#Huge
Hugelength <- subset(TableS1,TableS1$Cage.Size == "Huge")
Hugelength[order(Hugelength$Coexistence_periods, decreasing = TRUE),]
– The goal is to obtain the optimal embedding dimension E*
– Warning when executing simplex() is expected one so no problems
Simplex_Small_Y <- list() #list for assigning results
#Executing simplex()
for(i in 1:length(Small)){
Simplex_Small_Y[[i]] <- simplex(Small[[i]]$YotsumonTotalL, lib = c(1, nrow(Small[[i]])), pred = c (1, nrow(Small[[i]])), E = c(1:10))
}
#Picking up the optimal embedding dimension (E*)
Small_BestE_Y <- vector()
for(i in 1:length(Simplex_Small_Y)){
Small_BestE_Y[i] <- which.max(Simplex_Small_Y[[i]]$rho)
}
Small_BestE_Y
[1] 9 4 4 1 2 8 8 5 1 5 4 4 1 4 1 1
#Combining it with replicate ID
Small_ID_BestE_Y <- NULL
for(i in 1:length(Small)){
Small_ID_BestE_Y <- rbind.data.frame(Small_ID_BestE_Y, cbind(Small[[i]][1, ]$ID, Small_BestE_Y[i]))
}
#adding column names
colnames(Small_ID_BestE_Y) <- c("ID","C.m_E")
Small_ID_BestE_Y
Simplex_Small_A <- list() #list for assigning results
#Executing simplex()
for(i in 1:length(Small)){
Simplex_Small_A[[i]] <- simplex(Small[[i]]$AzukiTotalL, lib = c(1, nrow(Small[[i]])), pred = c (1, nrow(Small[[i]])), E = c(1:10))
}
#Picking up the optimal embedding dimension (E*)
Small_BestE_A <- vector()
for(i in 1:length(Simplex_Small_A)){
Small_BestE_A[i] <- which.max(Simplex_Small_A[[i]]$rho)
}
Small_BestE_A
[1] 5 6 4 2 1 5 5 3 3 5 5 1 1 2 1 3
#Combining it with replicate ID
Small_ID_BestE_A <- NULL
for(i in 1:length(Small)){
Small_ID_BestE_A <- rbind.data.frame(Small_ID_BestE_A, cbind(Small[[i]][1, ]$ID, Small_BestE_A[i]))
}
#adding column names
colnames(Small_ID_BestE_A) <- c("ID","C.c_E")
Small_ID_BestE_A
Simplex_Medium_Y <- list() #list for assigning results
#Executing simplex()
for(i in 1:length(Medium)){
Simplex_Medium_Y[[i]] <- simplex(Medium[[i]]$YotsumonTotalL, lib = c(1, nrow(Medium[[i]])), pred = c (1, nrow(Medium[[i]])), E = c(1:10))
}
#Picking up the optimal embedding dimension (E*)
Medium_BestE_Y <- vector()
for(i in 1:length(Simplex_Medium_Y)){
Medium_BestE_Y[i] <- which.max(Simplex_Medium_Y[[i]]$rho)
}
Medium_BestE_Y
[1] 4 4 1 4 10 7 5 4 5 7 1 1 4 6 4 1
for(i in 1:length(Medium)){
Medium_ID_BestE_Y <- rbind.data.frame(Medium_ID_BestE_Y, cbind(Medium[[i]][1, ]$ID, Medium_BestE_Y[i]))
}
#adding column names
colnames(Medium_ID_BestE_Y) <- c("ID","C.m_E")
Medium_ID_BestE_Y
Simplex_Medium_A <- list() #list for assigning results
#Executing simplex()
for(i in 1:length(Medium)){
Simplex_Medium_A[[i]] <- simplex(Medium[[i]]$AzukiTotalL, lib = c(1, nrow(Medium[[i]])), pred = c (1, nrow(Medium[[i]])), E = c(1:10))
}
#Picking up the optimal embedding dimension (E*)
Medium_BestE_A <- vector()
for(i in 1:length(Simplex_Medium_A)){
Medium_BestE_A[i] <- which.max(Simplex_Medium_A[[i]]$rho)
}
Medium_BestE_A
[1] 8 5 1 5 8 9 5 3 1 8 1 6 2 10 2 2
#Combining it with replicate ID
Medium_ID_BestE_A <- NULL
for(i in 1:length(Medium)){
Medium_ID_BestE_A <- rbind.data.frame(Medium_ID_BestE_A, cbind(Medium[[i]][1, ]$ID, Medium_BestE_A[i]))
}
#adding column names
colnames(Medium_ID_BestE_A) <- c("ID","C.c_E")
Medium_ID_BestE_A
Simplex_Large_Y <- list() #list for assigning results
#Executing simplex()
for(i in 1:length(Large)){
Simplex_Large_Y[[i]] <- simplex(Large[[i]]$YotsumonTotalL, lib = c(1, nrow(Large[[i]])), pred = c (1, nrow(Large[[i]])), E = c(1:10))
}
#Picking up the optimal embedding dimension (E*)
Large_BestE_Y <- vector()
for(i in 1:length(Simplex_Large_Y)){
Large_BestE_Y[i] <- which.max(Simplex_Large_Y[[i]]$rho)
}
Large_BestE_Y
[1] 8 5 4 9 1 7 6 7 5 4 4 8 8 4 4 4
#Combining it with replicate ID
Large_ID_BestE_Y <- NULL
for(i in 1:length(Large)){
Large_ID_BestE_Y <- rbind.data.frame(Large_ID_BestE_Y, cbind(Large[[i]][1, ]$ID, Large_BestE_Y[i]))
}
#adding column names
colnames(Large_ID_BestE_Y) <- c("ID","C.m_E")
Large_ID_BestE_Y
Simplex_Large_A <- list() #list for assigning results
#Executing simplex()
for(i in 1:length(Large)){
Simplex_Large_A[[i]] <- simplex(Large[[i]]$AzukiTotalL, lib = c(1, nrow(Large[[i]])), pred = c (1, nrow(Large[[i]])), E = c(1:10))
}
#Picking up the optimal embedding dimension (E*)
Large_BestE_A <- vector()
for(i in 1:length(Simplex_Large_A)){
Large_BestE_A[i] <- which.max(Simplex_Large_A[[i]]$rho)
}
Large_BestE_A
[1] 2 2 5 8 6 10 9 6 4 3 5 3 2 9 3 2
#Combining it with replicate ID
Large_ID_BestE_A <- NULL
for(i in 1:length(Large)){
Large_ID_BestE_A <- rbind.data.frame(Large_ID_BestE_A, cbind(Large[[i]][1, ]$ID, Large_BestE_A[i]))
}
#adding column names
colnames(Large_ID_BestE_A) <- c("ID","C.c_E")
Large_ID_BestE_A
Simplex_Huge_Y <- list() #list for assigning results
#Executing simplex()
for(i in 1:length(Huge)){
Simplex_Huge_Y[[i]] <- simplex(Huge[[i]]$YotsumonTotalL, lib = c(1, nrow(Huge[[i]])), pred = c (1, nrow(Huge[[i]])), E = c(1:10))
}
#Picking up the optimal embedding dimension (E*)
Huge_BestE_Y <- vector()
for(i in 1:length(Simplex_Huge_Y)){
Huge_BestE_Y[i] <- which.max(Simplex_Huge_Y[[i]]$rho)
}
Huge_BestE_Y
[1] 6 1 8 8 10 4 8 10 5 10 6 10 4 4 4 4
#Combining it with replicate ID
Huge_ID_BestE_Y <- NULL
for(i in 1:length(Huge)){
Huge_ID_BestE_Y <- rbind.data.frame(Huge_ID_BestE_Y, cbind(Huge[[i]][1, ]$ID, Huge_BestE_Y[i]))
}
colnames(Huge_ID_BestE_Y) <- c("ID","C.m_E")
Huge_ID_BestE_Y
Simplex_Huge_A <- list() #list for assigning results
#Executing simplex()
for(i in 1:length(Huge)){
Simplex_Huge_A[[i]] <- simplex(Huge[[i]]$AzukiTotalL, lib = c(1, nrow(Huge[[i]])), pred = c (1, nrow(Huge[[i]])), E = c(1:10))
}
#Picking up the optimal embedding dimension (E*)
Huge_BestE_A <- vector()
for(i in 1:length(Simplex_Huge_A)){
Huge_BestE_A
[1] 7 5 1 8 1 6 10 9 7 3 5 3 5 5 2 3
#Combining it with replicate ID
Huge_ID_BestE_A <- NULL
for(i in 1:length(Huge)){
Huge_ID_BestE_A <- rbind.data.frame(Huge_ID_BestE_A, cbind(Huge[[i]][1, ]$ID, Huge_BestE_A[i]))
}
#adding column names
colnames(Huge_ID_BestE_A) <- c("ID","C.c_E")
Huge_ID_BestE_A
#For Small
Small_ID_BestE <- merge(Small_ID_BestE_Y, Small_ID_BestE_A, by = "ID")
Small_ID_BestE
#For Medium
Medium_ID_BestE <- merge(Medium_ID_BestE_Y, Medium_ID_BestE_A, by = "ID")
Medium_ID_BestE
#For Large
Large_ID_BestE <- merge(Large_ID_BestE_Y, Large_ID_BestE_A, by = "ID")
Large_ID_BestE
#For Huge
Huge_ID_BestE <- merge(Huge_ID_BestE_Y, Huge_ID_BestE_A, by = "ID")
Huge_ID_BestE
NA
com_Small_E <- apply(Small_ID_BestE[, -1], 1, max)
com_Medium_E <- apply(Medium_ID_BestE[, -1], 1, max)
com_Large_E <- apply(Large_ID_BestE[, -1], 1, max)
com_Huge_E <- apply(Huge_ID_BestE[, -1], 1, max)
com_Small_E
[1] 9 6 4 2 2 8 8 5 3 5 5 4 1 4 1 3
com_Medium_E
[1] 8 5 1 5 10 9 5 4 5 8 1 6 4 10 4 2
com_Large_E
[1] 8 5 5 9 6 10 9 7 5 4 5 8 8 9 4 4
com_Huge_E
[1] 7 5 8 8 10 6 10 10 7 10 6 10 5 5 4 4
#Merging the optimal E and the number of NA
Small_ID_BestE_NA <- data.frame(Small_ID_BestE, No.NA = com_Small_E)
Medium_ID_BestE_NA <- data.frame(Medium_ID_BestE, No.NA = com_Medium_E)
Large_ID_BestE_NA <- data.frame(Large_ID_BestE, No.NA = com_Large_E)
Huge_ID_BestE_NA <- data.frame(Huge_ID_BestE, No.NA = com_Huge_E)
#example
Huge_ID_BestE_NA
#Merging it with the coexistence period
SmallTable1base <- merge(subset(TableS1,TableS1$Cage.Size == "Small", select = c(ID, Coexistence_periods)), Small_ID_BestE_NA, by = "ID")
MediumTable1base <- merge(subset(TableS1,TableS1$Cage.Size == "Medium", select = c(ID, Coexistence_periods)), Medium_ID_BestE_NA, by = "ID")
LargeTable1base <- merge(subset(TableS1,TableS1$Cage.Size == "Large", select = c(ID, Coexistence_periods)), Large_ID_BestE_NA, by = "ID")
HugeTable1base <- merge(subset(TableS1,TableS1$Cage.Size == "Huge", select = c(ID, Coexistence_periods)), Huge_ID_BestE_NA, by = "ID")
#example
LargeTable1base
– Noting that the length of each replicate is equal to the coexistence period + 1
– Therefore, 1) when the number of coexistence weeks should be greater than or equal to “34”, the combining with other replicates is not needed, but 2) for the number of the sum of the coexistence period for the combined replicate should be larger than or equal to “33” weeks
#Descending order & labeling replicates
SmallTable1base_m <- SmallTable1base[order(SmallTable1base$Coexistence_periods,decreasing = TRUE),]
SmallTable1base_m
ID_combined <- c(rep("S1", 2), rep("S2", 2), rep("S3", 2), rep("S4", 3), rep("S5", 3), rep("S6", 4))
SmallTable1base_m <- cbind.data.frame(ID_combined, SmallTable1base_m)
SmallTable1base_m
#Descending order & labeling replicates
MediumTable1base_m <- MediumTable1base[order(MediumTable1base$Coexistence_periods,decreasing = TRUE),]
MediumTable1base_m
ID_combined <- c(rep("M1", 1), rep("M2", 2), rep("M3", 2), rep("M4", 2), rep("M5", 2), rep("M6", 3), rep("M7", 4))
MediumTable1base_m <- cbind.data.frame(ID_combined, MediumTable1base_m)
MediumTable1base_m
#Descending order & labeling replicates
LargeTable1base_m <- LargeTable1base[order(LargeTable1base$Coexistence_periods,decreasing = TRUE),]
LargeTable1base_m
ID_combined <- c(rep("L1", 1), rep("L2", 1), rep("L3", 2), rep("L4", 2), rep("L5", 2), rep("L6", 2), rep("L7", 2), rep("L8", 4))
LargeTable1base_m <- cbind.data.frame(ID_combined, LargeTable1base_m)
LargeTable1base_m
#Descending order & labeling replicates
HugeTable1base_m <- HugeTable1base[order(HugeTable1base$Coexistence_periods,decreasing = TRUE),]
HugeTable1base_m
ID_combined <- c(rep("H1", 1), rep("H2", 2), rep("H3", 2), rep("H4", 2), rep("H5", 2), rep("H6", 2), rep("H7", 2), rep("H8", 3))
HugeTable1base_m <- cbind.data.frame(ID_combined, HugeTable1base_m)
HugeTable1base_m
S_combined_Y <- list()
S_combined_A <- list()
S_id_comb <- c("S1", "S2", "S3", "S4", "S5", "S6")
#S1-S3: combining two replicates
for(i in 1:3) {
#determining the number of NA to be added
no_NA_inserted <- max(subset(SmallTable1base_m, SmallTable1base_m$ID_combined == S_id_comb[i], select = No.NA))
#listing the replicate IDs used for i-th combined replicates
id_each <- subset(SmallTable1base_m, SmallTable1base_m$ID_combined == S_id_comb[i], select = ID)
S_combined_Y[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$YotsumonTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$YotsumonTotalL)
S_combined_A[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$AzukiTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$AzukiTotalL)
}
#S4-S5: combining three replicates
for(i in 4:5) {
no_NA_inserted <- max(subset(SmallTable1base_m, SmallTable1base_m$ID_combined == S_id_comb[i], select = No.NA))
id_each <- subset(SmallTable1base_m, SmallTable1base_m$ID_combined == S_id_comb[i], select = ID)
S_combined_Y[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$YotsumonTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$YotsumonTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[3,],]$YotsumonTotalL)
S_combined_A[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$AzukiTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$AzukiTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[3,],]$AzukiTotalL)
}
#S6: combining four replicates
no_NA_inserted <- max(subset(SmallTable1base_m, SmallTable1base_m$ID_combined == S_id_comb[i], select = No.NA))
id_each <- subset(SmallTable1base_m, SmallTable1base_m$ID_combined == S_id_comb[i], select = ID)
S_combined_Y[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$YotsumonTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$YotsumonTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[3,],]$YotsumonTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[4,],]$YotsumonTotalL)
S_combined_A[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$AzukiTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$AzukiTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[3,],]$AzukiTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[4,],]$AzukiTotalL)
#example
S_combined_Y
[[1]]
[1] 8 0 0 0 94 16 1 36 180 97 2 138 169 135 6 118 129 53
[19] 10 94 96 54 0 56 78 100 0 53 NA NA NA NA NA NA NA NA
[37] NA 8 0 0 0 110 4 12 37 192 21 11 133 55 9 98 197 76
[55] 30 150 141 88 100 91 119 76 130
[[2]]
[1] 8 0 0 0 110 39 1 41 233 58 3 141 179 82 8 112 106 68
[19] 24 70 120 114 48 100 NA NA NA 4
[37] 114 21 0 75 200 14 38 165 106 24 141 180 141 83 143 195 124 91
[[3]]
[1] 8 0 0 0 98 23 0 35 208 47 2 132 158 128 23 128 112 96
[19] 61 NA NA NA NA NA 8 0 0 0 111 62 5 51 246 62 4 158
[37] 177 105 72 133
[[4]]
[1] 8 0 0 16 117 0 0 171 62 18 125 197 139 113 NA NA NA NA
[19] NA NA 8 0 0 0 98 14 0 119 195 5 116 175 192 NA NA NA
[37] NA NA NA 8 0 0 0 139 33 1 95 214 128 34 199 138
[[5]]
[1] 8 0 0 0 135 24 0 53 215 87 20 402 NA NA NA NA 8 0
[19] 0 0 107 15 0 59 173 0 3 112 NA NA NA NA 8 0 0 2
[37] 101 10 0 67 151 6 8 246
[[6]]
[1] 8 0 0 0 88 21 3 109 213 11 NA NA NA 8 0 0 0 146
[19] 14 0 72 208 59 NA NA NA 8 0 0 0 114 16 1 123 179 12
[37] NA NA NA 8 0 0 0 139 27 0 129 238
M_combined_Y <- list()
M_combined_A <- list()
M_id_comb <- c("M1", "M2", "M3", "M4", "M5", "M6", "M7")
#M1: no need for combining
i <- 1
id_each <- subset(MediumTable1base_m, MediumTable1base_m$ID_combined == M_id_comb[i], select = ID)
M_combined_Y[[i]] <- shavedataframe[shavedataframe$ID == id_each[1,],]$YotsumonTotalL
M_combined_A[[i]] <- shavedataframe[shavedataframe$ID == id_each[1,],]$AzukiTotalL
#M2-M5: combining two replicates
for(i in 2:5) {
no_NA_inserted <- max(subset(MediumTable1base_m, MediumTable1base_m$ID_combined == M_id_comb[i], select = No.NA))
id_each <- subset(MediumTable1base_m, MediumTable1base_m$ID_combined == M_id_comb[i], select = ID)
M_combined_Y[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$YotsumonTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$YotsumonTotalL)
M_combined_A[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$AzukiTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$AzukiTotalL)
}
#M6: combining three replicates
for(i in 6:6) {
no_NA_inserted <- max(subset(MediumTable1base_m, MediumTable1base_m$ID_combined == M_id_comb[i], select = No.NA))
id_each <- subset(MediumTable1base_m, MediumTable1base_m$ID_combined == M_id_comb[i], select = ID)
M_combined_Y[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$YotsumonTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$YotsumonTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[3,],]$YotsumonTotalL)
M_combined_A[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$AzukiTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$AzukiTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[3,],]$AzukiTotalL)
}
#M7: combining four replicates
i <- 7
no_NA_inserted <- max(subset(MediumTable1base_m, MediumTable1base_m$ID_combined == M_id_comb[i], select = No.NA))
id_each <- subset(MediumTable1base_m, MediumTable1base_m$ID_combined == M_id_comb[i], select = ID)
M_combined_Y[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$YotsumonTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$YotsumonTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[3,],]$YotsumonTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[4,],]$YotsumonTotalL)
M_combined_A[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$AzukiTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$AzukiTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[3,],]$AzukiTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[4,],]$AzukiTotalL)
#example
M_combined_A
[[1]]
[1] 8 0 0 0 68 7 0 47 11 0 12 43 0 0 100 2 0 0
[19] 60 5 1 0 131 2 0 0 130 6 0 0 0 20 0 0 0 46
[37] 0 0 16 47 3 1 52 0 0 1 5
[[2]]
[1] 8 0 0 1 93 2 0 219 32 0 124 218 1 24 195 12 0 15
[19] 100 1 0 34 109 4 0 4 81 5 0 1 1 0 NA NA NA NA
[37] NA NA NA NA NA NA 8 0 0 3 109 4 0 231 42 0 176 150
[55] 1 22 201 65 0 6 79 12 0 0 147 3 0 0 100 4 0
[[3]]
[1] 8 0 0 4 91 0 0 171 10 0 225 119 0 41 197 17 0 15
[19] 18 1 0 0 12 0 0 0 7 1 0 NA NA NA NA NA NA NA
[37] NA NA NA 8 0 0 1 117 1 0 158 9 0 44 115 1 1 50
[55] 54 0 0 90 14 0 0 198 25 0
[[4]]
[1] 8 0 0 67 24 2 163 208 19 44 156 10 12 179 21 1 25 48
[19] 0 0 18 0 NA NA NA NA NA NA NA NA 8 0 0 0 126 0
[37] 0 195 25 0 219 92 1 13 180 8 0 10 12 0
[[5]]
[1] 8 0 0 4 70 11 0 43 30 10 101 22 1 5 72 0 0 0
[19] 11 0 NA NA NA NA NA 8 0 0 22 50 0 5 7 1 2 8
[37] 0 1 12 2 0 3 0
[[6]]
[1] 8 0 0 3 83 2 3 35 9 0 65 10 0 0 22 0 NA NA
[19] NA NA NA NA 8 0 0 2 80 22 0 28 9 1 16 7 0 0
[37] 2 0 NA NA NA NA NA NA 8 0 0 3 60 2 2 102 3 0
[55] 22 2 0 0 10 0
[[7]]
[1] 8 0 0 17 63 3 5 97 12 1 100 4 0 27 0 NA NA NA
[19] NA NA 8 0 0 5 15 0 4 8 0 0 29 0 NA NA NA NA
[37] NA 8 0 0 15 40 0 0 7 5 0 NA NA NA NA NA 8 0
[55] 0 1 24 1 0 0 1 0
L_combined_Y <- list()
L_combined_A <- list()
L_id_comb <- c("L1", "L2", "L3", "L4", "L5", "L6", "L7", "L8")
#L1-L2: no need for combining
for(i in 1:2){
id_each <- subset(LargeTable1base_m, LargeTable1base_m$ID_combined == L_id_comb[i], select = ID)
L_combined_Y[[i]] <- shavedataframe[shavedataframe$ID == id_each[1,],]$YotsumonTotalL
L_combined_A[[i]] <- shavedataframe[shavedataframe$ID == id_each[1,],]$AzukiTotalL
}
#L3-L7: combining two replicates
for(i in 3:7) {
no_NA_inserted <- max(subset(LargeTable1base_m, LargeTable1base_m$ID_combined == L_id_comb[i], select = No.NA))
id_each <- subset(LargeTable1base_m, LargeTable1base_m$ID_combined == L_id_comb[i], select = ID)
L_combined_Y[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$YotsumonTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$YotsumonTotalL)
L_combined_A[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$AzukiTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$AzukiTotalL)
}
#L8: combining four replicates
i <- 8
no_NA_inserted <- max(subset(LargeTable1base_m, LargeTable1base_m$ID_combined == L_id_comb[i], select = No.NA))
id_each <- subset(LargeTable1base_m, LargeTable1base_m$ID_combined == L_id_comb[i], select = ID)
L_combined_Y[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$YotsumonTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$YotsumonTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[3,],]$YotsumonTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[4,],]$YotsumonTotalL)
L_combined_A[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$AzukiTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$AzukiTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[3,],]$AzukiTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[4,],]$AzukiTotalL)
#example
L_combined_Y
[[1]]
[1] 8 0 0 0 104 9 0 31 186 9 0 59 136 20 24 49 137 44
[19] 4 34 113 73 0 50 86 85 26 0 0 0 59 0 0 26 138 19
[37] 0 142 152 5 94 150 204
[[2]]
[1] 8 0 0 0 77 19 0 36 207 22 4 56 178 8 2 118 140 6
[19] 22 155 67 15 71 97 137 5 110 85 56 5 11 48 1 0 58 52
[37] 96
[[3]]
[1] 8 0 0 19 111 36 2 2 167 66 8 186 37 15 0 133 52 9
[19] 119 126 16 5 156 106 32 13 112 107 86 25 63 82 97 65 NA NA
[37] NA NA NA NA NA NA NA NA 8 0 0 0 121 74 5 9 162 162
[55] 10 4 208 161 5 10 96 56 1 0 82 98 1 0 63 77 1 0
[73] 58
[[4]]
[1] 8 0 0 0 102 40 0 8 178 59 1 14 172 105 0 21 121 67
[19] 1 14 61 91 2 10 37 71 4 6 16 NA NA NA NA NA NA NA
[37] NA NA 8 0 0 109 72 0 72 157 36 0 12 231 92 2 14
[55] 136 76 8 8 126 90 1 0 93 64 4 0
[[5]]
[1] 8 0 0 0 107 29 0 44 212 20 0 68 110 35 10 143 155 126
[19] 15 130 125 116 5 129 NA NA NA NA NA NA NA NA 8 0 0
[37] 0 114 41 0 36 243 148 27 172 133 86 7 154 123 101 16 151
[[6]]
[1] 8 0 0 0 132 50 1 60 146 128 25 201 56 152 6 106 16 129
[19] 10 112 NA NA NA NA NA NA NA NA 8 0 0 0 142 74 1 31
[37] 223 134 19 119 197 185 13 157 81 86 28 92
[[7]]
[1] 8 0 0 10 162 31 7 126 176 53 102 163 185 164 113 144 247 148
[19] 140 NA NA NA NA NA 8 0 0 124 16 0 140 234 33 121 190
[37] 251 72 201 154 242 127 185
[[8]]
[1] 8 0 0 0 102 21 1 139 206 45 42 241 245 173 119 14 NA NA
[19] NA NA NA NA 8 0 0 0 123 50 3 28 231 128 10 133 192 147
[37] 11 165 NA NA NA NA NA NA 8 0 0 0 132 44 1 10 234 120
[55] 5 141 110 163 9 135 NA NA NA NA NA NA 8 0 0 0 137 27
[73] 0 116 201 25 0 155 141 69 57 106
H_combined_Y <- list()
H_combined_A <- list()
H_id_comb <- c("H1", "H2", "H3", "H4", "H5", "H6", "H7", "H8")
#H1: no need for combining
i <- 1
id_each <- subset(HugeTable1base_m, HugeTable1base_m$ID_combined == H_id_comb[i], select = ID)
H_combined_Y[[i]] <- shavedataframe[shavedataframe$ID == id_each[1,],]$YotsumonTotalL
H_combined_A[[i]] <- shavedataframe[shavedataframe$ID == id_each[1,],]$AzukiTotalL
#H2-H7: combining two replicates
for(i in 2:7) {
no_NA_inserted <- max(subset(HugeTable1base_m, HugeTable1base_m$ID_combined == H_id_comb[i], select = No.NA))
id_each <- subset(HugeTable1base_m, HugeTable1base_m$ID_combined == H_id_comb[i], select = ID)
H_combined_Y[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$YotsumonTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$YotsumonTotalL)
H_combined_A[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$AzukiTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$AzukiTotalL)
}
#H8: combining three replicates
for(i in 8:8) {
no_NA_inserted <- max(subset(HugeTable1base_m, HugeTable1base_m$ID_combined == H_id_comb[i], select = No.NA))
id_each <- subset(HugeTable1base_m, HugeTable1base_m$ID_combined == H_id_comb[i], select = ID)
H_combined_Y[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$YotsumonTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$YotsumonTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[3,],]$YotsumonTotalL)
H_combined_A[[i]] <- c(shavedataframe[shavedataframe$ID == id_each[1,],]$AzukiTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[2,],]$AzukiTotalL, rep(NA,no_NA_inserted), shavedataframe[shavedataframe$ID == id_each[3,],]$AzukiTotalL)
}
#example
H_combined_Y
[[1]]
[1] 8 0 0 1 127 24 1 4 173 59 7 94 189 61 2 157 76 15
[19] 17 165 37 0 40 130 44 0 56 125 70 2 65 86 98 3 69
[[2]]
[1] 8 0 0 0 98 54 0 9 197 115 6 27 216 187 0 78 133 113
[19] 9 99 99 111 0 74 92 63 4 74 101 41 5 2 0 NA NA NA
[37] NA NA NA NA NA NA NA 8 0 0 0 100 18 0 17 165 56 0
[55] 6 187 188 10 0 103 86 0 0 4 125 0 0 5 5 2 0 0
[73] 1 21 3
[[3]]
[1] 8 0 0 0 88 18 0 59 33 0 0 29 3 0 0 22 10 0
[19] 0 11 8 0 0 5 4 0 0 1 2 1 0 NA NA NA NA NA
[37] NA NA NA NA NA 8 0 0 0 93 61 1 19 216 103 11 7 188
[55] 15 1 3 121 23 0 3 85 6 0 0 23 6 0 0 2 1 0
[[4]]
[1] 8 0 0 0 68 44 0 6 97 22 3 1 258 129 0 3 129 98
[19] 1 0 106 73 2 0 76 66 3 0 79 NA NA NA NA NA NA NA
[37] NA NA NA 8 0 0 0 92 57 3 93 111 157 14 0 175 182 4
[55] 0 88 35 1 0 91 77 9 0 48 54 2 0 31
[[5]]
[1] 8 0 0 0 98 28 1 16 139 82 3 13 228 200 7 6 87 52
[19] 7 0 63 85 5 0 167 147 90 1 NA NA NA NA NA NA NA NA
[37] NA NA 8 0 0 0 135 75 2 15 152 202 7 83 227 104 0 144
[55] 90 89 8 107 124 108 3 53 31 52 4 40
[[6]]
[1] 8 0 0 0 103 38 0 39 225 117 6 177 212 165 7 139 118 109
[19] 30 42 125 117 8 31 72 40 3 11 NA NA NA NA NA NA NA NA
[37] NA NA 8 0 0 0 94 49 2 6 146 30 0 1 154 68 1 0
[55] 39 97 5 0 0 112 7 0
[[7]]
[1] 8 0 0 0 73 34 0 7 152 25 0 24 137 37 4 51 114 57
[19] 0 102 132 76 1 138 NA NA NA NA NA NA NA NA 8 0 0 0
[37] 111 57 0 52 204 157 11 173 171 156 39 129 117 111 116 82
[[8]]
[1] 8 0 0 0 77 12 0 20 199 112 14 64 121 141 9 69 108 NA
[19] NA NA NA NA 8 0 0 0 147 56 2 148 260 102 54 146 130 134
[37] 200 247 NA NA NA NA NA 8 0 0 0 113 65 1 64 233 109 12
[55] 140 124 147 42 143
S_combined_mean_coweek <- tapply(SmallTable1base_m$Coexistence_periods, SmallTable1base_m$ID_combined, mean)
M_combined_mean_coweek <- tapply(MediumTable1base_m$Coexistence_periods, MediumTable1base_m$ID_combined, mean)
L_combined_mean_coweek <- tapply(LargeTable1base_m$Coexistence_periods, LargeTable1base_m$ID_combined, mean)
H_combined_mean_coweek <- tapply(HugeTable1base_m$Coexistence_periods, HugeTable1base_m$ID_combined, mean)
S_combined_mean_coweek
S1 S2 S3 S4 S5 S6
26.00000 22.00000 16.50000 12.33333 11.00000 8.75000
M_combined_mean_coweek
M1 M2 M3 M4 M5 M6 M7
46.00 29.50 26.00 20.00 18.00 15.00 10.75
L_combined_mean_coweek
L1 L2 L3 L4 L5 L6 L7 L8
42.0 36.0 30.5 27.5 21.0 19.0 18.0 15.0
H_combined_mean_coweek
H1 H2 H3 H4 H5 H6 H7 H8
34.00000 31.50000 30.00000 28.00000 27.00000 25.00000 21.00000 15.33333
Small_combined_Y_simplex <- list()
Small_combined_Y_smap <- list()
length(S_id_comb)
[1] 6
#Executing simplex projection
for(i in 1:length(S_id_comb)){
Small_combined_Y_simplex[[i]] <- simplex(S_combined_Y[[i]],lib = c(1,length(S_combined_Y[[i]])), pred = c(1,length(S_combined_Y[[i]])), E = c(1:10), silent = TRUE)
}
#Calculating optimal E
Small_combined_E_Y <- vector()
for(i in 1:length(S_id_comb)){
Small_combined_E_Y[i] <- which.max(Small_combined_Y_simplex[[i]]$rho)
}
#Executing univariate S-map
for(i in 1:length(S_id_comb)){
Small_combined_Y_smap[[i]] <- s_map(S_combined_Y[[i]],lib = c(1,length(S_combined_Y[[i]])), pred = c(1,length(S_combined_Y[[i]])),E = Small_combined_E_Y[i], silent = TRUE)
}
#Calculating optimal theta
Small_combined_theta_Y <- vector()
for(i in 1:length(S_id_comb)){
Small_combined_theta_Y[i] <- Small_combined_Y_smap[[i]][which.max(Small_combined_Y_smap[[i]]$rho),]$theta
}
#Calculating delta-rho, maximum minus minimum
Small_combined_deltarho_Y <- vector()
for(i in 1:length(S_id_comb)){
Small_combined_deltarho_Y[i] <- max(Small_combined_Y_smap[[i]]$rho) - Small_combined_Y_smap[[i]]$rho[1]
}
Small_combined_E_Y
[1] 5 4 10 5 3 8
Small_combined_theta_Y
[1] 4.0 2.0 0.0 2.0 8.0 0.1
Small_combined_deltarho_Y
[1] 0.1077679 0.1088690 0.0000000 0.1371399 0.5610443 0.1443142
Small_combined_A_simplex <- list()
Small_combined_A_smap <- list()
#Executing simplex projection
for(i in 1:length(S_id_comb)){
Small_combined_A_simplex[[i]] <- simplex(S_combined_A[[i]],lib = c(1,length(S_combined_A[[i]])), pred = c(1,length(S_combined_A[[i]])), E = c(1:10), silent = TRUE)
}
#Calculating optimal E
Small_combined_E_A <- vector()
Small_combined_E_A[i] <- which.max(Small_combined_A_simplex[[i]]$rho)
}
#Executing univariate S-map
for(i in 1:length(S_id_comb)){
Small_combined_A_smap[[i]] <- s_map(S_combined_A[[i]],lib = c(1,length(S_combined_A[[i]])), pred = c(1,length(S_combined_A[[i]])), E = Small_combined_E_A[i], silent = TRUE)
}
#Calculating optimal theta
Small_combined_theta_A <- vector()
for(i in 1:length(S_id_comb)){
#Calculating delta-rho
Small_combined_deltarho_A <- vector()
for(i in 1:length(S_id_comb)){
Small_combined_deltarho_A[i] <- max(Small_combined_A_smap[[i]]$rho) - Small_combined_A_smap[[i]]$rho[1]
}
Small_combined_E_A
[1] 10 3 3 7 4 4
Small_combined_theta_A
[1] 1 8 0 0 3 3
Small_combined_deltarho_A
[1] 0.01861431 0.35407409 0.00000000 0.00000000 0.71102585 0.42537821
Medium_combined_Y_simplex <- list()
Medium_combined_Y_smap <- list()
length(M_id_comb)
[1] 7
#Executing simplex projection
for(i in 1:length(M_id_comb)){
Medium_combined_Y_simplex[[i]] <- simplex(M_combined_Y[[i]],lib = c(1,length(M_combined_Y[[i]])), pred = c(1,length(M_combined_Y[[i]])), E = c(1:10), silent = TRUE)
}
#Calculating optimal E
Medium_combined_E_Y <- vector()
for(i in 1:length(M_id_comb)){
Medium_combined_E_Y[i] <- which.max(Medium_combined_Y_simplex[[i]]$rho)
}
#Executing univariate S-map
for(i in 1:length(M_id_comb)){
Medium_combined_Y_smap[[i]] <- s_map(M_combined_Y[[i]],lib = c(1,length(M_combined_Y[[i]])), pred = c(1,length(M_combined_Y[[i]])), E = Medium_combined_E_Y[i], silent = TRUE)
}
#Calculating optimal theta
Medium_combined_theta_Y <- vector()
for(i in 1:length(M_id_comb)){
Medium_combined_theta_Y[i] <- Medium_combined_Y_smap[[i]][which.max(Medium_combined_Y_smap[[i]]$rho),]$theta
}
#Calculating delta-rho
Medium_combined_deltarho_Y <- vector()
for(i in 1:length(M_id_comb)){
Medium_combined_deltarho_Y[i] <- max(Medium_combined_Y_smap[[i]]$rho) - Medium_combined_Y_smap[[i]]$rho[1]
}
Medium_combined_E_Y
[1] 7 8 9 5 8 6 8
Medium_combined_theta_Y
[1] 3.00 2.00 1.50 2.00 0.75 3.00 1.00
Medium_combined_deltarho_Y
[1] 0.09693470 0.17438007 0.05197754 0.06097689 0.02553795 0.15239394
[7] 0.01075863
Medium_combined_A_simplex <- list()
Medium_combined_A_smap <- list()
#Executing simplex projection
for(i in 1:length(M_id_comb)){
Medium_combined_A_simplex[[i]] <- simplex(M_combined_A[[i]],lib = c(1,length(M_combined_A[[i]])), pred = c(1,length(M_combined_A[[i]])),E = c(1:10), silent = TRUE)
}
#Calculating optimal E
Medium_combined_E_A <- vector()
for(i in 1:length(M_id_comb)){
Medium_combined_E_A[i] <- which.max(Medium_combined_A_simplex[[i]]$rho)
}
#Executing univariate S-map
for(i in 1:length(M_id_comb)){
Medium_combined_A_smap[[i]] <- s_map(M_combined_A[[i]],lib = c(1,length(M_combined_A[[i]])), pred = c(1,length(M_combined_A[[i]])),E = Medium_combined_E_A[i], silent = TRUE)
}
#Calculating optimal theta
Medium_combined_theta_A <- vector()
Medium_combined_theta_A[i] <- Medium_combined_A_smap[[i]][which.max(Medium_combined_A_smap[[i]]$rho),]$theta
}
#Calculating delta-rho
Medium_combined_deltarho_A <- vector()
for(i in 1:length(M_id_comb)){
Medium_combined_deltarho_A[i] <- max(Medium_combined_A_smap[[i]]$rho) - Medium_combined_A_smap[[i]]$rho[1]
}
Medium_combined_E_A
[1] 8 8 7 10 5 5 5
Medium_combined_theta_A
[1] 0.00 3.00 4.00 0.75 3.00 3.00 8.00
Medium_combined_deltarho_A
[1] 0.00000000 0.14654892 0.06478939 0.02017538 0.18022398 0.08901785
[7] 0.21635228
Large_combined_Y_simplex <- list()
Large_combined_Y_smap <- list()
length(L_id_comb)
[1] 8
#Executing simplex projection
for(i in 1:length(L_id_comb)){
Large_combined_Y_simplex[[i]] <- simplex(L_combined_Y[[i]],lib = c(1,length(L_combined_Y[[i]])), pred = c(1,length(L_combined_Y[[i]])), E = c(1:10), silent = TRUE)
}
#Calculating optimal E
Large_combined_E_Y <- vector()
for(i in 1:length(L_id_comb)){
Large_combined_E_Y[i] <- which.max(Large_combined_Y_simplex[[i]]$rho)
}
#Executing univariate S-map
for(i in 1:length(L_id_comb)){
Large_combined_Y_smap[[i]] <- s_map(L_combined_Y[[i]],lib = c(1,length(L_combined_Y[[i]])), pred = c(1,length(L_combined_Y[[i]])), E = Large_combined_E_Y[i], silent = TRUE)
}
#Calculating optimal theta
Large_combined_theta_Y <- vector()
for(i in 1:length(L_id_comb)){
Large_combined_theta_Y[i] <- Large_combined_Y_smap[[i]][which.max(Large_combined_Y_smap[[i]]$rho),]$theta
}
#Calculating delta-rho
Large_combined_deltarho_Y <- vector()
for(i in 1:length(L_id_comb)){
Large_combined_deltarho_Y[i] <- max(Large_combined_Y_smap[[i]]$rho) - Large_combined_Y_smap[[i]]$rho[1]
}
Large_combined_E_Y
[1] 8 9 4 9 4 5 5 5
Large_combined_theta_Y
[1] 3.0 4.0 0.0 6.0 2.0 2.0 1.5 2.0
Large_combined_deltarho_Y
[1] 0.16517993 0.10936934 0.00000000 0.07956119 0.10937309 0.24821462
[7] 0.06062738 0.11924754
Large_combined_A_simplex <- list()
Large_combined_A_smap <- list()
#Executing simplex projection
for(i in 1:length(L_id_comb)){
Large_combined_A_simplex[[i]] <- simplex(L_combined_A[[i]],lib = c(1,length(L_combined_A[[i]])), pred = c(1,length(L_combined_A[[i]])),E = c(1:10), silent = TRUE)
}
#Calculating optimal E
Large_combined_E_A <- vector()
for(i in 1:length(L_id_comb)){
Large_combined_E_A[i] <- which.max(Large_combined_A_simplex[[i]]$rho)
}
#Executing univariate S-map
for(i in 1:length(L_id_comb)){
Large_combined_A_smap[[i]] <- s_map(L_combined_A[[i]],lib = c(1,length(L_combined_A[[i]])), pred = c(1,length(L_combined_A[[i]])),E = Large_combined_E_A[i], silent = TRUE)
}
#Calculating optimal theta
Large_combined_theta_A <- vector()
for(i in 1:length(L_id_comb)){
Large_combined_theta_A[i] <- Large_combined_A_smap[[i]][which.max(Large_combined_A_smap[[i]]$rho),]$theta
}
#Calculating delta-rho
Large_combined_deltarho_A <- vector()
for(i in 1:length(L_id_comb)){
Large_combined_deltarho_A[i] <- max(Large_combined_A_smap[[i]]$rho) - Large_combined_A_smap[[i]]$rho[1]
}
Large_combined_E_A
[1] 2 8 10 10 9 9 1 5
Large_combined_theta_A
[1] 6.0 0.0 0.0 1.5 8.0 1.5 8.0 4.0
Large_combined_deltarho_A
[1] 0.16542098 0.00000000 0.00000000 0.14433095 0.03180342 0.11819601
[7] 0.20139964 0.37881962
Huge_combined_Y_simplex <- list()
Huge_combined_Y_smap <- list()
length(H_id_comb)
[1] 8
#Executing simplex projection
for(i in 1:length(H_id_comb)){
Huge_combined_Y_simplex[[i]] <- simplex(H_combined_Y[[i]],lib = c(1,length(H_combined_Y[[i]])), pred = c(1,length(H_combined_Y[[i]])), E = c(1:10), silent = TRUE)
}
#Calculating optimal E
Huge_combined_E_Y <- vector()
for(i in 1:length(H_id_comb)){
Huge_combined_E_Y[i] <- which.max(Huge_combined_Y_simplex[[i]]$rho)
}
#Executing univariate S-map
for(i in 1:length(H_id_comb)){
Huge_combined_Y_smap[[i]] <- s_map(H_combined_Y[[i]],lib = c(1,length(H_combined_Y[[i]])), pred = c(1,length(H_combined_Y[[i]])), E = Huge_combined_E_Y[i], silent = TRUE)
}
#Calculating optimal theta
Huge_combined_theta_Y <- vector()
for(i in 1:length(H_id_comb)){
Huge_combined_theta_Y[i] <- Huge_combined_Y_smap[[i]][which.max(Huge_combined_Y_smap[[i]]$rho),]$theta
}
#Calculating delta-rho
Huge_combined_deltarho_Y <- vector()
for(i in 1:length(H_id_comb)){
Huge_combined_deltarho_Y[i] <- max(Huge_combined_Y_smap[[i]]$rho) - Huge_combined_Y_smap[[i]]$rho[1]
}
Huge_combined_E_Y
[1] 6 10 10 6 5 5 5 4
Huge_combined_theta_Y
[1] 6 2 6 0 2 1 1 2
Huge_combined_deltarho_Y
[1] 0.054355365 0.136859678 0.140825599 0.000000000 0.094501451 0.004716782
[7] 0.027814963 0.143490533
Huge_combined_A_simplex <- list()
Huge_combined_A_smap <- list()
#Executing simplex projection
for(i in 1:length(H_id_comb)){
Huge_combined_A_simplex[[i]] <- simplex(H_combined_A[[i]],lib = c(1,length(H_combined_A[[i]])), pred = c(1,length(H_combined_A[[i]])),E = c(1:10), silent = TRUE)
}
#Calculating optimal E
Huge_combined_E_A <- vector()
for(i in 1:length(H_id_comb)){
Huge_combined_E_A[i] <- which.max(Huge_combined_A_simplex[[i]]$rho)
}
#Executing univariate S-map
for(i in 1:length(H_id_comb)){
Huge_combined_A_smap[[i]] <- s_map(H_combined_A[[i]],lib = c(1,length(H_combined_A[[i]])), pred = c(1,length(H_combined_A[[i]])),E = Huge_combined_E_A[i], silent = TRUE)
}
#Calculating optimal theta
Huge_combined_theta_A <- vector()
for(i in 1:length(H_id_comb)){
Huge_combined_theta_A[i] <- Huge_combined_A_smap[[i]][which.max(Huge_combined_A_smap[[i]]$rho),]$theta
}
#Calculating delta-rho
Huge_combined_deltarho_A <- vector()
for(i in 1:length(H_id_comb)){
Huge_combined_deltarho_A[i] <- max(Huge_combined_A_smap[[i]]$rho) - Huge_combined_A_smap[[i]]$rho[1]
}
Huge_combined_E_A
[1] 7 4 9 6 8 6 9 5
Huge_combined_theta_A
[1] 1.5 8.0 3.0 3.0 0.0 1.0 6.0 2.0
Huge_combined_deltarho_A
[1] 0.04300393 0.02070305 0.02979581 0.15838825 0.00000000 0.01877942
[7] 0.08699150 0.03294874
Table1Small <- data.frame("ID" = S_id_comb, "average_coexistence_week" = S_combined_mean_coweek, "E_for_C.m" = Small_combined_E_Y, "E_for_C.c" = Small_combined_E_A, "theta_for_C.m" = Small_combined_theta_Y, "theta_for_C.c" = Small_combined_theta_A, "delta-rho_for_C.m" = Small_combined_deltarho_Y, "delta-rho_for_C.c" = Small_combined_deltarho_A)
Table1Medium <- data.frame("ID" = M_id_comb, "average_coexistence_week" = M_combined_mean_coweek, "E_for_C.m" = Medium_combined_E_Y, "E_for_C.c" = Medium_combined_E_A, "theta_for_C.m" = Medium_combined_theta_Y, "theta_for_C.c" = Medium_combined_theta_A, "delta-rho_for_C.m" = Medium_combined_deltarho_Y, "delta-rho_for_C.c" = Medium_combined_deltarho_A)
Table1Large <- data.frame("ID" = L_id_comb, "average_coexistence_week" = L_combined_mean_coweek, "E_for_C.m" = Large_combined_E_Y, "E_for_C.c" = Large_combined_E_A, "theta_for_C.m" = Large_combined_theta_Y, "theta_for_C.c" = Large_combined_theta_A, "delta-rho_for_C.m" = Large_combined_deltarho_Y, "delta-rho_for_C.c" = Large_combined_deltarho_A)
Table1Huge <- data.frame("ID" = H_id_comb, "average_coexistence_week" = H_combined_mean_coweek, "E_for_C.m" = Huge_combined_E_Y, "E_for_C.c" = Huge_combined_E_A, "theta_for_C.m" = Huge_combined_theta_Y, "theta_for_C.c" = Huge_combined_theta_A, "delta-rho_for_C.m" = Huge_combined_deltarho_Y, "delta-rho_for_C.c" = Huge_combined_deltarho_A)
#example
Table1Small
Table1 <- rbind.data.frame(Table1Small, Table1Medium, Table1Large, Table1Huge)
Table1
write.csv(Table1,file="Table1.csv")
write.csv(SmallTable1base_m,file="SmallTable1base_m.csv")
write.csv(MediumTable1base_m,file="MediumTable1base_m.csv")
write.csv(LargeTable1base_m,file="LargeTable1base_m.csv")
write.csv(HugeTable1base_m,file="HugeTable1base_m.csv")
#Extracting target the abundance of individuals as list
x_map_target <- Small
Small_Y_each <- lapply(1:length(x_map_target), function(i) x_map_target[[i]]$YotsumonTotalL)
Small_A_each <- lapply(1:length(x_map_target), function(i) x_map_target[[i]]$AzukiTotalL)
#Number of NAs to be added is the maximum of optimal E among all replicates
no_NAs_Small <- max(SmallTable1base_m$No.NA)
#Combining all replicate
Small_Y_allcombined <- unlist(lapply(1:16, function(i) c(Small_Y_each[[i]], rep(NA, no_NAs_Small))))
Small_A_allcombined <- unlist(lapply(1:16, function(i) c(Small_A_each[[i]], rep(NA, no_NAs_Small))))
Small_Y_allcombined
[1] 8 0 0 0 110 4 12 37 192 21 11 133 55 9 98 197 76 30
[19] 150 141 88 100 91 119 76 130 NA NA NA NA NA NA NA NA NA 8
[37] 0 0 16 117 0 0 171 62 18 125 197 139 113 NA NA NA NA NA
[55] NA NA NA NA 8 0 0 4 114 21 0 75 200 14 38 165 106 24
[73] 141 180 141 83 143 195 124 91 NA NA NA NA NA NA NA NA NA 8
[91] 0 0 0 88 21 3 109 213 11 NA NA NA NA NA NA NA NA NA
[109] 8 0 0 0 98 14 0 119 195 5 116 175 192 NA NA NA NA NA
[127] NA NA NA NA 8 0 0 0 94 16 1 36 180 97 2 138 169 135
[145] 6 118 129 53 10 94 96 54 0 56 78 100 0 53 NA NA NA NA
[163] NA NA NA NA NA 8 0 0 0 110 39 1 41 233 58 3 141 179
[181] 82 8 112 106 68 24 70 120 114 48 100 NA NA NA NA NA NA NA
[199] NA NA 8 0 0 0 98 23 0 35 208 47 2 132 158 128 23 128
[217] 112 96 61 NA NA NA NA NA NA NA NA NA 8 0 0 0 146 14
[235] 0 72 208 59 NA NA NA NA NA NA NA NA NA 8 0 0 0 111
[253] 62 5 51 246 62 4 158 177 105 72 133 NA NA NA NA NA NA NA
[271] NA NA 8 0 0 0 139 33 1 95 214 128 34 199 138 NA NA NA
[289] NA NA NA NA NA NA 8 0 0 0 135 24 0 53 215 87 20 402
[307] NA NA NA NA NA NA NA NA NA 8 0 0 0 107 15 0 59 173
[325] 0 3 112 NA NA NA NA NA NA NA NA NA 8 0 0 2 101 10
[343] 0 67 151 6 8 246 NA NA NA NA NA NA NA NA NA 8 0 0
[361] 0 139 27 0 129 238 NA NA NA NA NA NA NA NA NA 8 0 0
[379] 0 114 16 1 123 179 12 NA NA NA NA NA NA NA NA NA
#Visualization
SmallPopDynamics <- data.frame("Time_step" = rep(c(1:length(Small_Y_allcombined)),2), "Population" = c(Small_Y_allcombined, Small_A_allcombined), "Species" = c(rep("C.m",length(Small_Y_allcombined)), rep("C.c",length(Small_A_allcombined))))
PopulationSmall <- ggplot(data = SmallPopDynamics, aes(x = Time_step, y = Population, group = Species, color = Species)) + geom_line()
print(PopulationSmall)
#Extracting target the abundance of individuals as list
x_map_target <- Medium
Medium_Y_each <- lapply(1:length(x_map_target), function(i) x_map_target[[i]]$YotsumonTotalL)
Medium_A_each <- lapply(1:length(x_map_target), function(i) x_map_target[[i]]$AzukiTotalL)
#Number of NAs to be added is the maximum of optimal E among all replicates
no_NAs_Medium <- max(MediumTable1base_m$No.NA)
#Combining all replicate
Medium_Y_allcombined <- unlist(lapply(1:16, function(i) c(Medium_Y_each[[i]], rep(NA, no_NAs_Medium))))
Medium_Y_allcombined
[1] 8 0 0 1 93 9 0 17 175 22 18 97 169 5 124 153 115 36
[19] 133 157 125 81 NA NA NA NA NA NA NA NA NA NA 8 0 0 0
[37] 134 29 2 3 243 42 2 317 97 1 85 113 155 58 NA NA NA NA
[55] NA NA NA NA NA NA 8 0 0 4 110 22 1 179 215 38 NA NA
[73] NA NA NA NA NA NA NA NA 8 0 0 1 74 15 0 79 184 54
[91] 71 177 160 109 105 NA NA NA NA NA NA NA NA NA NA 8 0 0
[109] 0 98 51 0 2 199 23 0 10 125 72 1 14 119 69 0 6 76
[127] 113 4 7 59 7 0 0 78 NA NA NA NA NA NA NA NA NA NA
[145] 8 0 0 0 125 58 0 1 208 30 7 2 134 70 0 11 131 30
[163] 1 12 80 43 0 5 10 5 0 0 18 11 0 0 NA NA NA NA
[181] NA NA NA NA NA NA 8 0 0 0 106 26 1 23 207 60 1 53
[199] 117 75 3 64 112 101 12 123 130 126 6 33 9 59 4 51 1 NA
[217] NA NA NA NA NA NA NA NA NA 8 0 0 0 114 31 0 9 208
[235] 106 1 57 140 142 1 135 99 93 2 143 NA NA NA NA NA NA NA
[253] NA NA NA 8 0 0 0 127 38 0 41 232 139 9 153 180 162 33
[271] 105 NA NA NA NA NA NA 8 0 1 0 106 45 2
[289] 2 157 206 25 9 227 137 7 17 70 83 12 12 108 129 0 9 136
[307] 92 1 17 51 34 21 10 6 75 13 10 24 99 61 14 101 226 144 NA NA NA NA NA NA NA NA NA NA 8 0 0 0
[343] 90 62 1 31 205 95 16 148 181 158 10 157 102 76 13 78 NA NA
[361] NA NA NA NA NA NA NA NA 8 0 0 0 110 76 1 52 239 88
[379] 9 207 168 269 24 105 NA NA NA NA NA NA NA NA NA NA 8 0
[397] 0 0 99 23 0 41 167 21 13 107 101 61 64 128 NA NA NA NA
[415] NA NA NA NA NA NA 8 0 0 0 75 15 0 92 31 0 0 26
[433] 20 3 0 22 60 5 0 47 78 14 0 110 75 NA NA NA NA NA
[451] NA NA NA NA NA 8 0 0 0 128 20 1 84 191 43 9 153 NA
[469] NA NA NA NA NA NA NA NA NA 8 0 0 0 122 5 0 138 221
[487] 22 NA NA NA NA NA NA NA NA NA NA
#Visualization
MediumPopDynamics <- data.frame("Time_step" = rep(c(1:length(Medium_Y_allcombined)),2), "Population" = c(Medium_Y_allcombined, Medium_A_allcombined), "Species" = c(rep("C.m",length(Medium_Y_allcombined)), rep("C.c",length(Medium_A_allcombined))))
PopulationMedium <- ggplot(data = MediumPopDynamics, aes(x = Time_step, y = Population, group = Species, color = Species)) + geom_line()
print(PopulationMedium)
#Extracting target the abundance of individuals as list
x_map_target <- Large
Large_Y_each <- lapply(1:length(x_map_target), function(i) x_map_target[[i]]$YotsumonTotalL)
Large_A_each <- lapply(1:length(x_map_target), function(i) x_map_target[[i]]$AzukiTotalL)
no_NAs_Large <- max(LargeTable1base_m$No.NA)
#Combining all replicate
Large_Y_allcombined <- unlist(lapply(1:16, function(i) c(Large_Y_each[[i]], rep(NA, no_NAs_Large))))
Large_A_allcombined <- unlist(lapply(1:16, function(i) c(Large_A_each[[i]],rep(NA, no_NAs_Large))))
Large_Y_allcombined
[1] 8 0 0 19 111 36 2 2 167 66 8 186 37 15 0 133 52 9
[19] 119 126 16 5 156 106 32 13 112 107 86 25 63 82 97 65 NA NA
[37] NA NA NA NA NA NA NA NA 8 0 0 10 162 31 7 126 176 53
[55] 102 163 185 164 113 144 247 148 140 NA NA NA NA NA NA NA NA NA
[73] NA 8 0 0 0 124 16 0 140 234 33 121 190 251 72 201 154 242
[91] 127 185 NA NA NA NA NA NA NA NA NA NA 8 0 0 0 77 19
[109] 0 36 207 22 4 56 178 8 2 118 140 6 22 155 67 15 71 97
[127] 137 5 110 85 56 5 NA NA NA NA NA
[145] NA NA NA NA NA 8 0 0 0 102 21 1 139 206 45 42 241 245
[163] 173 119 14 NA NA NA NA NA NA NA NA NA NA 8 0 0 0 121
[181] 74 5 9 162 162 10 4 208 161 5 10 96 56 1 0 82 98 1
[199] 0 63 77 1 0 58 NA NA NA NA NA NA NA NA NA NA 8 0
[217] 0 0 109 72 0 72 157 36 0 12 231 92 2 14 136 76 8 8
[235] 126 90 1 0 93 64 4 0 NA NA NA NA NA NA NA NA NA NA
[253] 8 0 0 0 102 40 0 8 178 59 1 14 172 105 0 21 121 67
[271] 1 14 61 91 2 10 37 71 4 6 16 NA NA NA NA NA NA NA
[289] NA NA NA 8 0 0 0 123 50 3 28 231 128 10 133 192 147 11
[307] 165 NA NA NA NA NA NA NA NA NA NA 8 0 0 0 114 41 0
[325] 36 243 148 27 172 133 86 7 154 123 101 16 151 NA NA NA NA NA
[343] NA NA NA NA NA 8 0 0 0 132 50 1 60 146 128 25 201 56
[361] 152 6 106 16 129 10 112 NA NA NA NA NA NA NA NA NA NA 8
[379] 0 0 0 142 74 1 31 223 134 19 119 197 185 13 157 81 86 28
[397] 92 NA NA NA NA NA NA NA NA NA NA 8 0 0 0 104 9 0
[415] 31 186 9 0 59 136 20 24 49 137 44 4 34 113 73 0 50 86
[433] 85 26 0 0 0 59 0 0 26 138 19 0 142 152 5 94 150 NA NA NA NA NA NA NA NA 8 0 0 0 107 29 0 44
[469] 212 20 0 68 110 35 10 143 155 126 15 130 125 116 5 129 NA NA
[487] NA NA NA NA NA NA NA NA 8 0 0 0 132 44 1 10 234 120
[505] 5 141 110 163 9 135 NA NA NA NA NA NA NA NA NA NA 8 0
[523] 0 0 137 27 0 116 201 25 0 155 141 69 57 106 NA NA NA NA
[541] NA NA NA NA NA NA
#Visualization
PopulationLarge <- ggplot(data = LargePopDynamics, aes(x = Time_step, y = Population, group = Species, color = Species)) + geom_line()
print(PopulationLarge)
#Extracting target the abundance of individuals as list
x_map_target <- Huge
Huge_Y_each <- lapply(1:length(x_map_target), function(i) x_map_target[[i]]$YotsumonTotalL)
Huge_A_each <- lapply(1:length(x_map_target), function(i) x_map_target[[i]]$AzukiTotalL)
#Number of NAs to be added is the maximum of optimal E among all replicates
no_NAs_Huge <- max(HugeTable1base_m$No.NA)
#Combining all replicate
Huge_Y_allcombined <- unlist(lapply(1:16, function(i) c(Huge_Y_each[[i]], rep(NA, no_NAs_Huge))))
Huge_A_allcombined <- unlist(lapply(1:16, function(i) c(Huge_A_each[[i]],rep(NA, no_NAs_Huge))))
Huge_Y_allcombined
[1] 8 0 0 1 127 24 1 4 173 59 7 94 189 61 2 157 76 15
[19] 17 165 37 0 40 130 44 0 56 125 70 2 65 86 98 3 69 NA
[37] NA NA NA NA NA NA NA NA NA 8 0 0 0 147 56 2 148 260
[55] 102 54 146 130 134 200 247 NA NA NA NA NA NA NA NA NA NA 8
[73] 0 0 0 111 57 0 52 204 157 11 173 171 156 39 129 117 111 116
[91] 82 NA NA NA NA NA NA NA NA NA NA 8 0 0 0 88 18 0
[109] 59 33 0 0 29 3 0 0 22 10 0 0 11 8 0 0 5 4
[127] 0 0 1 2 1 0 NA NA NA NA NA NA NA NA NA NA 8 0
[145] 0 0 93 61 1 19 216 103 11 7 188 15 1 3 121 23 0 3
[163] 85 6 0 0 23 6 0 0 2 1 0 NA NA NA NA NA NA NA
[181] NA NA NA 8 0 0 0 68 44 0 6 97 22 3 1 258 129 0
[199] 3 129 98 1 0 106 73 2 0 76 66 3 0 79 NA NA NA NA
[217] NA NA NA NA NA NA 8 0 0 0 92 57 3 93 111 157 14 0
[235] 175 182 4 0 88 35 1 0 91 77 9 0 48 54 2 0 31 NA
[253] NA NA NA NA NA NA NA NA NA 8 0 0 0 98 28 1 16 139
[271] 82 3 13 228 200 7 6 87 52 7 0 63 85 5 0 167 147 90
[289] 1 NA NA NA NA NA NA NA NA NA NA 8 0 0 0 135 75 2
[307] 15 152 202 7 83 227 104 0 144 90 89 8 107 124 108 3 53 31
[325] 52 4 40 NA NA NA NA NA NA NA NA NA NA 8 0 0 0 103
[343] 38 0 39 225 117 6 177 212 165 7 139 118 109 30 42 125 117 8
[361] 31 72 40 3 11 NA NA NA NA NA NA NA NA NA NA 8 0 0
[379] 0 98 54 0 9 197 115 6 27 216 187 0 78 133 113 9 99 99
[397] 111 0 74 92 63 4 74 101 41 5 2 0 NA NA NA NA NA NA
[415] NA NA NA NA 8 0 0 0 100 18 0 17 165 56 0 6 187 188
[433] 10 0 103 86 0 0 4 125 0 0 5 5 2 0 0 1 21 3
[451] NA NA NA NA NA NA NA NA NA NA 8 0 0 0 94 49 2 6
[469] 146 30 0 1 154 68 1 0 39 97 5 0 0 112 7 0 NA NA
[487] NA NA NA NA NA NA NA NA 8 0 0 0 113 65 1 64 233 109
[505] 12 140 124 147 42 143 NA NA NA NA NA NA NA NA NA NA 8 0
[523] 0 0 73 34 0 7 152 25 0 24 137 37 4 51 114 57 0 102
[541] 132 76 1 138 NA NA NA NA NA NA NA NA NA NA 8 0 0 0
[559] 77 12 0 20 199 112 14 64 121 141 9 69 108 NA NA NA NA NA
[577] NA NA NA NA NA
#Visualization
HugePopDynamics <- data.frame("Time_step" = rep(c(1:length(Huge_Y_allcombined)),2), "Population" = c(Huge_Y_allcombined, Huge_A_allcombined), "Species" = c(rep("C.m",length(Huge_Y_allcombined)), rep("C.c",length(Huge_A_allcombined))))
PopulationHuge <- ggplot(data = HugePopDynamics, aes(x = Time_step, y = Population, group = Species, color = Species)) + geom_line()
print(PopulationHuge)
–The Warning message is an expected one since we set the library and prediction identical.
#data combined
Small_YA_allcombined <- data.frame("Small_Y" = Small_Y_allcombined, "Small_A" = Small_A_allcombined)
#simplex projection for Y and A, respectively
Small_simplex_output_Y <- simplex(Small_YA_allcombined$Small_Y,lib = c(1, nrow(Small_YA_allcombined)), pred = c(1, nrow(Small_YA_allcombined)))
Small_simplex_output_A <- simplex(Small_YA_allcombined$Small_A,lib = c(1, nrow(Small_YA_allcombined)), pred = c(1, nrow(Small_YA_allcombined)))
#optimal embedding dimension E*
Small_optE_Y <- which.max(Small_simplex_output_Y$rho)
Small_optE_A <- which.max(Small_simplex_output_A$rho)
#univariate S-map
Small_s_map_outputY <- s_map(Small_YA_allcombined$Small_Y,lib = c(1, nrow(Small_YA_allcombined)), pred = c(1, nrow(Small_YA_allcombined)), E = Small_optE_Y)
Small_s_map_outputA <- s_map(Small_YA_allcombined$Small_A,lib = c(1, nrow(Small_YA_allcombined)), pred = c(1, nrow(Small_YA_allcombined)), E = Small_optE_A)
#optimal theta
Small_opttheta_Y <- Small_s_map_outputY$theta[which.max(Small_s_map_outputY$rho)]
Small_opttheta_A <- Small_s_map_outputA$theta[which.max(Small_s_map_outputA$rho)]
#data combined
Medium_YA_allcombined <- data.frame("Medium_Y" = Medium_Y_allcombined, "Medium_A" = Medium_A_allcombined)
#simplex projection for Y and A, respectively
Medium_simplex_output_Y <- simplex(Medium_YA_allcombined$Medium_Y,lib = c(1, nrow(Medium_YA_allcombined)), pred = c(1, nrow(Medium_YA_allcombined)))
Medium_simplex_output_A <- simplex(Medium_YA_allcombined$Medium_A,lib = c(1, nrow(Medium_YA_allcombined)), pred = c(1, nrow(Medium_YA_allcombined)))
#optimal embedding dimension E*
Medium_optE_Y <- which.max(Medium_simplex_output_Y$rho)
Medium_optE_A <- which.max(Medium_simplex_output_A$rho)
#univariate S-map
Medium_s_map_outputY <- s_map(Medium_YA_allcombined$Medium_Y,lib = c(1, nrow(Medium_YA_allcombined)), pred = c(1, nrow(Medium_YA_allcombined)),E = Medium_optE_Y)
Medium_s_map_outputA <- s_map(Medium_YA_allcombined$Medium_A,lib = c(1, nrow(Medium_YA_allcombined)), pred = c(1, nrow(Medium_YA_allcombined)),E = Medium_optE_A)
#optimal theta
Medium_opttheta_Y <- Medium_s_map_outputY$theta[which.max(Medium_s_map_outputY$rho)]
Medium_opttheta_A <- Medium_s_map_outputA$theta[which.max(Medium_s_map_outputA$rho)]
#data combined
Large_YA_allcombined <- data.frame("Large_Y" = Large_Y_allcombined, "Large_A" = Large_A_allcombined)
#simplex projection for Y and A, respectively
Large_simplex_output_Y <- simplex(Large_YA_allcombined$Large_Y,lib = c(1, nrow(Large_YA_allcombined)), pred = c(1, nrow(Large_YA_allcombined)))
Large_simplex_output_A <- simplex(Large_YA_allcombined$Large_A,lib = c(1, nrow(Large_YA_allcombined)), pred = c(1, nrow(Large_YA_allcombined)))
#optimal embedding dimension E*
Large_optE_Y <- which.max(Large_simplex_output_Y$rho)
Large_optE_A <- which.max(Large_simplex_output_A$rho)
#univariate S-map
Large_s_map_outputY <- s_map(Large_YA_allcombined$Large_Y,lib = c(1, nrow(Large_YA_allcombined)), pred = c(1, nrow(Large_YA_allcombined)),E = Large_optE_Y)
Large_s_map_outputA <- s_map(Large_YA_allcombined$Large_A,lib = c(1, nrow(Large_YA_allcombined)), pred = c(1, nrow(Large_YA_allcombined)),E = Large_optE_A)
#optimal theta
Large_opttheta_Y <- Large_s_map_outputY$theta[which.max(Large_s_map_outputY$rho)]
#data combined
Huge_YA_allcombined <- data.frame("Huge_Y" = Huge_Y_allcombined, "Huge_A" = Huge_A_allcombined)
#simplex projection for Y and A, respectively
Huge_simplex_output_Y <- simplex(Huge_YA_allcombined$Huge_Y,lib = c(1, nrow(Huge_YA_allcombined)), pred = c(1, nrow(Huge_YA_allcombined)))
Huge_simplex_output_A <- simplex(Huge_YA_allcombined$Huge_A,lib = c(1, nrow(Huge_YA_allcombined)), pred = c(1, nrow(Huge_YA_allcombined)))
#optimal embedding dimension E*
Huge_optE_Y <- which.max(Huge_simplex_output_Y$rho)
Huge_optE_A <- which.max(Huge_simplex_output_A$rho)
#univariate S-map
Huge_s_map_outputY <- s_map(Huge_YA_allcombined$Huge_Y,lib = c(1, nrow(Huge_YA_allcombined)), pred = c(1, nrow(Huge_YA_allcombined)),E = Huge_optE_Y)
Huge_s_map_outputA <- s_map(Huge_YA_allcombined$Huge_A,lib = c(1, nrow(Huge_YA_allcombined)), pred = c(1, nrow(Huge_YA_allcombined)),E = Huge_optE_A)
#optimal theta
Huge_opttheta_Y <- Huge_s_map_outputY$theta[which.max(Huge_s_map_outputY$rho)]
Huge_opttheta_A <- Huge_s_map_outputA$theta[which.max(Huge_s_map_outputA$rho)]
#For Small
par(mfrow = c(1, 2), mar = c(4, 4, 1, 1), oma = c(3, 3, 3, 0))
plot(Small_simplex_output_Y$E, Small_simplex_output_Y$rho, type = "l", xlab = " ", ylab = "", xlim = c(0, 10), ylim = c(0, 1), cex.lab = 2, cex.axis = 1,main="(a) C.maculatus")
abline(h = 0, lty = 3, col = "black",lwd = 2)
plot(Small_simplex_output_A$E, Small_simplex_output_A$rho, type = "l", xlab = " ", ylab = "", xlim = c(0, 10), ylim = c(-0.1, 1), cex.axis = 2, cex.axis = 1, main = "(b) C.chinensis")
abline(h = 0, lty = 3, col ="black", lwd = 2)
mtext(side = 3, line = 1, outer = TRUE, text = "Optimal Embedding Dimension for Small cages")
#For Medium
par(mfrow = c(1, 2), mar = c(4, 4, 1, 1), oma = c(3, 3, 3, 0))
plot(Medium_simplex_output_Y$E, Medium_simplex_output_Y$rho, type = "l", xlab = " ", ylab = "", xlim = c(0, 10), ylim = c(0, 1), cex.lab = 2, cex.axis = 1,main="(c) C.maculatus")
abline(h = 0, lty = 3, col = "black", lwd = 2)
plot(Medium_simplex_output_A$E, Medium_simplex_output_A$rho, type = "l", xlab = " ", ylab = "", xlim = c(0, 10), ylim = c(0, 1), cex.axis = 2, cex.axis = 1, main = "(d) C.chinensis")
abline(h = 0, lty = 3, col ="black", lwd = 2)
mtext(side = 3, line = 1, outer = TRUE, text = "Optimal Embedding Dimension for Medium cages")
#For Large
par(mfrow = c(1, 2), mar = c(4, 4, 1, 1), oma = c(3, 3, 3, 0))
plot(Large_simplex_output_Y$E, Large_simplex_output_Y$rho, type = "l", xlab = " ", ylab = "", xlim = c(0, 10), ylim = c(0, 1), cex.lab = 2, cex.axis = 1,main="(e) C.maculatus")
abline(h = 0, lty = 3, col = "black", lwd = 2)
plot(Large_simplex_output_A$E, Large_simplex_output_A$rho, type = "l", xlab = " ", ylab = "", xlim = c(0, 10), ylim = c(-0.1, 1), cex.axis = 2, cex.axis = 1, main = "(f) C.chinensis")
abline(h = 0, lty = 3, col ="black", lwd = 2)
mtext(side = 3, line = 1, outer = TRUE, text = "Optimal Embedding Dimension for Large cages")
#For Huge
par(mfrow = c(1, 2), mar = c(4, 4, 1, 1), oma = c(3, 3, 3, 0))
plot(Huge_simplex_output_Y$E, Huge_simplex_output_Y$rho, type = "l", xlab = " ", ylab = "", xlim = c(0, 10), ylim = c(0, 1), cex.lab = 2, cex.axis = 1,main="(g) C.maculatus")
abline(h = 0, lty = 3, col = "black", lwd = 2)
plot(Huge_simplex_output_A$E, Huge_simplex_output_A$rho, type = "l", xlab = " ", ylab = "", xlim = c(0, 10), ylim = c(0, 1), cex.axis = 2, cex.axis = 1, main = "(h) C.chinensis")
abline(h = 0, lty = 3, col ="black", lwd = 2)
mtext(side = 3, line = 1, outer = TRUE, text = "Optimal Embedding Dimension for Huge cages")
#For Small
par(mfrow = c(1, 2), mar = c(4, 4, 1, 1), oma = c(3, 3, 3, 0))
plot(Small_s_map_outputY$theta, Small_s_map_outputY$rho, type = "l", xlab = "", ylab = "", cex.main = 1,ylim = c(0, 0.9), cex.lab = 2, cex.axis = 2, main = "(a) C.maculatus")
abline(h = 0,lty = 3, col = "black",lwd = 2)
plot(Small_s_map_outputA$theta, Small_s_map_outputA$rho, type = "l", xlab = "", ylab = "", cex.main = 1, ylim = c(0, 0.9), cex.lab = 2, cex.axis = 2, main = "(b) C.chinensis")
abline(h = 0,lty = 3, col = "black",lwd = 2)
mtext(side = 3, line = 1, outer = TRUE, text = "Optimal theta (nonlinearity) for Small cage")
#For Medium
par(mfrow = c(1, 2), mar = c(4, 4, 1, 1), oma = c(3, 3, 3, 0))
plot(Medium_s_map_outputY$theta, Medium_s_map_outputY$rho, type = "l", xlab = "", ylab = "", cex.main = 1,ylim = c(0, 0.9), cex.lab = 2, cex.axis = 2, main = "(c) C.maculatus")
abline(h = 0,lty = 3, col = "black",lwd = 2)
plot(Medium_s_map_outputA$theta, Medium_s_map_outputA$rho, type = "l", xlab = "", ylab = "", cex.main = 1, ylim = c(0, 0.9), cex.lab = 2, cex.axis = 2, main = "(d) C.chinensis")
abline(h = 0,lty = 3, col = "black",lwd = 2)
mtext(side = 3, line = 1, outer = TRUE, text = "Optimal theta (nonlinearity) for Medium cage")
#For Large
par(mfrow = c(1, 2), mar = c(4, 4, 1, 1), oma = c(3, 3, 3, 0))
plot(Large_s_map_outputY$theta, Large_s_map_outputY$rho, type = "l", xlab = "", ylab = "", cex.main = 1,ylim = c(0, 0.9), cex.lab = 2, cex.axis = 2, main = "(e) C.maculatus")
abline(h = 0,lty = 3, col = "black",lwd = 2)
plot(Large_s_map_outputA$theta, Large_s_map_outputA$rho, type = "l", xlab = "", ylab = "", cex.main = 1, ylim = c(0, 0.9), cex.lab = 2, cex.axis = 2, main = "(f) C.chinensis")
abline(h = 0,lty = 3, col = "black",lwd = 2)
mtext(side = 3, line = 1, outer = TRUE, text = "Optimal theta (nonlinearity) for Large cage")
#For Huge
par(mfrow = c(1, 2), mar = c(4, 4, 1, 1), oma = c(3, 3, 3, 0))
plot(Huge_s_map_outputY$theta, Huge_s_map_outputY$rho, type = "l", xlab = "", ylab = "", cex.main = 1,ylim = c(0, 0.9), cex.lab = 2, cex.axis = 2, main = "(g) C.maculatus")
abline(h = 0,lty = 3, col = "black",lwd = 2)
plot(Huge_s_map_outputA$theta, Huge_s_map_outputA$rho, type = "l", xlab = "", ylab = "", cex.main = 1, ylim = c(0, 0.9), cex.lab = 2, cex.axis = 2, main = "(h) C.chinensis")
abline(h = 0,lty = 3, col = "black",lwd = 2)
mtext(side = 3, line = 1, outer = TRUE, text = "Optimal theta (nonlinearity) for Huge cage")
– Y_xmap_A represents the causality from A (Azuki Cc) to Y (Yotsumon Cm) #F8766D (red)
– A_xmap_Y represents the causality from Y (Yotsumon Cm) to A (Azuki Cc) #00BFC4 (blue-green)
Small_lagY_xmap_A <- list()
Small_lagA_xmap_Y <- list()
#not showing warning
Small_lagY_xmap_A <- lapply(0:-4, function(i) ccm(Small_YA_allcombined, E = Small_optE_Y, lib_column = "Small_Y", target_column = "Small_A", lib_sizes = floor(seq(Small_optE_Y, nrow(Small_YA_allcombined), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
#example
Small_lagY_xmap_A[[2]]
Small_lagY_xmap_A_means <- list()
Small_lagA_xmap_Y_means <- list()
Small_lagY_xmap_A_means <- lapply(1:5, function(i) ccm_means(Small_lagY_xmap_A[[i]], na.rm = TRUE))
Small_lagA_xmap_Y_means <- lapply(1:5, function(i) ccm_means(Small_lagA_xmap_Y[[i]], na.rm = TRUE))
Small_lagA_xmap_Y_means <- lapply(1:5, function(i) data.frame(Small_lagA_xmap_Y_means[[i]], sd.rho = ccm_means(Small_lagA_xmap_Y[[i]], FUN = sd, na.rm = TRUE)$rho))
#With the maximum library size,since there is not replicate, so SD cannot be calculated and set as zero
for(i in 1:5){
}
for(i in 1:5){
Small_lagA_xmap_Y_means[[i]]$sd.rho[nrow(Small_lagA_xmap_Y_means[[i]])] <- 0
}
#example
–For Yotsumon (Cm)
#Generating Surrogate data; for each replicate of the time series, we generated 100 surrogate data
set.seed(12345)
Small_Y_each_s <- list()
Small_Y_each_s <- lapply(1:length(Small_Y_each), function(i) make_surrogate_data(Small_Y_each[[i]], method = c("ebisuzaki"), num_surr = 100))
Small_Y_each_s <- lapply(1:length(Small_Y_each_s), function(i) rbind.data.frame(Small_Y_each_s[[i]], matrix(rep(NA, 100*(no_NAs_Small)), nrow = no_NAs_Small)))
Small_Y_each_s2 <- NULL
for(i in 1:length(Small_Y_each_s)){
Small_Y_each_s2 <- rbind.data.frame(Small_Y_each_s2, Small_Y_each_s[[i]])
}
Small_Y_each_s2
#Surrogate Y (Ys) xmap A (CCM)
Small_AYs <- cbind.data.frame(Small_Y_each_s2, "Small_A" = Small_YA_allcombined$Small_A)
Small_lagYs_xmap_A <- list()
for(j in 1:5) {
Small_lagYs_xmap_A[[j]] <- lapply(1:100, function(i) ccm(Small_AYs, E = Small_optE_Y, lib = c(1,nrow(Small_AYs)), pred = c(1,nrow(Small_AYs)), lib_column = paste("V",i,sep=""), target_column = "Small_A", lib_sizes = floor(seq(Small_optE_Y, length(Small_AYs[, 1]),length = 10)),tp = (1 - j), num_sample = 100, replace = F, silent = T, RNGseed = 1234))
}
#example
Small_lagYs_xmap_A[[4]][[16]]
#mean and sd of surrogate Y xmap A
Small_lagYs_xmap_A_means_list <- list()
for(j in 1:5) {
Small_lagYs_xmap_A_means_list[[j]] <- lapply(1:length(Small_lagYs_xmap_A[[j]]), function(i) ccm_means(Small_lagYs_xmap_A[[j]][[i]], na.rm = TRUE))
}
#combining the list
Small_lagYs_xmap_A_means <- list()
for(j in 1:5) {
Small_lagYs_xmap_A_means[[j]] <- Small_lagYs_xmap_A_means_list[[j]][[1]]$rho
for(i in 2:length(Small_lagYs_xmap_A_means_list[[j]])){
Small_lagYs_xmap_A_means[[j]] <- cbind.data.frame(Small_lagYs_xmap_A_means[[j]], Small_lagYs_xmap_A_means_list[[j]][[i]]$rho)
}
}
#calculating means and sds of "means"
Small_lagYs_xmap_A_means_sds <- list()
for(j in 1:5) {
Small_lagYs_xmap_A_means_means[[j]] <- apply(Small_lagYs_xmap_A_means[[j]], 1, mean)
Small_lagYs_xmap_A_means_sds[[j]] <- apply(Small_lagYs_xmap_A_means[[j]], 1, sd)
}
#example
Small_lagYs_xmap_A_means_sds[[1]]
[1] 0.02096274 0.06128845 0.07976680 0.09243328 0.10352799 0.10601242
Small_lagYs_xmap_A_summary <- list()
for(j in 1:5) {
Small_lagYs_xmap_A_summary[[j]] <- data.frame("lib_size" = Small_lagYs_xmap_A_means_list[[j]][[1]]$lib_size, "mean" = Small_lagYs_xmap_A_means_means[[j]], "sd" = Small_lagYs_xmap_A_means_sds[[j]])
}
#example
Small_lagYs_xmap_A_summary[[3]]
–For Azuki (Cc)
#Generating Surrogate data
set.seed(12345)
Small_A_each_s <- list()
Small_A_each_s <- lapply(1:length(Small_A_each), function(i) make_surrogate_data(Small_A_each[[i]], method = c("ebisuzaki"), num_surr = 100))
Small_A_each_s <- lapply(1:length(Small_A_each_s), function(i) rbind.data.frame(Small_A_each_s[[i]], matrix(rep(NA,100*(no_NAs_Small)), nrow = no_NAs_Small)))
Small_A_each_s2 <- NULL
for(i in 1:length(Small_A_each_s)){
}
Small_A_each_s2
#Surrogate A (As) xmap Y (CCM)
Small_YAs <- cbind.data.frame(Small_A_each_s2, "Small_Y" = Small_YA_allcombined$Small_Y)
Small_lagAs_xmap_Y <- list()
for(j in 1:5) {
Small_lagAs_xmap_Y[[j]] <- lapply(1:100, function(i) ccm(Small_YAs, E = Small_optE_A, lib = c(1,nrow(Small_YAs)), pred = c(1,nrow(Small_YAs)), lib_column = paste("V", i, sep = ""), target_column = "Small_Y", lib_sizes = floor(seq(Small_optE_A, length(Small_YAs[, 1]),length = 10)),tp = (1 - j), num_sample = 100, replace = F, silent = T, RNGseed = 1234))
}
#example
Small_lagAs_xmap_Y[[4]][[16]]
#mean and sd of surrogate A xmap Y
Small_lagAs_xmap_Y_means_list <- list()
for(j in 1:5) {
Small_lagAs_xmap_Y_means_list[[j]] <- lapply(1:length(Small_lagAs_xmap_Y[[j]]), function(i) ccm_means(Small_lagAs_xmap_Y[[j]][[i]], na.rm = TRUE))
}
#combining the list
Small_lagAs_xmap_Y_means <- list()
for(j in 1:5) {
Small_lagAs_xmap_Y_means[[j]] <- Small_lagAs_xmap_Y_means_list[[j]][[1]]$rho
for(i in 2:length(Small_lagAs_xmap_Y_means_list[[j]])){
Small_lagAs_xmap_Y_means[[j]] <- cbind.data.frame(Small_lagAs_xmap_Y_means[[j]], Small_lagAs_xmap_Y_means_list[[j]][[i]]$rho)
}
}
#calculating means and sds of "means"
Small_lagAs_xmap_Y_means_means <- list()
Small_lagAs_xmap_Y_means_sds <- list()
for(j in 1:5) {
Small_lagAs_xmap_Y_means_means[[j]] <- apply(Small_lagAs_xmap_Y_means[[j]], 1, mean)
Small_lagAs_xmap_Y_means_sds[[j]] <- apply(Small_lagAs_xmap_Y_means[[j]], 1, sd)
}
#example
Small_lagAs_xmap_Y_means_sds[[1]]
[1] 0.03264759 0.07746794 0.10148603 0.10986380
Small_lagAs_xmap_Y_summary <- list()
for(j in 1:5) {
}
#example
Small_lagAs_xmap_Y_summary[[3]]
Small_ccmY_xmap_A_all <- list()
for(i in 1:5){
Small_ccmY_xmap_A_all[[i]] <- data.frame(
"Library_size" = Small_lagY_xmap_A_means[[i]]$lib_size,
"Y_xmap_A_rho" = Small_lagY_xmap_A_means[[i]]$rho,
"Ys_xmap_A_rho" = Small_lagYs_xmap_A_summary[[i]]$mean,
"Y_xmap_A_rho_psd" = Small_lagY_xmap_A_means[[i]]$rho + Small_lagY_xmap_A_means[[i]]$sd.rho,
"Y_xmap_A_rho_msd" = Small_lagY_xmap_A_means[[i]]$rho - Small_lagY_xmap_A_means[[i]]$sd.rho,
"Ys_xmap_A_rho_psd" = Small_lagYs_xmap_A_summary[[i]]$mean + Small_lagYs_xmap_A_summary[[i]]$sd,
"Ys_xmap_A_rho_msd" = Small_lagYs_xmap_A_summary[[i]]$mean - Small_lagYs_xmap_A_summary[[i]]$sd
)
}
#example
Small_ccmY_xmap_A_all[[2]]
Small_ccmA_xmap_Y_all <- list()
for(i in 1:5){
Small_ccmA_xmap_Y_all[[i]] <- data.frame(
"Library_size" = Small_lagA_xmap_Y_means[[i]]$lib_size,
"A_xmap_Y_rho" = Small_lagA_xmap_Y_means[[i]]$rho,
"As_xmap_Y_rho" = Small_lagAs_xmap_Y_summary[[i]]$mean,
"A_xmap_Y_rho_psd" = Small_lagA_xmap_Y_means[[i]]$rho + Small_lagA_xmap_Y_means[[i]]$sd.rho,
"A_xmap_Y_rho_msd" = Small_lagA_xmap_Y_means[[i]]$rho - Small_lagA_xmap_Y_means[[i]]$sd.rho,
"As_xmap_Y_rho_psd" = Small_lagAs_xmap_Y_summary[[i]]$mean + Small_lagAs_xmap_Y_summary[[i]]$sd,
)
#example
Small_ccmA_xmap_Y_all[[2]]
#CCM Smallに関するplot
Small_ccm_ggplot <- list()
Small_subtitle <- c("(a)","(b)","(c)","(d)","(e)")
for(i in 1:5){
#Y_xmap_A:Cm xmap Cc #F8766D (red)
Small_ccm_ggplot[[i]] <- ggplot(data = Small_ccmY_xmap_A_all[[i]], aes(x = Library_size, y = Y_xmap_A_rho)) +
geom_line(linetype = 1, color = "#F8766D") +
geom_line(aes(x = Library_size, y = Y_xmap_A_rho_psd), linetype = 2, color = "#F8766D") +
geom_line(aes(x = Library_size, y = Y_xmap_A_rho_msd), linetype = 2, color = "#F8766D") +
geom_point(shape = 16, size = 3, color = "#F8766D") +
#A_xmap_Y:Cc xmap Cm #00BFC4 (blue-green)
geom_line(data = Small_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = A_xmap_Y_rho_psd), linetype = 2, color = "#00BFC4") +
geom_line(data = Small_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = A_xmap_Y_rho_msd), linetype = 2, color = "#00BFC4") +
geom_point(data = Small_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = A_xmap_Y_rho), shape = 17, size = 3, color = "#00BFC4") +
geom_line(data = Small_ccmY_xmap_A_all[[i]], aes(x = Library_size, y = Ys_xmap_A_rho_psd), linetype = 2, color = "#F8766D") +
geom_line(data = Small_ccmY_xmap_A_all[[i]], aes(x = Library_size, y = Ys_xmap_A_rho_msd), linetype = 2, color = "#F8766D") +
geom_point(data = Small_ccmY_xmap_A_all[[i]], aes(x = Library_size, y = Ys_xmap_A_rho), shape = 1, size = 3, color = "#F8766D") +
geom_line(data = Small_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = As_xmap_Y_rho_psd), linetype = 2, color = "#00BFC4") +
geom_point(data = Small_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = As_xmap_Y_rho), shape = 2, size = 3, color = "#00BFC4") +
ylim(-0.1, 1) +
xlab("") +
ggtitle(paste(Small_subtitle[i], "Small", "tp =", 1 - i, sep = " "))
}
for(i in 1:5) {
print(Small_ccm_ggplot[[i]])
Medium_lagY_xmap_A <- list()
Medium_lagA_xmap_Y <- list()
#not showing warning
Medium_lagY_xmap_A <- lapply(0:-4, function(i) ccm(Medium_YA_allcombined, E = Medium_optE_Y, lib_column = "Medium_Y", target_column = "Medium_A", lib_sizes = floor(seq(Medium_optE_Y, nrow(Medium_YA_allcombined), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
Medium_lagA_xmap_Y <- lapply(0:-4, function(i) ccm(Medium_YA_allcombined, E = Medium_optE_A, lib_column = "Medium_A", target_column = "Medium_Y", lib_sizes = floor(seq(Medium_optE_A, nrow(Medium_YA_allcombined), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
#example
Medium_lagY_xmap_A[[2]]
Medium_lagY_xmap_A_means <- list()
Medium_lagA_xmap_Y_means <- list()
Medium_lagY_xmap_A_means <- lapply(1:5, function(i) ccm_means(Medium_lagY_xmap_A[[i]], na.rm = TRUE))
Medium_lagA_xmap_Y_means <- lapply(1:5, function(i) ccm_means(Medium_lagA_xmap_Y[[i]], na.rm = TRUE))
Medium_lagY_xmap_A_means <- lapply(1:5, function(i) data.frame(Medium_lagY_xmap_A_means[[i]], sd.rho = ccm_means(Medium_lagY_xmap_A[[i]], FUN = sd, na.rm = TRUE)$rho))
Medium_lagA_xmap_Y_means <- lapply(1:5, function(i) data.frame(Medium_lagA_xmap_Y_means[[i]], sd.rho = ccm_means(Medium_lagA_xmap_Y[[i]], FUN = sd, na.rm = TRUE)$rho))
#With the maximum library size,since there is not replicate, so SD cannot be calculated and set as zero
for(i in 1:5){
Medium_lagY_xmap_A_means[[i]]$sd.rho[nrow(Medium_lagY_xmap_A_means[[i]])] <- 0
for(i in 1:5){
}
#example
Medium_lagA_xmap_Y_means[[3]]
–For Yotsumon (Cm)
#Generating Surrogate data
set.seed(12345)
Medium_Y_each_s <- list()
Medium_Y_each_s <- lapply(1:length(Medium_Y_each), function(i) make_surrogate_data(Medium_Y_each[[i]], method = c("ebisuzaki"), num_surr = 100))
Medium_Y_each_s <- lapply(1:length(Medium_Y_each_s), function(i) rbind.data.frame(Medium_Y_each_s[[i]], matrix(rep(NA, 100*(no_NAs_Medium)), nrow = no_NAs_Medium)))
Medium_Y_each_s2 <- NULL
for(i in 1:length(Medium_Y_each_s)){
Medium_Y_each_s2 <- rbind.data.frame(Medium_Y_each_s2, Medium_Y_each_s[[i]])
}
Medium_Y_each_s2
#Surrogate Y (Ys) xmap A (CCM)
Medium_AYs <- cbind.data.frame(Medium_Y_each_s2, "Medium_A" = Medium_YA_allcombined$Medium_A)
Medium_lagYs_xmap_A <- list()
for(j in 1:5) {
Medium_lagYs_xmap_A[[j]] <- lapply(1:100, function(i) ccm(Medium_AYs, E = Medium_optE_Y, lib = c(1,nrow(Medium_AYs)), pred = c(1,nrow(Medium_AYs)), lib_column = paste("V",i,sep=""), target_column = "Medium_A", lib_sizes = floor(seq(Medium_optE_Y, length(Medium_AYs[, 1]),length = 10)),tp = (1 - j), num_sample = 100, replace = F, silent = T, RNGseed = 1234))
}
#example
Medium_lagYs_xmap_A[[4]][[16]]
#mean and sd of surrogate Y xmap A
Medium_lagYs_xmap_A_means_list <- list()
for(j in 1:5) {
Medium_lagYs_xmap_A_means_list[[j]] <- lapply(1:length(Medium_lagYs_xmap_A[[j]]), function(i) ccm_means(Medium_lagYs_xmap_A[[j]][[i]], na.rm = TRUE))
}
#combining the list
Medium_lagYs_xmap_A_means <- list()
for(j in 1:5) {
Medium_lagYs_xmap_A_means[[j]] <- Medium_lagYs_xmap_A_means_list[[j]][[1]]$rho
for(i in 2:length(Medium_lagYs_xmap_A_means_list[[j]])){
Medium_lagYs_xmap_A_means[[j]] <- cbind.data.frame(Medium_lagYs_xmap_A_means[[j]], Medium_lagYs_xmap_A_means_list[[j]][[i]]$rho)
}
}
#calculating means and sds of "means"
Medium_lagYs_xmap_A_means_means <- list()
Medium_lagYs_xmap_A_means_sds <- list()
for(j in 1:5) {
Medium_lagYs_xmap_A_means_means[[j]] <- apply(Medium_lagYs_xmap_A_means[[j]], 1, mean)
}
#example
Medium_lagYs_xmap_A_means_sds[[1]]
[1] 0.03032287 0.07089069 0.08802803 0.09887281 0.10925049
Medium_lagYs_xmap_A_summary <- list()
for(j in 1:5) {
Medium_lagYs_xmap_A_summary[[j]] <- data.frame("lib_size" = Medium_lagYs_xmap_A_means_list[[j]][[1]]$lib_size, "mean" = Medium_lagYs_xmap_A_means_means[[j]], "sd" = Medium_lagYs_xmap_A_means_sds[[j]])
}
#example
Medium_lagYs_xmap_A_summary[[3]]
–For Azuki (Cc)
#Generating Surrogate data
set.seed(12345)
Medium_A_each_s <- list()
Medium_A_each_s <- lapply(1:length(Medium_A_each), function(i) make_surrogate_data(Medium_A_each[[i]], method = c("ebisuzaki"), num_surr = 100))
Medium_A_each_s <- lapply(1:length(Medium_A_each_s), function(i) rbind.data.frame(Medium_A_each_s[[i]], matrix(rep(NA,100*(no_NAs_Medium)), nrow = no_NAs_Medium)))
for(i in 1:length(Medium_A_each_s)){
Medium_A_each_s2
#Surrogate A (As) xmap Y (CCM)
Medium_YAs <- cbind.data.frame(Medium_A_each_s2, "Medium_Y" = Medium_YA_allcombined$Medium_Y)
Medium_lagAs_xmap_Y <- list()
for(j in 1:5) {
Medium_lagAs_xmap_Y[[j]] <- lapply(1:100, function(i) ccm(Medium_YAs, E = Medium_optE_A, lib = c(1,nrow(Medium_YAs)), pred = c(1,nrow(Medium_YAs)), lib_column = paste("V", i, sep = ""), target_column = "Medium_Y", lib_sizes = floor(seq(Medium_optE_A, length(Medium_YAs[, 1]),length = 10)),tp = (1 - j), num_sample = 100, replace = F, silent = T, RNGseed = 1234))
}
#example
Medium_lagAs_xmap_Y[[4]][[16]]
#mean and sd of surrogate A xmap Y
Medium_lagAs_xmap_Y_means_list <- list()
for(j in 1:5) {
Medium_lagAs_xmap_Y_means_list[[j]] <- lapply(1:length(Medium_lagAs_xmap_Y[[j]]), function(i) ccm_means(Medium_lagAs_xmap_Y[[j]][[i]], na.rm = TRUE))
}
#combining the list
Medium_lagAs_xmap_Y_means <- list()
for(j in 1:5) {
Medium_lagAs_xmap_Y_means[[j]] <- Medium_lagAs_xmap_Y_means_list[[j]][[1]]$rho
for(i in 2:length(Medium_lagAs_xmap_Y_means_list[[j]])){
Medium_lagAs_xmap_Y_means[[j]] <- cbind.data.frame(Medium_lagAs_xmap_Y_means[[j]], Medium_lagAs_xmap_Y_means_list[[j]][[i]]$rho)
}
}
#calculating means and sds of "means"
Medium_lagAs_xmap_Y_means_means <- list()
Medium_lagAs_xmap_Y_means_sds <- list()
for(j in 1:5) {
Medium_lagAs_xmap_Y_means_means[[j]] <- apply(Medium_lagAs_xmap_Y_means[[j]], 1, mean)
Medium_lagAs_xmap_Y_means_sds[[j]] <- apply(Medium_lagAs_xmap_Y_means[[j]], 1, sd)
}
#example
Medium_lagAs_xmap_Y_means_sds[[1]]
[1] 0.03593354 0.06558470 0.07813127 0.08933782 0.09360111
Medium_lagAs_xmap_Y_summary <- list()
for(j in 1:5) {
Medium_lagAs_xmap_Y_summary[[j]] <- data.frame("lib_size" = Medium_lagAs_xmap_Y_means_list[[j]][[1]]$lib_size, "mean" = Medium_lagAs_xmap_Y_means_means[[j]], "sd" = Medium_lagAs_xmap_Y_means_sds[[j]])
}
#example
Medium_lagAs_xmap_Y_summary[[3]]
Medium_ccmY_xmap_A_all <- list()
for(i in 1:5){
Medium_ccmY_xmap_A_all[[i]] <- data.frame(
"Library_size" = Medium_lagY_xmap_A_means[[i]]$lib_size,
"Y_xmap_A_rho" = Medium_lagY_xmap_A_means[[i]]$rho,
"Ys_xmap_A_rho" = Medium_lagYs_xmap_A_summary[[i]]$mean,
"Y_xmap_A_rho_psd" = Medium_lagY_xmap_A_means[[i]]$rho + Medium_lagY_xmap_A_means[[i]]$sd.rho,
"Y_xmap_A_rho_msd" = Medium_lagY_xmap_A_means[[i]]$rho - Medium_lagY_xmap_A_means[[i]]$sd.rho,
"Ys_xmap_A_rho_msd" = Medium_lagYs_xmap_A_summary[[i]]$mean - Medium_lagYs_xmap_A_summary[[i]]$sd
}
#example
Medium_ccmY_xmap_A_all[[2]]
Medium_ccmA_xmap_Y_all <- list()
for(i in 1:5){
Medium_ccmA_xmap_Y_all[[i]] <- data.frame(
"Library_size" = Medium_lagA_xmap_Y_means[[i]]$lib_size,
"A_xmap_Y_rho" = Medium_lagA_xmap_Y_means[[i]]$rho,
"As_xmap_Y_rho" = Medium_lagAs_xmap_Y_summary[[i]]$mean,
"A_xmap_Y_rho_psd" = Medium_lagA_xmap_Y_means[[i]]$rho + Medium_lagA_xmap_Y_means[[i]]$sd.rho,
"A_xmap_Y_rho_msd" = Medium_lagA_xmap_Y_means[[i]]$rho - Medium_lagA_xmap_Y_means[[i]]$sd.rho,
"As_xmap_Y_rho_psd" = Medium_lagAs_xmap_Y_summary[[i]]$mean + Medium_lagAs_xmap_Y_summary[[i]]$sd,
"As_xmap_Y_rho_msd" = Medium_lagAs_xmap_Y_summary[[i]]$mean - Medium_lagAs_xmap_Y_summary[[i]]$sd
}
#example
#CCM Mediumに関するplot
Medium_ccm_ggplot <- list()
Medium_subtitle <- c("(f)","(g)","(h)","(i)","(j)")
for(i in 1:5){
#Y_xmap_A:Cm xmap Cc #F8766D (red)
Medium_ccm_ggplot[[i]] <- ggplot(data = Medium_ccmY_xmap_A_all[[i]], aes(x = Library_size, y = Y_xmap_A_rho)) +
geom_line(linetype = 1, color = "#F8766D") +
geom_line(aes(x = Library_size, y = Y_xmap_A_rho_psd), linetype = 2, color = "#F8766D") +
geom_line(aes(x = Library_size, y = Y_xmap_A_rho_msd), linetype = 2, color = "#F8766D") +
geom_point(shape = 16, size = 3, color = "#F8766D") +
#A_xmap_Y:Cc xmap Cm #00BFC4 (blue-green)
geom_line(data = Medium_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = A_xmap_Y_rho), linetype = 1, color = "#00BFC4") +
geom_line(data = Medium_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = A_xmap_Y_rho_psd), linetype = 2, color = "#00BFC4") +
geom_line(data = Medium_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = A_xmap_Y_rho_msd), linetype = 2, color = "#00BFC4") +
geom_point(data = Medium_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = A_xmap_Y_rho), shape = 17, size = 3, color = "#00BFC4") +
geom_line(data = Medium_ccmY_xmap_A_all[[i]], aes(x = Library_size, y = Ys_xmap_A_rho), linetype = 1, color = "#F8766D") +
geom_line(data = Medium_ccmY_xmap_A_all[[i]], aes(x = Library_size, y = Ys_xmap_A_rho_psd), linetype = 2, color = "#F8766D") +
geom_line(data = Medium_ccmY_xmap_A_all[[i]], aes(x = Library_size, y = Ys_xmap_A_rho_msd), linetype = 2, color = "#F8766D") +
geom_point(data = Medium_ccmY_xmap_A_all[[i]], aes(x = Library_size, y = Ys_xmap_A_rho), shape = 1, size = 3, color = "#F8766D") +
geom_line(data = Medium_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = As_xmap_Y_rho), linetype = 1, color = "#00BFC4") +
geom_line(data = Medium_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = As_xmap_Y_rho_psd), linetype = 2, color = "#00BFC4") +
geom_line(data = Medium_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = As_xmap_Y_rho_msd), linetype = 2, color = "#00BFC4") +
geom_point(data = Medium_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = As_xmap_Y_rho), shape = 2, size = 3, color = "#00BFC4") +
ylim(-0.1, 1) +
ylab("") +
xlab("") +
ggtitle(paste(Medium_subtitle[i], "Medium", "tp =", 1 - i, sep = " "))
}
for(i in 1:5) {
print(Medium_ccm_ggplot[[i]])
}
Large_lagY_xmap_A <- list()
Large_lagA_xmap_Y <- list()
#not showing warning
Large_lagY_xmap_A <- lapply(0:-4, function(i) ccm(Large_YA_allcombined, E = Large_optE_Y, lib_column = "Large_Y", target_column = "Large_A", lib_sizes = floor(seq(Large_optE_Y, nrow(Large_YA_allcombined), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
Large_lagA_xmap_Y <- lapply(0:-4, function(i) ccm(Large_YA_allcombined, E = Large_optE_A, lib_column = "Large_A", target_column = "Large_Y", lib_sizes = floor(seq(Large_optE_A, nrow(Large_YA_allcombined), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
#example
Large_lagY_xmap_A[[2]]
Large_lagY_xmap_A_means <- list()
Large_lagA_xmap_Y_means <- list()
Large_lagY_xmap_A_means <- lapply(1:5, function(i) ccm_means(Large_lagY_xmap_A[[i]], na.rm = TRUE))
Large_lagA_xmap_Y_means <- lapply(1:5, function(i) ccm_means(Large_lagA_xmap_Y[[i]], na.rm = TRUE))
Large_lagY_xmap_A_means <- lapply(1:5, function(i) data.frame(Large_lagY_xmap_A_means[[i]], sd.rho = ccm_means(Large_lagY_xmap_A[[i]], FUN = sd, na.rm = TRUE)$rho))
Large_lagA_xmap_Y_means <- lapply(1:5, function(i) data.frame(Large_lagA_xmap_Y_means[[i]], sd.rho = ccm_means(Large_lagA_xmap_Y[[i]], FUN = sd, na.rm = TRUE)$rho))
#With the maximum library size,since there is not replicate, so SD cannot be calculated and set as zero
Large_lagY_xmap_A_means[[i]]$sd.rho[nrow(Large_lagY_xmap_A_means[[i]])] <- 0
}
for(i in 1:5){
Large_lagA_xmap_Y_means[[i]]$sd.rho[nrow(Large_lagA_xmap_Y_means[[i]])] <- 0
}
#example
–For Yotsumon (Cm)
#Generating Surrogate data
set.seed(12345)
Large_Y_each_s <- list()
Large_Y_each_s <- lapply(1:length(Large_Y_each), function(i) make_surrogate_data(Large_Y_each[[i]], method = c("ebisuzaki"), num_surr = 100))
Large_Y_each_s <- lapply(1:length(Large_Y_each_s), function(i) rbind.data.frame(Large_Y_each_s[[i]], matrix(rep(NA, 100*(no_NAs_Large)), nrow = no_NAs_Large)))
Large_Y_each_s2 <- NULL
Large_Y_each_s2 <- rbind.data.frame(Large_Y_each_s2, Large_Y_each_s[[i]])
}
Large_Y_each_s2
#Surrogate Y (Ys) xmap A (CCM)
Large_AYs <- cbind.data.frame(Large_Y_each_s2, "Large_A" = Large_YA_allcombined$Large_A)
Large_lagYs_xmap_A <- list()
for(j in 1:5) {
}
#example
Large_lagYs_xmap_A[[4]][[16]]
#mean and sd of surrogate Y xmap A
Large_lagYs_xmap_A_means_list <- list()
for(j in 1:5) {
Large_lagYs_xmap_A_means_list[[j]] <- lapply(1:length(Large_lagYs_xmap_A[[j]]), function(i) ccm_means(Large_lagYs_xmap_A[[j]][[i]], na.rm = TRUE))
}
#combining the list
Large_lagYs_xmap_A_means <- list()
for(j in 1:5) {
Large_lagYs_xmap_A_means[[j]] <- Large_lagYs_xmap_A_means_list[[j]][[1]]$rho
for(i in 2:length(Large_lagYs_xmap_A_means_list[[j]])){
Large_lagYs_xmap_A_means[[j]] <- cbind.data.frame(Large_lagYs_xmap_A_means[[j]], Large_lagYs_xmap_A_means_list[[j]][[i]]$rho)
}
}
#calculating means and sds of "means"
Large_lagYs_xmap_A_means_means <- list()
Large_lagYs_xmap_A_means_sds <- list()
for(j in 1:5) {
Large_lagYs_xmap_A_means_means[[j]] <- apply(Large_lagYs_xmap_A_means[[j]], 1, mean)
Large_lagYs_xmap_A_means_sds[[j]] <- apply(Large_lagYs_xmap_A_means[[j]], 1, sd)
}
#example
Large_lagYs_xmap_A_means_sds[[1]]
[1] 0.01427833 0.04441981 0.06011655 0.07361516 0.08296351 0.09109343
Large_lagYs_xmap_A_summary <- list()
for(j in 1:5) {
Large_lagYs_xmap_A_summary[[j]] <- data.frame("lib_size" = Large_lagYs_xmap_A_means_list[[j]][[1]]$lib_size, "mean" = Large_lagYs_xmap_A_means_means[[j]], "sd" = Large_lagYs_xmap_A_means_sds[[j]])
#example
–For Azuki (Cc)
#Generating Surrogate data
set.seed(12345)
Large_A_each_s <- list()
Large_A_each_s <- lapply(1:length(Large_A_each), function(i) make_surrogate_data(Large_A_each[[i]], method = c("ebisuzaki"), num_surr = 100))
Large_A_each_s <- lapply(1:length(Large_A_each_s), function(i) rbind.data.frame(Large_A_each_s[[i]], matrix(rep(NA,100*(no_NAs_Large)), nrow = no_NAs_Large)))
Large_A_each_s2 <- NULL
for(i in 1:length(Large_A_each_s)){
Large_A_each_s2 <- rbind.data.frame(Large_A_each_s2, Large_A_each_s[[i]])
}
Large_A_each_s2
#Surrogate A (As) xmap Y (CCM)
Large_YAs <- cbind.data.frame(Large_A_each_s2, "Large_Y" = Large_YA_allcombined$Large_Y)
Large_lagAs_xmap_Y <- list()
#example
Large_lagAs_xmap_Y[[4]][[16]]
#mean and sd of surrogate A xmap Y
Large_lagAs_xmap_Y_means_list <- list()
for(j in 1:5) {
Large_lagAs_xmap_Y_means_list[[j]] <- lapply(1:length(Large_lagAs_xmap_Y[[j]]), function(i) ccm_means(Large_lagAs_xmap_Y[[j]][[i]], na.rm = TRUE))
}
#combining the list
Large_lagAs_xmap_Y_means <- list()
for(j in 1:5) {
Large_lagAs_xmap_Y_means[[j]] <- Large_lagAs_xmap_Y_means_list[[j]][[1]]$rho
for(i in 2:length(Large_lagAs_xmap_Y_means_list[[j]])){
Large_lagAs_xmap_Y_means[[j]] <- cbind.data.frame(Large_lagAs_xmap_Y_means[[j]], Large_lagAs_xmap_Y_means_list[[j]][[i]]$rho)
}
}
#calculating means and sds of "means"
Large_lagAs_xmap_Y_means_means <- list()
Large_lagAs_xmap_Y_means_sds <- list()
for(j in 1:5) {
Large_lagAs_xmap_Y_means_means[[j]] <- apply(Large_lagAs_xmap_Y_means[[j]], 1, mean)
Large_lagAs_xmap_Y_means_sds[[j]] <- apply(Large_lagAs_xmap_Y_means[[j]], 1, sd)
}
#example
Large_lagAs_xmap_Y_means_sds[[1]]
[1] 0.02397913 0.03634388 0.04533728 0.04988835 0.05444542 0.05519462
Large_lagAs_xmap_Y_summary <- list()
for(j in 1:5) {
Large_lagAs_xmap_Y_summary[[j]] <- data.frame("lib_size" = Large_lagAs_xmap_Y_means_list[[j]][[1]]$lib_size, "mean" = Large_lagAs_xmap_Y_means_means[[j]], "sd" = Large_lagAs_xmap_Y_means_sds[[j]])
}
#example
Large_lagAs_xmap_Y_summary[[3]]
Large_ccmY_xmap_A_all <- list()
for(i in 1:5){
Large_ccmY_xmap_A_all[[i]] <- data.frame(
"Library_size" = Large_lagY_xmap_A_means[[i]]$lib_size,
"Y_xmap_A_rho" = Large_lagY_xmap_A_means[[i]]$rho,
"Ys_xmap_A_rho" = Large_lagYs_xmap_A_summary[[i]]$mean,
"Y_xmap_A_rho_psd" = Large_lagY_xmap_A_means[[i]]$rho + Large_lagY_xmap_A_means[[i]]$sd.rho,
"Y_xmap_A_rho_msd" = Large_lagY_xmap_A_means[[i]]$rho - Large_lagY_xmap_A_means[[i]]$sd.rho,
"Ys_xmap_A_rho_psd" = Large_lagYs_xmap_A_summary[[i]]$mean + Large_lagYs_xmap_A_summary[[i]]$sd,
"Ys_xmap_A_rho_msd" = Large_lagYs_xmap_A_summary[[i]]$mean - Large_lagYs_xmap_A_summary[[i]]$sd
)
}
#example
Large_ccmY_xmap_A_all[[2]]
Large_ccmA_xmap_Y_all <- list()
for(i in 1:5){
"Library_size" = Large_lagA_xmap_Y_means[[i]]$lib_size,
"A_xmap_Y_rho" = Large_lagA_xmap_Y_means[[i]]$rho,
"As_xmap_Y_rho" = Large_lagAs_xmap_Y_summary[[i]]$mean,
"A_xmap_Y_rho_psd" = Large_lagA_xmap_Y_means[[i]]$rho + Large_lagA_xmap_Y_means[[i]]$sd.rho,
"A_xmap_Y_rho_msd" = Large_lagA_xmap_Y_means[[i]]$rho - Large_lagA_xmap_Y_means[[i]]$sd.rho,
"As_xmap_Y_rho_psd" = Large_lagAs_xmap_Y_summary[[i]]$mean + Large_lagAs_xmap_Y_summary[[i]]$sd,
)
}
#example
Large_ccmA_xmap_Y_all[[2]]
#CCM Largeに関するplot
Large_ccm_ggplot <- list()
Large_subtitle <- c("(k)","(l)","(m)","(n)","(o)")
for(i in 1:5){
#Y_xmap_A:Cm xmap Cc #F8766D (red)
Large_ccm_ggplot[[i]] <- ggplot(data = Large_ccmY_xmap_A_all[[i]], aes(x = Library_size, y = Y_xmap_A_rho)) +
geom_line(linetype = 1, color = "#F8766D") +
geom_line(aes(x = Library_size, y = Y_xmap_A_rho_psd), linetype = 2, color = "#F8766D") +
geom_line(aes(x = Library_size, y = Y_xmap_A_rho_msd), linetype = 2, color = "#F8766D") +
geom_point(shape = 16, size = 3, color = "#F8766D") +
#A_xmap_Y:Cc xmap Cm #00BFC4 (blue-green)
geom_line(data = Large_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = A_xmap_Y_rho), linetype = 1, color = "#00BFC4") +
geom_line(data = Large_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = A_xmap_Y_rho_psd), linetype = 2, color = "#00BFC4") +
geom_line(data = Large_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = A_xmap_Y_rho_msd), linetype = 2, color = "#00BFC4") +
geom_point(data = Large_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = A_xmap_Y_rho), shape = 17, size = 3, color = "#00BFC4") +
geom_line(data = Large_ccmY_xmap_A_all[[i]], aes(x = Library_size, y = Ys_xmap_A_rho), linetype = 1, color = "#F8766D") +
geom_line(data = Large_ccmY_xmap_A_all[[i]], aes(x = Library_size, y = Ys_xmap_A_rho_psd), linetype = 2, color = "#F8766D") +
geom_line(data = Large_ccmY_xmap_A_all[[i]], aes(x = Library_size, y = Ys_xmap_A_rho_msd), linetype = 2, color = "#F8766D") +
geom_point(data = Large_ccmY_xmap_A_all[[i]], aes(x = Library_size, y = Ys_xmap_A_rho), shape = 1, size = 3, color = "#F8766D") +
geom_line(data = Large_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = As_xmap_Y_rho), linetype = 1, color = "#00BFC4") +
geom_line(data = Large_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = As_xmap_Y_rho_psd), linetype = 2, color = "#00BFC4") +
geom_line(data = Large_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = As_xmap_Y_rho_msd), linetype = 2, color = "#00BFC4") +
geom_point(data = Large_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = As_xmap_Y_rho), shape = 2, size = 3, color = "#00BFC4") +
ylim(-0.1, 1) +
ylab("") +
xlab("") +
ggtitle(paste(Large_subtitle[i], "Large", "tp =", 1 - i, sep = " "))
for(i in 1:5) {
print(Large_ccm_ggplot[[i]])
}
Huge_lagY_xmap_A <- list()
Huge_lagA_xmap_Y <- list()
#not showing warning
Huge_lagY_xmap_A <- lapply(0:-4, function(i) ccm(Huge_YA_allcombined, E = Huge_optE_Y, lib_column = "Huge_Y", target_column = "Huge_A", lib_sizes = floor(seq(Huge_optE_Y, nrow(Huge_YA_allcombined), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
Huge_lagA_xmap_Y <- lapply(0:-4, function(i) ccm(Huge_YA_allcombined, E = Huge_optE_A, lib_column = "Huge_A", target_column = "Huge_Y", lib_sizes = floor(seq(Huge_optE_A, nrow(Huge_YA_allcombined), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
#example
Huge_lagY_xmap_A[[2]]
Huge_lagY_xmap_A_means <- list()
Huge_lagA_xmap_Y_means <- list()
Huge_lagY_xmap_A_means <- lapply(1:5, function(i) ccm_means(Huge_lagY_xmap_A[[i]], na.rm = TRUE))
Huge_lagA_xmap_Y_means <- lapply(1:5, function(i) ccm_means(Huge_lagA_xmap_Y[[i]], na.rm = TRUE))
Huge_lagY_xmap_A_means <- lapply(1:5, function(i) data.frame(Huge_lagY_xmap_A_means[[i]], sd.rho = ccm_means(Huge_lagY_xmap_A[[i]], FUN = sd, na.rm = TRUE)$rho))
Huge_lagA_xmap_Y_means <- lapply(1:5, function(i) data.frame(Huge_lagA_xmap_Y_means[[i]], sd.rho = ccm_means(Huge_lagA_xmap_Y[[i]], FUN = sd, na.rm = TRUE)$rho))
#With the maximum library size,since there is not replicate, so SD cannot be calculated and set as zero
for(i in 1:5){
Huge_lagY_xmap_A_means[[i]]$sd.rho[nrow(Huge_lagY_xmap_A_means[[i]])] <- 0
}
for(i in 1:5){
Huge_lagA_xmap_Y_means[[i]]$sd.rho[nrow(Huge_lagA_xmap_Y_means[[i]])] <- 0
}
#example
Huge_lagA_xmap_Y_means[[3]]
–For Yotsumon (Cm)
#Generating Surrogate data
set.seed(12345)
Huge_Y_each_s <- list()
Huge_Y_each_s <- lapply(1:length(Huge_Y_each), function(i) make_surrogate_data(Huge_Y_each[[i]], method = c("ebisuzaki"), num_surr = 100))
Huge_Y_each_s <- lapply(1:length(Huge_Y_each_s), function(i) rbind.data.frame(Huge_Y_each_s[[i]], matrix(rep(NA, 100*(no_NAs_Huge)), nrow = no_NAs_Huge)))
Huge_Y_each_s2 <- NULL
for(i in 1:length(Huge_Y_each_s)){
Huge_Y_each_s2 <- rbind.data.frame(Huge_Y_each_s2, Huge_Y_each_s[[i]])
}
Huge_Y_each_s2
#Surrogate Y (Ys) xmap A (CCM)
Huge_AYs <- cbind.data.frame(Huge_Y_each_s2, "Huge_A" = Huge_YA_allcombined$Huge_A)
Huge_lagYs_xmap_A <- list()
for(j in 1:5) {
}
#example
Huge_lagYs_xmap_A[[4]][[16]]
#mean and sd of surrogate Y xmap A
Huge_lagYs_xmap_A_means_list <- list()
for(j in 1:5) {
Huge_lagYs_xmap_A_means_list[[j]] <- lapply(1:length(Huge_lagYs_xmap_A[[j]]), function(i) ccm_means(Huge_lagYs_xmap_A[[j]][[i]], na.rm = TRUE))
}
#combining the list
Huge_lagYs_xmap_A_means <- list()
for(j in 1:5) {
Huge_lagYs_xmap_A_means[[j]] <- Huge_lagYs_xmap_A_means_list[[j]][[1]]$rho
for(i in 2:length(Huge_lagYs_xmap_A_means_list[[j]])){
Huge_lagYs_xmap_A_means[[j]] <- cbind.data.frame(Huge_lagYs_xmap_A_means[[j]], Huge_lagYs_xmap_A_means_list[[j]][[i]]$rho)
}
}
#calculating means and sds of "means"
Huge_lagYs_xmap_A_means_means <- list()
Huge_lagYs_xmap_A_means_sds <- list()
for(j in 1:5) {
Huge_lagYs_xmap_A_means_means[[j]] <- apply(Huge_lagYs_xmap_A_means[[j]], 1, mean)
Huge_lagYs_xmap_A_means_sds[[j]] <- apply(Huge_lagYs_xmap_A_means[[j]], 1, sd)
}
#example
Huge_lagYs_xmap_A_means_sds[[1]]
0.02726424 0.05057614 0.05527516 0.05871576 0.06292198 0.06850279
Huge_lagYs_xmap_A_summary <- list()
for(j in 1:5) {
Huge_lagYs_xmap_A_summary[[j]] <- data.frame("lib_size" = Huge_lagYs_xmap_A_means_list[[j]][[1]]$lib_size, "mean" = Huge_lagYs_xmap_A_means_means[[j]], "sd" = Huge_lagYs_xmap_A_means_sds[[j]])
}
#example
Huge_lagYs_xmap_A_summary[[3]]
–For Azuki (Cc)
#Generating Surrogate data
set.seed(12345)
Huge_A_each_s <- list()
Huge_A_each_s <- lapply(1:length(Huge_A_each), function(i) make_surrogate_data(Huge_A_each[[i]], method = c("ebisuzaki"), num_surr = 100))
Huge_A_each_s <- lapply(1:length(Huge_A_each_s), function(i) rbind.data.frame(Huge_A_each_s[[i]], matrix(rep(NA,100*(no_NAs_Huge)), nrow = no_NAs_Huge)))
Huge_A_each_s2 <- NULL
for(i in 1:length(Huge_A_each_s)){
}
Huge_A_each_s2
#Surrogate A (As) xmap Y (CCM)
Huge_YAs <- cbind.data.frame(Huge_A_each_s2, "Huge_Y" = Huge_YA_allcombined$Huge_Y)
Huge_lagAs_xmap_Y <- list()
for(j in 1:5) {
}
#example
Huge_lagAs_xmap_Y[[4]][[16]]
#mean and sd of surrogate A xmap Y
Huge_lagAs_xmap_Y_means_list <- list()
for(j in 1:5) {
Huge_lagAs_xmap_Y_means_list[[j]] <- lapply(1:length(Huge_lagAs_xmap_Y[[j]]), function(i) ccm_means(Huge_lagAs_xmap_Y[[j]][[i]], na.rm = TRUE))
}
#combining the list
Huge_lagAs_xmap_Y_means <- list()
for(j in 1:5) {
Huge_lagAs_xmap_Y_means[[j]] <- Huge_lagAs_xmap_Y_means_list[[j]][[1]]$rho
for(i in 2:length(Huge_lagAs_xmap_Y_means_list[[j]])){
Huge_lagAs_xmap_Y_means[[j]] <- cbind.data.frame(Huge_lagAs_xmap_Y_means[[j]], Huge_lagAs_xmap_Y_means_list[[j]][[i]]$rho)
}
}
#calculating means and sds of "means"
Huge_lagAs_xmap_Y_means_means <- list()
Huge_lagAs_xmap_Y_means_sds <- list()
for(j in 1:5) {
Huge_lagAs_xmap_Y_means_means[[j]] <- apply(Huge_lagAs_xmap_Y_means[[j]], 1, mean)
Huge_lagAs_xmap_Y_means_sds[[j]] <- apply(Huge_lagAs_xmap_Y_means[[j]], 1, sd)
}
#example
[1] 0.03000238 0.04001159 0.04301337 0.04740172 0.05334495 0.05620980
Huge_lagAs_xmap_Y_summary <- list()
for(j in 1:5) {
Huge_lagAs_xmap_Y_summary[[j]] <- data.frame("lib_size" = Huge_lagAs_xmap_Y_means_list[[j]][[1]]$lib_size, "mean" = Huge_lagAs_xmap_Y_means_means[[j]], "sd" = Huge_lagAs_xmap_Y_means_sds[[j]])
}
#example
Huge_ccmY_xmap_A_all <- list()
for(i in 1:5){
Huge_ccmY_xmap_A_all[[i]] <- data.frame(
"Library_size" = Huge_lagY_xmap_A_means[[i]]$lib_size,
"Y_xmap_A_rho" = Huge_lagY_xmap_A_means[[i]]$rho,
"Ys_xmap_A_rho" = Huge_lagYs_xmap_A_summary[[i]]$mean,
"Y_xmap_A_rho_psd" = Huge_lagY_xmap_A_means[[i]]$rho + Huge_lagY_xmap_A_means[[i]]$sd.rho,
"Y_xmap_A_rho_msd" = Huge_lagY_xmap_A_means[[i]]$rho - Huge_lagY_xmap_A_means[[i]]$sd.rho,
"Ys_xmap_A_rho_psd" = Huge_lagYs_xmap_A_summary[[i]]$mean + Huge_lagYs_xmap_A_summary[[i]]$sd,
"Ys_xmap_A_rho_msd" = Huge_lagYs_xmap_A_summary[[i]]$mean - Huge_lagYs_xmap_A_summary[[i]]$sd
)
}
#example
Huge_ccmY_xmap_A_all[[2]]
Huge_ccmA_xmap_Y_all <- list()
for(i in 1:5){
Huge_ccmA_xmap_Y_all[[i]] <- data.frame(
"Library_size" = Huge_lagA_xmap_Y_means[[i]]$lib_size,
"A_xmap_Y_rho" = Huge_lagA_xmap_Y_means[[i]]$rho,
"As_xmap_Y_rho" = Huge_lagAs_xmap_Y_summary[[i]]$mean,
"A_xmap_Y_rho_psd" = Huge_lagA_xmap_Y_means[[i]]$rho + Huge_lagA_xmap_Y_means[[i]]$sd.rho,
"A_xmap_Y_rho_msd" = Huge_lagA_xmap_Y_means[[i]]$rho - Huge_lagA_xmap_Y_means[[i]]$sd.rho,
"As_xmap_Y_rho_psd" = Huge_lagAs_xmap_Y_summary[[i]]$mean + Huge_lagAs_xmap_Y_summary[[i]]$sd,
"As_xmap_Y_rho_msd" = Huge_lagAs_xmap_Y_summary[[i]]$mean - Huge_lagAs_xmap_Y_summary[[i]]$sd
)
}
Huge_ccmA_xmap_Y_all[[2]]
#CCM Hugeに関するplot
Huge_ccm_ggplot <- list()
Huge_subtitle <- c("(p)","(q)","(r)","(s)","(t)")
for(i in 1:5){
#Y_xmap_A:Cm xmap Cc #F8766D (red)
Huge_ccm_ggplot[[i]] <- ggplot(data = Huge_ccmY_xmap_A_all[[i]], aes(x = Library_size, y = Y_xmap_A_rho)) +
geom_line(linetype = 1, color = "#F8766D") +
geom_line(aes(x = Library_size, y = Y_xmap_A_rho_psd), linetype = 2, color = "#F8766D") +
geom_line(aes(x = Library_size, y = Y_xmap_A_rho_msd), linetype = 2, color = "#F8766D") +
geom_point(shape = 16, size = 3, color = "#F8766D") +
#A_xmap_Y:Cc xmap Cm #00BFC4 (blue-green)
geom_line(data = Huge_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = A_xmap_Y_rho), linetype = 1, color = "#00BFC4") +
geom_line(data = Huge_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = A_xmap_Y_rho_psd), linetype = 2, color = "#00BFC4") +
geom_line(data = Huge_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = A_xmap_Y_rho_msd), linetype = 2, color = "#00BFC4") +
geom_point(data = Huge_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = A_xmap_Y_rho), shape = 17, size = 3, color = "#00BFC4") +
geom_line(data = Huge_ccmY_xmap_A_all[[i]], aes(x = Library_size, y = Ys_xmap_A_rho), linetype = 1, color = "#F8766D") +
geom_line(data = Huge_ccmY_xmap_A_all[[i]], aes(x = Library_size, y = Ys_xmap_A_rho_psd), linetype = 2, color = "#F8766D") +
geom_line(data = Huge_ccmY_xmap_A_all[[i]], aes(x = Library_size, y = Ys_xmap_A_rho_msd), linetype = 2, color = "#F8766D") +
geom_point(data = Huge_ccmY_xmap_A_all[[i]], aes(x = Library_size, y = Ys_xmap_A_rho), shape = 1, size = 3, color = "#F8766D") +
geom_line(data = Huge_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = As_xmap_Y_rho), linetype = 1, color = "#00BFC4") +
geom_line(data = Huge_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = As_xmap_Y_rho_psd), linetype = 2, color = "#00BFC4") +
geom_line(data = Huge_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = As_xmap_Y_rho_msd), linetype = 2, color = "#00BFC4") +
geom_point(data = Huge_ccmA_xmap_Y_all[[i]], aes(x = Library_size, y = As_xmap_Y_rho), shape = 2, size = 3, color = "#00BFC4") +
ylim(-0.1, 1) +
ylab("") +
xlab("") +
ggtitle(paste(Huge_subtitle[i], "Huge", "tp =", 1 - i, sep = " "))
}
for(i in 1:5) {
print(Huge_ccm_ggplot[[i]])
}
#When either C.m or C.c showed linear dynamics, we did not apply CCM to them
subset(Table1, (theta_for_C.m == 0 | theta_for_C.c == 0))
subset(Table1, (theta_for_C.m == 0 | theta_for_C.c == 0))$ID
[1] "S3" "S4" "M1" "L2" "L3" "H4" "H5"
Sk <- c(1, 2, 5, 6) #index for sample ID used for CCM (S1, S2, S5, & S6)
Small_ccmdata_combined <- list()
Small_ccmY_xmap_A_combined <- list()
Small_ccmA_xmap_Y_combined <- list()
Small_ccmY_xmap_Y_combined <- list()
Small_ccmA_xmap_A_combined <- list()
for(j in 1:6) Small_ccmdata_combined[[j]] <- data.frame("Y" = S_combined_Y[[j]], "A" = S_combined_A[[j]])
#Cm (Y) xmap Cc (A)
for(k in 1:4) {
Small_ccmY_xmap_A_combined[[k]] <- lapply(0:-4, function(i) ccm(Small_ccmdata_combined[[Sk[k]]], E = Table1$E_for_C.m[Sk[k]], lib_column = "Y", target_column = "A", lib_sizes = floor(seq(Table1$E_for_C.m[Sk[k]], length(Small_ccmdata_combined[[Sk[k]]]$Y), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
}
#Cc (A) xmap Cm (Y)
for(k in 1:4) {
Small_ccmA_xmap_Y_combined[[k]] <- lapply(0:-4, function(i) ccm(Small_ccmdata_combined[[Sk[k]]], E = Table1$E_for_C.c[Sk[k]], lib_column = "A", target_column = "Y", lib_sizes = floor(seq(Table1$E_for_C.c[Sk[k]], length(Small_ccmdata_combined[[Sk[k]]]$A), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
}
#example
Small_ccmA_xmap_Y_combined[[3]][[1]]
#Cm (Y) xmap Cm (Y) for intraspecific density dependence
for(k in 1:4) {
Small_ccmY_xmap_Y_combined[[k]] <- lapply(0:-4, function(i) ccm(Small_ccmdata_combined[[Sk[k]]], E = Table1$E_for_C.m[Sk[k]], lib_column = "Y", target_column = "Y", lib_sizes = floor(seq(Table1$E_for_C.m[Sk[k]], length(Small_ccmdata_combined[[Sk[k]]]$Y), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
}
#Cc (A) xmap Cc (A)
Small_ccmA_xmap_A_combined[[k]] <- lapply(0:-4, function(i) ccm(Small_ccmdata_combined[[Sk[k]]], E = Table1$E_for_C.c[Sk[k]], lib_column = "A", target_column = "A", lib_sizes = floor(seq(Table1$E_for_C.c[Sk[k]], length(Small_ccmdata_combined[[Sk[k]]]$A), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
}
#example
Small_ccmA_xmap_A_combined[[3]][[1]]
Mk <- c(2, 3, 4, 5, 6, 7) #index for sample ID used for CCM (M2 to M7)
Medium_ccmdata_combined <- list()
Medium_ccmY_xmap_A_combined <- list()
Medium_ccmA_xmap_Y_combined <- list()
Medium_ccmY_xmap_Y_combined <- list()
Medium_ccmA_xmap_A_combined <- list()
for(j in 1:7) Medium_ccmdata_combined[[j]] <- data.frame("Y" = M_combined_Y[[j]], "A" = M_combined_A[[j]])
#Cm (Y) xmap Cc (A)
for(k in 1:6) {
Medium_ccmY_xmap_A_combined[[k]] <- lapply(0:-4, function(i) ccm(Medium_ccmdata_combined[[Mk[k]]], E = Table1$E_for_C.m[Mk[k]], lib_column = "Y", target_column = "A", lib_sizes = floor(seq(Table1$E_for_C.m[Mk[k]], length(Medium_ccmdata_combined[[Mk[k]]]$Y), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
}
#Cc (A) xmap Cm (Y)
Medium_ccmA_xmap_Y_combined[[k]] <- lapply(0:-4, function(i) ccm(Medium_ccmdata_combined[[Mk[k]]], E = Table1$E_for_C.c[Mk[k]], lib_column = "A", target_column = "Y", lib_sizes = floor(seq(Table1$E_for_C.c[Mk[k]], length(Medium_ccmdata_combined[[Mk[k]]]$A), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
}
#example
Medium_ccmA_xmap_Y_combined[[3]][[1]]
#Cm (Y) xmap Cm (Y) for intraspecific density dependence
for(k in 1:6) {
Medium_ccmY_xmap_Y_combined[[k]] <- lapply(0:-4, function(i) ccm(Medium_ccmdata_combined[[Mk[k]]], E = Table1$E_for_C.m[Mk[k]], lib_column = "Y", target_column = "Y", lib_sizes = floor(seq(Table1$E_for_C.m[Mk[k]], length(Medium_ccmdata_combined[[Mk[k]]]$Y), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
}
#Cc (A) xmap Cc (A)
for(k in 1:6) {
Medium_ccmA_xmap_A_combined[[k]] <- lapply(0:-4, function(i) ccm(Medium_ccmdata_combined[[Mk[k]]], E = Table1$E_for_C.c[Mk[k]], lib_column = "A", target_column = "A", lib_sizes = floor(seq(Table1$E_for_C.c[Mk[k]], length(Medium_ccmdata_combined[[Mk[k]]]$A), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
}
#example
Medium_ccmA_xmap_A_combined[[3]][[1]]
Lk <- c(1, 4, 5, 6, 7, 8) #index for sample ID used for CCM (L1, L4 to L8)
Large_ccmdata_combined <- list()
Large_ccmY_xmap_A_combined <- list()
Large_ccmA_xmap_Y_combined <- list()
Large_ccmY_xmap_Y_combined <- list()
Large_ccmA_xmap_A_combined <- list()
for(j in 1:8) Large_ccmdata_combined[[j]] <- data.frame("Y" = L_combined_Y[[j]], "A" = L_combined_A[[j]])
#Cm (Y) xmap Cc (A)
for(k in 1:6) {
Large_ccmY_xmap_A_combined[[k]] <- lapply(0:-4, function(i) ccm(Large_ccmdata_combined[[Lk[k]]], E = Table1$E_for_C.m[Lk[k]], lib_column = "Y", target_column = "A", lib_sizes = floor(seq(Table1$E_for_C.m[Lk[k]], length(Large_ccmdata_combined[[Lk[k]]]$Y), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
}
#Cc (A) xmap Cm (Y)
Large_ccmA_xmap_Y_combined[[k]] <- lapply(0:-4, function(i) ccm(Large_ccmdata_combined[[Lk[k]]], E = Table1$E_for_C.c[Lk[k]], lib_column = "A", target_column = "Y", lib_sizes = floor(seq(Table1$E_for_C.c[Lk[k]], length(Large_ccmdata_combined[[Lk[k]]]$A), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
}
#example
Large_ccmA_xmap_Y_combined[[3]][[1]]
#Cm (Y) xmap Cm (Y) for intraspecific density dependence
for(k in 1:6) {
Large_ccmY_xmap_Y_combined[[k]] <- lapply(0:-4, function(i) ccm(Large_ccmdata_combined[[Lk[k]]], E = Table1$E_for_C.m[Lk[k]], lib_column = "Y", target_column = "Y", lib_sizes = floor(seq(Table1$E_for_C.m[Lk[k]], length(Large_ccmdata_combined[[Lk[k]]]$Y), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
}
#Cc (A) xmap Cc (A)
for(k in 1:6) {
Large_ccmA_xmap_A_combined[[k]] <- lapply(0:-4, function(i) ccm(Large_ccmdata_combined[[Lk[k]]], E = Table1$E_for_C.c[Lk[k]], lib_column = "A", target_column = "A", lib_sizes = floor(seq(Table1$E_for_C.c[Lk[k]], length(Large_ccmdata_combined[[Lk[k]]]$A), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
}
#example
Large_ccmA_xmap_A_combined[[3]][[1]]
Hk <- c(1, 2, 3, 6, 7, 8) #index for sample ID used for CCM (H1, H2, H3, H6, H7, H8)
Huge_ccmdata_combined <- list()
Huge_ccmY_xmap_A_combined <- list()
Huge_ccmA_xmap_Y_combined <- list()
Huge_ccmY_xmap_Y_combined <- list()
Huge_ccmA_xmap_A_combined <- list()
for(j in 1:8) Huge_ccmdata_combined[[j]] <- data.frame("Y" = H_combined_Y[[j]], "A" = H_combined_A[[j]])
for(k in 1:6) {
Huge_ccmY_xmap_A_combined[[k]] <- lapply(0:-4, function(i) ccm(Huge_ccmdata_combined[[Hk[k]]], E = Table1$E_for_C.m[Hk[k]], lib_column = "Y", target_column = "A", lib_sizes = floor(seq(Table1$E_for_C.m[Hk[k]], length(Huge_ccmdata_combined[[Hk[k]]]$Y), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
}
#Cc (A) xmap Cm (Y)
for(k in 1:6) {
Huge_ccmA_xmap_Y_combined[[k]] <- lapply(0:-4, function(i) ccm(Huge_ccmdata_combined[[Hk[k]]], E = Table1$E_for_C.c[Hk[k]], lib_column = "A", target_column = "Y", lib_sizes = floor(seq(Table1$E_for_C.c[Hk[k]], length(Huge_ccmdata_combined[[Hk[k]]]$A), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
}
#example
Huge_ccmA_xmap_Y_combined[[3]][[1]]
#Cm (Y) xmap Cm (Y) for intraspecific density dependence
for(k in 1:6) {
Huge_ccmY_xmap_Y_combined[[k]] <- lapply(0:-4, function(i) ccm(Huge_ccmdata_combined[[Hk[k]]], E = Table1$E_for_C.m[Hk[k]], lib_column = "Y", target_column = "Y", lib_sizes = floor(seq(Table1$E_for_C.m[Hk[k]], length(Huge_ccmdata_combined[[Hk[k]]]$Y), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
}
#Cc (A) xmap Cc (A)
for(k in 1:6) {
Huge_ccmA_xmap_A_combined[[k]] <- lapply(0:-4, function(i) ccm(Huge_ccmdata_combined[[Hk[k]]], E = Table1$E_for_C.c[Hk[k]], lib_column = "A", target_column = "A", lib_sizes = floor(seq(Table1$E_for_C.c[Hk[k]], length(Huge_ccmdata_combined[[Hk[k]]]$A), length = 10)), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
}
#example
Huge_ccmA_xmap_A_combined[[3]][[1]]
saveRDS(Small_ccmY_xmap_A_combined, "Small_ccmY_xmap_A_combined.obj")
saveRDS(Small_ccmY_xmap_Y_combined, "Small_ccmY_xmap_Y_combined.obj")
saveRDS(Small_ccmA_xmap_Y_combined, "Small_ccmA_xmap_Y_combined.obj")
saveRDS(Small_ccmA_xmap_A_combined, "Small_ccmA_xmap_A_combined.obj")
saveRDS(Medium_ccmY_xmap_A_combined, "Medium_ccmY_xmap_A_combined.obj")
saveRDS(Medium_ccmY_xmap_Y_combined, "Medium_ccmY_xmap_Y_combined.obj")
saveRDS(Medium_ccmA_xmap_Y_combined, "Medium_ccmA_xmap_Y_combined.obj")
saveRDS(Medium_ccmA_xmap_A_combined, "Medium_ccmA_xmap_A_combined.obj")
saveRDS(Large_ccmY_xmap_A_combined, "Large_ccmY_xmap_A_combined.obj")
saveRDS(Large_ccmY_xmap_Y_combined, "Large_ccmY_xmap_Y_combined.obj")
saveRDS(Large_ccmA_xmap_Y_combined, "Large_ccmA_xmap_Y_combined.obj")
saveRDS(Large_ccmA_xmap_A_combined, "Large_ccmA_xmap_A_combined.obj")
saveRDS(Huge_ccmY_xmap_A_combined, "Huge_ccmY_xmap_A_combined.obj")
saveRDS(Huge_ccmY_xmap_Y_combined, "Huge_ccmY_xmap_Y_combined.obj")
saveRDS(Huge_ccmA_xmap_Y_combined, "Huge_ccmA_xmap_Y_combined.obj")
saveRDS(Huge_ccmA_xmap_A_combined, "Huge_ccmA_xmap_A_combined.obj")
saveRDS(Table1, "Table1.obj")
#Cage volume
Smallvolume <- 8.5*6*3.5
Mediumvolume <- 13.5*8.5*4.5
Largevolume <- 18.5*14*6.5
Hugevolume <- 22*15.5*7.5
#Table, which include the coexistence week information
Table1 <- readRDS("Table1.obj")
#CCM results for combined time series >= 35 points
Small_ccmY_xmap_A_combined <- readRDS("Small_ccmY_xmap_A_combined.obj")
Small_ccmY_xmap_Y_combined <- readRDS("Small_ccmY_xmap_Y_combined.obj")
Small_ccmA_xmap_Y_combined <- readRDS("Small_ccmA_xmap_Y_combined.obj")
Small_ccmA_xmap_A_combined <- readRDS("Small_ccmA_xmap_A_combined.obj")
Medium_ccmY_xmap_A_combined <- readRDS("Medium_ccmY_xmap_A_combined.obj")
Medium_ccmA_xmap_Y_combined <- readRDS("Medium_ccmA_xmap_Y_combined.obj")
Medium_ccmA_xmap_A_combined <- readRDS("Medium_ccmA_xmap_A_combined.obj")
Large_ccmY_xmap_A_combined <- readRDS("Large_ccmY_xmap_A_combined.obj")
Large_ccmY_xmap_Y_combined <- readRDS("Large_ccmY_xmap_Y_combined.obj")
Large_ccmA_xmap_Y_combined <- readRDS("Large_ccmA_xmap_Y_combined.obj")
Large_ccmA_xmap_A_combined <- readRDS("Large_ccmA_xmap_A_combined.obj")
Huge_ccmY_xmap_A_combined <- readRDS("Huge_ccmY_xmap_A_combined.obj")
Huge_ccmY_xmap_Y_combined <- readRDS("Huge_ccmY_xmap_Y_combined.obj")
Huge_ccmA_xmap_Y_combined <- readRDS("Huge_ccmA_xmap_Y_combined.obj")
Huge_ccmA_xmap_A_combined <- readRDS("Huge_ccmA_xmap_A_combined.obj")
– There are multiple possibilities to calculate the index of interspecific interaction strength. We finally decided to use the rho value of interspecific CCM with the maximum library length standardized by using the interspecific CCM rho with the minimum library length, and the intraspecific CCM rho with the maximum and minimum library lengths.
– Note that SIZE_ccmU_xmap_V_combined[[k]][[j]] represents the result from the combined time series ID “SIZEk” with delay tp (1 - j). (j from 1 to 5)
Small_Y_xmap_A_interaction <- list()
Small_Y_xmap_Y_interaction <- list()
Small_A_xmap_Y_interaction <- list()
Small_A_xmap_A_interaction <- list()
for(k in 1:length(Small_ccmY_xmap_A_combined)){
Small_Y_xmap_A_interaction[[k]] <- list()
Small_Y_xmap_Y_interaction[[k]] <- list()
Small_A_xmap_Y_interaction[[k]] <- list()
Small_A_xmap_A_interaction[[k]] <- list()
for(j in 1:5) {
Small_Y_xmap_A_interaction[[k]][[j]] <- subset(Small_ccmY_xmap_A_combined[[k]][[j]]$rho, Small_ccmY_xmap_A_combined[[k]][[j]]$lib_size == max(Small_ccmY_xmap_A_combined[[k]][[j]]$lib_size)) - subset(Small_ccmY_xmap_A_combined[[k]][[j]]$rho, Small_ccmY_xmap_A_combined[[k]][[j]]$lib_size == min(Small_ccmY_xmap_A_combined[[k]][[j]]$lib_size))
Small_Y_xmap_Y_interaction[[k]][[j]] <- subset(Small_ccmY_xmap_Y_combined[[k]][[j]]$rho, Small_ccmY_xmap_Y_combined[[k]][[j]]$lib_size == max(Small_ccmY_xmap_Y_combined[[k]][[j]]$lib_size)) - subset(Small_ccmY_xmap_Y_combined[[k]][[j]]$rho, Small_ccmY_xmap_Y_combined[[k]][[j]]$lib_size == min(Small_ccmY_xmap_Y_combined[[k]][[j]]$lib_size))
Small_A_xmap_Y_interaction[[k]][[j]] <- subset(Small_ccmA_xmap_Y_combined[[k]][[j]]$rho, Small_ccmA_xmap_Y_combined[[k]][[j]]$lib_size == max(Small_ccmA_xmap_Y_combined[[k]][[j]]$lib_size)) - subset(Small_ccmA_xmap_Y_combined[[k]][[j]]$rho, Small_ccmA_xmap_Y_combined[[k]][[j]]$lib_size == min(Small_ccmA_xmap_Y_combined[[k]][[j]]$lib_size))
Small_A_xmap_A_interaction[[k]][[j]] <- subset(Small_ccmA_xmap_A_combined[[k]][[j]]$rho, Small_ccmA_xmap_A_combined[[k]][[j]]$lib_size == max(Small_ccmA_xmap_A_combined[[k]][[j]]$lib_size)) - subset(Small_ccmA_xmap_A_combined[[k]][[j]]$rho, Small_ccmA_xmap_A_combined[[k]][[j]]$lib_size == min(Small_ccmA_xmap_A_combined[[k]][[j]]$lib_size))
#convert negative values into zero
Small_Y_xmap_A_interaction[[k]][[j]][Small_Y_xmap_A_interaction[[k]][[j]] < 0] <- 0
Small_Y_xmap_Y_interaction[[k]][[j]][Small_Y_xmap_Y_interaction[[k]][[j]] < 0] <- 0
Small_A_xmap_Y_interaction[[k]][[j]][Small_A_xmap_Y_interaction[[k]][[j]] < 0] <- 0
Small_A_xmap_A_interaction[[k]][[j]][Small_A_xmap_A_interaction[[k]][[j]] < 0] <- 0
}
}
Medium_Y_xmap_A_interaction <- list()
Medium_Y_xmap_Y_interaction <- list()
Medium_A_xmap_Y_interaction <- list()
Medium_A_xmap_A_interaction <- list()
for(k in 1:length(Medium_ccmY_xmap_A_combined)){
Medium_Y_xmap_A_interaction[[k]] <- list()
Medium_Y_xmap_Y_interaction[[k]] <- list()
Medium_A_xmap_Y_interaction[[k]] <- list()
Medium_A_xmap_A_interaction[[k]] <- list()
for(j in 1:5) {
Medium_Y_xmap_A_interaction[[k]][[j]] <- subset(Medium_ccmY_xmap_A_combined[[k]][[j]]$rho, Medium_ccmY_xmap_A_combined[[k]][[j]]$lib_size == max(Medium_ccmY_xmap_A_combined[[k]][[j]]$lib_size)) - subset(Medium_ccmY_xmap_A_combined[[k]][[j]]$rho, Medium_ccmY_xmap_A_combined[[k]][[j]]$lib_size == min(Medium_ccmY_xmap_A_combined[[k]][[j]]$lib_size))
Medium_Y_xmap_Y_interaction[[k]][[j]] <- subset(Medium_ccmY_xmap_Y_combined[[k]][[j]]$rho, Medium_ccmY_xmap_Y_combined[[k]][[j]]$lib_size == max(Medium_ccmY_xmap_Y_combined[[k]][[j]]$lib_size)) - subset(Medium_ccmY_xmap_Y_combined[[k]][[j]]$rho, Medium_ccmY_xmap_Y_combined[[k]][[j]]$lib_size == min(Medium_ccmY_xmap_Y_combined[[k]][[j]]$lib_size))
Medium_A_xmap_Y_interaction[[k]][[j]] <- subset(Medium_ccmA_xmap_Y_combined[[k]][[j]]$rho, Medium_ccmA_xmap_Y_combined[[k]][[j]]$lib_size == max(Medium_ccmA_xmap_Y_combined[[k]][[j]]$lib_size)) - subset(Medium_ccmA_xmap_Y_combined[[k]][[j]]$rho, Medium_ccmA_xmap_Y_combined[[k]][[j]]$lib_size == min(Medium_ccmA_xmap_Y_combined[[k]][[j]]$lib_size))
Medium_A_xmap_A_interaction[[k]][[j]] <- subset(Medium_ccmA_xmap_A_combined[[k]][[j]]$rho, Medium_ccmA_xmap_A_combined[[k]][[j]]$lib_size == max(Medium_ccmA_xmap_A_combined[[k]][[j]]$lib_size)) - subset(Medium_ccmA_xmap_A_combined[[k]][[j]]$rho, Medium_ccmA_xmap_A_combined[[k]][[j]]$lib_size == min(Medium_ccmA_xmap_A_combined[[k]][[j]]$lib_size))
#convert negative values into zero
Medium_Y_xmap_A_interaction[[k]][[j]][Medium_Y_xmap_A_interaction[[k]][[j]] < 0] <- 0
Medium_Y_xmap_Y_interaction[[k]][[j]][Medium_Y_xmap_Y_interaction[[k]][[j]] < 0] <- 0
Medium_A_xmap_Y_interaction[[k]][[j]][Medium_A_xmap_Y_interaction[[k]][[j]] < 0] <- 0
Medium_A_xmap_A_interaction[[k]][[j]][Medium_A_xmap_A_interaction[[k]][[j]] < 0] <- 0
}
}
Large_Y_xmap_A_interaction <- list()
Large_Y_xmap_Y_interaction <- list()
Large_A_xmap_Y_interaction <- list()
Large_A_xmap_A_interaction <- list()
for(k in 1:length(Large_ccmY_xmap_A_combined)){
Large_Y_xmap_A_interaction[[k]] <- list()
Large_Y_xmap_Y_interaction[[k]] <- list()
Large_A_xmap_Y_interaction[[k]] <- list()
Large_A_xmap_A_interaction[[k]] <- list()
for(j in 1:5) {
Large_Y_xmap_A_interaction[[k]][[j]] <- subset(Large_ccmY_xmap_A_combined[[k]][[j]]$rho, Large_ccmY_xmap_A_combined[[k]][[j]]$lib_size == max(Large_ccmY_xmap_A_combined[[k]][[j]]$lib_size)) - subset(Large_ccmY_xmap_A_combined[[k]][[j]]$rho, Large_ccmY_xmap_A_combined[[k]][[j]]$lib_size == min(Large_ccmY_xmap_A_combined[[k]][[j]]$lib_size))
Large_Y_xmap_Y_interaction[[k]][[j]] <- subset(Large_ccmY_xmap_Y_combined[[k]][[j]]$rho, Large_ccmY_xmap_Y_combined[[k]][[j]]$lib_size == max(Large_ccmY_xmap_Y_combined[[k]][[j]]$lib_size)) - subset(Large_ccmY_xmap_Y_combined[[k]][[j]]$rho, Large_ccmY_xmap_Y_combined[[k]][[j]]$lib_size == min(Large_ccmY_xmap_Y_combined[[k]][[j]]$lib_size))
Large_A_xmap_Y_interaction[[k]][[j]] <- subset(Large_ccmA_xmap_Y_combined[[k]][[j]]$rho, Large_ccmA_xmap_Y_combined[[k]][[j]]$lib_size == max(Large_ccmA_xmap_Y_combined[[k]][[j]]$lib_size)) - subset(Large_ccmA_xmap_Y_combined[[k]][[j]]$rho, Large_ccmA_xmap_Y_combined[[k]][[j]]$lib_size == min(Large_ccmA_xmap_Y_combined[[k]][[j]]$lib_size))
}
}
Huge_Y_xmap_A_interaction <- list()
Huge_Y_xmap_Y_interaction <- list()
Huge_A_xmap_Y_interaction <- list()
Huge_A_xmap_A_interaction <- list()
for(k in 1:length(Huge_ccmY_xmap_A_combined)){
Huge_Y_xmap_A_interaction[[k]] <- list()
Huge_Y_xmap_Y_interaction[[k]] <- list()
Huge_A_xmap_Y_interaction[[k]] <- list()
Huge_A_xmap_A_interaction[[k]] <- list()
for(j in 1:5) {
Huge_Y_xmap_A_interaction[[k]][[j]] <- subset(Huge_ccmY_xmap_A_combined[[k]][[j]]$rho, Huge_ccmY_xmap_A_combined[[k]][[j]]$lib_size == max(Huge_ccmY_xmap_A_combined[[k]][[j]]$lib_size)) - subset(Huge_ccmY_xmap_A_combined[[k]][[j]]$rho, Huge_ccmY_xmap_A_combined[[k]][[j]]$lib_size == min(Huge_ccmY_xmap_A_combined[[k]][[j]]$lib_size))
Huge_Y_xmap_Y_interaction[[k]][[j]] <- subset(Huge_ccmY_xmap_Y_combined[[k]][[j]]$rho, Huge_ccmY_xmap_Y_combined[[k]][[j]]$lib_size == max(Huge_ccmY_xmap_Y_combined[[k]][[j]]$lib_size)) - subset(Huge_ccmY_xmap_Y_combined[[k]][[j]]$rho, Huge_ccmY_xmap_Y_combined[[k]][[j]]$lib_size == min(Huge_ccmY_xmap_Y_combined[[k]][[j]]$lib_size))
Huge_A_xmap_Y_interaction[[k]][[j]] <- subset(Huge_ccmA_xmap_Y_combined[[k]][[j]]$rho, Huge_ccmA_xmap_Y_combined[[k]][[j]]$lib_size == max(Huge_ccmA_xmap_Y_combined[[k]][[j]]$lib_size)) - subset(Huge_ccmA_xmap_Y_combined[[k]][[j]]$rho, Huge_ccmA_xmap_Y_combined[[k]][[j]]$lib_size == min(Huge_ccmA_xmap_Y_combined[[k]][[j]]$lib_size))
Huge_A_xmap_A_interaction[[k]][[j]] <- subset(Huge_ccmA_xmap_A_combined[[k]][[j]]$rho, Huge_ccmA_xmap_A_combined[[k]][[j]]$lib_size == max(Huge_ccmA_xmap_A_combined[[k]][[j]]$lib_size)) - subset(Huge_ccmA_xmap_A_combined[[k]][[j]]$rho, Huge_ccmA_xmap_A_combined[[k]][[j]]$lib_size == min(Huge_ccmA_xmap_A_combined[[k]][[j]]$lib_size))
#convert negative values into zero
Huge_Y_xmap_Y_interaction[[k]][[j]][Huge_Y_xmap_Y_interaction[[k]][[j]] < 0] <- 0
Huge_A_xmap_A_interaction[[k]][[j]][Huge_A_xmap_A_interaction[[k]][[j]] < 0] <- 0
}
}
– Note the direction of CCM and that of interaction are opposite. – Relative interaction from U to V is V_xmap_U_interaction/U_xmap_U_interaction ##### For Small
Small_fromYtoA_interaction <- list()
Small_fromAtoY_interaction <- list()
for(k in 1:length(Small_ccmY_xmap_A_combined)){
Small_fromYtoA_interaction[[k]] <- list()
Small_fromAtoY_interaction[[k]] <- list()
for(j in 1:5) {
Small_fromYtoA_interaction[[k]][[j]] <- Small_A_xmap_Y_interaction[[k]][[j]]/Small_Y_xmap_Y_interaction[[k]][[j]]
Small_fromAtoY_interaction[[k]][[j]] <- Small_Y_xmap_A_interaction[[k]][[j]]/Small_A_xmap_A_interaction[[k]][[j]]
}
}
Medium_fromYtoA_interaction <- list()
Medium_fromAtoY_interaction <- list()
for(k in 1:length(Medium_ccmY_xmap_A_combined)){
Medium_fromYtoA_interaction[[k]] <- list()
Medium_fromAtoY_interaction[[k]] <- list()
for(j in 1:5) {
Medium_fromYtoA_interaction[[k]][[j]] <- Medium_A_xmap_Y_interaction[[k]][[j]]/Medium_Y_xmap_Y_interaction[[k]][[j]]
Medium_fromAtoY_interaction[[k]][[j]] <- Medium_Y_xmap_A_interaction[[k]][[j]]/Medium_A_xmap_A_interaction[[k]][[j]]
}
}
Large_fromYtoA_interaction <- list()
Large_fromAtoY_interaction <- list()
for(k in 1:length(Large_ccmY_xmap_A_combined)){
Large_fromYtoA_interaction[[k]] <- list()
Large_fromAtoY_interaction[[k]] <- list()
for(j in 1:5) {
Large_fromYtoA_interaction[[k]][[j]] <- Large_A_xmap_Y_interaction[[k]][[j]]/Large_Y_xmap_Y_interaction[[k]][[j]]
Large_fromAtoY_interaction[[k]][[j]] <- Large_Y_xmap_A_interaction[[k]][[j]]/Large_A_xmap_A_interaction[[k]][[j]]
}
}
Huge_fromYtoA_interaction <- list()
Huge_fromAtoY_interaction <- list()
for(k in 1:length(Huge_ccmY_xmap_A_combined)){
Huge_fromYtoA_interaction[[k]] <- list()
Huge_fromAtoY_interaction[[k]] <- list()
for(j in 1:5) {
Huge_fromYtoA_interaction[[k]][[j]] <- Huge_A_xmap_Y_interaction[[k]][[j]]/Huge_Y_xmap_Y_interaction[[k]][[j]]
Huge_fromAtoY_interaction[[k]][[j]] <- Huge_Y_xmap_A_interaction[[k]][[j]]/Huge_A_xmap_A_interaction[[k]][[j]]
}
}
#Defining the function of standard error
se <- function(x, na.inf.rm = T) {
if(na.inf.rm == T) {
y <- x[!is.na(x)]
y <- y[!is.infinite(y)]
} else {
y <- x
}
sd(y)/sqrt(length(y))
}
#For Small
Small_fromYtoA_interaction_mean <- c()
Small_fromAtoY_interaction_mean <- c()
Small_fromYtoA_interaction_se <- c()
Small_fromAtoY_interaction_se <- c()
for(j in 1:5){
Small_fromYtoA_interaction_mean[[j]] <- list()
Small_fromAtoY_interaction_mean[[j]] <- list()
Small_fromYtoA_interaction_se[[j]] <- list()
Small_fromAtoY_interaction_se[[j]] <- list()
for(k in 1:length(Small_ccmY_xmap_A_combined)) {
Small_fromYtoA_interaction_mean[[j]][[k]] <- mean(Small_fromYtoA_interaction[[k]][[j]][!is.infinite(Small_fromYtoA_interaction[[k]][[j]])], na.rm = T)
Small_fromAtoY_interaction_mean[[j]][[k]] <- mean(Small_fromAtoY_interaction[[k]][[j]][!is.infinite(Small_fromAtoY_interaction[[k]][[j]])], na.rm = T)
Small_fromYtoA_interaction_se[[j]][[k]] <- se(Small_fromYtoA_interaction[[k]][[j]], na.inf.rm = T)
}
}
#For Medium
Medium_fromYtoA_interaction_mean <- c()
Medium_fromYtoA_interaction_se <- c()
Medium_fromAtoY_interaction_se <- c()
for(j in 1:5){
Medium_fromYtoA_interaction_mean[[j]] <- list()
Medium_fromAtoY_interaction_mean[[j]] <- list()
Medium_fromYtoA_interaction_se[[j]] <- list()
Medium_fromAtoY_interaction_se[[j]] <- list()
for(k in 1:length(Medium_ccmY_xmap_A_combined)) {
Medium_fromYtoA_interaction_mean[[j]][[k]] <- mean(Medium_fromYtoA_interaction[[k]][[j]][!is.infinite(Medium_fromYtoA_interaction[[k]][[j]])], na.rm = T)
Medium_fromAtoY_interaction_mean[[j]][[k]] <- mean(Medium_fromAtoY_interaction[[k]][[j]][!is.infinite(Medium_fromAtoY_interaction[[k]][[j]])], na.rm = T)
Medium_fromYtoA_interaction_se[[j]][[k]] <- se(Medium_fromYtoA_interaction[[k]][[j]], na.inf.rm = T)
Medium_fromAtoY_interaction_se[[j]][[k]] <- se(Medium_fromAtoY_interaction[[k]][[j]], na.inf.rm = T)
}
#For Large
Large_fromYtoA_interaction_mean <- c()
Large_fromAtoY_interaction_mean <- c()
Large_fromYtoA_interaction_se <- c()
Large_fromAtoY_interaction_se <- c()
for(j in 1:5){
Large_fromYtoA_interaction_mean[[j]] <- list()
Large_fromYtoA_interaction_se[[j]] <- list()
Large_fromAtoY_interaction_se[[j]] <- list()
for(k in 1:length(Large_ccmY_xmap_A_combined)) {
Large_fromYtoA_interaction_mean[[j]][[k]] <- mean(Large_fromYtoA_interaction[[k]][[j]][!is.infinite(Large_fromYtoA_interaction[[k]][[j]])], na.rm = T)
Large_fromAtoY_interaction_mean[[j]][[k]] <- mean(Large_fromAtoY_interaction[[k]][[j]][!is.infinite(Large_fromAtoY_interaction[[k]][[j]])], na.rm = T)
Large_fromYtoA_interaction_se[[j]][[k]] <- se(Large_fromYtoA_interaction[[k]][[j]], na.inf.rm = T)
}
}
#For Huge
Huge_fromAtoY_interaction_mean <- c()
Huge_fromYtoA_interaction_se <- c()
Huge_fromAtoY_interaction_se <- c()
for(j in 1:5){
Huge_fromYtoA_interaction_mean[[j]] <- list()
Huge_fromAtoY_interaction_mean[[j]] <- list()
Huge_fromYtoA_interaction_se[[j]] <- list()
Huge_fromAtoY_interaction_se[[j]] <- list()
for(k in 1:length(Huge_ccmY_xmap_A_combined)) {
Huge_fromYtoA_interaction_mean[[j]][[k]] <- mean(Huge_fromYtoA_interaction[[k]][[j]][!is.infinite(Huge_fromYtoA_interaction[[k]][[j]])], na.rm = T)
Huge_fromAtoY_interaction_mean[[j]][[k]] <- mean(Huge_fromAtoY_interaction[[k]][[j]][!is.infinite(Huge_fromAtoY_interaction[[k]][[j]])], na.rm = T)
Huge_fromAtoY_interaction_se[[j]][[k]] <- se(Huge_fromAtoY_interaction[[k]][[j]], na.inf.rm = T)
}
}
Total_interaction <- list()
Small_interaction <- list()
Medium_interaction <- list()
Large_interaction <- list()
Huge_interaction <- list()
#Small combined (S1, S2, S5, & S6)
#Medium combined (M2 to M7)
#Large combined (L1, L4 to L8)
#Huge combined (H1, H2, H3, H6, H7, H8)
total_average_coexistence_week <- subset(Table1, (theta_for_C.m > 0 & theta_for_C.c > 0))$average_coexistence_week
#Generating dataframe for each cage size
for(j in 1:5){
Small_interaction[[j]] <- data.frame(volume = rep(Smallvolume, length(Small_ccmY_xmap_A_combined)),
volume_class = rep("Small", length(Small_ccmY_xmap_A_combined)),
YtoA_mean = unlist(Small_fromYtoA_interaction_mean[[j]]),
AtoY_mean = unlist(Small_fromAtoY_interaction_mean[[j]]),
AtoY_se = unlist(Small_fromAtoY_interaction_se[[j]])
Medium_interaction[[j]] <- data.frame(volume = rep(Mediumvolume, length(Medium_ccmY_xmap_A_combined)),
volume_class = rep("Medium", length(Medium_ccmY_xmap_A_combined)),
YtoA_mean = unlist(Medium_fromYtoA_interaction_mean[[j]]),
AtoY_mean = unlist(Medium_fromAtoY_interaction_mean[[j]]),
YtoA_se = unlist(Medium_fromYtoA_interaction_se[[j]]),
AtoY_se = unlist(Medium_fromAtoY_interaction_se[[j]])
)
Large_interaction[[j]] <- data.frame(volume = rep(Largevolume, length(Large_ccmY_xmap_A_combined)),
volume_class = rep("Large", length(Large_ccmY_xmap_A_combined)),
YtoA_mean = unlist(Large_fromYtoA_interaction_mean[[j]]),
AtoY_mean = unlist(Large_fromAtoY_interaction_mean[[j]]),
YtoA_se = unlist(Large_fromYtoA_interaction_se[[j]]),
AtoY_se = unlist(Large_fromAtoY_interaction_se[[j]])
Huge_interaction[[j]] <- data.frame(volume = rep(Hugevolume, length(Huge_ccmY_xmap_A_combined)),
volume_class = rep("Huge", length(Huge_ccmY_xmap_A_combined)),
YtoA_mean = unlist(Huge_fromYtoA_interaction_mean[[j]]),
AtoY_mean = unlist(Huge_fromAtoY_interaction_mean[[j]]),
AtoY_se = unlist(Huge_fromAtoY_interaction_se[[j]])
)
}
#Combining all cage size results
for(j in 1:5) {
Total_interaction[[j]] <- rbind.data.frame(Small_interaction[[j]], Medium_interaction[[j]])
Total_interaction[[j]] <- rbind.data.frame(Total_interaction[[j]], Large_interaction[[j]])
Total_interaction[[j]] <- data.frame(coexistence_week = total_average_coexistence_week, Total_interaction[[j]])
}
#from tp = 0 (j = 1) to tp = -4 (j = 5)
for(j in 1:5) {
#OLS
cat("tp = ", 1 - j, "\n")
print(summary(lm(YtoA_mean ~ volume, data = Total_interaction[[j]]))$coefficients)
print(summary(lm(YtoA_mean ~ volume, data = Total_interaction[[j]]))$adj.r.squared)
}
tp = 0
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.956516740 2.471470177 2.005493 0.05863042
volume -0.001867062 0.001522433 -1.226368 0.23430373
[1] 0.02343647
tp = -1
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.488500152 4.685638659 1.811599 0.08509321
volume -0.003556826 0.002886367 -1.232285 0.23213374
[1] 0.02409671
tp = -2
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.0990408516 1.4175473480 2.186199 0.04086026
volume -0.0009601758 0.0008732134 -1.099589 0.28457147
[1] 0.009858782
tp = -3
Estimate Std. Error t value Pr(>|t|)
(Intercept) 62.04503631 37.84346091 1.639518 0.1167408
volume -0.02944008 0.02331168 -1.262889 0.2211548
[1] 0.02754771
tp = -4
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.299108574 1.3504268264 2.443012 0.02396969
volume -0.001147594 0.0008318669 -1.379540 0.18295893
[1] 0.04123299
for(j in 1:5) {
#WLS
cat("tp = ", 1 - j, "\n")
print(summary(lm(YtoA_mean ~ volume, weights = 1/YtoA_se, data = Total_interaction[[j]]))$coefficients)
print(summary(lm(YtoA_mean ~ volume, weights = 1/YtoA_se, data = Total_interaction[[j]]))$adj.r.squared)
}
tp = 0
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.394569e-01 0.2576228569 3.646636 0.001604627
volume -4.843997e-05 0.0001643582 -0.294722 0.771243725
[1] -0.04545952
tp = -1
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.674990e-01 0.353851071 2.4515935 0.0235376
volume 5.861501e-05 0.000232539 0.2520653 0.8035628
[1] -0.04667488
tp = -2
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.944328e-01 0.3122800624 2.5439755 0.01932516
volume 5.950965e-05 0.0001975242 0.3012777 0.76631277
[1] -0.0452562
tp = -3
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.0426898814 0.6050111500 1.7234226 0.1002383
volume -0.0002314866 0.0003547476 -0.6525388 0.5214840
[1] -0.02811116
tp = -4
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.0884033754 0.279558373 3.893296 0.0009028016
volume -0.0001672836 0.000162685 -1.028267 0.3161005535
[1] 0.002722675
#from tp = 0 (j = 1) to tp = -4 (j = 5)
for(j in 1:5) {
#OLS
cat("tp = ", 1 - j, "\n")
print(summary(lm(AtoY_mean ~ volume, data = Total_interaction[[j]]))$coefficients)
print(summary(lm(AtoY_mean ~ volume, data = Total_interaction[[j]]))$adj.r.squared)
}
tp = 0
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.528796e-01 1.510284e-01 5.6471479 1.583017e-05
volume -5.946242e-05 9.303393e-05 -0.6391477 5.299798e-01
[1] -0.02898253
tp = -1
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.511947e-01 0.1762732625 5.3961374 2.784154e-05
volume -5.523822e-05 0.0001085849 -0.5087102 6.165214e-01
[1] -0.03658729
tp = -2
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.134340e-01 0.1695187643 4.7984896 0.0001093882
volume 5.672491e-05 0.0001044241 0.5432168 0.5929839809
[1] -0.03473331
tp = -3
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.2912565024 0.9593215587 2.3884134 0.02689643
volume -0.0003962964 0.0005909449 -0.6706149 0.51013661
[1] -0.02690871
tp = -4
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.208008954 0.4174694407 2.893646 0.008980455
volume 0.000263145 0.0002571624 1.023264 0.318401643
[1] 0.002236346
for(j in 1:5) {
#WLS
cat("tp = ", 1 - j, "\n")
print(summary(lm(AtoY_mean ~ volume, weights = 1/AtoY_se, data = Total_interaction[[j]]))$coefficients)
print(summary(lm(AtoY_mean ~ volume, weights = 1/AtoY_se, data = Total_interaction[[j]]))$adj.r.squared)
}
tp = 0
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.397951e-01 1.492400e-01 3.6169613 0.001719239
volume 5.806017e-05 9.540103e-05 0.6085906 0.549646249
[1] -0.03090848
tp = -1
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.4820636968 0.1760168864 2.7387355 0.01265596
volume 0.0000979778 0.0001143731 0.8566505 0.40178610
[1] -0.01283649
tp = -2
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.4987991265 0.1571672501 3.173684 0.004774312
volume 0.0001227377 0.0001029975 1.191657 0.247345226
[1] 0.01960996
tp = -3
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.489166765 0.2096514891 2.333238 0.03018845
volume 0.000219598 0.0001442595 1.522243 0.14360308
[1] 0.05902278
tp = -4
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.7330156259 0.247194431 2.9653404 0.007649225
volume 0.0001433075 0.000162822 0.8801483 0.389232774
[1] -0.01084682
for(j in 1:5) {
Fig3 <- ggplot(data = Total_interaction[[j]], aes(x = volume, y = YtoA_mean))
Fig3 <- Fig3 + geom_point(size = 1 + 0.2/Total_interaction[[j]]$YtoA_se, shape = 21)
Fig3 <- Fig3 + geom_smooth(method = 'lm', se = F, color = 'black', linetype = "dashed")
Fig3 <- Fig3 + geom_smooth(method = 'lm', mapping = aes(weight = 1/Total_interaction[[j]]$YtoA_se), se = F, color = 'red', linetype = "dashed")
Fig3 <- Fig3 + labs(title = paste("From Cm to Cc: tp = ", 1 - j, sep = ""), y = "Relative interaction strength", x = "Volume (cm3)")
print(Fig3)
}
for(j in 1:5) {
Fig4 <- ggplot(data = Total_interaction[[j]], aes(x = volume, y = AtoY_mean))
Fig4 <- Fig4 + geom_point(size = 1 + 0.2/Total_interaction[[j]]$AtoY_se, shape = 21)
Fig4 <- Fig4 + geom_smooth(method = 'lm', se = F, color = 'black', linetype = "dashed")
Fig4 <- Fig4 + geom_smooth(method = 'lm', mapping = aes(weight = 1/Total_interaction[[j]]$AtoY_se), se = F, color = 'red', linetype = "dashed")
Fig4 <- Fig4 + labs(title = paste("From Cc to Cm: tp = ", 1 - j, sep = ""), y = "Relative interaction strength", x = "Volume (cm3)")
print(Fig4)
}
– This is because the optimal embedding dimention (E*) is relatively high (= 8) but the minimum library size is 8 and the maximum library size is 11. The difference between them is too short.
#For S6 CCM results
Small_ccmY_xmap_A_combined[[4]][[1]]
for(j in 1:5) {
#effect of Y (Cm) on A (Cc)
cat("tp = ", 1 - j, "\n")
cat("OLS\n")
print(anova(lm(YtoA_mean ~ volume_class, data = Total_interaction[[j]])))
cat("\n")
}
tp = 0
OLS
Analysis of Variance Table
Response: YtoA_mean
Df Sum Sq Mean Sq F value Pr(>F)
volume_class 3 200.74 66.913 1.6069 0.2228
Residuals 18 749.55 41.642
tp = -1
OLS
Analysis of Variance Table
Response: YtoA_mean
Df Sum Sq Mean Sq F value Pr(>F)
volume_class 3 718.55 239.52 1.5971 0.225
Residuals 18 2699.48 149.97
tp = -2
OLS
Analysis of Variance Table
Response: YtoA_mean
Df Sum Sq Mean Sq F value Pr(>F)
volume_class 3 65.625 21.875 1.6223 0.2193
Residuals 18 242.711 13.484
tp = -3
OLS
Analysis of Variance Table
Response: YtoA_mean
Df Sum Sq Mean Sq F value Pr(>F)
volume_class 3 48231 16077.0 1.6488 0.2135
Residuals 18 175517 9750.9
tp = -4
OLS
Analysis of Variance Table
Response: YtoA_mean
Df Sum Sq Mean Sq Pr(>F)
volume_class 3 66.60 22.200 1.7969 0.1839
Residuals 18 222.38 12.355
for(j in 1:5) {
#effect of Y (Cm) on A (Cc)
cat("WLS\n")
print(anova(lm(YtoA_mean ~ volume_class, weights = 1/YtoA_se, data = Total_interaction[[j]])))
cat("\n")
}
WLS
Analysis of Variance Table
Response: YtoA_mean
Df Sum Sq Mean Sq F value Pr(>F)
volume_class 3 7.538 2.5127 0.4034 0.7523
Residuals 18 112.105 6.2281
WLS
Analysis of Variance Table
Response: YtoA_mean
Df Sum Sq Mean Sq F value Pr(>F)
volume_class 3 7.346 2.4486 0.1641 0.9191
Residuals 18 268.504 14.9169
WLS
Analysis of Variance Table
Response: YtoA_mean
Df Sum Sq Mean Sq F value Pr(>F)
volume_class 3 10.319 3.4395 0.3136 0.8153
Residuals 18 197.435 10.9686
WLS
Analysis of Variance Table
Response: YtoA_mean
Df Sum Sq Mean Sq F value Pr(>F)
volume_class 3 17.70 5.901 0.1603 0.9217
Residuals 18 662.69 36.816
WLS
Analysis of Variance Table
Response: YtoA_mean
Df Sum Sq Mean Sq F value Pr(>F)
volume_class 3 12.975 4.3249 0.6923 0.5686
Residuals 18 112.454 6.2474
for(j in 1:5) {
#effect of A (Cc) on Y (Cm)
cat("OLS\n")
print(anova(lm(AtoY_mean ~ volume_class, data = Total_interaction[[j]])))
cat("\n")
}
OLS
Analysis of Variance Table
Df Sum Sq Mean Sq F value Pr(>F)
volume_class 3 0.50459 0.16820 1.0574 0.3916
Residuals 18 2.86328 0.15907
OLS
Analysis of Variance Table
Response: AtoY_mean
Df Sum Sq Mean Sq F value Pr(>F)
volume_class 3 0.7377 0.24590 1.1597 0.3524
Residuals 18 3.8165 0.21203
OLS
Analysis of Variance Table
Response: AtoY_mean
Df Sum Sq Mean Sq F value Pr(>F)
volume_class 3 0.4752 0.15840 0.7615 0.5303
Residuals 18 3.7442 0.20801
OLS
Analysis of Variance Table
Response: AtoY_mean
Df Sum Sq Mean Sq F value Pr(>F)
volume_class 3 13.578 4.5261 0.6646 0.5846
Residuals 18 122.579 6.8100
OLS
Analysis of Variance Table
Response: AtoY_mean
Df Sum Sq Mean Sq F value Pr(>F)
volume_class 1.965 0.65499 0.4798 0.7004
Residuals 18 24.573 1.36516
for(j in 1:5) {
#effect of A (Cc) on Y (Cm)
cat("WLS\n")
print(anova(lm(AtoY_mean ~ volume_class, weights = 1/AtoY_se, data = Total_interaction[[j]])))
cat("\n")
}
WLS
Analysis of Variance Table
Response: AtoY_mean
Df Sum Sq Mean Sq F value Pr(>F)
volume_class 3 9.116 3.0388 0.8753 0.4722
Residuals 18 62.489 3.4716
WLS
Analysis of Variance Table
Response: AtoY_mean
Df Sum Sq Mean Sq F value Pr(>F)
volume_class 3 20.940 6.9799 1.4846 0.2524
Residuals 18 84.625 4.7014
WLS
Analysis of Variance Table
Response: AtoY_mean
Df Sum Sq Mean Sq F value Pr(>F)
volume_class 3 11.746 3.9152 1.1233 0.3659
Residuals 18 62.741 3.4856
WLS
Analysis of Variance Table
Response: AtoY_mean
Df Sum Sq Mean Sq F value Pr(>F)
volume_class 3 20.678 6.8925 0.9943 0.4179
Residuals 18 124.771 6.9317
WLS
Analysis of Variance Table
Response: AtoY_mean
Df Sum Sq Mean Sq F value Pr(>F)
volume_class 3 16.188 5.3960 1.1479 0.3567
Residuals 18 84.612 4.7007
– The regression relationship (and differences between different sizes, indicated by ANOVA) was not statistically significant.
Why we used the Gamma distribution? – see the following review paper https://www.jstage.jst.go.jp/article/weed/55/4/55_4_268/_pdf
The workflow: We started with checking the individual effect of a specific interaction (Cm (Y) -> Cc (A) or Cc (A) -> Cm (Y)) for each tp (in this section) by applying stepAIC() function. Then, we included all of the statistically significant factors into a single model and selected the best model based on AIC/BIC (in the next section).
#Although we know the Gamma distribution is more appropriate than gaussian, we also tried to check if there is no substantial difference between them
#GLM with Gamma
glm_YtoA_tp0 <- glm(coexistence_week ~ volume + YtoA_mean + volume:YtoA_mean, data = Total_interaction[[1]], family = Gamma(link = identity))
#Get GLM regression lines
glm_YtoA_tp0_predicted <- ggpredict(glm_YtoA_tp0, term = c("YtoA_mean","volume"))
#define linetypes & linewidth
linewidths <- c(0.5, 0.75, 1.0, 1.25)
#ggplot
glm_YtoA_tp0_fig <- ggplot() + geom_point(data = Total_interaction[[1]], aes(x = YtoA_mean, y = coexistence_week, shape = as.factor(volume_class), size = 3), alpha = 0.5)
glm_YtoA_tp0_fig <- glm_YtoA_tp0_fig + geom_line(data = glm_YtoA_tp0_predicted, aes(x = x, y = predicted, linewidth = group)) + scale_linewidth_manual(values = linewidths)
plot(glm_YtoA_tp0_fig)
##With gaussian distribution
glm_YtoA_tp0g <- glm(coexistence_week ~ volume + YtoA_mean + volume:YtoA_mean, data = Total_interaction[[1]], family = gaussian(link = identity))
##Compare Gamma and gaussian
summary(glm_YtoA_tp0)
Call:
glm(formula = coexistence_week ~ volume + YtoA_mean + volume:YtoA_mean,
family = Gamma(link = identity), data = Total_interaction[[1]])
Deviance Residuals:
Min 1Q Median 3Q Max
-0.56953 -0.24133 -0.02554 0.18817 0.65057
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.1487253 2.9087137 6.583 3.49e-06 ***
volume 0.0024261 0.0029965 0.810 0.429
YtoA_mean -0.4143180 0.4060919 -1.020 0.321
volume:YtoA_mean 0.0004502 0.0021277 0.212 0.835
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for Gamma family taken to be 0.1202773)
Null deviance: 3.1418 on 21 degrees of freedom
Residual deviance: 2.1771 on 18 degrees of freedom
AIC: 154.98
Number of Fisher Scoring iterations: 8
summary(glm_YtoA_tp0g)
Call:
glm(formula = coexistence_week ~ volume + YtoA_mean + volume:YtoA_mean,
family = gaussian(link = identity), data = Total_interaction[[1]])
Deviance Residuals:
Min 1Q Median 3Q Max
-10.6109 -5.1950 -0.6275 5.2976 18.8217
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.1900327 3.2285755 5.944 1.26e-05 ***
volume 0.0022484 0.0028646 0.785 0.443
YtoA_mean -0.4477818 0.4535497 -0.987 0.337
volume:YtoA_mean 0.0006182 0.0019413 0.318 0.754
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 60.77681)
Null deviance: 1438.1 on 21 degrees of freedom
Residual deviance: 1094.0 on 18 degrees of freedom
AIC: 158.38
Number of Fisher Scoring iterations: 2
##Confirm the best model and significant factors
summary(stepAIC(glm_YtoA_tp0))
Start: AIC=154.98
coexistence_week ~ volume + YtoA_mean + volume:YtoA_mean
Df Deviance AIC
- volume:YtoA_mean 1 2.1842 153.04
<none> 2.1771 154.99
Step: AIC=153.06
coexistence_week ~ volume + YtoA_mean
Df Deviance AIC
<none> 2.1842 153.06
- volume 1 2.4857 153.72
- YtoA_mean 1 2.6300 155.00
Call:
glm(formula = coexistence_week ~ volume + YtoA_mean, family = Gamma(link = identity),
data = Total_interaction[[1]])
Deviance Residuals:
Min 1Q Median 3Q Max
-0.56700 -0.25474 -0.02682 0.18579 0.62444
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.013530 2.767102 6.871 1.48e-06 ***
volume 0.002943 0.001823 1.615 0.1229
YtoA_mean -0.331938 0.121350 -2.735 0.0131 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for Gamma family taken to be 0.1132062)
Null deviance: 3.1418 on 21 degrees of freedom
Residual deviance: 2.1842 on 19 degrees of freedom
AIC: 153.06
Number of Fisher Scoring iterations: 4
#Prepare objects to save results
glm_YtoA_tp <- list()
glm_YtoA_tp_predicted <- list()
glm_YtoA_tp_fig <- list()
#assign results from tp = 0
glm_YtoA_tp[[1]] <- glm_YtoA_tp0
glm_YtoA_tp_predicted[[1]] <- glm_YtoA_tp0_predicted
glm_YtoA_tp_fig[[1]] <- glm_YtoA_tp0_fig
#GLM with Gamma and getting GLM regression lines
for(j in 2:5) {
glm_YtoA_tp[[j]] <- glm(coexistence_week ~ volume + YtoA_mean + volume:YtoA_mean, data = Total_interaction[[j]], family = Gamma(link = identity))
glm_YtoA_tp_predicted[[j]] <- ggpredict(glm_YtoA_tp[[j]], term = c("YtoA_mean","volume"))
}
##Checking glm results
for(j in 1:5) {
cat("tp = ", 1 - j, "\n")
print(summary(glm_YtoA_tp[[j]]))
}
tp = 0
Call:
glm(formula = coexistence_week ~ volume + YtoA_mean + volume:YtoA_mean,
family = Gamma(link = identity), data = Total_interaction[[1]])
Deviance Residuals:
Min 1Q Median 3Q Max
-0.56953 -0.24133 -0.02554 0.18817 0.65057
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.1487253 2.9087137 6.583 3.49e-06 ***
volume 0.0024261 0.0029965 0.810 0.429
YtoA_mean -0.4143180 0.4060919 -1.020 0.321
volume:YtoA_mean 0.0004502 0.0021277 0.212 0.835
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for Gamma family taken to be 0.1202773)
Null deviance: 3.1418 on 21 degrees of freedom
Residual deviance: 2.1771 on 18 degrees of freedom
AIC: 154.98
Number of Fisher Scoring iterations: 8
tp = -1
Call:
glm(formula = coexistence_week ~ volume + YtoA_mean + volume:YtoA_mean,
family = Gamma(link = identity), data = Total_interaction[[j]])
Deviance Residuals:
Min 1Q Median 3Q Max
-0.56270 -0.25173 -0.03036 0.18338 0.62510
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 18.8269291 2.8831247 6.530 3.88e-06 ***
volume 0.0030717 0.0035599 0.863 0.400
YtoA_mean -0.1538509 0.4839894 -0.318 0.754
volume:YtoA_mean -0.0001205 0.0026681 -0.045 0.964
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for Gamma family taken to be 0.1187723)
Null deviance: 3.1418 on 21 degrees of freedom
Residual deviance: 2.1718 on 18 degrees of freedom
AIC: 154.93
Number of Fisher Scoring iterations: 7
tp = -2
Call:
glm(formula = coexistence_week ~ volume + YtoA_mean + volume:YtoA_mean,
family = Gamma(link = identity), data = Total_interaction[[j]])
Deviance Residuals:
Min 1Q Median 3Q Max
-0.56753 -0.27702 -0.01154 0.22011 0.61204
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.0546231 2.9806117 6.393 5.09e-06 ***
volume 0.0036129 0.0032318 1.118 0.278
YtoA_mean -0.4831568 0.4590060 -1.053 0.306
volume:YtoA_mean -0.0005512 0.0021267 -0.259 0.798
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for Gamma family taken to be 0.1187819)
Null deviance: 3.1418 on 21 degrees of freedom
Residual deviance: 2.1650 on 18 degrees of freedom
AIC: 154.86
Number of Fisher Scoring iterations: 6
tp = -3
Call:
glm(formula = coexistence_week ~ volume + YtoA_mean + volume:YtoA_mean,
family = Gamma(link = identity), data = Total_interaction[[j]])
Deviance Residuals:
Min 1Q Median 3Q Max
-0.56743 -0.26607 -0.05633 0.22742 0.61632
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 18.713520 2.776101 6.741 2.56e-06 ***
volume 0.003616 0.002644 1.368 0.188
YtoA_mean 0.171660 0.492285 0.349 0.731
volume:YtoA_mean -0.001084 0.002761 -0.393 0.699
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for Gamma family taken to be 0.1193296)
Null deviance: 3.1418 on 21 degrees of freedom
Residual deviance: 2.1659 on 18 degrees of freedom
AIC: 154.87
Number of Fisher Scoring iterations: 4
tp = -4
Call:
glm(formula = coexistence_week ~ volume + YtoA_mean + volume:YtoA_mean,
family = Gamma(link = identity), data = Total_interaction[[j]])
Deviance Residuals:
Min 1Q Median 3Q Max
-0.57584 -0.23760 -0.02773 0.18441 0.68569
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.556062 3.013417 6.490 4.2e-06 ***
volume 0.001660 0.003033 0.547 0.591
YtoA_mean -0.846160 0.542587 -1.559 0.136
volume:YtoA_mean 0.001356 0.002692 0.504 0.621
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for Gamma family taken to be 0.122981)
Null deviance: 3.1418 on 21 degrees of freedom
Residual deviance: 2.1864 on 18 degrees of freedom
AIC: 155.08
Number of Fisher Scoring iterations: 9
##Confirm the best model and significant factors
for(j in 1:5) {
cat("tp = ", 1 - j, "\n")
print(summary(stepAIC(glm_YtoA_tp[[j]]))$coefficients)
cat("\n")
}
tp = 0
Start: AIC=154.98
coexistence_week ~ volume + YtoA_mean + volume:YtoA_mean
Df Deviance AIC
- volume:YtoA_mean 1 2.1842 153.04
<none> 2.1771 154.99
Step: AIC=153.06
coexistence_week ~ volume + YtoA_mean
Df Deviance AIC
<none> 2.1842 153.06
- volume 1 2.4857 153.72
- YtoA_mean 1 2.6300 155.00
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.013529677 2.767101953 6.871279 1.483905e-06
volume 0.002943276 0.001822849 1.614657 1.228696e-01
YtoA_mean -0.331938169 0.121350031 -2.735378 1.314486e-02
tp = -1
Start: AIC=154.93
coexistence_week ~ volume + YtoA_mean + volume:YtoA_mean
Df Deviance AIC
- volume:YtoA_mean 1 2.1721 152.93
<none> 2.1718 154.93
Step: AIC=152.93
coexistence_week ~ volume + YtoA_mean
Df Deviance AIC
<none> 2.1721 152.93
- volume 1 2.4699 153.57
- YtoA_mean 1 2.6300 154.99
Estimate Std. Error t value Pr(>|t|)
(Intercept) 18.868489371 2.725136529 6.923869 1.335421e-06
volume 0.002927838 0.001820163 1.608558 1.242032e-01
YtoA_mean -0.175652588 0.063553574 -2.763851 1.235798e-02
tp = -2
Start: AIC=154.86
coexistence_week ~ volume + YtoA_mean + volume:YtoA_mean
Df Deviance AIC
- volume:YtoA_mean 1 2.1735 152.93
<none> 2.1650 154.86
Step: AIC=152.95
coexistence_week ~ volume + YtoA_mean
Df Deviance AIC
<none> 2.1735 152.95
- volume 1 2.4822 153.68
- YtoA_mean 1 2.6300 154.99
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.222023454 2.806067088 6.850165 1.548227e-06
volume 0.002968237 0.001809226 1.640611 1.173288e-01
YtoA_mean -0.586320536 0.212239304 -2.762545 1.239307e-02
tp = -3
Start: AIC=154.87
coexistence_week ~ volume + YtoA_mean + volume:YtoA_mean
Df Deviance AIC
- volume:YtoA_mean 1 2.1837 153.02
<none> 2.1659 154.87
Step: AIC=153.05
coexistence_week ~ volume + YtoA_mean
Df Deviance AIC
<none> 2.1837 153.05
- volume 1 2.4767 153.63
- YtoA_mean 1 2.6300 154.98
Estimate Std. Error t value Pr(>|t|)
(Intercept) 18.724736327 2.700453036 6.933924 1.308828e-06
volume 0.002914395 0.001828742 1.593661 1.275120e-01
YtoA_mean -0.021647919 0.007910371 -2.736650 1.310870e-02
tp = -4
Start: AIC=155.08
coexistence_week ~ volume + YtoA_mean + volume:YtoA_mean
Df Deviance AIC
- volume:YtoA_mean 1 2.2300 153.44
<none> 2.1864 155.08
Step: AIC=153.52
coexistence_week ~ volume + YtoA_mean
Df Deviance AIC
<none> 2.2300 153.52
- volume 1 2.5178 154.02
- YtoA_mean 1 2.6300 154.99
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.331803605 2.880271831 6.711798 2.047614e-06
volume 0.002895789 0.001852905 1.562837 1.345936e-01
YtoA_mean -0.598473622 0.228182515 -2.622785 1.674965e-02
– No statistically significant interactions of cage size and interaction strength
– for all tp, the interaction strength (YtoA_mean) has negative relationship with coexistence
– for all tp, the volume was not statistically significant
#define linetypes & linewidth
linetypes <- c("dashed", "dotted", "dotdash", "longdash")
#ggplot
for(j in 1:5) {
glm_YtoA_tp_fig[[j]] <- ggplot() + geom_point(data = Total_interaction[[j]], aes(x = YtoA_mean, y = coexistence_week, shape = as.factor(volume_class), size = 3), alpha = 0.5)
glm_YtoA_tp_fig[[j]] <- glm_YtoA_tp_fig[[j]] + geom_line(data = glm_YtoA_tp_predicted[[j]], aes(x = x, y = predicted, linetype = group)) + scale_linetype_manual(values = linetypes)
glm_YtoA_tp_fig[[j]] <- glm_YtoA_tp_fig[[j]] + geom_smooth(data = Total_interaction[[j]], aes(x = YtoA_mean, y = coexistence_week), method = 'glm', method.args = list(family = Gamma(link = identity)), se = F, color = 'black')
}
#plot all
for(j in 1:5) {
glm_YtoA_tp_fig[[j]] <- glm_YtoA_tp_fig[[j]] + ggtitle(paste("Effect of Y (Cm) on A (Cc), tp = ", 1 - j, sep = ""))
plot(glm_YtoA_tp_fig[[j]])
}
#plot again for tp = -3 for narrower y axis
glm_YtoA_tp_fig[[4]] <- glm_YtoA_tp_fig[[4]] + ylim(c(0, NA))
plot(glm_YtoA_tp_fig[[4]])
#plot again for tp = -4 for narrower y axis
glm_YtoA_tp_fig[[5]] <- glm_YtoA_tp_fig[[5]] + ylim(c(0, 50))
plot(glm_YtoA_tp_fig[[5]])
#Prepare objects to save results
glm_AtoY_tp <- list()
glm_AtoY_tp_predicted <- list()
glm_AtoY_tp_fig <- list()
#GLM with Gamma * getting GLM regression lines
for(j in 1:5) {
glm_AtoY_tp[[j]] <- glm(coexistence_week ~ volume + AtoY_mean + volume:AtoY_mean, data = Total_interaction[[j]], family = Gamma(link = identity))
glm_AtoY_tp_predicted[[j]] <- ggpredict(glm_AtoY_tp[[j]], term = c("AtoY_mean","volume"))
}
##Checking glm results
for(j in 1:5) {
cat("tp = ", 1 - j, "\n")
print(summary(glm_AtoY_tp[[j]]))
}
tp = 0
Call:
glm(formula = coexistence_week ~ volume + AtoY_mean + volume:AtoY_mean,
family = Gamma(link = identity), data = Total_interaction[[j]])
Deviance Residuals:
Min 1Q Median 3Q Max
-0.48332 -0.22948 -0.08491 0.16138 0.64320
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.963321 3.610635 1.929 0.0697 .
volume 0.009841 0.004156 2.368 0.0293 *
AtoY_mean 11.453596 4.651596 2.462 0.0241 *
volume:AtoY_mean -0.006891 0.005098 -1.352 0.1932
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for Gamma family taken to be 0.1204898)
Null deviance: 3.1418 on 21 degrees of freedom
Residual deviance: 2.0584 on 18 degrees of freedom
AIC: 153.73
Number of Fisher Scoring iterations: 6
tp = -1
Call:
glm(formula = coexistence_week ~ volume + AtoY_mean + volume:AtoY_mean,
family = Gamma(link = identity), data = Total_interaction[[j]])
Deviance Residuals:
Min 1Q Median 3Q Max
-0.53760 -0.18667 -0.03739 0.08615 0.75734
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.039627 3.187765 2.208 0.0404 *
volume 0.007681 0.003606 2.130 0.0472 *
AtoY_mean 10.029670 3.725952 2.692 0.0149 *
volume:AtoY_mean -0.003494 0.004002 -0.873 0.3941
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for Gamma family taken to be 0.1131852)
Null deviance: 3.1418 on 21 degrees of freedom
Residual deviance: 1.8382 on 18 degrees of freedom
AIC: 151.21
Number of Fisher Scoring iterations: 6
tp = -2
Call:
glm(formula = coexistence_week ~ volume + AtoY_mean + volume:AtoY_mean,
family = Gamma(link = identity), data = Total_interaction[[j]])
Deviance Residuals:
Min 1Q Median 3Q Max
-0.51990 -0.22059 -0.00403 0.10655 0.71976
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.790505 3.295742 2.060 0.0541 .
volume 0.008559 0.003386 2.528 0.0211 *
AtoY_mean 11.319107 4.208559 2.690 0.0150 *
volume:AtoY_mean -0.004951 0.003324 -1.489 0.1537
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for Gamma family taken to be 0.11513)
Null deviance: 3.1418 on 21 degrees of freedom
Residual deviance: 1.9266 on 18 degrees of freedom
AIC: 152.25
Number of Fisher Scoring iterations: 9
tp = -3
Call:
glm(formula = coexistence_week ~ volume + AtoY_mean + volume:AtoY_mean,
family = Gamma(link = identity), data = Total_interaction[[j]])
Deviance Residuals:
Min 1Q Median 3Q Max
-0.62348 -0.26013 -0.04892 0.17585 0.61966
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.121e+01 3.034e+00 3.696 0.00165 **
volume 4.570e-03 3.001e-03 1.523 0.14517
AtoY_mean 3.059e+00 1.762e+00 1.736 0.09973 .
volume:AtoY_mean -6.568e-06 1.978e-03 -0.003 0.99739
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for Gamma family taken to be 0.1116182)
Null deviance: 3.1418 on 21 degrees of freedom
Residual deviance: 2.0209 on 18 degrees of freedom
AIC: 153.32
Number of Fisher Scoring iterations: 22
tp = -4
Call:
glm(formula = coexistence_week ~ volume + AtoY_mean + volume:AtoY_mean,
family = Gamma(link = identity), data = Total_interaction[[j]])
Deviance Residuals:
Min 1Q Median 3Q Max
-0.54459 -0.33084 -0.04793 0.26828 0.59890
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.906263 4.466563 2.218 0.0397 *
volume 0.008610 0.003550 2.425 0.0260 *
AtoY_mean 5.403578 3.488728 1.549 0.1388
volume:AtoY_mean -0.003289 0.002127 -1.546 0.1395
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for Gamma family taken to be 0.1342207)
Null deviance: 3.1418 on 21 degrees of freedom
Residual deviance: 2.3987 on 18 degrees of freedom
AIC: 157.15
Number of Fisher Scoring iterations: 10
##Confirm the best model and significant factors
for(j in 1:5) {
cat("tp = ", 1 - j, "\n")
print(summary(stepAIC(glm_AtoY_tp[[j]]))$coefficients)
cat("\n")
}
tp = 0
Start: AIC=153.73
coexistence_week ~ volume + AtoY_mean + volume:AtoY_mean
Df Deviance AIC
- volume:AtoY_mean 1 2.2895 153.65
<none> 2.0584 153.73
Step: AIC=154.11
coexistence_week ~ volume + AtoY_mean
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.570563492 3.665854102 2.883520 0.009516926
volume 0.004896986 0.001886922 2.595224 0.017765295
AtoY_mean 6.676695608 3.965930505 1.683513 0.108634858
tp = -1
Start: AIC=151.21
coexistence_week ~ volume + AtoY_mean + volume:AtoY_mean
Df Deviance AIC
- volume:AtoY_mean 1 1.9308 150.02
<none> 1.8382 151.21
Step: AIC=150.3
coexistence_week ~ volume + AtoY_mean
Df Deviance AIC
<none> 1.9308 150.30
- AtoY_mean 1 2.6300 154.12
- volume 1 2.7761 155.34
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.823516067 2.97971854 2.961191 0.008021308
volume 0.004994063 0.00175177 2.850866 0.010222869
AtoY_mean 7.724096295 3.00029654 2.574444 0.018569295
tp = -2
Start: AIC=152.25
coexistence_week ~ volume + AtoY_mean + volume:AtoY_mean
Df Deviance AIC
<none> 1.9266 152.25
- volume:AtoY_mean 1 2.1990 152.62
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.790505430 3.295742353 2.060387 0.05411742
volume 0.008559064 0.003386359 2.527512 0.02106609
AtoY_mean 11.319107188 4.208558974 2.689545 0.01498047
volume:AtoY_mean -0.004951192 0.003324305 -1.489392 0.15369740
tp = -3
Start: AIC=153.32
coexistence_week ~ volume + AtoY_mean + volume:AtoY_mean
Df Deviance AIC
- volume:AtoY_mean 1 2.0209 151.32
<none> 2.0209 153.32
Step: AIC=151.32
coexistence_week ~ volume + AtoY_mean
Df Deviance AIC
<none> 2.0209 151.32
- AtoY_mean 1 2.6300 155.08
- volume 1 2.8015 156.70
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.220754735 2.450685303 4.578619 0.0002051581
volume 0.004561095 0.001615429 2.823457 0.0108540142
AtoY_mean 3.052954403 1.156210339 2.640484 0.0161266887
tp = -4
Start: AIC=157.15
coexistence_week ~ volume + AtoY_mean + volume:AtoY_mean
Df Deviance AIC
- volume:AtoY_mean 1 2.6292 156.87
<none> 2.3987 157.15
Step: AIC=157.21
coexistence_week ~ volume + AtoY_mean
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.032196529 3.255919795 5.23114745 4.756549e-05
volume 0.003653972 0.001960856 1.86345754 7.792917e-02
AtoY_mean 0.145022180 1.673824307 0.08664122 9.318633e-01
– No statistically significant interactions of cage size and interaction strength
– tp = 0 & -4: interaction strength was not significant
– tp = -1, -2, & -3: interaction strength was statistically significant with positive slope
– for all tp (except for tp = -4), the volume was statistically significant with positive slope
#define linetypes & linewidth
linetypes <- c("dashed", "dotted", "dotdash", "longdash")
#ggplot
for(j in 1:5) {
glm_AtoY_tp_fig[[j]] <- ggplot() + geom_point(data = Total_interaction[[j]], aes(x = AtoY_mean, y = coexistence_week, shape = as.factor(volume_class), size = 3), alpha = 0.5)
glm_AtoY_tp_fig[[j]] <- glm_AtoY_tp_fig[[j]] + geom_line(data = glm_AtoY_tp_predicted[[j]], aes(x = x, y = predicted, linetype = group)) + scale_linetype_manual(values = linetypes)
}
#adding the overall slopes
glm_AtoY_tp_fig[[1]] <- glm_AtoY_tp_fig[[1]] + geom_smooth(data = Total_interaction[[j]], aes(x = AtoY_mean, y = coexistence_week), method = 'glm', method.args = list(family = Gamma(link = identity)), se = F, color = 'black', linetype = "dashed")
glm_AtoY_tp_fig[[5]] <- glm_AtoY_tp_fig[[5]] + geom_smooth(data = Total_interaction[[j]], aes(x = AtoY_mean, y = coexistence_week), method = 'glm', method.args = list(family = Gamma(link = identity)), se = F, color = 'black', linetype = "dashed")
for(j in 2:4) {
glm_AtoY_tp_fig[[j]] <- glm_AtoY_tp_fig[[j]] + geom_smooth(data = Total_interaction[[j]], aes(x = AtoY_mean, y = coexistence_week), method = 'glm', method.args = list(family = Gamma(link = identity)), se = F, color = 'black')
}
#plot all
for(j in 1:5) {
glm_AtoY_tp_fig[[j]] <- glm_AtoY_tp_fig[[j]] + ggtitle(paste("Effect of A (Cc) on Y (Cm), tp = ", 1 - j, sep = ""))
plot(glm_AtoY_tp_fig[[j]])
}
– For the effect from Y to A: there was negative effect for all tp.
– For the effect from A to Y: tp = 0 & -4: interaction strength from A to Y was not significan, while tp = -1, -2, & -3: interaction strength from A to Y was statistically significant with positive slope.
– for all tp (except for tp = -4), the volume was statistically significant with positive slope.
selected_data_all <- data.frame(coweek = Total_interaction[[1]]$coexistence_week, volume = Total_interaction[[1]]$volume, fromYtoA_tp0 = Total_interaction[[1]]$YtoA_mean, fromYtoA_tp1 = Total_interaction[[2]]$YtoA_mean, fromYtoA_tp2 = Total_interaction[[3]]$YtoA_mean, fromYtoA_tp3 = Total_interaction[[4]]$YtoA_mean, fromYtoA_tp4 = Total_interaction[[5]]$YtoA_mean, fromAtoY_tp0 = Total_interaction[[1]]$AtoY_mean, fromAtoY_tp1 = Total_interaction[[2]]$AtoY_mean, fromAtoY_tp2 = Total_interaction[[3]]$AtoY_mean, fromAtoY_tp3 = Total_interaction[[4]]$AtoY_mean, fromAtoY_tp4 = Total_interaction[[5]]$AtoY_mean)
initial_model <- glm(coweek ~ volume, data = selected_data_all, family = Gamma(link = identity))
selected_model <- stepAIC(initial_model, direction = "both", scope = list(lower = ~ 1, upper = ~ volume + fromYtoA_tp0 + fromYtoA_tp1 + fromYtoA_tp2 + fromYtoA_tp3 + fromYtoA_tp4 + fromAtoY_tp1 + fromAtoY_tp2 + fromAtoY_tp3))
Start: AIC=155.22
coweek ~ volume
Df Deviance AIC
+ fromAtoY_tp1 1 1.9308 151.79
+ fromAtoY_tp3 1 2.0209 152.49
+ fromYtoA_tp1 1 2.1721 153.67
+ fromYtoA_tp2 1 2.1735 153.68
+ fromYtoA_tp3 1 2.1837 153.76
+ fromYtoA_tp0 1 2.1842 153.76
+ fromAtoY_tp2 1 2.1990 153.88
+ fromYtoA_tp4 1 2.2300 154.12
<none> 2.6300 155.22
- volume 1 3.1418 157.19
Step: AIC=150.3
coweek ~ volume + fromAtoY_tp1
Df Deviance AIC
+ fromAtoY_tp3 1 1.5537 149.16
<none> 1.9308 150.30
+ fromAtoY_tp2 1 1.8425 151.57
+ fromYtoA_tp0 1 1.8728 151.82
+ fromYtoA_tp2 1 1.8733 151.82
+ fromYtoA_tp1 1 1.8761 151.85
+ fromYtoA_tp3 1 1.8818 151.89
+ fromYtoA_tp4 1 1.8843 151.91
- fromAtoY_tp1 1 2.6300 154.12
- volume 1 2.7761 155.34
Step: AIC=147.46
coweek ~ volume + fromAtoY_tp1 + fromAtoY_tp3
Df Deviance AIC
<none> 1.5537 147.46
+ fromAtoY_tp2 1 1.5237 149.16
- fromAtoY_tp3 1 1.9308 149.18
+ fromYtoA_tp2 1 1.5474 149.40
+ fromYtoA_tp4 1 1.5489 149.41
+ fromYtoA_tp0 1 1.5495 149.42
+ fromYtoA_tp1 1 1.5501 149.42
+ fromYtoA_tp3 1 1.5513 149.44
- fromAtoY_tp1 1 2.0209 150.07
- volume 1 2.5149 154.94
full_model <- glm(coweek ~ volume + fromYtoA_tp0 + fromYtoA_tp1 + fromYtoA_tp2 + fromYtoA_tp3 + fromYtoA_tp4 + fromAtoY_tp1 + fromAtoY_tp2 + fromAtoY_tp3, data = selected_data_all, family = Gamma(link = identity))
summary(full_model)$coefficients
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.270319682 6.803034136 1.95064723 0.07484329
volume 0.005628476 0.002240214 2.51247296 0.02728368
fromYtoA_tp0 1.903466285 7.712847793 0.24679163 0.80924089
fromYtoA_tp1 -3.606477068 5.719661012 -0.63054035 0.54016374
fromYtoA_tp2 -3.428055623 5.239472859 -0.65427491 0.52527601
fromYtoA_tp3 0.439453895 0.476246352 0.92274491 0.37432115
fromYtoA_tp4 0.172362236 7.428915983 0.02320153 0.98187086
fromAtoY_tp1 10.147933294 9.144140193 1.10977447 0.28884058
fromAtoY_tp2 -5.458304483 9.045245319 -0.60344460 0.55744738
fromAtoY_tp3 1.178939026 1.096043760 1.07563135 0.30324330
summary(selected_model)$coefficients
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.668590059 2.534802765 3.025320 0.007273366
volume 0.005097605 0.001565156 3.256930 0.004378047
fromAtoY_tp1 6.142044968 2.786337749 2.204343 0.040754197
fromAtoY_tp3 1.359815404 0.847675698 1.604169 0.126078459
model02 <- glm(coweek~ volume + fromYtoA_tp0 + fromYtoA_tp1 + fromYtoA_tp2 + fromYtoA_tp3 + fromYtoA_tp4 + fromAtoY_tp1 + fromAtoY_tp2 + fromAtoY_tp3, data = selected_data_all, family = Gamma(link = identity), na.action = "na.fail")
dd2 <- dredge(model02, rank = "AIC")
Fixed term is "(Intercept)"
top.models <- get.models(dd2, cumsum(weight) <= .95)
#topmost model:
summary(top.models[[1]])
Call:
glm(formula = coweek ~ fromAtoY_tp1 + fromAtoY_tp3 + volume +
1, family = Gamma(link = identity), data = selected_data_all,
na.action = "na.fail")
Deviance Residuals:
Min 1Q Median 3Q Max
-0.56013 -0.15749 -0.00414 0.07185 0.80370
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.668590 2.534803 3.025 0.00727 **
fromAtoY_tp1 6.142045 2.786338 2.204 0.04075 *
fromAtoY_tp3 1.359815 0.847676 1.604 0.12608
volume 0.005098 0.001565 3.257 0.00438 **
(Dispersion parameter for Gamma family taken to be 0.1014021)
Null deviance: 3.1418 on 21 degrees of freedom
Residual deviance: 1.5537 on 18 degrees of freedom
AIC: 147.46
Number of Fisher Scoring iterations: 7
model03 <- glm(coweek~ volume + fromYtoA_tp0 + fromYtoA_tp1 + fromYtoA_tp2 + fromYtoA_tp3 + fromYtoA_tp4 + fromAtoY_tp0 + fromAtoY_tp1 + fromAtoY_tp2 + fromAtoY_tp3 + fromAtoY_tp4 + fromAtoY_tp2:volume, data = selected_data_all, family = Gamma(link = identity), na.action = "na.fail")
dd3 <- dredge(model03, rank = "AIC")
Fixed term is "(Intercept)"
top.models <- get.models(dd3, cumsum(weight) <= .95)
#topmost model:
summary(top.models[[1]])
Call:
glm(formula = coweek ~ fromAtoY_tp0 + fromAtoY_tp1 + fromAtoY_tp3 +
fromYtoA_tp2 + fromYtoA_tp3 + volume + 1, family = Gamma(link = identity),
data = selected_data_all, na.action = "na.fail")
Deviance Residuals:
Min 1Q Median 3Q Max
-0.38254 -0.16539 -0.02881 0.11844 0.62597
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.291575 4.678264 3.269 0.00518 **
fromAtoY_tp0 -20.673901 9.452274 -2.187 0.04498 *
fromAtoY_tp1 22.629879 7.997786 2.830 0.01268 *
fromAtoY_tp3 1.059965 0.697067 1.521 0.14916
fromYtoA_tp2 -5.091734 2.788474 -1.826 0.08782 .
fromYtoA_tp3 0.183203 0.103253 1.774 0.09631 .
volume 0.004647 0.001623 2.863 0.01186 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for Gamma family taken to be 0.07830766)
Null deviance: 3.1418 on 21 degrees of freedom
Residual deviance: 1.0493 on 15 degrees of freedom
AIC: 144.74
Number of Fisher Scoring iterations: 8
– When all replicates were combined for each cage size, the minimum of the maximum library length of CCM’s results are 106. Examples from the small cage only.
Small_Y_xmap_A_supp <- list()
Small_A_xmap_Y_supp <- list()
Small_Y_xmap_Y_supp <- list()
Small_A_xmap_A_supp <- list()
min_max_lib_length <- min(max(Small_lagA_xmap_Y[[1]]$lib_size), max(Small_lagY_xmap_A[[1]]$lib_size))
#not showing warning
Small_Y_xmap_A_supp <- lapply(0:-4, function(i) ccm(Small_YA_allcombined, E = Small_optE_Y, lib_column = "Small_Y", target_column = "Small_A", lib_sizes = seq(Small_optE_Y, min_max_lib_length, length = 2), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
Small_A_xmap_Y_supp <- lapply(0:-4, function(i) ccm(Small_YA_allcombined, E = Small_optE_A, lib_column = "Small_A", target_column = "Small_Y", lib_sizes = seq(Small_optE_A, min_max_lib_length, length = 2), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
Small_Y_xmap_Y_supp <- lapply(0:-4, function(i) ccm(Small_YA_allcombined, E = Small_optE_Y, lib_column = "Small_Y", target_column = "Small_Y", lib_sizes = seq(Small_optE_Y, min_max_lib_length, length = 2), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
Small_A_xmap_A_supp <- lapply(0:-4, function(i) ccm(Small_YA_allcombined, E = Small_optE_A, lib_column = "Small_A", target_column = "Small_A", lib_sizes = seq(Small_optE_A, min_max_lib_length, length = 2), num_samples = 100, replace = F, tp = i, RNGseed = 1234))
#example
Small_A_xmap_Y_supp[[2]]
NA
Small_Y_xmap_A_interaction_supp <- list()
Small_Y_xmap_Y_interaction_supp <- list()
Small_A_xmap_Y_interaction_supp <- list()
Small_A_xmap_A_interaction_supp <- list()
Small_fromYtoA_interaction_supp <- list()
Small_fromAtoY_interaction_supp <- list()
for(j in 1:5) {
Small_Y_xmap_A_interaction_supp[[j]] <- subset(Small_Y_xmap_A_supp[[j]]$rho, Small_Y_xmap_A_supp[[j]]$lib_size == max(Small_Y_xmap_A_supp[[j]]$lib_size)) - subset(Small_Y_xmap_A_supp[[j]]$rho, Small_Y_xmap_A_supp[[j]]$lib_size == min(Small_Y_xmap_A_supp[[j]]$lib_size))
Small_A_xmap_Y_interaction_supp[[j]] <- subset(Small_A_xmap_Y_supp[[j]]$rho, Small_A_xmap_Y_supp[[j]]$lib_size == max(Small_A_xmap_Y_supp[[j]]$lib_size)) - subset(Small_A_xmap_Y_supp[[j]]$rho, Small_A_xmap_Y_supp[[j]]$lib_size == min(Small_A_xmap_Y_supp[[j]]$lib_size))
Small_Y_xmap_Y_interaction_supp[[j]] <- subset(Small_Y_xmap_Y_supp[[j]]$rho, Small_Y_xmap_Y_supp[[j]]$lib_size == max(Small_Y_xmap_Y_supp[[j]]$lib_size)) - subset(Small_Y_xmap_Y_supp[[j]]$rho, Small_Y_xmap_Y_supp[[j]]$lib_size == min(Small_Y_xmap_Y_supp[[j]]$lib_size))
Small_A_xmap_A_interaction_supp[[j]] <- subset(Small_A_xmap_A_supp[[j]]$rho, Small_A_xmap_A_supp[[j]]$lib_size == max(Small_A_xmap_A_supp[[j]]$lib_size)) - subset(Small_A_xmap_A_supp[[j]]$rho, Small_A_xmap_A_supp[[j]]$lib_size == min(Small_A_xmap_A_supp[[j]]$lib_size))
#convert negative values into zero
Small_Y_xmap_A_interaction_supp[[j]][Small_Y_xmap_A_interaction_supp[[j]] < 0] <- 0
Small_A_xmap_Y_interaction_supp[[j]][Small_A_xmap_Y_interaction_supp[[j]] < 0] <- 0
Small_Y_xmap_Y_interaction_supp[[j]][Small_Y_xmap_Y_interaction_supp[[j]] < 0] <- 0
Small_A_xmap_A_interaction_supp[[j]][Small_A_xmap_A_interaction_supp[[j]] < 0] <- 0
#Standardizing intraspecific effect
Small_fromYtoA_interaction_supp[[j]] <- Small_A_xmap_Y_interaction_supp[[j]]/Small_Y_xmap_Y_interaction_supp[[j]]
Small_fromAtoY_interaction_supp[[j]] <- Small_Y_xmap_A_interaction_supp[[j]]/Small_A_xmap_A_interaction_supp[[j]]
}
small_interaction_supp <- data.frame(tp = c(rep(0, 100), rep(-1, 100), rep(-2, 100), rep(-3, 100), rep(-4, 100)), fromYtoA = unlist(Small_fromYtoA_interaction_supp), fromAtoY = unlist(Small_fromAtoY_interaction_supp))
boxplot(
fromYtoA ~ tp, data = small_interaction_supp,
outline = TRUE,
main = "from Cm (Y) to Cc (A)"
)
boxplot(
fromAtoY ~ tp, data = small_interaction_supp,
outline = TRUE,
main = "from Cc (A) to Cm (Y)"
)
volume_combined <- c(rep(Smallvolume, 6), rep(Mediumvolume, 7), rep(Largevolume, 8), rep(Hugevolume, 8))
forFigS5_data <- data.frame(volume = volume_combined, Table1)
#GLM with Poisson for FigS5a
summary(glm(E_for_C.m ~ volume, data = forFigS5_data, family = poisson(link = "log")))
Call:
glm(formula = E_for_C.m ~ volume, family = poisson(link = "log"),
data = forFigS5_data)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.5403 -0.6068 -0.1918 0.5749 1.3519
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.876e+00 1.264e-01 14.838 <2e-16 ***
volume -1.310e-05 7.802e-05 -0.168 0.867
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 19.365 on 28 degrees of freedom
Residual deviance: 19.337 on 27 degrees of freedom
AIC: 129.82
Number of Fisher Scoring iterations: 4
#GLM with Poisson for FigS5b
summary(glm(E_for_C.c ~ volume, data = forFigS5_data, family = poisson(link = "log")))
Call:
glm(formula = E_for_C.c ~ volume, family = poisson(link = "log"),
data = forFigS5_data)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.71888 -0.77883 0.01921 0.74371 1.49328
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.779e+00 1.302e-01 13.664 <2e-16 ***
volume 6.244e-05 7.767e-05 0.804 0.421
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 32.855 on 28 degrees of freedom
Residual deviance: 32.209 on 27 degrees of freedom
AIC: 141.34
Number of Fisher Scoring iterations: 4
#GLM with Poisson for FigS5c
summary(glm(delta.rho_for_C.m ~ volume, data = forFigS5_data, family = gaussian))
Call:
glm(formula = delta.rho_for_C.m ~ volume, family = gaussian,
data = forFigS5_data)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.13638 -0.06703 -0.01951 0.02439 0.42466
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.408e-01 3.417e-02 4.121 0.000321 ***
volume -2.479e-05 2.096e-05 -1.183 0.247278
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 0.01126655)
Null deviance: 0.31995 on 28 degrees of freedom
Residual deviance: 0.30420 on 27 degrees of freedom
AIC: -43.865
Number of Fisher Scoring iterations: 2
#GLM with Poisson for FigS5d
summary(glm(delta.rho_for_C.c ~ volume, data = forFigS5_data, family = gaussian))
Call:
glm(formula = delta.rho_for_C.c ~ volume, family = gaussian,
data = forFigS5_data)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.19262 -0.10579 -0.02556 0.04322 0.51840
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.029e-01 5.046e-02 4.021 0.000418 ***
volume -5.770e-05 3.096e-05 -1.864 0.073260 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 0.02457561)
Null deviance: 0.74891 on 28 degrees of freedom
Residual deviance: 0.66354 on 27 degrees of freedom
AIC: -21.248
Number of Fisher Scoring iterations: 2
#GLM with Poisson for FigS5e
summary(glm(theta_for_C.m ~ volume, data = forFigS5_data, family = gaussian))
Call:
glm(formula = theta_for_C.m ~ volume, family = gaussian, data = forFigS5_data)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.5220 -1.3332 -0.4412 0.6668 5.6980
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.2854572 0.6515364 3.508 0.0016 **
volume 0.0000925 0.0003997 0.231 0.8187
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 4.096701)
Null deviance: 110.83 on 28 degrees of freedom
Residual deviance: 110.61 on 27 degrees of freedom
AIC: 127.12
Number of Fisher Scoring iterations: 2
#GLM with Poisson for FigS5f
summary(glm(theta_for_C.c ~ volume, data = forFigS5_data, family = gaussian))
Call:
glm(formula = theta_for_C.c ~ volume, family = gaussian, data = forFigS5_data)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.3617 -2.1961 -0.3617 1.0539 5.1227
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.8409920 0.9154855 3.103 0.00445 **
volume 0.0002036 0.0005616 0.362 0.71980
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 8.088349)
Null deviance: 219.45 on 28 degrees of freedom
Residual deviance: 218.39 on 27 degrees of freedom
AIC: 146.85
Number of Fisher Scoring iterations: 2
#Simple plot with GLM line for FigS5a
FigS5a <- ggplot(data = forFigS5_data, aes(x = volume, y= E_for_C.m))
FigS5a <- FigS5a + geom_point(size = 3, shape = 21)
FigS5a <- FigS5a + geom_smooth(method = 'glm', method.args = list(family = 'poisson'), se = T, color = 'black', linetype = "dashed")
FigS5a <- FigS5a + labs(y = "E*", x = "Volume (cm3)") + ggtitle("Fig. S5a C.m")
#Simple plot with GLM line for FigS5b
FigS5b <- ggplot(data = forFigS5_data, aes(x = volume, y= E_for_C.c))
FigS5b <- FigS5b + geom_point(size = 3, shape = 21)
FigS5b <- FigS5b + geom_smooth(method = 'glm', method.args = list(family = 'poisson'), se = T, color = 'black', linetype = "dashed")
#Simple plot with GLM line for FigS5c
FigS5c <- ggplot(data = forFigS5_data, aes(x = volume, y= delta.rho_for_C.m))
FigS5c <- FigS5c + geom_smooth(method = 'lm', se = T, color = 'black', linetype = "dashed")
FigS5c <- FigS5c + labs(y = "delta-rho", x = "Volume (cm3)") + ggtitle("Fig. S5c C.m")
#Simple plot with GLM line for FigS5d
FigS5d <- ggplot(data = forFigS5_data, aes(x = volume, y= delta.rho_for_C.c))
FigS5d <- FigS5d + geom_smooth(method = 'lm', se = T, color = 'black', linetype = "dashed")
#Simple plot with GLM line for FigS5e
FigS5e <- ggplot(data = forFigS5_data, aes(x = volume, y= theta_for_C.m))
FigS5e <- FigS5e + geom_point(size = 3, shape = 21)
FigS5e <- FigS5e + geom_smooth(method = 'lm', se = T, color = 'black', linetype = "dashed")
FigS5e <- FigS5e + labs(y = "theta*", x = "Volume (cm3)") + ggtitle("Fig. S5e C.m")
#Simple plot with GLM line for FigS5f
FigS5f <- FigS5f + geom_point(size = 3, shape = 21)
FigS5f <- FigS5f + geom_smooth(method = 'lm', se = T, color = 'black', linetype = "dashed")
FigS5f <- FigS5f + labs(y = "theta*", x = "Volume (cm3)") + ggtitle("Fig. S5f C.c")
print(FigS5a)
print(FigS5b)
print(FigS5c)
print(FigS5d)
print(FigS5e)
print(FigS5f)
#Data summary
Huge_rep <- list()
for(i in 1:length(Huge)){
Huge_rep[[i]] <- data.frame("Time_step" = c(1:length(Huge[[i]]$YotsumonTotalL)),
"Population_size_C.m" = Huge[[i]]$YotsumonTotalL,
"Population_size_C.c" = Huge[[i]]$AzukiTotalL
)
}
#Graphics
FigS6panellabs <- c("(a)", "(b)", "(c)", "(d)", "(e)", "(f)", "(g)", "(h)", "(i)", "(j)", "(k)", "(l)", "(m)", "(n)", "(o)", "(p)")
FigS6 <- list()
for(i in 1:length(Huge_rep)){
FigS6[[i]] <- ggplot(data = Huge_rep[[i]], aes(x = Time_step, y = Population_size_C.m)) +
geom_line(lty = 1) +
geom_line(data = Huge_rep[[i]], aes(x = Time_step, y = Population_size_C.c), lty = 2) +
ylab("Population size") +
ggtitle(FigS6panellabs[i])
print(FigS6[[i]])
}
Small_cv_Y <- c()
Small_cv_A <- c()
Medium_cv_Y <- c()
Medium_cv_A <- c()
Large_cv_Y <- c()
Large_cv_A <- c()
Huge_cv_Y <- c()
Huge_cv_A <- c()
for(i in 1:length(Small)){
Small_cv_Y[i] <- sd(Small[[i]]$YotsumonTotalL, na.rm = T)/mean(Small[[i]]$YotsumonTotalL, na.rm = T)
Small_cv_A[i] <- sd(Small[[i]]$AzukiTotalL, na.rm = T)/mean(Small[[i]]$AzukiTotalL, na.rm = T)
}
for(i in 1:length(Medium)){
Medium_cv_Y[i] <- sd(Medium[[i]]$YotsumonTotalL, na.rm = T)/mean(Medium[[i]]$YotsumonTotalL, na.rm = T)
Medium_cv_A[i] <- sd(Medium[[i]]$AzukiTotalL, na.rm = T)/mean(Medium[[i]]$AzukiTotalL, na.rm = T)
}
for(i in 1:length(Large)){
Large_cv_Y[i] <- sd(Large[[i]]$YotsumonTotalL, na.rm = T)/mean(Large[[i]]$YotsumonTotalL, na.rm = T)
Large_cv_A[i] <- sd(Large[[i]]$AzukiTotalL, na.rm = T)/mean(Large[[i]]$AzukiTotalL, na.rm = T)
}
for(i in 1:length(Huge)){
Huge_cv_Y[i] <- sd(Huge[[i]]$YotsumonTotalL, na.rm = T)/mean(Huge[[i]]$YotsumonTotalL, na.rm = T)
Huge_cv_A[i] <- sd(Huge[[i]]$AzukiTotalL, na.rm = T)/mean(Huge[[i]]$AzukiTotalL, na.rm = T)
}
–For C.m (Y), slope was significant when Small, Large, and Huge
–For C.c (A), slope was significant when Small and Large
#GLM for each cage size
#for C.m (Y)
summary(glm(subset(TableS1$Coexistence_periods, TableS1$Cage.Size == "Small") ~ Small_cv_Y, family = "poisson"))$coefficients
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.165355 0.2865541 14.53601 7.164409e-48
Small_cv_Y -1.283970 0.2516943 -5.10131 3.373112e-07
summary(glm(subset(TableS1$Coexistence_periods, TableS1$Cage.Size == "Medium") ~ Medium_cv_Y, family = "poisson"))$coefficients
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.1117072 0.3194051 9.7421975 1.991994e-22
Medium_cv_Y -0.0956626 0.2670559 -0.3582119 7.201847e-01
summary(glm(subset(TableS1$Coexistence_periods, TableS1$Cage.Size == "Large") ~ Large_cv_Y, family = "poisson"))$coefficients
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.2033482 0.4027249 5.471099 4.472522e-08
Large_cv_Y 0.9040723 0.3819684 2.366877 1.793889e-02
summary(glm(subset(TableS1$Coexistence_periods, TableS1$Cage.Size == "Huge") ~ Huge_cv_Y, family = "poisson"))$coefficients
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.8271170 0.1856177 15.230858 2.20678e-52
Huge_cv_Y 0.3243216 0.1412992 2.295283 2.17169e-02
#for C.c (A)
summary(glm(subset(TableS1$Coexistence_periods, TableS1$Cage.Size == "Small") ~ Small_cv_A, family = "poisson"))$coefficients
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.0018796 0.4561198 8.773747 1.728173e-18
Small_cv_A -0.7369437 0.2556813 -2.882274 3.948157e-03
summary(glm(subset(TableS1$Coexistence_periods, TableS1$Cage.Size == "Medium") ~ Medium_cv_A, family = "poisson"))$coefficients
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.10902401 0.4563374 6.8129934 9.558855e-12
Medium_cv_A -0.06392303 0.2630055 -0.2430483 8.079680e-01
summary(glm(subset(TableS1$Coexistence_periods, TableS1$Cage.Size == "Large") ~ Large_cv_A, family = "poisson"))$coefficients
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.1418655 0.3670250 5.835749 5.354957e-09
Large_cv_A 0.5798425 0.2087438 2.777771 5.473315e-03
summary(glm(subset(TableS1$Coexistence_periods, TableS1$Cage.Size == "Huge") ~ Huge_cv_A, family = "poisson"))$coefficients
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.4606243 0.3208369 10.7862421 3.998052e-27
Huge_cv_A -0.1433893 0.1987414 -0.7214868 4.706100e-01
#GLM for all cage sizes and prediction for each cage
cv_Y_all <- c(Small_cv_Y, Medium_cv_Y, Large_cv_Y, Huge_cv_Y)
cv_A_all <- c(Small_cv_A, Medium_cv_A, Large_cv_A, Huge_cv_A)
cv_data_all <- data.frame(Coexistence_periods = TableS1$Coexistence_periods, volume = TableS1$Cage.Size, cv_Y = cv_Y_all, cv_A = cv_A_all)
cv_coexistence_model_Y <- glm(Coexistence_periods ~ volume*cv_Y, data = cv_data_all, family = "poisson")
cv_coexistence_model_Y_predicted <- ggpredict(cv_coexistence_model_Y, term = c("cv_Y","volume"))
cv_coexistence_model_A <- glm(Coexistence_periods ~ volume*cv_A, data = cv_data_all, family = "poisson")
cv_coexistence_model_A_predicted <- ggpredict(cv_coexistence_model_A, term = c("cv_A","volume"))
#define linetypes & linewidth
linetypes_Y <- c("solid", "dashed", "solid", "solid")
linewidths_Y <- c(0.5, 0.5, 0.75, 1)
linetypes_A <- c("solid", "dashed", "solid", "dotted")
linewidths_A <- c(0.5, 0.5, 1, 0.5)
#ggplot
cv_coexist_Y_fig <- ggplot() + geom_point(data = cv_data_all, aes(x = cv_Y, y = Coexistence_periods, shape = as.factor(volume)), size = 3, alpha = 0.5)
cv_coexist_Y_fig <- cv_coexist_Y_fig + geom_line(data = cv_coexistence_model_Y_predicted, aes(x = x, y = predicted, linetype = group, linewidth = group)) + scale_linetype_manual(values = linetypes_Y)
cv_coexist_Y_fig <- cv_coexist_Y_fig + scale_linewidth_manual(values = linewidths_Y)
cv_coexist_Y_fig <- cv_coexist_Y_fig + ggtitle("CV vs coexistence for C.m (Y)")
plot(cv_coexist_Y_fig)
cv_coexist_A_fig <- ggplot() + geom_point(data = cv_data_all, aes(x = cv_A, y = Coexistence_periods, shape = as.factor(volume)), size = 3, alpha = 0.5)
cv_coexist_A_fig <- cv_coexist_A_fig + geom_line(data = cv_coexistence_model_A_predicted, aes(x = x, y = predicted, linetype = group, linewidth = group)) + scale_linetype_manual(values = linetypes_A)
cv_coexist_A_fig <- cv_coexist_A_fig + scale_linewidth_manual(values = linewidths_A)
cv_coexist_A_fig <- cv_coexist_A_fig + ggtitle("CV vs coexistence for C.c (A)")
plot(cv_coexist_A_fig)
NA
Small_5th_Y <- c()
Small_5th_A <- c()
Medium_5th_Y <- c()
Medium_5th_A <- c()
Large_5th_Y <- c()
Large_5th_A <- c()
Huge_5th_Y <- c()
Huge_5th_A <- c()
for(i in 1:16) {
Small_5th_Y[i] <- Small[[i]]$YotsumonTotalL[5]
Small_5th_A[i] <- Small[[i]]$AzukiTotalL[5]
Medium_5th_Y[i] <- Medium[[i]]$YotsumonTotalL[5]
Medium_5th_A[i] <- Medium[[i]]$AzukiTotalL[5]
Large_5th_Y[i] <- Large[[i]]$YotsumonTotalL[5]
Large_5th_A[i] <- Large[[i]]$AzukiTotalL[5]
Huge_5th_Y[i] <- Huge[[i]]$YotsumonTotalL[5]
Huge_5th_A[i] <- Huge[[i]]$AzukiTotalL[5]
}
Small_5th_coexist <- data.frame(Coexistence_periods = subset(TableS1$Coexistence_periods, TableS1$Cage.Size == "Small"), Y_5th = Small_5th_Y, A_5th = Small_5th_A)
Medium_5th_coexist <- data.frame(Coexistence_periods = subset(TableS1$Coexistence_periods, TableS1$Cage.Size == "Medium"), Y_5th = Medium_5th_Y, A_5th = Medium_5th_A)
Large_5th_coexist <- data.frame(Coexistence_periods = subset(TableS1$Coexistence_periods, TableS1$Cage.Size == "Large"), Y_5th = Large_5th_Y, A_5th = Large_5th_A)
Huge_5th_coexist <- data.frame(Coexistence_periods = subset(TableS1$Coexistence_periods, TableS1$Cage.Size == "Huge"), Y_5th = Huge_5th_Y, A_5th = Huge_5th_A)
summary(glm(Coexistence_periods ~ Y_5th, data = Small_5th_coexist, family = "poisson"))
Call:
glm(formula = Coexistence_periods ~ Y_5th, family = "poisson",
data = Small_5th_coexist)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.4496 -1.1968 -0.2544 0.6242 2.3607
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.755182 0.455821 8.238 <2e-16 ***
Y_5th -0.009534 0.004056 -2.350 0.0188 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 36.944 on 15 degrees of freedom
Residual deviance: 31.229 on 14 degrees of freedom
AIC: 106.54
Number of Fisher Scoring iterations: 4
summary(glm(Coexistence_periods ~ A_5th, data = Small_5th_coexist, family = "poisson"))
Call:
glm(formula = Coexistence_periods ~ A_5th, family = "poisson",
data = Small_5th_coexist)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.2116 -0.9607 -0.3643 1.1848 2.2719
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.023925 0.202592 9.990 < 2e-16 ***
A_5th 0.014068 0.003931 3.579 0.000345 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 36.944 on 15 degrees of freedom
Residual deviance: 24.094 on 14 degrees of freedom
AIC: 99.403
Number of Fisher Scoring iterations: 4
summary(glm(Coexistence_periods ~ Y_5th, data = Medium_5th_coexist, family = "poisson"))
glm(formula = Coexistence_periods ~ Y_5th, family = "poisson",
data = Medium_5th_coexist)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.7181 -1.4484 -0.3840 0.5927 4.9394
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.408069 0.337243 10.106 <2e-16 ***
Y_5th -0.003848 0.003145 -1.223 0.221
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 63.081 on 15 degrees of freedom
Residual deviance: 61.594 on 14 degrees of freedom
AIC: 141.59
Number of Fisher Scoring iterations: 4
summary(glm(Coexistence_periods ~ A_5th, data = Medium_5th_coexist, family = "poisson"))
Call:
glm(formula = Coexistence_periods ~ A_5th, family = "poisson",
data = Medium_5th_coexist)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.9681 -1.4000 -0.6540 0.6791 5.1097
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.533221 0.144071 17.583 < 2e-16 ***
A_5th 0.006387 0.001743 3.665 0.000247 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 63.081 on 15 degrees of freedom
Residual deviance: 49.496 on 14 degrees of freedom
AIC: 129.49
Number of Fisher Scoring iterations: 4
summary(glm(Coexistence_periods ~ Y_5th, data = Large_5th_coexist, family = "poisson"))
Call:
glm(formula = Coexistence_periods ~ Y_5th, family = "poisson",
data = Large_5th_coexist)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.56714 -0.86804 -0.05889 0.55035 2.74155
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.451380 0.316252 14.075 < 2e-16 ***
Y_5th -0.011238 0.002723 -4.127 3.68e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 42.624 on 15 degrees of freedom
Residual deviance: 25.466 on 14 degrees of freedom
AIC: 108.36
Number of Fisher Scoring iterations: 4
summary(glm(Coexistence_periods ~ A_5th, data = Large_5th_coexist, family = "poisson"))
Call:
glm(formula = Coexistence_periods ~ A_5th, family = "poisson",
data = Large_5th_coexist)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.0574 -0.9449 -0.5074 0.3003 3.9059
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.889068 0.125974 22.934 <2e-16 ***
A_5th 0.003887 0.001724 2.254 0.0242 *
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
family taken to be 1)
Null deviance: 42.624 on 15 degrees of freedom
Residual deviance: 37.584 on 14 degrees of freedom
AIC: 120.48
Number of Fisher Scoring iterations: 4
summary(glm(Coexistence_periods ~ Y_5th, data = Huge_5th_coexist, family = "poisson"))
Call:
glm(formula = Coexistence_periods ~ Y_5th, family = "poisson",
data = Huge_5th_coexist)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.1849 -0.8404 0.3302 0.7996 1.8855
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.409885 0.245346 13.90 <2e-16 ***
Y_5th -0.001774 0.002396 -0.74 0.459
---
Signif. codes: poisson family taken to be 1)
Null deviance: 23.991 on 15 degrees of freedom
Residual deviance: 23.439 on 14 degrees of freedom
AIC: 108.14
Number of Fisher Scoring iterations: 4
summary(glm(Coexistence_periods ~ A_5th, data = Huge_5th_coexist, family = "poisson"))
Call:
glm(formula = Coexistence_periods ~ A_5th, family = "poisson",
data = Huge_5th_coexist)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.06108 -0.46498 -0.04061 0.83194 1.87350
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.009259 0.151495 19.864 <2e-16 ***
A_5th 0.003172 0.002017 1.573 0.116
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1Residual deviance: 21.550 on 14 degrees of freedom
AIC: 106.25
Number of Fisher Scoring iterations: 4
#Small
coexist_5th_Small_fig <- ggplot() + geom_point(data = Small_5th_coexist, aes(x = Y_5th, y = Coexistence_periods), shape = 1, size = 3, alpha = 0.5) + geom_point(data = Small_5th_coexist, aes(x = A_5th, y = Coexistence_periods), shape = 2, size = 3, alpha = 0.5) + geom_smooth(data = Small_5th_coexist, aes(x = Y_5th, y = Coexistence_periods), method = 'glm', method.args = list(family = 'poisson'), color = "black", se = T) + geom_smooth(data = Small_5th_coexist, aes(x = A_5th, y = Coexistence_periods), method = 'glm', method.args = list(family = 'poisson'), color = "black", se = T)
coexist_5th_Small_fig <- coexist_5th_Small_fig + xlab("population size") + ggtitle("Early poplation size vs coexistence (Small)")
#Medium
coexist_5th_Medium_fig <- ggplot() + geom_point(data = Medium_5th_coexist, aes(x = Y_5th, y = Coexistence_periods), shape = 1, size = 3, alpha = 0.5) + geom_point(data = Medium_5th_coexist, aes(x = A_5th, y = Coexistence_periods), shape = 2, size = 3, alpha = 0.5) + geom_smooth(data = Medium_5th_coexist, aes(x = Y_5th, y = Coexistence_periods), method = 'glm', method.args = list(family = 'poisson'), color = "black", linetype = "dashed", se = T) + geom_smooth(data = Medium_5th_coexist, aes(x = A_5th, y = Coexistence_periods), method = 'glm', method.args = list(family = 'poisson'), color = "black", se = T)
coexist_5th_Medium_fig <- coexist_5th_Medium_fig + xlab("population size") + ggtitle("Early poplation size vs coexistence (Medium)")
#Large
coexist_5th_Large_fig <- ggplot() + geom_point(data = Large_5th_coexist, aes(x = Y_5th, y = Coexistence_periods), shape = 1, size = 3, alpha = 0.5) + geom_point(data = Large_5th_coexist, aes(x = A_5th, y = Coexistence_periods), shape = 2, size = 3, alpha = 0.5) + geom_smooth(data = Large_5th_coexist, aes(x = Y_5th, y = Coexistence_periods), method = 'glm', method.args = list(family = 'poisson'), color = "black", se = T) + geom_smooth(data = Large_5th_coexist, aes(x = A_5th, y = Coexistence_periods), method = 'glm', method.args = list(family = 'poisson'), color = "black", se = T)
coexist_5th_Large_fig <- coexist_5th_Large_fig + xlab("population size") + ggtitle("Early poplation size vs coexistence (Large)")
#Huge
coexist_5th_Huge_fig <- ggplot() + geom_point(data = Huge_5th_coexist, aes(x = Y_5th, y = Coexistence_periods), shape = 1, size = 3, alpha = 0.5) + geom_point(data = Huge_5th_coexist, aes(x = A_5th, y = Coexistence_periods), shape = 2, size = 3, alpha = 0.5) + geom_smooth(data = Huge_5th_coexist, aes(x = Y_5th, y = Coexistence_periods), method = 'glm', method.args = list(family = 'poisson'), color = "black", linetype = "dashed", se = T) + geom_smooth(data = Huge_5th_coexist, aes(x = A_5th, y = Coexistence_periods), method = 'glm', method.args = list(family = 'poisson'), color = "black", linetype = "dashed", se = T)
coexist_5th_Huge_fig <- coexist_5th_Huge_fig + xlab("population size") + ggtitle("Early poplation size vs coexistence (Huge)")
print(coexist_5th_Small_fig)
print(coexist_5th_Medium_fig)
print(coexist_5th_Large_fig)
print(coexist_5th_Huge_fig)