require(igraph) # builds bipartites and calculates modularity require(MASS) # fits distributions require(fitdistrplus) # fits distributions #source("./requiredFuns_V1.R") source("./TriSim_RequiredFuns.R") #---------------------------------------------------------------# #--------------------- EXAMPLE 1: APPLIED EXAMPLE --------------# #---------------------------------------------------------------# #--- # General constraints nIndivids <- 50 nPatches <- 30 nTimesteps <- 40 hrSize <- ceiling(rgamma(n = nIndivids, shape = 5, scale = 2)) hrSize <- ifelse(hrSize > nPatches, nPatches, hrSize) # Scenario 1 constraints (overdispersed space use) patchQual <- rgamma(nPatches, shape = .5, scale = 50) # Scenario 2 constraints (underdispersed space use) hrSize2 <- ceiling(rgamma(n = nIndivids, shape = 2, scale = 2)) hrSize2 <- ifelse(hrSize2 > nPatches, nPatches, hrSize2) patchQual2 <- rgamma(nPatches, shape = 5, scale = 10) # specify number of sims to run nreps <- 1000 # build storage objects individModOut <- individModOut2 <- GSOut <- GSOut2 <- rep(NA, nreps) individModOutNull <- GSOutNull <- rep(NA, nreps) outFull <- outFull2 <- outFullNull <- vector("list", nreps) patchOccEstdShapeOut <- patchOccEstdScaleOut <- patchOccEstdShapeOut2 <- patchOccEstdScaleOut2 <- rep(NA, nreps) # loop through sims for(r in 1:nreps){ k <- try(tripSim_DataExample(nIndivids, nPatches, nTimesteps, hrSize, patchQual)) k.out <- class(k) while(k.out == "try-error"){ k <- try(tripSim_DataExample(nIndivids, nPatches, nTimesteps, hrSize, patchQual)) k.out <- class(k) } outFull[[r]] <- k individModOut[r] <- outFull[[r]]$individMod GSOut[r] <- outFull[[r]]$meanGS patchOccEstdShapeOut[r] <- outFull[[r]]$patchQualEstdShape patchOccEstdScaleOut[r] <- 1/outFull[[r]]$patchQualEstdScale k2 <- try(tripSim_DataExample(nIndivids, nPatches, nTimesteps, hrSize2, patchQual2)) k2.out <- class(k2) while(k2.out == "try-error"){ k2 <- try(tripSim_DataExample(nIndivids, nPatches, nTimesteps, hrSize2, patchQual2)) k2.out <- class(k2) } outFull2[[r]] <- k2 individModOut2[r] <- outFull2[[r]]$individMod GSOut2[r] <- outFull2[[r]]$meanGS patchOccEstdShapeOut2[r] <- outFull2[[r]]$patchQualEstdShape patchOccEstdScaleOut2[r] <- 1/outFull2[[r]]$patchQualEstdScale outFullNull[[r]] <- tripSim_Null(nIndivids, nPatches, nTimesteps) individModOutNull[r] <- outFullNull[[r]]$individMod GSOutNull[r] <- outFullNull[[r]]$meanGS } #-- Figure 4A: modularity comparison # get full range for x-axis of Figure 4. modsFull <- c(individModOutNull, individModOut, individModOut2) # store appropriate breaks k <- hist(modsFull, breaks = 40)$breaks layout(matrix(c(1, 1, 1, 1, 2, 3), byrow = F, nrow = 2)) par(mar = c(4, 4, 1, 1)) hist(individModOutNull, breaks = k, main = "", xlab = "Modularity", ylab = "Frequency", las = 1, cex.axis = .8, border = rgb(0, 0, 0, 1), ylim = c(0, 500)) hist(individModOut, col = rgb(0, 0, 0, alpha = .2), breaks = k,add = T) hist(individModOut2, col = rgb(0, 0, 0, alpha = .6), breaks = k,add = T) leg.text <- c("Null model -- no constraints", "Constrained toward varied patch quality", "Constrained toward consistent patch quality") legend("topright", leg.text, fill = c("white", rgb(0, 0, 0, alpha = .2), rgb(0, 0, 0, alpha = .6)), bty = "n", cex = .8) #-- median and range of refit patch qualities shapeEsts <- quantile(patchOccEstdShapeOut, c(0.025, .5, 0.975)) scaleEsts <- 1/quantile(patchOccEstdScaleOut, c(0.025, .5, 0.975)) shapeEsts2 <- quantile(patchOccEstdShapeOut2, c(0.025, .5, 0.975)) scaleEsts2 <- 1/quantile(patchOccEstdScaleOut2, c(0.025, .5, 0.975)) #-- Figure 4B/4C: Gamma fit/refit figure # unconstrained sim x <- seq(0, 100, length.out = 1001) y.lb <- dgamma(x, shape = shapeEsts[1], scale = scaleEsts[1]) y.ub <- dgamma(x, shape = shapeEsts[3], scale = scaleEsts[3]) y.lb2 <- dgamma(x, shape = shapeEsts2[1], scale = scaleEsts2[1]) y.ub2 <- dgamma(x, shape = shapeEsts2[3], scale = scaleEsts2[3]) yOrig <- dgamma(x, shape = .5, scale = 50) yOrig2 <- dgamma(x, shape = 5, scale = 10) #par(mfrow = c(1, 2)) plot(yOrig ~ x, type = "l", lwd = 1.5, ylim = c(0, 0.03), #log = "y", las = 1, ylab = "Relative frequency", xlab = "Patch quality", xlim = c(0, 100), cex.axis = .8) lines(y.lb ~ x, col = "grey60", lwd = 1, lty = 3) lines(y.ub ~ x, col = "grey60", lwd = 1, lty = 3) leg.text <- c("Originally specified", "Emergent") legend("topright", leg.text, col = c("black", "grey60"), lwd = c(1.5, 1), lty = c(1, 3), bty = "n", cex = .8, title = "Parameters") plot(yOrig2 ~ x, type = "l", lwd = 1.5, #log = "y", las = 1, ylab = "Relative frequency", xlab = "Patch quality", xlim = c(0, 100), ylim = c(0, .03), cex.axis = .8) lines(y.lb2 ~ x, col = "grey60", lwd = 1, lty = 3) lines(y.ub2 ~ x, col = "grey60", lwd = 1, lty = 3) leg.text <- c("Originally specified", "Emergent") legend("topright", leg.text, col = c("black", "grey60"), lwd = c(1.5, 1), lty = c(1, 3), bty = "n", cex = .8, title = "Parameters") #------------------------------------------------------------# #--------------------- EXAMPLE 2: CORRELATIONS --------------# #------------------------------------------------------------# nIndivids <- 50 nPatches <- 30 nTimesteps <- 20 gsShape <- exp(seq(log(.5), log(50), length.out = 10)) gsScale <- c(2:10) patchQuality <- "Yes" gridIn <- expand.grid(nIndivids, nPatches, nTimesteps, gsShape, gsScale, patchQuality ) names(gridIn) <- c("nIndivids", "nPatches", "nTimesteps", "gsShape", "gsScale", "patchQuality") metaRepTest <- repSimMeta(gridIn = gridIn, nreps = 10) metaRepTestMat <- as.data.frame( do.call("rbind", metaRepTest$metaRepOutMats) ) names(metaRepTestMat) <- c("nIndivids", "nPatches", "nTimesteps", "gsShape", "gsScale", "patchQuality", "individMod", "daysOccdFit_meanlog", "daysOccdFit_sdlog", "gsFit_meanlog", "gsFit_sdlog", "hrSizeFit_meanlog", "hrSizeFit_sdlog", "meanGS", "meanDaysOccd", "meanHrSize") head(metaRepTestMat) # specify colors cols.in <- rainbow(10, alpha = .7) # build Figure 5. par(mfrow = c(1, 3), mar = c(5, 5, 2, 2)) plot(as.numeric(as.character(metaRepTestMat$individMod)) ~ as.numeric(as.character(metaRepTestMat$meanGS)), type = "p", log = "x", las = 1, xlab = "Mean group size", ylab = "Individual association \nmodularity", col = cols.in[factor(metaRepTestMat$gsShape)]) plot(as.numeric(as.character(metaRepTestMat$individMod)) ~ as.numeric(as.character(metaRepTestMat$meanHrSize)), type = "p", xlab = "Mean home range size", ylab = "Individual association \nmodularity", las = 1, col = cols.in[factor(metaRepTestMat$gsShape)]) plot(as.numeric(as.character(metaRepTestMat$meanHrSize)) ~ as.numeric(as.character(metaRepTestMat$meanGS)), type = "p", log = "x", las = 1, xlab = "Mean group size", ylab = "Mean home range size", col = cols.in[factor(metaRepTestMat$gsShape)])