Non-native (exotic) single moderator models
Lia Ossanna
2022-10-17
Analysis for Ossanna & Gornish (2022), “Efficacy of labile carbon addition to reduce fast-growing, invasive non-native plants: A review and meta-analysis”, Journal of Applied Ecology.
Introduction
- This is the analysis of the summary effect and effect of a single
moderator on non-native (exotic) plant abundance.
- See “Data cleaning and effect size” document for code for data
cleaning and calculation of effect size
(
code_data-cleaning.html).
- See “Data cleaning and effect size” document for code for data
cleaning and calculation of effect size
(
- Note that the word “exotic” is used to refer to the non-native dataset and models in the code.
- For each model there are the model results, proportion of heterogeneity explained by the model (\(\frac{Q_M}{Q_E+Q_M}\)), and a forest plot.
- When there are missing data for a moderator, a new model is run with
the subset of data that the single-moderator model uses; this is
reflected in the forest plot.
- Missing data will return
TRUEfor the first line of code,any(is.na(dat$variable)).
- Missing data will return
- Analysis includes a model run with all the data, and a model run
with outliers removed (see “Non-native (exotic) outliers for models”
document for details;
code_nonnative-outliers.html).- When only one model is presented, no outliers were removed.
Setup
library(tidyverse)
library(metafor)
# Load data ---------------------------------------------------------------
ex <- read.csv("data/cleaned/exotic-cleaned.csv") # data cleaning and effect size calculation in separate script
# Functions ---------------------------------------------------------------
single.mod.cat <- function(dat, var, dat.mv, dat.mv.var, h) {
print(dat.mv.var)
print(noquote(paste("QM/QT =", dat.mv.var[["QM"]] / sum(dat.mv.var[["QE"]], dat.mv.var[["QM"]]))))
n.forest <- count(dat, !!sym(var))
colnames(n.forest) <- c("moderator", "n")
n.forest$end <- rep(")", dim(n.forest)[1])
n.forest$n <- paste(n.forest$n, n.forest$end, sep = "")
n.forest$start <- rep("(n = ", dim(n.forest)[1])
n.forest$n <- paste(n.forest$start, n.forest$n)
n.forest$var_col <- paste(n.forest$moderator, n.forest$n)
n.forest <- n.forest %>%
filter(!is.na(moderator))
overall <- paste("Overall (n = ", dat.mv$k, ")", sep = "")
forest(x = c(dat.mv[[1]], dat.mv.var[[1]]),
ci.lb = c(dat.mv[[6]], dat.mv.var[[6]]),
ci.ub = c(dat.mv[[7]], dat.mv.var[[7]]),
slab = c(overall, n.forest$var_col),
header = h,
cex = 0.7)
outliers <- function(dat, out) {
dat %>%
filter(!obs_ID %in% out)
}
# Outliers ----------------------------------------------------------------
out <- c(457, 601)
out.reg <- c(601, 457)
out.dfc <- c(457)
out.crat <- c(54, 120, 132)
out.cratc <- c(437)
out.capm <- c(257, 247, 601)
out.pgfs <- c(257, 258)No moderators (summary effect)
All data
ex.mv <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex)print(ex.mv)##
## Multivariate Meta-Analysis Model (k = 488; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.4018 0.6339 50 no exp_ID
## sigma^2.2 0.1701 0.4124 488 no exp_ID/obs_ID
##
## Test for Heterogeneity:
## Q(df = 487) = 1042.8256, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.5831 0.1068 -5.4613 <.0001 -0.7923 -0.3738 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1predict(ex.mv)##
## pred se ci.lb ci.ub pi.lb pi.ub
## -0.5831 0.1068 -0.7923 -0.3738 -2.0800 0.9138# QM/QT
ex.mv[["QM"]] / sum(ex.mv[["QE"]], ex.mv[["QM"]])## [1] 0.0278053# I^2
W <- diag(1 / ex$vi)
X <- model.matrix(ex.mv)
P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W
100 * sum(ex.mv$sigma2) / (sum(ex.mv$sigma2) + (ex.mv$k-ex.mv$p)/sum(diag(P)))## [1] 59.59988# I^2 across nested experiment and within nested experiment
100 * ex.mv$sigma2 / (sum(ex.mv$sigma2) + (ex.mv$k-ex.mv$p)/sum(diag(P)))## [1] 41.87386 17.72602# Probability based on PI (IntHout 2016)
D <- 0 # threshold D (null effect)
M <- ex.mv[["b"]] # summary treatment effect
SD.PI <- predict(ex.mv)$se * sqrt(ex.mv[["k"]]) # SD of prediction interval, calculated from SE
t <- as.numeric((D - M) / (SD.PI))
pt(t, ex.mv[["k"]] - 1) # left-tail cumulative t-distribution with k-1 df## [1] 0.5975788Outliers removed
out <- c(457, 601)
ex.o <- outliers(ex, out)
ex.mv.o <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex.o)print(ex.mv.o)##
## Multivariate Meta-Analysis Model (k = 486; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.3457 0.5880 50 no exp_ID
## sigma^2.2 0.1554 0.3942 486 no exp_ID/obs_ID
##
## Test for Heterogeneity:
## Q(df = 485) = 1012.8128, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.5768 0.1009 -5.7168 <.0001 -0.7746 -0.3791 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1predict(ex.mv.o)##
## pred se ci.lb ci.ub pi.lb pi.ub
## -0.5768 0.1009 -0.7746 -0.3791 -1.9782 0.8246# QM/QT
ex.mv.o[["QM"]] / sum(ex.mv.o[["QE"]], ex.mv.o[["QM"]])## [1] 0.03125952# I^2 overall
W <- diag(1 / ex.o$vi)
X <- model.matrix(ex.mv.o)
P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W
100 * sum(ex.mv.o$sigma2) / (sum(ex.mv.o$sigma2) + (ex.mv.o$k-ex.mv.o$p)/sum(diag(P)))## [1] 56.4348# I^2 across nested experiment and within nested experiment
100 * ex.mv.o$sigma2 / (sum(ex.mv.o$sigma2) + (ex.mv.o$k-ex.mv.o$p)/sum(diag(P)))## [1] 38.93651 17.49828# Probability based on PI (IntHout 2016)
D <- 0 # threshold D (null effect)
M <- ex.mv.o[["b"]] # summary treatment effect
SD.PI <- predict(ex.mv.o)$se * sqrt(ex.mv.o[["k"]]) # SD of prediction interval, calculated from SE
t <- as.numeric((D - M) / (SD.PI))
pt(t, ex.mv.o[["k"]] - 1) # left-tail cumulative t-distribution with k-1 df## [1] 0.6022503Comparison
- Removing outliers does not change significance.
- Removing outliers increases QM/QT, but decreases I2.
Biome
any(is.na(ex$biome)) # if TRUE, run no-moderators model on subset with no missing data for that moderator (i.e. ex.mv.bio.na)## [1] TRUEex.bio.na <- ex %>%
filter(!is.na(biome))
ex.mv.bio <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex.bio.na,
mods = ~ factor(biome) - 1)
ex.mv.bio.na <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex.bio.na)h <- c("Biome", "Estimate [95% CI]")
single.mod.cat(ex, "biome", ex.mv.bio.na, ex.mv.bio, h)##
## Multivariate Meta-Analysis Model (k = 487; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.2301 0.4797 49 no exp_ID
## sigma^2.2 0.1755 0.4189 487 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 470) = 871.1333, p-val < .0001
##
## Test of Moderators (coefficients 1:17):
## QM(df = 17) = 80.1483, p-val < .0001
##
## Model Results:
##
## estimate se zval pval
## factor(biome)coastal grassland -0.1447 0.2488 -0.5817 0.5608
## factor(biome)coastal sage scrub -0.1019 0.7200 -0.1415 0.8875
## factor(biome)coastal sandplain grassland 0.2509 0.5270 0.4762 0.6339
## factor(biome)cold desert -0.5461 0.5585 -0.9777 0.3282
## factor(biome)creosote scrub -1.4725 0.4098 -3.5928 0.0003
## factor(biome)fescue prairie -0.2472 0.4013 -0.6160 0.5379
## factor(biome)fynbos 0.1676 0.7309 0.2293 0.8186
## factor(biome)grassy woodland -1.4107 0.2553 -5.5258 <.0001
## factor(biome)mixed-grass prairie -0.6727 0.2779 -2.4207 0.0155
## factor(biome)Puget lowland prairie 0.0675 0.4742 0.1424 0.8868
## factor(biome)sagebrush steppe -0.1758 0.2348 -0.7488 0.4540
## factor(biome)sedge meadow -1.6249 0.9214 -1.7635 0.0778
## factor(biome)shadscale-bunchgrass -0.8943 0.7910 -1.1307 0.2582
## factor(biome)shortgrass steppe -0.2212 0.2796 -0.7910 0.4289
## factor(biome)tallgrass prairie -0.6662 0.3613 -1.8439 0.0652
## factor(biome)upland prairie -1.2566 0.2838 -4.4283 <.0001
## factor(biome)wetland prairie -0.2731 0.5288 -0.5163 0.6056
## ci.lb ci.ub
## factor(biome)coastal grassland -0.6324 0.3429
## factor(biome)coastal sage scrub -1.5131 1.3093
## factor(biome)coastal sandplain grassland -0.7820 1.2838
## factor(biome)cold desert -1.6407 0.5486
## factor(biome)creosote scrub -2.2758 -0.6692 ***
## factor(biome)fescue prairie -1.0338 0.5394
## factor(biome)fynbos -1.2649 1.6002
## factor(biome)grassy woodland -1.9110 -0.9103 ***
## factor(biome)mixed-grass prairie -1.2173 -0.1280 *
## factor(biome)Puget lowland prairie -0.8618 0.9968
## factor(biome)sagebrush steppe -0.6360 0.2844
## factor(biome)sedge meadow -3.4308 0.1810 .
## factor(biome)shadscale-bunchgrass -2.4445 0.6559
## factor(biome)shortgrass steppe -0.7693 0.3269
## factor(biome)tallgrass prairie -1.3743 0.0419 .
## factor(biome)upland prairie -1.8128 -0.7004 ***
## factor(biome)wetland prairie -1.3096 0.7634
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0842529853004472
Region
All data
any(is.na(ex$region)) # if FALSE, use ex.mv for no-moderator model## [1] FALSE # but still requires new model for outliers removed and levels n > 3ex.mv.reg <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex,
mods = ~ factor(region) - 1)h <- c("Region (all)", "Estimate [95% CI]")
single.mod.cat(ex, "region", ex.mv, ex.mv.reg, h)##
## Multivariate Meta-Analysis Model (k = 488; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.2720 0.5216 50 no exp_ID
## sigma^2.2 0.1770 0.4207 488 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 477) = 901.5078, p-val < .0001
##
## Test of Moderators (coefficients 1:11):
## QM(df = 11) = 61.9513, p-val < .0001
##
## Model Results:
##
## estimate se zval pval
## factor(region)central Canada -0.0708 0.3372 -0.2101 0.8336
## factor(region)coastal CA, USA -0.1428 0.2507 -0.5696 0.5689
## factor(region)Great Basin, USA -0.3926 0.2663 -1.4744 0.1404
## factor(region)Great Plains, USA -0.4209 0.2790 -1.5083 0.1315
## factor(region)Midwest, USA -0.6098 0.3308 -1.8435 0.0653
## factor(region)northeast USA 0.2510 0.5655 0.4439 0.6571
## factor(region)Pacific Northwest, USA -0.7366 0.2176 -3.3850 0.0007
## factor(region)Rocky Mountains, USA -0.4615 0.3111 -1.4835 0.1379
## factor(region)South Africa 0.1677 0.7595 0.2207 0.8253
## factor(region)southeast Australia -1.4402 0.2705 -5.3249 <.0001
## factor(region)southwest USA -1.4731 0.4348 -3.3876 0.0007
## ci.lb ci.ub
## factor(region)central Canada -0.7317 0.5900
## factor(region)coastal CA, USA -0.6341 0.3485
## factor(region)Great Basin, USA -0.9146 0.1293
## factor(region)Great Plains, USA -0.9678 0.1260
## factor(region)Midwest, USA -1.2582 0.0385 .
## factor(region)northeast USA -0.8573 1.3592
## factor(region)Pacific Northwest, USA -1.1631 -0.3101 ***
## factor(region)Rocky Mountains, USA -1.0712 0.1482
## factor(region)South Africa -1.3210 1.6563
## factor(region)southeast Australia -1.9703 -0.9101 ***
## factor(region)southwest USA -2.3254 -0.6208 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0643009439262971
Outliers removed
out.reg <- c(601, 457)
ex.reg.o <- outliers(ex, out.reg.1)
ex.mv.reg.o <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex.reg.o,
mods = ~ factor(region) - 1)
ex.mv.reg.na.o <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex.reg.o)h <- c("Region (outliers removed)", "Estimate [95% CI]")
single.mod.cat(ex.reg.o, "region", ex.mv.reg.na.o, ex.mv.reg.o, h)##
## Multivariate Meta-Analysis Model (k = 486; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.2076 0.4557 50 no exp_ID
## sigma^2.2 0.1626 0.4032 486 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 475) = 871.3672, p-val < .0001
##
## Test of Moderators (coefficients 1:11):
## QM(df = 11) = 71.9946, p-val < .0001
##
## Model Results:
##
## estimate se zval pval
## factor(region)central Canada -0.0613 0.3087 -0.1985 0.8427
## factor(region)coastal CA, USA -0.1389 0.2254 -0.6163 0.5377
## factor(region)Great Basin, USA -0.3869 0.2383 -1.6237 0.1044
## factor(region)Great Plains, USA -0.4088 0.2548 -1.6043 0.1087
## factor(region)Midwest, USA -0.5812 0.3029 -1.9192 0.0550
## factor(region)northeast USA 0.2507 0.5045 0.4968 0.6193
## factor(region)Pacific Northwest, USA -0.7396 0.1969 -3.7561 0.0002
## factor(region)Rocky Mountains, USA -0.4694 0.3000 -1.5650 0.1176
## factor(region)South Africa 0.1672 0.7108 0.2352 0.8140
## factor(region)southeast Australia -1.3910 0.2460 -5.6537 <.0001
## factor(region)southwest USA -1.4678 0.3943 -3.7227 0.0002
## ci.lb ci.ub
## factor(region)central Canada -0.6662 0.5437
## factor(region)coastal CA, USA -0.5807 0.3029
## factor(region)Great Basin, USA -0.8539 0.0801
## factor(region)Great Plains, USA -0.9082 0.0906
## factor(region)Midwest, USA -1.1749 0.0124 .
## factor(region)northeast USA -0.7382 1.2395
## factor(region)Pacific Northwest, USA -1.1255 -0.3536 ***
## factor(region)Rocky Mountains, USA -1.0573 0.1185
## factor(region)South Africa -1.2260 1.5604
## factor(region)southeast Australia -1.8733 -0.9088 ***
## factor(region)southwest USA -2.2406 -0.6950 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0763170564949973
Comparison
- Removing outliers does not change significance.
- Removing outliers increases QM/QT.
Soil suborder
any(is.na(ex$soil_suborder))## [1] TRUEex.soil.na <- ex %>%
filter(!is.na(soil_suborder))
ex.mv.soil <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex.soil.na,
mods = ~ factor(soil_suborder) - 1)
ex.mv.soil.na <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex.soil.na)h <- c("Soil suborder", "Estimate [95% CI]")
single.mod.cat(ex, "soil_suborder", ex.mv.soil.na, ex.mv.soil, h)##
## Multivariate Meta-Analysis Model (k = 410; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.2156 0.4644 40 no exp_ID
## sigma^2.2 0.1646 0.4057 410 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 391) = 718.7237, p-val < .0001
##
## Test of Moderators (coefficients 1:19):
## QM(df = 19) = 80.3878, p-val < .0001
##
## Model Results:
##
## estimate se zval pval
## factor(soil_suborder)Aqualfs -0.2728 0.5138 -0.5310 0.5955
## factor(soil_suborder)Aquolls -1.1822 0.4929 -2.3984 0.0165
## factor(soil_suborder)Argids -0.4990 0.3010 -1.6579 0.0973
## factor(soil_suborder)Blacktown -0.9221 0.3970 -2.3226 0.0202
## factor(soil_suborder)Calcids -1.3218 0.5652 -2.3385 0.0194
## factor(soil_suborder)Cambids 0.0960 0.3388 0.2834 0.7769
## factor(soil_suborder)Chernozems -0.0627 0.3124 -0.2006 0.8410
## factor(soil_suborder)Humults -0.8103 0.3902 -2.0768 0.0378
## factor(soil_suborder)Orthents -0.4614 0.5219 -0.8841 0.3766
## factor(soil_suborder)Psamments, Udepts 0.2507 0.5125 0.4892 0.6247
## factor(soil_suborder)Udolls -0.0999 0.4701 -0.2126 0.8316
## factor(soil_suborder)Udults 0.7790 0.8596 0.9062 0.3648
## factor(soil_suborder)Ustalfs -0.0695 0.3353 -0.2073 0.8357
## factor(soil_suborder)Ustepts -0.8221 0.5408 -1.5200 0.1285
## factor(soil_suborder)Ustolls -0.4746 0.2572 -1.8451 0.0650
## factor(soil_suborder)Vertisols -7.0545 1.2253 -5.7573 <.0001
## factor(soil_suborder)Xerands 0.1273 0.5805 0.2193 0.8264
## factor(soil_suborder)Xerepts -0.0057 0.6278 -0.0091 0.9928
## factor(soil_suborder)Xerolls -1.2368 0.3145 -3.9328 <.0001
## ci.lb ci.ub
## factor(soil_suborder)Aqualfs -1.2798 0.7342
## factor(soil_suborder)Aquolls -2.1484 -0.2161 *
## factor(soil_suborder)Argids -1.0890 0.0909 .
## factor(soil_suborder)Blacktown -1.7003 -0.1440 *
## factor(soil_suborder)Calcids -2.4296 -0.2139 *
## factor(soil_suborder)Cambids -0.5680 0.7601
## factor(soil_suborder)Chernozems -0.6750 0.5497
## factor(soil_suborder)Humults -1.5749 -0.0456 *
## factor(soil_suborder)Orthents -1.4842 0.5614
## factor(soil_suborder)Psamments, Udepts -0.7538 1.2552
## factor(soil_suborder)Udolls -1.0212 0.8214
## factor(soil_suborder)Udults -0.9058 2.4637
## factor(soil_suborder)Ustalfs -0.7267 0.5876
## factor(soil_suborder)Ustepts -1.8821 0.2380
## factor(soil_suborder)Ustolls -0.9788 0.0296 .
## factor(soil_suborder)Vertisols -9.4561 -4.6529 ***
## factor(soil_suborder)Xerands -1.0105 1.2651
## factor(soil_suborder)Xerepts -1.2361 1.2248
## factor(soil_suborder)Xerolls -1.8531 -0.6204 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.100596475688904
Duration since first C app (categorical)
All data
# Create categories
count(ex, duration_first)##
## duration_first n
## 1 3.0 88
## 2 5.0 8
## 3 6.0 3
## 4 7.0 42
## 5 8.0 44
## 6 9.0 11
## 7 11.0 27
## 8 12.0 3
## 9 15.0 75
## 10 16.0 5
## 11 17.0 18
## 12 18.0 6
## 13 18.5 10
## 14 20.0 8
## 15 22.0 6
## 16 23.0 2
## 17 24.0 6
## 18 26.0 1
## 19 27.0 1
## 20 30.0 9
## 21 39.0 29
## 22 45.0 45
## 23 46.0 19
## 24 50.0 8
## 25 112.0 1
## 26 116.0 2
## 27 120.0 2
## 28 123.0 1
## 29 207.0 6
## 30 266.0 2ex$dfc <- NA
for(i in 1:nrow(ex)) {
if(ex$duration_first[i] == 3) {
ex$dfc[i] <- "3"
} else if(between(ex$duration_first[i], 5, 6)) {
ex$dfc[i] <- "5-6"
} else if(between(ex$duration_first[i], 7, 12)) {
ex$dfc[i] <- "7-12"
} else if(between(ex$duration_first[i], 13, 18)) {
ex$dfc[i] <- "13-18"
} else if(between(ex$duration_first[i], 19, 24)) {
ex$dfc[i] <- "19-24"
} else if(between(ex$duration_first[i], 25, 36)) {
ex$dfc[i] <- "25-36"
} else if(between(ex$duration_first[i], 37, 50)) {
ex$dfc[i] <- "37-50"
} else if(between(ex$duration_first[i], 100, 200)) {
ex$dfc[i] <- "100-200"
} else {
ex$dfc[i] <- ">200"
}
}
ex$dfc <- factor(ex$df, levels = c("3", "5-6", "7-12", "12-18", "19-24", "25-36", "37-50", "100-200", ">200"))count(ex, dfc)##
## dfc n
## 1 3 88
## 2 5-6 11
## 3 7-12 127
## 4 13-18 104
## 5 19-24 22
## 6 25-36 11
## 7 37-50 101
## 8 100-200 6
## 9 >200 18any(is.na(ex$dfc))## [1] FALSEex.mv.dfc <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex,
mods = ~ factor(dfc) - 1)h <- c("Duration since first C app (all)", "Estimate [95% CI]")
single.mod.cat(ex, "dfc", ex.mv, ex.mv.dfc, h)##
## Multivariate Meta-Analysis Model (k = 488; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.4942 0.7030 50 no exp_ID
## sigma^2.2 0.1452 0.3810 488 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 479) = 976.4876, p-val < .0001
##
## Test of Moderators (coefficients 1:9):
## QM(df = 9) = 56.1709, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## factor(dfc)3 -1.5026 0.2559 -5.8722 <.0001 -2.0042 -1.0011 ***
## factor(dfc)5-6 -0.2381 0.3760 -0.6334 0.5265 -0.9750 0.4988
## factor(dfc)7-12 -0.7388 0.2932 -2.5195 0.0118 -1.3135 -0.1641 *
## factor(dfc)13-18 -0.8729 0.2369 -3.6848 0.0002 -1.3372 -0.4086 ***
## factor(dfc)19-24 -0.3404 0.2852 -1.1937 0.2326 -0.8993 0.2185
## factor(dfc)25-36 -0.9957 0.4008 -2.4840 0.0130 -1.7813 -0.2101 *
## factor(dfc)37-50 -0.2822 0.2378 -1.1868 0.2353 -0.7482 0.1839
## factor(dfc)100-200 0.1022 0.5297 0.1930 0.8470 -0.9360 1.1405
## factor(dfc)>200 -0.3740 0.3395 -1.1017 0.2706 -1.0394 0.2914
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0543944281307179
Outliers removed
out.dfc <- c(457)
ex.dfc.o <- outliers(ex, out.dfc)
ex.mv.dfc.o <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex.dfc.o,
mods = ~ factor(dfc) - 1)
ex.mv.dfc.na.o <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex.dfc.o)h <- c("Duration since first C app (outliers removed)", "Estimate [95% CI]")
single.mod.cat(ex.dfc.o, "dfc", ex.mv.dfc.na.o, ex.mv.dfc.o, h)##
## Multivariate Meta-Analysis Model (k = 487; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.4202 0.6482 50 no exp_ID
## sigma^2.2 0.1441 0.3796 487 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 478) = 959.4345, p-val < .0001
##
## Test of Moderators (coefficients 1:9):
## QM(df = 9) = 60.6108, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## factor(dfc)3 -1.4735 0.2425 -6.0764 <.0001 -1.9487 -0.9982 ***
## factor(dfc)5-6 -0.2699 0.3590 -0.7518 0.4522 -0.9735 0.4337
## factor(dfc)7-12 -0.7167 0.2751 -2.6052 0.0092 -1.2560 -0.1775 **
## factor(dfc)13-18 -0.8532 0.2225 -3.8352 0.0001 -1.2892 -0.4172 ***
## factor(dfc)19-24 -0.3350 0.2698 -1.2417 0.2144 -0.8638 0.1938
## factor(dfc)25-36 -1.0084 0.3808 -2.6484 0.0081 -1.7547 -0.2621 **
## factor(dfc)37-50 -0.2850 0.2239 -1.2728 0.2031 -0.7238 0.1539
## factor(dfc)100-200 0.0127 0.5102 0.0248 0.9802 -0.9874 1.0127
## factor(dfc)>200 -0.5610 0.3281 -1.7100 0.0873 -1.2040 0.0820 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0594197570042484
Comparison
- Removing outliers does not change significance.
- Removing outliers increases QM/QT.
Duration since last C app (categorical)
# Create categories
count(ex, duration_last)##
## duration_last n
## 1 0.50 10
## 2 0.75 32
## 3 1.00 11
## 4 2.00 64
## 5 3.00 92
## 6 4.00 48
## 7 5.00 4
## 8 6.00 2
## 9 7.00 42
## 10 8.00 44
## 11 9.00 8
## 12 12.00 1
## 13 14.00 7
## 14 15.00 65
## 15 16.00 8
## 16 22.00 5
## 17 23.00 2
## 18 24.00 6
## 19 26.00 1
## 20 36.00 5
## 21 45.00 2
## 22 46.00 19
## 23 49.00 8
## 24 118.00 2ex$dlc <- NA
for(i in 1:nrow(ex)) {
if(between(ex$duration_last[i], 0, 1.5)) {
ex$dlc[i] <- "0-1.5"
} else if(ex$duration_last[i] == 2) {
ex$dlc[i] <- "2"
} else if(between(ex$duration_last[i], 3, 3.5)) {
ex$dlc[i] <- "3-3.5"
} else if(between(ex$duration_last[i], 4, 6)) {
ex$dlc[i] <- "4-6"
} else if(between(ex$duration_last[i], 7, 12)) {
ex$dlc[i] <- "7-12"
} else if(between(ex$duration_last[i], 13, 18)) {
ex$dlc[i] <- "13-18"
} else if(between(ex$duration_last[i], 19, 24)) {
ex$dlc[i] <- "19-24"
} else if(between(ex$duration_last[i], 36, 49)) {
ex$dlc[i] <- "36-49"
} else {
ex$dlc[i] <- ">100"
}
}
ex$dlc <- factor(ex$dl, levels = c("0-1.5", "2", "3-3.5", "4-6", "7-12", "13-18", "19-24",
"36-49", ">100"))count(ex, dlc)##
## dlc n
## 1 0-1.5 53
## 2 2 64
## 3 3-3.5 92
## 4 4-6 54
## 5 7-12 95
## 6 13-18 80
## 7 19-24 13
## 8 36-49 34
## 9 >100 3any(is.na(ex$dlc))## [1] FALSEex.mv.dlc <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex,
mods = ~ factor(dlc) - 1)h <- c("Duration since last C app", "Estimate [95% CI]")
single.mod.cat(ex, "dlc", ex.mv, ex.mv.dlc, h)##
## Multivariate Meta-Analysis Model (k = 488; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.3964 0.6296 50 no exp_ID
## sigma^2.2 0.1461 0.3822 488 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 479) = 954.5516, p-val < .0001
##
## Test of Moderators (coefficients 1:9):
## QM(df = 9) = 64.1130, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## factor(dlc)0-1.5 -1.0594 0.2301 -4.6038 <.0001 -1.5104 -0.6084 ***
## factor(dlc)2 -0.7199 0.2262 -3.1822 0.0015 -1.1633 -0.2765 **
## factor(dlc)3-3.5 -0.9071 0.2384 -3.8055 0.0001 -1.3743 -0.4399 ***
## factor(dlc)4-6 -0.4289 0.2847 -1.5063 0.1320 -0.9870 0.1292
## factor(dlc)7-12 -0.5226 0.2544 -2.0539 0.0400 -1.0212 -0.0239 *
## factor(dlc)13-18 -0.2168 0.2311 -0.9384 0.3480 -0.6697 0.2360
## factor(dlc)19-24 -0.3596 0.3061 -1.1750 0.2400 -0.9595 0.2402
## factor(dlc)36-49 -0.1759 0.3313 -0.5309 0.5955 -0.8253 0.4735
## factor(dlc)>100 0.0654 0.6332 0.1033 0.9177 -1.1756 1.3064
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0629382636315438
C type
any(is.na(ex$C_type))## [1] FALSEex.mv.ctyp <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex,
mods = ~ factor(C_type) - 1)h <- c("C_type", "Estimate [95% CI]")
single.mod.cat(ex, "C_type", ex.mv, ex.mv.ctyp, h)##
## Multivariate Meta-Analysis Model (k = 488; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.3840 0.6196 50 no exp_ID
## sigma^2.2 0.1726 0.4155 488 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 476) = 948.4239, p-val < .0001
##
## Test of Moderators (coefficients 1:12):
## QM(df = 12) = 42.8679, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb
## factor(C_type)dextrose -2.0034 0.8939 -2.2411 0.0250 -3.7555
## factor(C_type)lignin -0.0859 0.7962 -0.1078 0.9141 -1.6465
## factor(C_type)lignin AND sucrose 0.3410 0.7999 0.4263 0.6699 -1.2268
## factor(C_type)mulch 0.3041 0.4842 0.6281 0.5300 -0.6449
## factor(C_type)sawdust -0.5629 0.2564 -2.1951 0.0282 -1.0654
## factor(C_type)sawdust AND woodchip -0.5100 0.9786 -0.5211 0.6023 -2.4279
## factor(C_type)straw -0.3747 0.3973 -0.9430 0.3457 -1.1534
## factor(C_type)sucrose -0.7078 0.1360 -5.2056 <.0001 -0.9743
## factor(C_type)sucrose AND mulch -1.0689 0.7910 -1.3513 0.1766 -2.6194
## factor(C_type)sucrose AND sawdust -0.3614 0.2637 -1.3703 0.1706 -0.8783
## factor(C_type)woodchip -0.1018 0.5041 -0.2020 0.8399 -1.0898
## factor(C_type)woodchip mulch 0.1675 0.8286 0.2022 0.8398 -1.4565
## ci.ub
## factor(C_type)dextrose -0.2513 *
## factor(C_type)lignin 1.4747
## factor(C_type)lignin AND sucrose 1.9087
## factor(C_type)mulch 1.2532
## factor(C_type)sawdust -0.0603 *
## factor(C_type)sawdust AND woodchip 1.4080
## factor(C_type)straw 0.4041
## factor(C_type)sucrose -0.4413 ***
## factor(C_type)sucrose AND mulch 0.4815
## factor(C_type)sucrose AND sawdust 0.1555
## factor(C_type)woodchip 0.8861
## factor(C_type)woodchip mulch 1.7916
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0432444925112445
C rate (categorical)
All data
# Create categories
count(ex, C_rate)##
## C_rate n
## 1 30.00000 2
## 2 50.00000 10
## 3 63.00000 1
## 4 76.33719 10
## 5 84.00000 10
## 6 93.60000 2
## 7 100.00000 10
## 8 133.00000 10
## 9 150.00000 16
## 10 160.00000 30
## 11 173.50000 2
## 12 183.94505 8
## 13 195.00000 43
## 14 200.00000 7
## 15 210.00000 12
## 16 229.01158 4
## 17 250.00000 2
## 18 260.40000 1
## 19 275.91757 6
## 20 300.00000 7
## 21 333.00000 10
## 22 350.00000 5
## 23 400.00000 8
## 24 420.00000 5
## 25 459.86261 3
## 26 500.00000 4
## 27 505.84887 3
## 28 529.00000 10
## 29 597.82140 3
## 30 600.00000 8
## 31 630.00000 28
## 32 666.00000 10
## 33 689.79392 4
## 34 700.00000 2
## 35 714.31207 1
## 36 800.00000 9
## 37 839.00000 10
## 38 840.00000 27
## 39 875.25000 3
## 40 1000.00000 9
## 41 1057.00000 10
## 42 1200.00000 8
## 43 1263.00000 8
## 44 1330.00000 10
## 45 1600.00000 8
## 46 1675.00000 10
## 47 1810.00000 5
## 48 2110.00000 10
## 49 2400.00000 8
## 50 2657.00000 10
## 51 3000.00000 1
## 52 3346.00000 10
## 53 3500.00000 2
## 54 3620.00000 5
## 55 4191.60000 1
## 56 5430.00000 5
## 57 NA 22ex$cratc <- NA
for(i in 1:nrow(ex)) {
if(is.na(ex$C_rate[i])) {
ex$cratc[i] <- "idk"
} else if(between(ex$C_rate[i], 0, 100)) {
ex$cratc[i] <- "30-100"
} else if(between(ex$C_rate[i], 101, 160)) {
ex$cratc[i] <- "133-160"
} else if(between(ex$C_rate[i], 161, 200)) {
ex$cratc[i] <- "174-200"
} else if(between(ex$C_rate[i], 201, 300)) {
ex$cratc[i] <- "210-300"
} else if(between(ex$C_rate[i], 301, 400)) {
ex$cratc[i] <- "330-400"
} else if(between(ex$C_rate[i], 401, 500)) {
ex$cratc[i] <- "420-500"
} else if(between(ex$C_rate[i], 501, 600)) {
ex$cratc[i] <- "506-600"
} else if(between(ex$C_rate[i], 630, 700)) {
ex$cratc[i] <- "620-700"
} else if(between(ex$C_rate[i], 714, 999)) {
ex$cratc[i] <- "714-999"
} else if(between(ex$C_rate[i], 1000, 1330)) {
ex$cratc[i] <- "1000-1330"
} else if(between(ex$C_rate[i], 1600, 2000)) {
ex$cratc[i] <- "1600-2000"
} else if(between(ex$C_rate[i], 2001, 3000)) {
ex$cratc[i] <- "2110-3000"
} else if(between(ex$C_rate[i], 3001, 5000)) {
ex$cratc[i] <- "3346-5000"
} else {
ex$cratc[i] <- ">5000"
}
}
ex$cratc <- factor(ex$cratc, levels = c("30-100", "133-160", "174-200", "210-300", "330-400", "420-500",
"506-600", "620-700", "714-999", "1000-1330", "1600-2000",
"2110-3000", "3346-5000", ">5000"))count(ex, cratc)##
## cratc n
## 1 30-100 45
## 2 133-160 56
## 3 174-200 60
## 4 210-300 32
## 5 330-400 23
## 6 420-500 12
## 7 506-600 24
## 8 620-700 44
## 9 714-999 50
## 10 1000-1330 45
## 11 1600-2000 23
## 12 2110-3000 29
## 13 3346-5000 18
## 14 >5000 5
## 15 <NA> 22any(is.na(ex$cratc))## [1] TRUEex.cratc.na <- ex %>%
filter(!is.na(cratc))
ex.mv.cratc <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex.cratc.na,
mods = ~ factor(cratc) - 1)
ex.mv.cratc.na <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex.cratc.na)h <- c("C rate (all)", "Estimate [95% CI]")
single.mod.cat(ex.cratc.na, "cratc", ex.mv.cratc.na, ex.mv.cratc, h)##
## Multivariate Meta-Analysis Model (k = 466; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.2248 0.4741 44 no exp_ID
## sigma^2.2 0.1923 0.4386 466 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 452) = 891.0899, p-val < .0001
##
## Test of Moderators (coefficients 1:14):
## QM(df = 14) = 66.9383, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## factor(cratc)30-100 -0.3758 0.1693 -2.2193 0.0265 -0.7076 -0.0439
## factor(cratc)133-160 -0.2434 0.1730 -1.4065 0.1596 -0.5825 0.0958
## factor(cratc)174-200 -0.0934 0.2160 -0.4323 0.6656 -0.5167 0.3300
## factor(cratc)210-300 -0.4915 0.1839 -2.6723 0.0075 -0.8521 -0.1310
## factor(cratc)330-400 -0.5412 0.2210 -2.4490 0.0143 -0.9743 -0.1081
## factor(cratc)420-500 -0.8747 0.3096 -2.8254 0.0047 -1.4815 -0.2679
## factor(cratc)506-600 -0.4821 0.2198 -2.1940 0.0282 -0.9129 -0.0514
## factor(cratc)620-700 -0.9629 0.1846 -5.2159 <.0001 -1.3247 -0.6011
## factor(cratc)714-999 -0.7834 0.2004 -3.9094 <.0001 -1.1761 -0.3906
## factor(cratc)1000-1330 -0.6346 0.1752 -3.6226 0.0003 -0.9780 -0.2913
## factor(cratc)1600-2000 -0.8816 0.2202 -4.0029 <.0001 -1.3133 -0.4499
## factor(cratc)2110-3000 -0.7311 0.2014 -3.6302 0.0003 -1.1259 -0.3364
## factor(cratc)3346-5000 -0.8978 0.2416 -3.7170 0.0002 -1.3713 -0.4244
## factor(cratc)>5000 -1.6735 0.4591 -3.6449 0.0003 -2.5734 -0.7736
##
## factor(cratc)30-100 *
## factor(cratc)133-160
## factor(cratc)174-200
## factor(cratc)210-300 **
## factor(cratc)330-400 *
## factor(cratc)420-500 **
## factor(cratc)506-600 *
## factor(cratc)620-700 ***
## factor(cratc)714-999 ***
## factor(cratc)1000-1330 ***
## factor(cratc)1600-2000 ***
## factor(cratc)2110-3000 ***
## factor(cratc)3346-5000 ***
## factor(cratc)>5000 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0698708976533511
Outliers removed
out.cratc <- c(437)
ex.cratc.o <- outliers(ex.cratc.na, out.cratc)
ex.mv.cratc.o <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex.cratc.o,
mods = ~ factor(cratc) - 1)
ex.mv.cratc.na.o <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex.cratc.o)h <- c("C rate (outliers removed)", "Estimate [95% CI]")
single.mod.cat(ex.cratc.o, "cratc", ex.mv.cratc.na.o, ex.mv.cratc.o, h)##
## Multivariate Meta-Analysis Model (k = 465; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.2293 0.4789 44 no exp_ID
## sigma^2.2 0.1785 0.4225 465 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 451) = 879.8370, p-val < .0001
##
## Test of Moderators (coefficients 1:14):
## QM(df = 14) = 72.9270, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## factor(cratc)30-100 -0.3756 0.1685 -2.2289 0.0258 -0.7060 -0.0453
## factor(cratc)133-160 -0.2455 0.1726 -1.4226 0.1548 -0.5837 0.0927
## factor(cratc)174-200 -0.0940 0.2161 -0.4351 0.6635 -0.5176 0.3295
## factor(cratc)210-300 -0.4929 0.1829 -2.6945 0.0070 -0.8514 -0.1344
## factor(cratc)330-400 -0.5423 0.2198 -2.4674 0.0136 -0.9730 -0.1115
## factor(cratc)420-500 -0.8706 0.3085 -2.8219 0.0048 -1.4753 -0.2659
## factor(cratc)506-600 -0.4849 0.2185 -2.2189 0.0265 -0.9131 -0.0566
## factor(cratc)620-700 -0.9554 0.1842 -5.1863 <.0001 -1.3165 -0.5944
## factor(cratc)714-999 -0.7820 0.1996 -3.9182 <.0001 -1.1731 -0.3908
## factor(cratc)1000-1330 -0.6337 0.1744 -3.6335 0.0003 -0.9756 -0.2919
## factor(cratc)1600-2000 -0.8794 0.2189 -4.0172 <.0001 -1.3084 -0.4503
## factor(cratc)2110-3000 -0.7295 0.2001 -3.6452 0.0003 -1.1217 -0.3373
## factor(cratc)3346-5000 -0.8951 0.2400 -3.7292 0.0002 -1.3655 -0.4246
## factor(cratc)>5000 -2.3387 0.5103 -4.5826 <.0001 -3.3390 -1.3385
##
## factor(cratc)30-100 *
## factor(cratc)133-160
## factor(cratc)174-200
## factor(cratc)210-300 **
## factor(cratc)330-400 *
## factor(cratc)420-500 **
## factor(cratc)506-600 *
## factor(cratc)620-700 ***
## factor(cratc)714-999 ***
## factor(cratc)1000-1330 ***
## factor(cratc)1600-2000 ***
## factor(cratc)2110-3000 ***
## factor(cratc)3346-5000 ***
## factor(cratc)>5000 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0765425838259561
Comparison
- Removing outliers does not change significance.
- Removing outliers increases QM/QT.
C applications total (categorical)
# Create categories
count(ex, C_app)##
## C_app n
## 1 1.0 283
## 2 2.0 20
## 3 3.0 19
## 4 3.5 1
## 5 4.0 27
## 6 5.0 2
## 7 6.0 36
## 8 7.0 3
## 9 8.0 3
## 10 9.0 43
## 11 10.0 8
## 12 15.0 1
## 13 18.0 1
## 14 22.0 1
## 15 32.0 28
## 16 36.0 6
## 17 39.5 2
## 18 66.0 1
## 19 71.0 1
## 20 104.0 2ex$capc <- NA
for(i in 1:nrow(ex)) {
if(ex$C_app[i] == 1) {
ex$capc[i] <- "1"
} else if(ex$C_app[i] == 2) {
ex$capc[i] <- "2"
} else if(between(ex$C_app[i], 3, 6)) {
ex$capc[i] <- "3-6"
} else if(between(ex$C_app[i], 7, 10)) {
ex$capc[i] <- "7-10"
} else if(between(ex$C_app[i], 15, 22)) {
ex$capc[i] <- "15-22"
} else if(between(ex$C_app[i], 32, 40)) {
ex$capc[i] <- "32-40"
} else {
ex$capc[i] <- ">40"
}
}
ex$capc <- factor(ex$capc, levels = c("1", "2", "3-6", "7-10", "15-22", "32-40", ">40")))count(ex, capc)##
## capc n
## 1 1 283
## 2 2 20
## 3 3-6 85
## 4 7-10 57
## 5 15-22 3
## 6 32-40 36
## 7 >40 4any(is.na(ex$capc))## [1] FALSEex.mv.capc <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex,
mods = ~ factor(capc) - 1)h <- c("C applications", "Estimate [95% CI]")
single.mod.cat(ex, "capc", ex.mv, ex.mv.capc, h)##
## Multivariate Meta-Analysis Model (k = 488; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.3538 0.5948 50 no exp_ID
## sigma^2.2 0.1639 0.4048 488 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 481) = 956.6665, p-val < .0001
##
## Test of Moderators (coefficients 1:7):
## QM(df = 7) = 49.7203, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## factor(capc)1 -0.4563 0.1508 -3.0265 0.0025 -0.7518 -0.1608 **
## factor(capc)2 -0.5667 0.2748 -2.0618 0.0392 -1.1054 -0.0280 *
## factor(capc)3-6 -1.0618 0.1774 -5.9865 <.0001 -1.4095 -0.7142 ***
## factor(capc)7-10 -0.6289 0.2568 -2.4493 0.0143 -1.1321 -0.1256 *
## factor(capc)15-22 -0.1230 0.4666 -0.2637 0.7920 -1.0375 0.7914
## factor(capc)32-40 -0.1734 0.2678 -0.6475 0.5173 -0.6983 0.3515
## factor(capc)>40 0.5305 0.4757 1.1151 0.2648 -0.4019 1.4629
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0494047972209333
Months applying C
count(ex, C_app_tm)##
## C_app_tm n
## 1 0.00 283
## 2 1.00 8
## 3 2.00 1
## 4 4.00 15
## 5 4.25 4
## 6 5.00 1
## 7 6.00 5
## 8 7.00 7
## 9 9.00 27
## 10 10.00 2
## 11 12.00 5
## 12 15.00 23
## 13 16.00 1
## 14 18.00 10
## 15 20.00 1
## 16 26.00 5
## 17 29.00 4
## 18 37.00 1
## 19 38.25 28
## 20 41.00 43
## 21 76.00 1
## 22 80.00 2
## 23 84.00 2
## 24 121.00 1
## 25 148.00 2
## 26 205.00 6ex$capt <- NA
for(i in 1:nrow(ex)) {
if(ex$C_app_tm[i] == 0) {
ex$capt[i] <- "1 application total"
} else if(between(ex$C_app_tm[i], 0.001, 6)) {
ex$capt[i] <- "1-6"
} else if(between(ex$C_app_tm[i], 7, 12)) {
ex$capt[i] <- "7-12"
} else if(between(ex$C_app_tm[i], 13, 24)) {
ex$capt[i] <- "13-24"
} else if(between(ex$C_app_tm[i], 25, 36)) {
ex$capt[i] <- "25-36"
} else if(between(ex$C_app_tm[i], 37, 48)) {
ex$capt[i] <- "37-48"
} else {
ex$capt[i] <- ">48"
}
}
ex$capt <- factor(ex$capt, levels = c("1 application total", "1-6", "7-12", "13-24",
"25-36", "37-48", ">48"))count(ex, capt)##
## capt n
## 1 1 application total 283
## 2 1-6 34
## 3 7-12 41
## 4 13-24 35
## 5 25-36 9
## 6 37-48 72
## 7 >48 14any(is.na(ex$capt))## [1] FALSEex.mv.capt <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex,
mods = ~ factor(capt) - 1)h <- c("Months applying C", "Estimate [95% CI]")
single.mod.cat(ex, "capt", ex.mv, ex.mv.capt, h)##
## Multivariate Meta-Analysis Model (k = 488; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1225 0.3500 50 no exp_ID
## sigma^2.2 0.1805 0.4248 488 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 481) = 890.9178, p-val < .0001
##
## Test of Moderators (coefficients 1:7):
## QM(df = 7) = 98.5005, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb
## factor(capt)1 application total -0.4234 0.1117 -3.7900 0.0002 -0.6423
## factor(capt)1-6 -0.3997 0.1799 -2.2219 0.0263 -0.7523
## factor(capt)7-12 -0.7820 0.2071 -3.7769 0.0002 -1.1879
## factor(capt)13-24 -1.3965 0.1974 -7.0739 <.0001 -1.7835
## factor(capt)25-36 -1.3745 0.3338 -4.1171 <.0001 -2.0288
## factor(capt)37-48 -0.1728 0.1789 -0.9662 0.3339 -0.5234
## factor(capt)>48 0.3709 0.2782 1.3330 0.1825 -0.1744
## ci.ub
## factor(capt)1 application total -0.2044 ***
## factor(capt)1-6 -0.0471 *
## factor(capt)7-12 -0.3762 ***
## factor(capt)13-24 -1.0096 ***
## factor(capt)25-36 -0.7202 ***
## factor(capt)37-48 0.1778
## factor(capt)>48 0.9163
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0995539865238045
Annual/perennial
any(is.na(ex$plant_anper))## [1] FALSEex.mv.panp <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex,
mods = ~ factor(plant_anper) - 1)h <- c("Annual/perennial", "Estimate [95% CI]")
single.mod.cat(ex, "plant_anper", ex.mv, ex.mv.panp, h)##
## Multivariate Meta-Analysis Model (k = 488; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.3738 0.6114 50 no exp_ID
## sigma^2.2 0.1692 0.4114 488 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 484) = 1003.4371, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## QM(df = 4) = 48.8764, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb
## factor(plant_anper)annual -0.6748 0.1172 -5.7561 <.0001 -0.9045
## factor(plant_anper)biennial -0.5771 0.3632 -1.5890 0.1121 -1.2889
## factor(plant_anper)perennial -0.2597 0.1325 -1.9609 0.0499 -0.5193
## factor(plant_anper)unknown -0.6758 0.1504 -4.4937 <.0001 -0.9706
## ci.ub
## factor(plant_anper)annual -0.4450 ***
## factor(plant_anper)biennial 0.1347
## factor(plant_anper)perennial -0.0001 *
## factor(plant_anper)unknown -0.3811 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0464465767305553
Grass/forb/shrub
All data
any(is.na(ex$plant_gfs))## [1] FALSEex.mv.pgfs <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex,
mods = ~ factor(plant_gfs) - 1)h <- c("Grass/forb/shrub (all)", "Estimate [95% CI]")
single.mod.cat(ex, "plant_gfs", ex.mv, ex.mv.pgfs, h)##
## Multivariate Meta-Analysis Model (k = 488; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.4246 0.6516 50 no exp_ID
## sigma^2.2 0.1714 0.4140 488 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 484) = 1034.4447, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## QM(df = 4) = 33.0747, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## factor(plant_gfs)forb -0.4828 0.1256 -3.8448 0.0001 -0.7289 -0.2367
## factor(plant_gfs)graminoid -0.6413 0.1196 -5.3606 <.0001 -0.8757 -0.4068
## factor(plant_gfs)shrub -0.1644 0.4495 -0.3657 0.7146 -1.0454 0.7167
## factor(plant_gfs)unknown -0.7345 0.2830 -2.5955 0.0094 -1.2891 -0.1798
##
## factor(plant_gfs)forb ***
## factor(plant_gfs)graminoid ***
## factor(plant_gfs)shrub
## factor(plant_gfs)unknown **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0309827235397018
Outliers removed
out.pgfs <- c(257, 258)
ex.pgfs.o <- outliers(ex, out.pgfs)
ex.mv.pgfs.o <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex.pgfs.o,
mods = ~ factor(plant_gfs) - 1)
ex.mv.pgfs.na.o <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex.pgfs.o)h <- c("Grass/forb/shrub (outliers removed)", "Estimate [95% CI]")
single.mod.cat(ex.pgfs.o, "plant_gfs", ex.mv.pgfs.na.o, ex.mv.pgfs.o, h)##
## Multivariate Meta-Analysis Model (k = 486; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.3287 0.5733 49 no exp_ID
## sigma^2.2 0.1610 0.4012 486 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 482) = 998.8359, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## QM(df = 4) = 33.1288, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## factor(plant_gfs)forb -0.4683 0.1156 -4.0501 <.0001 -0.6949 -0.2417
## factor(plant_gfs)graminoid -0.6169 0.1096 -5.6299 <.0001 -0.8317 -0.4021
## factor(plant_gfs)shrub -0.1526 0.4395 -0.3472 0.7285 -1.0140 0.7088
## factor(plant_gfs)unknown -0.4416 0.2712 -1.6285 0.1034 -0.9731 0.0899
##
## factor(plant_gfs)forb ***
## factor(plant_gfs)graminoid ***
## factor(plant_gfs)shrub
## factor(plant_gfs)unknown
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0321026485785864
Comparison
- Removing outliers changes unknown from significant to not significant.
- Removing outliers increases QM/QT.
Annual/perennial and grass/forb/shrub
any(is.na(ex$plant_apgfs))## [1] FALSEex.mv.papgfs <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex,
mods = ~ factor(plant_apgfs) - 1)h <- c("Annual/perennial and grass/forb/shrub", "Estimate [95% CI]")
single.mod.cat(ex, "plant_apgfs", ex.mv, ex.mv.papgfs, h)##
## Multivariate Meta-Analysis Model (k = 488; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.4111 0.6412 50 no exp_ID
## sigma^2.2 0.1464 0.3826 488 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 478) = 973.9460, p-val < .0001
##
## Test of Moderators (coefficients 1:10):
## QM(df = 10) = 68.7655, p-val < .0001
##
## Model Results:
##
## estimate se zval pval
## factor(plant_apgfs)annual forb -0.6646 0.1376 -4.8287 <.0001
## factor(plant_apgfs)annual graminoid -0.5792 0.1304 -4.4407 <.0001
## factor(plant_apgfs)biennial forb -0.4980 0.3578 -1.3919 0.1640
## factor(plant_apgfs)perennial forb -0.0556 0.1583 -0.3512 0.7254
## factor(plant_apgfs)perennial graminoid -0.4768 0.1636 -2.9146 0.0036
## factor(plant_apgfs)perennial shrub -0.4567 0.7628 -0.5987 0.5494
## factor(plant_apgfs)unknown forb -0.4049 0.1893 -2.1386 0.0325
## factor(plant_apgfs)unknown graminoid -1.3667 0.2634 -5.1889 <.0001
## factor(plant_apgfs)unknown shrub -0.1915 0.5486 -0.3491 0.7270
## factor(plant_apgfs)unknown unknown -0.7030 0.2772 -2.5361 0.0112
## ci.lb ci.ub
## factor(plant_apgfs)annual forb -0.9344 -0.3948 ***
## factor(plant_apgfs)annual graminoid -0.8348 -0.3235 ***
## factor(plant_apgfs)biennial forb -1.1993 0.2033
## factor(plant_apgfs)perennial forb -0.3659 0.2547
## factor(plant_apgfs)perennial graminoid -0.7974 -0.1562 **
## factor(plant_apgfs)perennial shrub -1.9518 1.0384
## factor(plant_apgfs)unknown forb -0.7760 -0.0338 *
## factor(plant_apgfs)unknown graminoid -1.8830 -0.8505 ***
## factor(plant_apgfs)unknown shrub -1.2667 0.8837
## factor(plant_apgfs)unknown unknown -1.2462 -0.1597 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0659487637848055
Plot size (categorical)
count(ex, plot)##
## plot n
## 1 7.853982e-03 1
## 2 9.000000e-02 1
## 3 2.500000e-01 9
## 4 1.000000e+00 55
## 5 1.300000e+00 3
## 6 1.500000e+00 3
## 7 2.000000e+00 4
## 8 2.250000e+00 11
## 9 3.750000e+00 8
## 10 4.000000e+00 136
## 11 6.000000e+00 18
## 12 9.000000e+00 10
## 13 1.200000e+01 86
## 14 1.250000e+01 1
## 15 1.600000e+01 24
## 16 2.500000e+01 69
## 17 5.000000e+01 6
## 18 7.500000e+01 7
## 19 1.000000e+02 28
## 20 1.600000e+02 2
## 21 NA 6ex$plotc <- NA
for(i in 1:nrow(ex)) {
if(is.na(ex$plot[i])) {
ex$plotc[i] <- "idk"
} else if(ex$plot[i] < 1) {
ex$plotc[i] <- "<1"
} else if(ex$plot[i] == 1) {
ex$plotc[i] <- "1"
} else if(between(ex$plot[i], 1.3, 3.75)) {
ex$plotc[i] <- "1.3-3.75"
} else if(ex$plot[i] == 4) {
ex$plotc[i] <- "4"
} else if(between(ex$plot[i], 6, 9)) {
ex$plotc[i] <- "6-9"
} else if(ex$plot[i] == 12) {
ex$plotc[i] <- "12"
} else if(between(ex$plot[i], 12.5, 16)) {
ex$plotc[i] <- "12.5-16"
} else if(ex$plot[i] == 25) {
ex$plotc[i] <- "25"
} else if(between(ex$plot[i], 28, 50)) {
ex$plotc[i] <- "28-50"
} else if(ex$plot[i] == 75) {
ex$plotc[i] <- "75"
} else {
ex$plotc[i] <- ">100"
}
}
ex$plotc <- factor(ex$plotc,
levels = c("<1", "1", "1.3-3.75", "4", "6-9", "12", "12.5-16",
"25", "28-50", "75", ">100"))count(ex, plotc)##
## plotc n
## 1 <1 11
## 2 1 55
## 3 1.3-3.75 29
## 4 4 136
## 5 6-9 28
## 6 12 86
## 7 12.5-16 25
## 8 25 69
## 9 28-50 6
## 10 75 7
## 11 >100 30
## 12 <NA> 6any(is.na(ex$plotc))## [1] TRUEex.plotc.na <- ex %>%
filter(!is.na(plotc))
ex.mv.plotc <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex.plotc.na,
mods = ~ factor(plotc) - 1)
ex.mv.plotc.na <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex.plotc.na)h <- c("Plot size", "Estimate [95% CI]")
single.mod.cat(ex.plotc.na, "plotc", ex.mv.plotc.na, ex.mv.plotc, h)##
## Multivariate Meta-Analysis Model (k = 482; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.4717 0.6868 47 no exp_ID
## sigma^2.2 0.1777 0.4216 482 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 471) = 937.3447, p-val < .0001
##
## Test of Moderators (coefficients 1:11):
## QM(df = 11) = 37.8372, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## factor(plotc)<1 -0.4909 0.4858 -1.0106 0.3122 -1.4430 0.4612
## factor(plotc)1 -0.9030 0.2287 -3.9487 <.0001 -1.3513 -0.4548 ***
## factor(plotc)1.3-3.75 -0.8331 0.2928 -2.8450 0.0044 -1.4070 -0.2591 **
## factor(plotc)4 -0.8002 0.4212 -1.8998 0.0575 -1.6257 0.0253 .
## factor(plotc)6-9 -0.7799 0.3535 -2.2060 0.0274 -1.4728 -0.0870 *
## factor(plotc)12 -0.4064 0.4966 -0.8183 0.4132 -1.3798 0.5670
## factor(plotc)12.5-16 -0.7359 0.4679 -1.5729 0.1157 -1.6529 0.1811
## factor(plotc)25 -0.1443 0.3708 -0.3892 0.6971 -0.8710 0.5824
## factor(plotc)28-50 -0.5463 0.7443 -0.7340 0.4630 -2.0050 0.9124
## factor(plotc)75 0.3421 0.5852 0.5846 0.5588 -0.8048 1.4891
## factor(plotc)>100 0.0751 0.3764 0.1996 0.8418 -0.6626 0.8128
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0388001791698277
Seeding of native
All data
count(ex, seeding_native)##
## seeding_native n
## 1 0.0 208
## 2 150.0 4
## 3 262.5 4
## 4 584.0 1
## 5 700.0 1
## 6 1480.0 2
## 7 2250.0 1
## 8 3230.0 1
## 9 4000.0 6
## 10 8000.0 6
## 11 10000.0 9
## 12 NA 245ex <- ex %>%
mutate(seedn = rep(NA, nrow(ex)))
for(i in 1:nrow(ex)) {
if(is.na(ex$seeding_native[i])) {
ex$seedn[i] <- "native seeded"
} else if(ex$seeding_native[i] == 0) {
ex$seedn[i] <- "native not seeded"
} else {
ex$seedn[i] <- "native seeded"
}
}count(ex, seedn)##
## seedn n
## 1 native not seeded 208
## 2 native seeded 280any(is.na(ex$seedn))## [1] FALSEex.mv.seedn <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex,
mods = ~ factor(seedn) - 1)h <- c("Seeding of native (all)", "Estimate [95% CI]")
single.mod.cat(ex, "seedn", ex.mv, ex.mv.seedn, h)##
## Multivariate Meta-Analysis Model (k = 488; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.4133 0.6429 50 no exp_ID
## sigma^2.2 0.1710 0.4135 488 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 486) = 1017.5916, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## QM(df = 2) = 29.4290, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb
## factor(seedn)native not seeded -0.6118 0.1311 -4.6652 <.0001 -0.8689
## factor(seedn)native seeded -0.5483 0.1446 -3.7911 0.0002 -0.8318
## ci.ub
## factor(seedn)native not seeded -0.3548 ***
## factor(seedn)native seeded -0.2648 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0281074128742978
Outliers removed
out.seedn <- NA
ex.seedn.o <- outliers(ex, out.seedn)
ex.mv.seedn.o <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex.seedn.o,
mods = ~ factor(seedn) - 1)
ex.mv.seedn.na.o <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = ex.seedn.o)h <- c("Seeding of native (outliers removed)", "Estimate [95% CI]")
single.mod.cat(ex.seedn.o, "seedn", ex.mv.seedn.na.o, ex.mv.seedn.o, h)##
## Multivariate Meta-Analysis Model (k = 487; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.3668 0.6057 50 no exp_ID
## sigma^2.2 0.1668 0.4084 487 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 485) = 1004.8059, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## QM(df = 2) = 33.4249, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb
## factor(seedn)native not seeded -0.6120 0.1263 -4.8469 <.0001 -0.8594
## factor(seedn)native seeded -0.5771 0.1396 -4.1354 <.0001 -0.8506
## ci.ub
## factor(seedn)native not seeded -0.3645 ***
## factor(seedn)native seeded -0.3036 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0321940894475952