# META-ANALYSES OF SCALE-DEPENDENCE IN THE NATIVE-EXOTIC RICHNESS RELATIONSHIP # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # LOAD PACKAGES, INITIALIZE --------------------------------------------------- # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # library(metafor) library(plyr) library(polycor) # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # LOAD DATA, CLEAN ------------------------------------------------------------ # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # meta-analysis subset path <- "/Users/nicolekinlock/Documents/ScaleDepNER_MA/Docs/Peng_et_al_Resubmission/" # path to data file dat <- read.csv(paste(path, "DataS2.csv", sep = "")) dat$abs_Latitude <- abs(dat$Latitude) # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # META-ANALYSES, REGRESSIONS -------------------------------------------------- # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # random effects: case nested within study, used in all models random <- list(~ 1 | Study / Case) # meta-regression of z as a function of ln grain size metar <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain, random = random, data = dat, method = "ML") # meta-regressions with other covariates # # covariate diagnostics # are the continuous covariates correlated? pairs(~ abs_Latitude + Longitude + ln_Grain, data = dat) # grain and extent are highly correlated # hetcor in the polycor allows for the calculation of r between numeric, polyserial correlations between numeric and ordinal, polychoric correlations between ordinal hetcor(data = dat[, c("ln_Grain", "abs_Latitude", "Longitude", "Extent", "Habitat_type")], std.err = FALSE) # no other covariates highly correlated # single covariate meta-regressions of z as a function of latitude, longitude (linear, with quadratic term, circular), extent (as a categorical and as an ordered covariate), and habitat type lat <- rma.mv(yi = z, V = var_z, mods = ~ abs_Latitude, random = random, data = dat, method = "ML") lon2 <- rma.mv(yi = z, V = var_z, mods = ~ Longitude + I(Longitude^2), random = random, data = dat, method = "ML") ext <- rma.mv(yi = z, V = var_z, mods = ~ -1 + Extent, random = random, data = dat, method = "ML") ext.ord <- rma.mv(yi = z, V = var_z, mods = ~ ordered(Extent), random = random, data = dat, method = "ML") subset.hab <- dat[which(dat$Habitat_type != "many"), ] hab.full <- rma.mv(yi = z, V = var_z, mods = ~ -1 + Habitat_type, random = random, data = dat, method = "ML") hab <- rma.mv(yi = z, V = var_z, mods = ~ -1 + Habitat_type, random = random, data = subset.hab, method = "ML") # fit separate regressions to the forest and grassland studies (have different grains, may have different slopes) forest.only <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain, random = random, data = dat[which(dat$Habitat_type == "forest"), ], method = "ML") grassland.only <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain, random = random, data = dat[which(dat$Habitat_type == "grassland"), ], method = "ML") # including intercept to get an accurate Q-test for heterogeneity ext.Q <- rma.mv(yi = z, V = var_z, mods = ~ Extent, random = random, data = dat, method = "ML") hab.full.Q <- rma.mv(yi = z, V = var_z, mods = ~ Habitat_type, random = random, data = dat, method = "ML") hab.Q <- rma.mv(yi = z, V = var_z, mods = ~ Habitat_type, random = random, data = subset.hab, method = "ML") # multiple meta-regressions including interaction terms # these interaction terms are hypothesized to be important in explaining variation in z metar.ext.int <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain * ordered(Extent), random = random, data = dat, method = "ML") metar.hab.int <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain * Habitat_type, random = random, data = subset.hab, method = "ML") # test each factor as a whole # Wald-type chi-square test of the null hypothesis that coefficient A = coefficient B = 0 # test of extent * grain interaction anova(metar.ext.int, btt = 9:14) # test of habitat * grain interaction anova(metar.hab.int, btt = 8:12) # multiple meta-regressions of z as a function of grain, including latitude, longitude (linear and quadratic term), extent, and habitat type as covariates # all combinations of these covariates with ln grain always in model metar.lat <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain + abs_Latitude, random = random, data = dat, method = "ML") metar.lon2 <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain + Longitude + I(Longitude^2), random = random, data = dat, method = "ML") metar.ext <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain + Extent, random = random, data = dat, method = "ML") metar.hab <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain + Habitat_type, random = random, data = dat, method = "ML") metar.latlon2 <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain + abs_Latitude + Longitude + I(Longitude^2), random = random, data = dat, method = "ML") metar.latext <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain + abs_Latitude + Extent, random = random, data = dat, method = "ML") metar.lathab <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain + abs_Latitude + Habitat_type, random = random, data = dat, method = "ML") metar.lon2ext <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain + Longitude + I(Longitude^2) + Extent, random = random, data = dat, method = "ML") metar.lon2hab <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain + Longitude + I(Longitude^2) + Habitat_type, random = random, data = dat, method = "ML") metar.exthab <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain + Extent + Habitat_type, random = random, data = dat, method = "ML") metar.exthablat <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain + Extent + Habitat_type + abs_Latitude, random = random, data = dat, method = "ML") metar.exthablon2 <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain + Extent + Habitat_type + Longitude + I(Longitude^2), random = random, data = dat, method = "ML") metar.hablatlon2 <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain + Habitat_type + abs_Latitude + Longitude + I(Longitude^2), random = random, data = dat, method = "ML") metar.extlatlon2 <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain + Extent + abs_Latitude + Longitude + I(Longitude^2), random = random, data = dat, method = "ML") metar.full <- rma.mv(yi = z, V = var_z, mods = ~ ln_Grain + abs(Latitude) + Longitude + I(Longitude^2) + Extent + Habitat_type, random = random, data = dat, method = "ML") # random effects only, no fixed effects (no covariates), for pseudo R2 calculation randomonly <- rma.mv(yi = z, V = var_z, random = random, data = dat, method = "ML") models <- list(metar, metar.lat, metar.lon2, metar.ext, metar.hab, metar.latlon2, metar.latext, metar.lathab, metar.lon2ext, metar.lon2hab, metar.exthab, metar.exthablat, metar.exthablon2, metar.hablatlon2, metar.extlatlon2, metar.full) # model fit diagnostics: pseudo R2, marginal/conditional R2, AICc, and I2 pseudo.r2 <- c() marg.r2 <- c() cond.r2 <- c() aicc <- c() I2 <- c() for (i in 1:length(models)) { # pseudo R2 of proportional reduction of variance explained when fitting reduced model relative to model with random effects only pseudo.r2[i] <- (sum(randomonly$sigma2) - sum(models[[i]]$sigma2)) / sum(randomonly$sigma2) # marginal and conditional R2 for linear mixed-models from Nakagawa and Scheilzeth 2013 # variance from fixed effects fe.total <- c() for (j in 2:ncol(model.matrix(models[[i]]))) { fe <- models[[i]]$b[j] * model.matrix(models[[i]])[, j] fe.total <- c(fe.total, fe) } v.fix <- var(fe.total) # variance from random effects v.rand <- sum(models[[i]]$sigma2) # variance from residuals v.resid <- var(residuals(models[[i]])) # marginal R2 - total variance explained from fixed effects marg.r2[i] <- v.fix / (v.fix + v.rand + v.resid) # conditional R2 - total variance explained from fixed and random effects cond.r2[i] <- (v.fix + v.rand) / (v.fix + v.rand + v.resid) # from Wolfgang Viechtbauer (http://www.metafor-project.org/doku.php/tips:i2_multilevel_multivariate) # I2 statistic - amount of heterogeneity relative to the total amount of variance in observed effects (how much heterogeneity contributes to total variance) W <- diag(1/dat$var_z) X <- model.matrix(models[[i]]) P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W I2[i] <- 100 * sum(models[[i]]$sigma2) / (sum(models[[i]]$sigma2) + (models[[i]]$k - models[[i]]$p) / sum(diag(P))) aicc[i] <- fitstats(models[[i]])[5] }