#################################################################### ## Predicting the effect of habitat modification on networks of ## interacting species ## ## Please cite as: ## Staniczenko, P.P.A., Lewis, O.T., Tylianakis, J.M., Albrecht, M., Coudrain, V., Klein, A.-M. & Reed-Tsochas, F. 2017 ## Predicting the effect of habitat modification on networks of interacting species ## Nature Communications, 8, 792 ## ## Phillip P.A. Staniczenko (pstaniczenko@sesync.org) #################################################################### require(bipartite) require(limSolve) rm(list=ls(all=TRUE)) # Clear workspace # FUNCTIONS getaggregatenetwork <- function(G){ # return an aggregate network of the networks in G for (i in 1:length(G)){ network <- G[[i]] if (i == 1){ A <- network } else{ A <- A + network } } return(A) } getswitches <- function(Network_aggregate_G1, Network_aggregate_G2){ # return binary incidence matrix of switches in G2 Network_aggregate_G1_binary <- Network_aggregate_G1 Network_aggregate_G1_binary[Network_aggregate_G1_binary > 0] <- 1 Networkaggregate_G2_binary <- Network_aggregate_G2 Networkaggregate_G2_binary[Networkaggregate_G2_binary > 0] <- 1 Switches_G2 <- Networkaggregate_G2_binary - Network_aggregate_G1_binary Switches_G2[Switches_G2 > 0] <- 1 Switches_G2[Switches_G2 < 0] <- 0 return(Switches_G2) } preparedata <- function(G){ # return data in list form for (i in 1:length(G)){ network <- G[[i]] if (i == 1){ y <- network[which(network > 0)] y_network <- rep(i, length(y)) xs_index <- data.frame(which(network > 0, arr.ind=TRUE, useNames=FALSE)[, 1], which(network > 0, arr.ind=TRUE, useNames=FALSE)[, 2]) colnames(xs_index) <- c("Host", "Parasitoid") } else{ yy <- network[which(network > 0)] yy_network <- rep(i, length(yy)) xsxs_index <- data.frame(which(network > 0, arr.ind=TRUE, useNames=FALSE)[, 1], which(network > 0, arr.ind=TRUE, useNames=FALSE)[, 2]) colnames(xsxs_index) <- c("Host", "Parasitoid") y <- c(y, yy) y_network <- c(y_network, yy_network) xs_index <- rbind(xs_index, xsxs_index) } prepG <- new.env() prepG$y <- y prepG$y_network <- y_network prepG$xs_index <- xs_index prepG <- as.list(prepG) } return(prepG) } # DATA # Example networks a1 <- matrix(c(3,0,5,2,9,8), nrow=2, ncol=3) a2 <- matrix(c(5,0,0,2,7,7), nrow=2, ncol=3) a3 <- matrix(c(1,0,3,4,10,8), nrow=2, ncol=3) a4 <- matrix(c(1,0,6,2,9,0), nrow=2, ncol=3) Ga <- list(a1,a2,a3,a4) b1 <- matrix(c(8,2,0,4,20,8), nrow=2, ncol=3) b2 <- matrix(c(7,1,0,3,21,6), nrow=2, ncol=3) b3 <- matrix(c(8,2,0,5,18,6), nrow=2, ncol=3) b4 <- matrix(c(8,2,0,5,15,5), nrow=2, ncol=3) Gb <- list(b1,b2,b3,b4) # VARIABLES # For Ga alpha_Ga <- 1 beta_Ga <- 1 # For Gb alpha_Gb <- 1 beta_Gb <- 1 # For correlated preferences model delta <- 0.6 # MAIN # get aggregate webs Network_aggregate_Ga <- getaggregatenetwork(Ga) Network_aggregate_Gb <- getaggregatenetwork(Gb) # get switches Switches_Gb <- getswitches(Network_aggregate_Ga, Network_aggregate_Gb) Switches_Ga <- getswitches(Network_aggregate_Gb, Network_aggregate_Ga) # prepare data prep_Ga <- preparedata(Ga) prep_Gb <- preparedata(Gb) # get log-likelihood of mass action for Gb K <- mat.or.vec(nr=length(prep_Gb$y), nc=(dim(Network_aggregate_Gb)[1] + dim(Network_aggregate_Gb)[2])) y_log <- mat.or.vec(length(prep_Gb$y), nc=1) for (i in 1:length(prep_Gb$y)){ y_log[i] <- log(prep_Gb$y[i]) K[i, prep_Gb$xs_index[i, 1]] <- alpha_Gb K[i, prep_Gb$xs_index[i, 2] + dim(Network_aggregate_Gb)[1]] <- beta_Gb } inverseoutput <- lsei(K, y_log, verbose=FALSE) # solve the inverse problem estden_Gb <- exp(inverseoutput$X) probvector <- mat.or.vec(nr=length(prep_Gb$y), nc=1) for (i in 1:length(probvector)){ probvector[i] <- estden_Gb[prep_Gb$xs_index[i, 1]] * estden_Gb[dim(Network_aggregate_Gb)[1] + prep_Gb$xs_index[i, 2]] # mass action model } loglikelihood_Gb_massaction <- dmultinom(prep_Gb$y, size=NULL, probvector, log=TRUE) # get preference matrix for Gb K <- mat.or.vec(nr=length(prep_Gb$y), nc=(dim(Network_aggregate_Gb)[1] + dim(Network_aggregate_Gb)[2])) y_log <- mat.or.vec(length(prep_Gb$y), nc=1) for (i in 1:length(prep_Gb$y)){ y_log[i] <- log(prep_Gb$y[i]) K[i, prep_Gb$xs_index[i, 1]] <- alpha_Gb K[i, prep_Gb$xs_index[i, 2] + dim(Network_aggregate_Gb)[1]] <- beta_Gb } inverseoutput <- lsei(K, y_log, verbose=FALSE) # solve the inverse problem estden_Gb <- exp(inverseoutput$X) probvector <- mat.or.vec(nr=length(prep_Gb$y), nc=1) for (i in 1:length(probvector)){ probvector[i] <- estden_Gb[prep_Gb$xs_index[i, 1]] * estden_Gb[dim(Network_aggregate_Gb)[1] + prep_Gb$xs_index[i, 2]] # mass action model } factor <- sum(prep_Gb$y) / sum(probvector) probvector <- probvector * factor # scale estimated densities to match empirical counts prefs_networklevel <- prep_Gb$y/probvector prefs_grouplevel_matrix <- mat.or.vec(dim(Network_aggregate_Gb)[1], nc=dim(Network_aggregate_Gb)[2]) occurrenceinnetworks_matrix <- mat.or.vec(dim(Network_aggregate_Gb)[1], nc=dim(Network_aggregate_Gb)[2]) for (i in 1:length(prefs_networklevel)){ prefs_grouplevel_matrix[prep_Gb$xs_index[i, 1], prep_Gb$xs_index[i, 2]] <- prefs_grouplevel_matrix[prep_Gb$xs_index[i, 1], prep_Gb$xs_index[i, 2]] + prefs_networklevel[i] occurrenceinnetworks_matrix[prep_Gb$xs_index[i, 1], prep_Gb$xs_index[i, 2]] <- occurrenceinnetworks_matrix[prep_Gb$xs_index[i, 1], prep_Gb$xs_index[i, 2]] + 1 } prefs_grouplevel_matrix <- prefs_grouplevel_matrix / occurrenceinnetworks_matrix prefs_grouplevel_matrix[is.nan(prefs_grouplevel_matrix)] <- 0 prefs_grouplevel_matrix_Gb <- prefs_grouplevel_matrix # get log-likelihood of complete characterisation model for Gb probvector <- mat.or.vec(nr=length(prep_Gb$y), nc=1) for (i in 1:length(probvector)){ probvector[i] <- estden_Gb[prep_Gb$xs_index[i, 1]] * estden_Gb[dim(Network_aggregate_Gb)[1] + prep_Gb$xs_index[i, 2]] * prefs_grouplevel_matrix_Gb[prep_Gb$xs_index[i, 1], prep_Gb$xs_index[i, 2]] } loglikelihood_Gb_complete <- dmultinom(prep_Gb$y, size=NULL, probvector, log=TRUE) # get log-likelihood of null model for Gb loglikelihood_Gb_null <- dmultinom(prep_Gb$y, size=NULL, rep(1, length(prep_Gb$y)), log=TRUE) print("Group B---log-likelihood for null model:") print(loglikelihood_Gb_null) # get log-likelihood of aggregate counts model for Gb probvector <- mat.or.vec(nr=length(prep_Gb$y), nc=1) for (i in 1:length(probvector)){ probvector[i] <- Network_aggregate_Ga[prep_Gb$xs_index[i, 1], prep_Gb$xs_index[i, 2]] } for (i in 1:length(probvector)){ if (Network_aggregate_Ga[prep_Gb$xs_index[i, 1], prep_Gb$xs_index[i, 2]] == 0){ probvector[i] <- Switches_Gb[prep_Gb$xs_index[i, 1], prep_Gb$xs_index[i, 2]] # mass action switches } } loglikelihood_Gb_aggcounts <- dmultinom(prep_Gb$y, size=NULL, probvector, log=TRUE) print("Group B---log-likelihood for aggregate counts model:") print(loglikelihood_Gb_aggcounts) # print log-likelihood for random encounter model (which is the same as mass action, above) print("Group B---log-likelihood for random encounter model (mass action):") print(loglikelihood_Gb_massaction) # get log-likelihood for alternative preferences model with mass action switches # get preference matrix for Ga K <- mat.or.vec(nr=length(prep_Ga$y), nc=(dim(Network_aggregate_Ga)[1] + dim(Network_aggregate_Ga)[2])) y_log <- mat.or.vec(length(prep_Ga$y), nc=1) for (i in 1:length(prep_Ga$y)){ y_log[i] <- log(prep_Ga$y[i]) K[i, prep_Ga$xs_index[i, 1]] <- alpha_Ga K[i, prep_Ga$xs_index[i, 2] + dim(Network_aggregate_Ga)[1]] <- beta_Ga } inverseoutput <- lsei(K, y_log, verbose=FALSE) # solve the inverse problem estden_Ga <- exp(inverseoutput$X) probvector <- mat.or.vec(nr=length(prep_Ga$y), nc=1) for (i in 1:length(probvector)){ probvector[i] <- estden_Ga[prep_Ga$xs_index[i, 1]] * estden_Ga[dim(Network_aggregate_Ga)[1] + prep_Ga$xs_index[i, 2]] # mass action model } factor <- sum(prep_Ga$y) / sum(probvector) probvector <- probvector * factor # scale estimated densities to match empirical counts prefs_networklevel <- prep_Ga$y/probvector prefs_grouplevel_matrix <- mat.or.vec(dim(Network_aggregate_Ga)[1], nc=dim(Network_aggregate_Ga)[2]) occurrenceinnetworks_matrix <- mat.or.vec(dim(Network_aggregate_Ga)[1], nc=dim(Network_aggregate_Ga)[2]) for (i in 1:length(prefs_networklevel)){ prefs_grouplevel_matrix[prep_Ga$xs_index[i, 1], prep_Ga$xs_index[i, 2]] <- prefs_grouplevel_matrix[prep_Ga$xs_index[i, 1], prep_Ga$xs_index[i, 2]] + prefs_networklevel[i] occurrenceinnetworks_matrix[prep_Ga$xs_index[i, 1], prep_Ga$xs_index[i, 2]] <- occurrenceinnetworks_matrix[prep_Ga$xs_index[i, 1], prep_Ga$xs_index[i, 2]] + 1 } prefs_grouplevel_matrix <- prefs_grouplevel_matrix / occurrenceinnetworks_matrix prefs_grouplevel_matrix[is.nan(prefs_grouplevel_matrix)] <- 0 prefs_grouplevel_matrix_Ga <- prefs_grouplevel_matrix # probvector <- mat.or.vec(nr=length(prep_Gb$y), nc=1) for (i in 1:length(probvector)){ probvector[i] <- estden_Gb[prep_Gb$xs_index[i, 1]] * estden_Gb[dim(Network_aggregate_Gb)[1] + prep_Gb$xs_index[i, 2]] * prefs_grouplevel_matrix_Ga[prep_Gb$xs_index[i, 1], prep_Gb$xs_index[i, 2]] } for (i in 1:length(probvector)){ if (prefs_grouplevel_matrix_Ga[prep_Gb$xs_index[i, 1], prep_Gb$xs_index[i, 2]] == 0){ probvector[i] <- estden_Gb[prep_Gb$xs_index[i, 1]] * estden_Gb[dim(Network_aggregate_Gb)[1] + prep_Gb$xs_index[i, 2]] * Switches_Gb[prep_Gb$xs_index[i, 1], prep_Gb$xs_index[i, 2]] # mass action switches } } loglikelihood_Gb_alternativeprefs <- dmultinom(prep_Gb$y, size=NULL, probvector, log=TRUE) print("Group B---log-likelihood for alternative preferences model:") print(loglikelihood_Gb_alternativeprefs) # get log-likelihood for correlated preferences model with mass action switches # get rank-order preferences for Gb prefs_grouplevel_matrix_Gb_temp <- prefs_grouplevel_matrix_Gb + Switches_Ga prefs_grouplevel_matrix_Gb_rankorder <- prefs_grouplevel_matrix_Gb_temp for (j in 1:dim(prefs_grouplevel_matrix_Gb_temp)[2]){ prefsset <- prefs_grouplevel_matrix_Gb_temp[prefs_grouplevel_matrix_Gb_temp[, j] > 0, j] presset <- unique(prefsset) prefsset <- sort(prefsset, decreasing=TRUE) for (i in 1:dim(prefs_grouplevel_matrix_Gb_rankorder)[1]){ if (prefs_grouplevel_matrix_Gb_rankorder[i, j] > 0){ for (k in 1:length(prefsset)){ if (prefs_grouplevel_matrix_Gb_temp[i, j]==prefsset[k]){ prefs_grouplevel_matrix_Gb_rankorder[i, j] <- k } } } } } # reorder preference matrix for Ga prefs_grouplevel_matrix_Ga_temp <- prefs_grouplevel_matrix_Ga + Switches_Gb prefs_grouplevel_matrix_Ga_reorder <- prefs_grouplevel_matrix_Ga_temp * 0.0 for (j in 1:dim(prefs_grouplevel_matrix_Ga_temp)[2]){ prefsset <- prefs_grouplevel_matrix_Ga_temp[prefs_grouplevel_matrix_Ga_temp[, j] > 0, j] prefsset <- sort(prefsset, decreasing=TRUE) prefs_grouplevel_matrix_Gb_rankorder_j <- prefs_grouplevel_matrix_Gb_rankorder[prefs_grouplevel_matrix_Gb_rankorder[, j] > 0, j] prefs_grouplevel_matrix_Gb_rankorder_j <- sort(prefs_grouplevel_matrix_Gb_rankorder_j, decreasing=FALSE) prefs_grouplevel_matrix_Gb_rankorder_j <- prefs_grouplevel_matrix_Gb_rankorder_j[1:length(prefsset)] prefs_grouplevel_matrix_Gb_rankorder_j_max <- max(prefs_grouplevel_matrix_Gb_rankorder_j) possible_elements <- which(prefs_grouplevel_matrix_Ga_temp[, j] > 0) for (i in 1:dim(prefs_grouplevel_matrix_Ga_temp)[1]){ if(prefs_grouplevel_matrix_Gb_rankorder[i, j] > 0){ myrankorder <- prefs_grouplevel_matrix_Gb_rankorder[i, j] if (myrankorder <= prefs_grouplevel_matrix_Gb_rankorder_j_max){ prefoptions <- prefsset[which(prefs_grouplevel_matrix_Gb_rankorder_j == myrankorder)] mypick <- sample(1:length(prefoptions), 1) prefs_grouplevel_matrix_Ga_reorder[i, j] <- prefoptions[mypick] setzero <- which(names(prefsset)==names(prefoptions[mypick])) prefsset[setzero[1]] <- 0 prefsset <- prefsset[prefsset > 0] prefs_grouplevel_matrix_Gb_rankorder_j[setzero[1]] <- 0 prefs_grouplevel_matrix_Gb_rankorder_j <- prefs_grouplevel_matrix_Gb_rankorder_j[prefs_grouplevel_matrix_Gb_rankorder_j > 0] } } } } # Commented-out code below is for finding mle \delta #delta_trial <- seq(0, 5, 0.1) #for (loop in 1:length(delta_trial)){ #delta <- delta_trial[loop] probvector <- mat.or.vec(nr=length(prep_Gb$y), nc=1) for (i in 1:length(probvector)){ probvector[i] <- estden_Gb[prep_Gb$xs_index[i, 1]] * estden_Gb[dim(Network_aggregate_Gb)[1] + prep_Gb$xs_index[i, 2]] * (prefs_grouplevel_matrix_Ga_reorder[prep_Gb$xs_index[i, 1], prep_Gb$xs_index[i, 2]])^delta } loglikelihood_Gb_correlatedprefs <- dmultinom(prep_Gb$y, size=NULL, probvector, log=TRUE) print("Group B---log-likelihood for correlated preferences model:") print(loglikelihood_Gb_correlatedprefs) print("delta:") print(delta) #} # get log-likelihood for specified preferences model with mass action switches # calculate log-likelihood contribution for each interaction in Gb LL_matrix_Gb <- prefs_grouplevel_matrix_Gb LL_matrix_Gb <- LL_matrix_Gb * 0.0 for (i in 1:dim(prefs_grouplevel_matrix_Gb)[1]){ for (j in 1:dim(prefs_grouplevel_matrix_Gb)[2]){ if (prefs_grouplevel_matrix_Gb[i, j] > 0){ prefs_grouplevel_matrix_Gb_temp <- prefs_grouplevel_matrix_Gb prefs_grouplevel_matrix_Gb_temp[prefs_grouplevel_matrix_Gb_temp > 0] <- 1 prefs_grouplevel_matrix_Gb_temp[i, j] <- prefs_grouplevel_matrix_Gb[i, j] probvector <- mat.or.vec(nr=length(prep_Gb$y), nc=1) for (k in 1:length(probvector)){ probvector[k] <- estden_Gb[prep_Gb$xs_index[k, 1]] * estden_Gb[dim(Network_aggregate_Gb)[1] + prep_Gb$xs_index[k, 2]] * prefs_grouplevel_matrix_Gb_temp[prep_Gb$xs_index[k, 1], prep_Gb$xs_index[k, 2]] } LL_matrix_Gb[i, j] <- dmultinom(prep_Gb$y, size=NULL, probvector, log=TRUE) } } } # hardcode four preferences prefs_grouplevel_matrix_Ga_temp <- prefs_grouplevel_matrix_Ga prefs_grouplevel_matrix_Ga_temp[1, 1] <- prefs_grouplevel_matrix_Gb[1, 1] prefs_grouplevel_matrix_Ga_temp[1, 3] <- prefs_grouplevel_matrix_Gb[1, 3] prefs_grouplevel_matrix_Ga_temp[2, 3] <- prefs_grouplevel_matrix_Gb[2, 3] prefs_grouplevel_matrix_Ga_temp[2, 2] <- prefs_grouplevel_matrix_Gb[2, 2] # probvector <- mat.or.vec(nr=length(prep_Gb$y), nc=1) for (i in 1:length(probvector)){ probvector[i] <- estden_Gb[prep_Gb$xs_index[i, 1]] * estden_Gb[dim(Network_aggregate_Gb)[1] + prep_Gb$xs_index[i, 2]] * prefs_grouplevel_matrix_Ga_temp[prep_Gb$xs_index[i, 1], prep_Gb$xs_index[i, 2]] } for (i in 1:length(probvector)){ if (prefs_grouplevel_matrix_Ga_temp[prep_Gb$xs_index[i, 1], prep_Gb$xs_index[i, 2]] == 0){ probvector[i] <- estden_Gb[prep_Gb$xs_index[i, 1]] * estden_Gb[dim(Network_aggregate_Gb)[1] + prep_Gb$xs_index[i, 2]] * Switches_Gb[prep_Gb$xs_index[i, 1], prep_Gb$xs_index[i, 2]] # mass action switches } } loglikelihood_Gb_specifiedprefs <- dmultinom(prep_Gb$y, size=NULL, probvector, log=TRUE) print("Group B---log-likelihood for specified preferences model:") print(loglikelihood_Gb_specifiedprefs) # print log-likelihood for complete characterisation model (calculated above) print("Group B---log-likelihood for complete characterisation model:") print(loglikelihood_Gb_complete) print("Model performance at the group level, predicting network structure in Group B using data from Group A") print("Alternative preferences model:") print((loglikelihood_Gb_massaction - loglikelihood_Gb_alternativeprefs) / (loglikelihood_Gb_massaction - loglikelihood_Gb_complete)) print("Correlated preferences model:") print((loglikelihood_Gb_massaction - loglikelihood_Gb_correlatedprefs) / (loglikelihood_Gb_massaction - loglikelihood_Gb_complete)) print("Specified preferences model:") print((loglikelihood_Gb_massaction - loglikelihood_Gb_specifiedprefs) / (loglikelihood_Gb_massaction - loglikelihood_Gb_complete))