Native 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 native 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
(
- 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 “Native outliers for models” document for
details;
code_native-outliers.html).- When only one model is presented, no outliers were removed.
Setup
library(tidyverse)
library(metafor)
# Load data ---------------------------------------------------------------
nt <- read.csv("data/cleaned/native-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.cratc <- c(440)No moderators (summary effect)
nt.mv <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt)print(nt.mv)##
## Multivariate Meta-Analysis Model (k = 655; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1155 0.3398 53 no exp_ID
## sigma^2.2 0.0000 0.0001 655 no exp_ID/obs_ID
##
## Test for Heterogeneity:
## Q(df = 654) = 925.7573, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.0742 0.0592 -1.2543 0.2097 -0.1902 0.0417
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1predict(nt.mv)##
## pred se ci.lb ci.ub pi.lb pi.ub
## -0.0742 0.0592 -0.1902 0.0417 -0.7502 0.6018# QM/QT
nt.mv[["QM"]] / sum(nt.mv[["QE"]], nt.mv[["QM"]])## [1] 0.001696612# I^2
W <- diag(1 / nt$vi)
X <- model.matrix(nt.mv)
P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W
100 * sum(nt.mv$sigma2) / (sum(nt.mv$sigma2) + (nt.mv$k-nt.mv$p)/sum(diag(P)))## [1] 20.53777# I^2 across nested experiment and within nested experiment
100 * nt.mv$sigma2 / (sum(nt.mv$sigma2) + (nt.mv$k-nt.mv$p)/sum(diag(P)))## [1] 2.053777e+01 1.511454e-06# Probability based on PI (IntHout 2016)
D <- 0 # threshold D (null effect)
M <- nt.mv[["b"]] # summary treatment effect
SD.PI <- predict(nt.mv)$se * sqrt(nt.mv[["k"]]) # SD of prediction interval, calculated from SE
t <- as.numeric((D - M) / (SD.PI))
pt(t, nt.mv[["k"]] - 1, lower.tail = FALSE) # right-tail cumulative t-distribution with k-1 df## [1] 0.480463Biome
any(is.na(nt$biome)) # if FALSE, use nt.mv for no-moderator model## [1] FALSE # but still requires new model for outliers removed, if applicable nt.mv.bio <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt,
mods = ~ factor(biome) - 1)h <- c("Biome", "Estimate [95% CI]")
single.mod.cat(nt, "biome", nt.mv, nt.mv.bio, h)##
## Multivariate Meta-Analysis Model (k = 655; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0510 0.2259 53 no exp_ID
## sigma^2.2 0.0000 0.0001 655 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 641) = 801.1230, p-val < .0001
##
## Test of Moderators (coefficients 1:14):
## QM(df = 14) = 48.7937, p-val < .0001
##
## Model Results:
##
## estimate se zval pval
## factor(biome)coastal grassland -0.0032 0.1389 -0.0232 0.9815
## factor(biome)coastal sandplain grassland 0.2487 0.6321 0.3934 0.6940
## factor(biome)cold desert -0.5495 0.3230 -1.7014 0.0889
## factor(biome)creosote scrub -1.5178 0.3980 -3.8137 0.0001
## factor(biome)fescue prairie 0.9606 0.2770 3.4677 0.0005
## factor(biome)grassy woodland -0.0295 0.1613 -0.1830 0.8548
## factor(biome)mixed-grass prairie -0.3743 0.1694 -2.2097 0.0271
## factor(biome)Puget lowland prairie -0.0784 0.3740 -0.2096 0.8339
## factor(biome)sagebrush steppe -0.0762 0.0830 -0.9178 0.3587
## factor(biome)sedge meadow 0.0572 0.2397 0.2384 0.8116
## factor(biome)shortgrass steppe -0.1856 0.1230 -1.5089 0.1313
## factor(biome)tallgrass prairie 0.3279 0.1634 2.0068 0.0448
## factor(biome)upland prairie -0.2855 0.1736 -1.6445 0.1001
## factor(biome)wetland prairie 0.5439 0.2621 2.0754 0.0379
## ci.lb ci.ub
## factor(biome)coastal grassland -0.2754 0.2690
## factor(biome)coastal sandplain grassland -0.9903 1.4876
## factor(biome)cold desert -1.1825 0.0835 .
## factor(biome)creosote scrub -2.2979 -0.7378 ***
## factor(biome)fescue prairie 0.4176 1.5035 ***
## factor(biome)grassy woodland -0.3457 0.2867
## factor(biome)mixed-grass prairie -0.7064 -0.0423 *
## factor(biome)Puget lowland prairie -0.8114 0.6546
## factor(biome)sagebrush steppe -0.2389 0.0865
## factor(biome)sedge meadow -0.4127 0.5270
## factor(biome)shortgrass steppe -0.4267 0.0555
## factor(biome)tallgrass prairie 0.0076 0.6482 *
## factor(biome)upland prairie -0.6259 0.0548
## factor(biome)wetland prairie 0.0303 1.0576 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0574100216171447
Region
any(is.na(nt$region))## [1] FALSEnt.mv.reg <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt,
mods = ~ factor(region) - 1)h <- c("Region", "Estimate [95% CI]")
single.mod.cat(nt, "region", nt.mv, nt.mv.reg, h)##
## Multivariate Meta-Analysis Model (k = 655; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0726 0.2695 53 no exp_ID
## sigma^2.2 0.0000 0.0001 655 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 645) = 841.7813, p-val < .0001
##
## Test of Moderators (coefficients 1:10):
## QM(df = 10) = 32.4791, p-val = 0.0003
##
## Model Results:
##
## estimate se zval pval
## factor(region)central Canada 0.5886 0.2138 2.7534 0.0059
## factor(region)coastal CA, USA -0.0119 0.1556 -0.0763 0.9392
## factor(region)Great Basin, USA -0.1259 0.0985 -1.2789 0.2009
## factor(region)Great Plains, USA -0.2504 0.1379 -1.8161 0.0694
## factor(region)Midwest, USA 0.2374 0.1523 1.5589 0.1190
## factor(region)northeast, USA 0.2487 0.6490 0.3832 0.7016
## factor(region)Pacific Northwest, USA -0.0091 0.1375 -0.0661 0.9473
## factor(region)Rocky Mountains, USA -0.3713 0.1752 -2.1190 0.0341
## factor(region)southeast Australia -0.0526 0.1771 -0.2970 0.7665
## factor(region)southwest, USA -1.5178 0.4243 -3.5775 0.0003
## ci.lb ci.ub
## factor(region)central Canada 0.1696 1.0076 **
## factor(region)coastal CA, USA -0.3169 0.2932
## factor(region)Great Basin, USA -0.3189 0.0671
## factor(region)Great Plains, USA -0.5206 0.0198 .
## factor(region)Midwest, USA -0.0611 0.5358
## factor(region)northeast, USA -1.0233 1.5207
## factor(region)Pacific Northwest, USA -0.2787 0.2605
## factor(region)Rocky Mountains, USA -0.7146 -0.0279 *
## factor(region)southeast Australia -0.3996 0.2944
## factor(region)southwest, USA -2.3494 -0.6863 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0371503181401346
Soil suborder
any(is.na(nt$soil_suborder)) # if TRUE, run no-moderators model on subset with no missing data for that moderator (i.e. nt.mv.soil.na)## [1] TRUEnt.soil.na <- nt %>%
filter(!is.na(soil_suborder))
nt.mv.soil <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt.soil.na,
mods = ~ factor(soil_suborder) - 1)
nt.mv.soil.na <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt.soil.na)h <- c("Soil suborder", "Estimate [95% CI]")
single.mod.cat(nt, "soil_suborder", nt.mv.soil.na, nt.mv.soil, h)##
## Multivariate Meta-Analysis Model (k = 480; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0635 0.2519 39 no exp_ID
## sigma^2.2 0.0165 0.1286 480 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 463) = 619.4585, p-val < .0001
##
## Test of Moderators (coefficients 1:17):
## QM(df = 17) = 32.0858, p-val = 0.0147
##
## Model Results:
##
## estimate se zval pval
## factor(soil_suborder)Aqualfs 0.5456 0.2867 1.9027 0.0571
## factor(soil_suborder)Aquolls -0.2513 0.1976 -1.2720 0.2034
## factor(soil_suborder)Argids -0.3060 0.1951 -1.5680 0.1169
## factor(soil_suborder)Blacktown -0.0096 0.2461 -0.0389 0.9689
## factor(soil_suborder)Cambids -0.1669 0.1508 -1.1071 0.2682
## factor(soil_suborder)Chernozems 0.5745 0.2099 2.7366 0.0062
## factor(soil_suborder)Durids -0.2951 0.3027 -0.9749 0.3296
## factor(soil_suborder)Humults -0.3437 0.2590 -1.3271 0.1845
## factor(soil_suborder)Psamments, Udepts 0.2491 0.6483 0.3843 0.7008
## factor(soil_suborder)Udolls 0.6425 0.2274 2.8259 0.0047
## factor(soil_suborder)Ustalfs 0.0042 0.1936 0.0219 0.9825
## factor(soil_suborder)Ustepts -0.6665 0.4604 -1.4477 0.1477
## factor(soil_suborder)Ustolls -0.2258 0.1452 -1.5555 0.1198
## factor(soil_suborder)Vertisols -0.2006 0.4930 -0.4068 0.6842
## factor(soil_suborder)Xerands -0.1586 0.5374 -0.2951 0.7679
## factor(soil_suborder)Xerepts -0.0065 0.5146 -0.0126 0.9900
## factor(soil_suborder)Xerolls -0.0523 0.2153 -0.2430 0.8080
## ci.lb ci.ub
## factor(soil_suborder)Aqualfs -0.0164 1.1076 .
## factor(soil_suborder)Aquolls -0.6386 0.1359
## factor(soil_suborder)Argids -0.6884 0.0765
## factor(soil_suborder)Blacktown -0.4920 0.4728
## factor(soil_suborder)Cambids -0.4624 0.1286
## factor(soil_suborder)Chernozems 0.1631 0.9860 **
## factor(soil_suborder)Durids -0.8883 0.2982
## factor(soil_suborder)Humults -0.8513 0.1639
## factor(soil_suborder)Psamments, Udepts -1.0215 1.5198
## factor(soil_suborder)Udolls 0.1969 1.0882 **
## factor(soil_suborder)Ustalfs -0.3753 0.3838
## factor(soil_suborder)Ustepts -1.5688 0.2358
## factor(soil_suborder)Ustolls -0.5104 0.0587
## factor(soil_suborder)Vertisols -1.1669 0.7658
## factor(soil_suborder)Xerands -1.2119 0.8947
## factor(soil_suborder)Xerepts -1.0150 1.0021
## factor(soil_suborder)Xerolls -0.4744 0.3697
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0492457784373502
Duration since first C app (categorical)
# Create categories
count(nt, duration_first)##
## duration_first n
## 1 3.0 5
## 2 5.0 3
## 3 6.0 2
## 4 8.0 5
## 5 9.0 20
## 6 11.0 39
## 7 12.0 3
## 8 15.0 89
## 9 16.0 5
## 10 17.5 5
## 11 18.5 176
## 12 20.0 1
## 13 22.0 60
## 14 24.0 6
## 15 27.0 2
## 16 29.0 12
## 17 30.0 12
## 18 30.5 47
## 19 39.0 109
## 20 45.0 23
## 21 46.0 5
## 22 50.0 6
## 23 112.0 2
## 24 116.0 3
## 25 120.0 4
## 26 207.0 3
## 27 266.0 8nt$dfc <- NA
for(i in 1:nrow(nt)) {
if(nt$duration_first[i] == 3) {
nt$dfc[i] <- "3"
} else if(between(nt$duration_first[i], 5, 6)) {
nt$dfc[i] <- "5-6"
} else if(between(nt$duration_first[i], 7, 12)) {
nt$dfc[i] <- "7-12"
} else if(between(nt$duration_first[i], 13, 18)) {
nt$dfc[i] <- "13-18"
} else if(between(nt$duration_first[i], 19, 24)) {
nt$dfc[i] <- "19-24"
} else if(between(nt$duration_first[i], 25, 36)) {
nt$dfc[i] <- "25-36"
} else if(between(nt$duration_first[i], 37, 50)) {
nt$dfc[i] <- "37-50"
} else if(between(nt$duration_first[i], 100, 200)) {
nt$dfc[i] <- "100-200"
} else {
nt$dfc[i] <- ">200"
}
}
nt$dfc <- factor(nt$df, levels = c("3", "5-6", "7-12", "12-18", "19-24", "25-36", "37-50", "100-200", ">200"))count(nt, dfc)##
## dfc n
## 1 3 5
## 2 5-6 5
## 3 7-12 67
## 4 13-18 99
## 5 19-24 67
## 6 25-36 73
## 7 37-50 143
## 8 100-200 9
## 9 >200 187any(is.na(nt$dfc))## [1] FALSEnt.mv.dfc <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt,
mods = ~ factor(dfc) - 1)h <- c("Duration since first C app", "Estimate [95% CI]")
single.mod.cat(nt, "dfc", nt.mv, nt.mv.dfc, h)##
## Multivariate Meta-Analysis Model (k = 655; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1347 0.3670 53 no exp_ID
## sigma^2.2 0.0000 0.0001 655 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 646) = 884.8217, p-val < .0001
##
## Test of Moderators (coefficients 1:9):
## QM(df = 9) = 13.6220, p-val = 0.1364
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## factor(dfc)3 0.0971 0.3489 0.2783 0.7808 -0.5867 0.7808
## factor(dfc)5-6 -0.5350 0.2945 -1.8165 0.0693 -1.1122 0.0423 .
## factor(dfc)7-12 0.1162 0.1833 0.6338 0.5262 -0.2431 0.4755
## factor(dfc)13-18 -0.2736 0.1559 -1.7544 0.0794 -0.5793 0.0321 .
## factor(dfc)19-24 0.2901 0.1766 1.6423 0.1005 -0.0561 0.6362
## factor(dfc)25-36 -0.1930 0.1408 -1.3713 0.1703 -0.4689 0.0829
## factor(dfc)37-50 0.0041 0.1327 0.0307 0.9755 -0.2561 0.2642
## factor(dfc)100-200 0.0616 0.4163 0.1480 0.8823 -0.7544 0.8776
## factor(dfc)>200 -0.2101 0.1116 -1.8825 0.0598 -0.4289 0.0086 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.015161816260635
Duration since last C app (categorical)
# Create categories
count(nt, duration_last)##
## duration_last n
## 1 0.50 6
## 2 0.75 108
## 3 1.00 12
## 4 1.50 1
## 5 2.00 78
## 6 3.00 9
## 7 3.50 1
## 8 4.00 31
## 9 5.00 2
## 10 6.00 2
## 11 9.00 13
## 12 12.00 1
## 13 14.50 217
## 14 15.00 76
## 15 16.00 1
## 16 22.00 58
## 17 24.00 6
## 18 30.00 5
## 19 36.00 9
## 20 46.00 5
## 21 49.00 6
## 22 118.00 8nt$dlc <- NA
for(i in 1:nrow(nt)) {
if(between(nt$duration_last[i], 0, 1.5)) {
nt$dlc[i] <- "0-1.5"
} else if(nt$duration_last[i] == 2) {
nt$dlc[i] <- "2"
} else if(between(nt$duration_last[i], 3, 3.5)) {
nt$dlc[i] <- "3-3.5"
} else if(between(nt$duration_last[i], 4, 6)) {
nt$dlc[i] <- "4-6"
} else if(between(nt$duration_last[i], 7, 12)) {
nt$dlc[i] <- "7-12"
} else if(between(nt$duration_last[i], 13, 18)) {
nt$dlc[i] <- "13-18"
} else if(between(nt$duration_last[i], 19, 24)) {
nt$dlc[i] <- "19-24"
} else if(between(nt$duration_last[i], 36, 49)) {
nt$dlc[i] <- "36-49"
} else {
nt$dlc[i] <- ">100"
}
}
nt$dlc <- factor(nt$dl, levels = c("0-1.5", "2", "3-3.5", "4-6", "7-12", "13-18", "19-24",
"36-49", ">100"))count(nt, dlc)##
## dlc n
## 1 0-1.5 127
## 2 2 78
## 3 3-3.5 10
## 4 4-6 35
## 5 7-12 14
## 6 13-18 294
## 7 19-24 64
## 8 36-49 20
## 9 >100 13any(is.na(nt$dlc))## [1] FALSEnt.mv.dlc <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt,
mods = ~ factor(dlc) - 1)h <- c("Duration since last C app", "Estimate [95% CI]")
single.mod.cat(nt, "dlc", nt.mv, nt.mv.dlc, h)##
## Multivariate Meta-Analysis Model (k = 655; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1126 0.3355 53 no exp_ID
## sigma^2.2 0.0000 0.0001 655 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 646) = 883.2405, p-val < .0001
##
## Test of Moderators (coefficients 1:9):
## QM(df = 9) = 23.0335, p-val = 0.0061
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## factor(dlc)0-1.5 -0.4035 0.1256 -3.2132 0.0013 -0.6496 -0.1574 **
## factor(dlc)2 0.1358 0.1447 0.9387 0.3479 -0.1477 0.4193
## factor(dlc)3-3.5 -0.0622 0.2416 -0.2576 0.7967 -0.5357 0.4113
## factor(dlc)4-6 -0.4023 0.1666 -2.4143 0.0158 -0.7289 -0.0757 *
## factor(dlc)7-12 -0.1896 0.1958 -0.9684 0.3329 -0.5732 0.1941
## factor(dlc)13-18 -0.0072 0.1112 -0.0645 0.9486 -0.2251 0.2108
## factor(dlc)19-24 0.3342 0.1892 1.7668 0.0773 -0.0365 0.7050 .
## factor(dlc)36-49 0.2199 0.2364 0.9299 0.3524 -0.2436 0.6833
## factor(dlc)>100 0.1702 0.3157 0.5393 0.5897 -0.4485 0.7890
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0254156098123526
C type
any(is.na(nt$C_type))## [1] FALSEnt.mv.ctyp <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt,
mods = ~ factor(C_type) - 1)h <- c("C_type", "Estimate [95% CI]")
single.mod.cat(nt, "C_type", nt.mv, nt.mv.ctyp, h)##
## Multivariate Meta-Analysis Model (k = 655; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0500 0.2235 53 no exp_ID
## sigma^2.2 0.0000 0.0001 655 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 644) = 813.9969, p-val < .0001
##
## Test of Moderators (coefficients 1:11):
## QM(df = 11) = 56.4027, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb
## factor(C_type)dextrose -2.8141 0.8548 -3.2922 0.0010 -4.4895
## factor(C_type)lignin -0.3549 0.4583 -0.7744 0.4387 -1.2532
## factor(C_type)lignin AND sucrose -0.8719 0.6199 -1.4066 0.1596 -2.0869
## factor(C_type)mulch -1.2537 0.4557 -2.7512 0.0059 -2.1468
## factor(C_type)sawdust 0.0079 0.1145 0.0690 0.9450 -0.2165
## factor(C_type)sawdust AND woodchip 1.0550 0.7909 1.3340 0.1822 -0.4950
## factor(C_type)straw 0.7935 0.2486 3.1916 0.0014 0.3062
## factor(C_type)sucrose -0.1769 0.0569 -3.1107 0.0019 -0.2884
## factor(C_type)sucrose AND mulch -1.5518 0.6380 -2.4321 0.0150 -2.8023
## factor(C_type)sucrose AND sawdust 0.2737 0.1123 2.4377 0.0148 0.0536
## factor(C_type)woodchip 0.8793 0.5826 1.5092 0.1312 -0.2626
## ci.ub
## factor(C_type)dextrose -1.1388 ***
## factor(C_type)lignin 0.5434
## factor(C_type)lignin AND sucrose 0.3431
## factor(C_type)mulch -0.3605 **
## factor(C_type)sawdust 0.2323
## factor(C_type)sawdust AND woodchip 2.6051
## factor(C_type)straw 1.2807 **
## factor(C_type)sucrose -0.0654 **
## factor(C_type)sucrose AND mulch -0.3013 *
## factor(C_type)sucrose AND sawdust 0.4938 *
## factor(C_type)woodchip 2.0212
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0648009480291643
C rate (categorical)
All data
count(nt, C_rate)##
## C_rate n
## 1 30.00000 2
## 2 50.00000 1
## 3 76.33719 1
## 4 84.00000 4
## 5 93.60000 1
## 6 100.00000 1
## 7 133.00000 5
## 8 150.00000 223
## 9 160.00000 115
## 10 173.50000 1
## 11 183.94505 8
## 12 195.00000 23
## 13 200.00000 1
## 14 210.00000 7
## 15 229.93131 4
## 16 250.00000 1
## 17 260.40000 2
## 18 333.00000 4
## 19 350.00000 9
## 20 420.00000 3
## 21 459.86261 4
## 22 500.00000 1
## 23 505.84887 5
## 24 529.00000 5
## 25 597.82140 6
## 26 630.00000 18
## 27 666.00000 5
## 28 689.79392 4
## 29 700.00000 3
## 30 714.31207 1
## 31 800.00000 2
## 32 839.00000 5
## 33 840.00000 39
## 34 875.25000 2
## 35 1000.00000 6
## 36 1057.00000 7
## 37 1330.00000 5
## 38 1675.00000 6
## 39 1810.00000 4
## 40 2110.00000 7
## 41 2529.24437 5
## 42 2657.00000 7
## 43 3000.00000 4
## 44 3346.00000 7
## 45 3500.00000 2
## 46 3620.00000 5
## 47 4191.60000 54
## 48 5430.00000 5
## 49 NA 15nt$cratc <- NA
for(i in 1:nrow(nt)) {
if(is.na(nt$C_rate[i])) {
nt$cratc[i] <- "idk"
} else if(between(nt$C_rate[i], 0, 100)) {
nt$cratc[i] <- "30-100"
} else if(between(nt$C_rate[i], 101, 160)) {
nt$cratc[i] <- "133-160"
} else if(between(nt$C_rate[i], 161, 200)) {
nt$cratc[i] <- "174-200"
} else if(between(nt$C_rate[i], 201, 300)) {
nt$cratc[i] <- "210-300"
} else if(between(nt$C_rate[i], 301, 400)) {
nt$cratc[i] <- "330-400"
} else if(between(nt$C_rate[i], 401, 500)) {
nt$cratc[i] <- "420-500"
} else if(between(nt$C_rate[i], 501, 600)) {
nt$cratc[i] <- "506-600"
} else if(between(nt$C_rate[i], 630, 700)) {
nt$cratc[i] <- "620-700"
} else if(between(nt$C_rate[i], 714, 999)) {
nt$cratc[i] <- "714-999"
} else if(between(nt$C_rate[i], 1000, 1330)) {
nt$cratc[i] <- "1000-1330"
} else if(between(nt$C_rate[i], 1600, 2000)) {
nt$cratc[i] <- "1600-2000"
} else if(between(nt$C_rate[i], 2001, 3000)) {
nt$cratc[i] <- "2110-3000"
} else if(between(nt$C_rate[i], 3001, 5000)) {
nt$cratc[i] <- "3346-5000"
} else {
nt$cratc[i] <- ">5000"
}
}
nt$cratc <- factor(nt$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(nt, cratc)##
## cratc n
## 1 30-100 10
## 2 133-160 343
## 3 174-200 33
## 4 210-300 14
## 5 330-400 13
## 6 420-500 8
## 7 506-600 16
## 8 620-700 30
## 9 714-999 49
## 10 1000-1330 18
## 11 1600-2000 10
## 12 2110-3000 23
## 13 3346-5000 68
## 14 >5000 5
## 15 <NA> 15any(is.na(nt$cratc))## [1] TRUEnt.cratc.na <- nt %>%
filter(!is.na(cratc))
nt.mv.cratc <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt.cratc.na,
mods = ~ factor(cratc) - 1)
nt.mv.cratc.na <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt.cratc.na)h <- c("C rate (all)", "Estimate [95% CI]")
single.mod.cat(nt.cratc.na, "cratc", nt.mv.cratc.na, nt.mv.cratc, h)##
## Multivariate Meta-Analysis Model (k = 640; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1420 0.3768 49 no exp_ID
## sigma^2.2 0.0000 0.0000 640 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 626) = 809.0913, p-val < .0001
##
## Test of Moderators (coefficients 1:14):
## QM(df = 14) = 38.6142, p-val = 0.0004
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## factor(cratc)30-100 -0.4572 0.2236 -2.0448 0.0409 -0.8955 -0.0190
## factor(cratc)133-160 -0.1989 0.1019 -1.9515 0.0510 -0.3987 0.0009
## factor(cratc)174-200 -0.1101 0.2052 -0.5366 0.5915 -0.5122 0.2920
## factor(cratc)210-300 -0.4470 0.1984 -2.2535 0.0242 -0.8358 -0.0582
## factor(cratc)330-400 -0.3818 0.2823 -1.3526 0.1762 -0.9350 0.1714
## factor(cratc)420-500 -0.2592 0.2529 -1.0252 0.3053 -0.7548 0.2364
## factor(cratc)506-600 0.0328 0.2530 0.1294 0.8970 -0.4632 0.5287
## factor(cratc)620-700 -0.1735 0.1715 -1.0114 0.3118 -0.5097 0.1627
## factor(cratc)714-999 -0.1063 0.2108 -0.5043 0.6141 -0.5194 0.3069
## factor(cratc)1000-1330 0.2738 0.2060 1.3296 0.1837 -0.1298 0.6775
## factor(cratc)1600-2000 0.4009 0.2474 1.6206 0.1051 -0.0840 0.8858
## factor(cratc)2110-3000 0.6591 0.1920 3.4322 0.0006 0.2827 1.0354
## factor(cratc)3346-5000 0.2382 0.1980 1.2031 0.2289 -0.1498 0.6261
## factor(cratc)>5000 -0.0924 0.3826 -0.2414 0.8092 -0.8422 0.6575
##
## 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.0455513965054985
Outliers removed
any(is.na(nt$cratc))## [1] TRUEnt.cratc.na <- nt %>%
filter(!is.na(cratc))
out.cratc <- c(452)
nt.cratc.o <- outliers(nt.cratc.na, nt.out.cratc)
nt.mv.cratc.o <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt.cratc.o,
mods = ~ factor(cratc) - 1)
nt.mv.cratc.na.o <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt.cratc.o)h <- c("C rate (out rm)", "Estimate [95% CI]")
single.mod.cat(nt.cratc.o, "cratc", nt.mv.cratc.na.o, nt.mv.cratc.o, h)##
## Multivariate Meta-Analysis Model (k = 639; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1412 0.3757 49 no exp_ID
## sigma^2.2 0.0000 0.0000 639 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 625) = 799.4106, p-val < .0001
##
## Test of Moderators (coefficients 1:14):
## QM(df = 14) = 41.8690, p-val = 0.0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## factor(cratc)30-100 -0.4565 0.2235 -2.0428 0.0411 -0.8944 -0.0185
## factor(cratc)133-160 -0.1986 0.1017 -1.9523 0.0509 -0.3980 0.0008
## factor(cratc)174-200 -0.1095 0.2047 -0.5351 0.5926 -0.5108 0.2917
## factor(cratc)210-300 -0.4463 0.1982 -2.2519 0.0243 -0.8347 -0.0579
## factor(cratc)330-400 -0.3804 0.2820 -1.3493 0.1772 -0.9331 0.1722
## factor(cratc)420-500 -0.2588 0.2526 -1.0248 0.3055 -0.7539 0.2362
## factor(cratc)506-600 0.0339 0.2528 0.1340 0.8934 -0.4615 0.5293
## factor(cratc)620-700 -0.1730 0.1713 -1.0102 0.3124 -0.5087 0.1627
## factor(cratc)714-999 -0.1057 0.2105 -0.5022 0.6155 -0.5183 0.3069
## factor(cratc)1000-1330 0.2744 0.2059 1.3331 0.1825 -0.1291 0.6779
## factor(cratc)1600-2000 0.4007 0.2473 1.6206 0.1051 -0.0839 0.8854
## factor(cratc)2110-3000 0.6590 0.1919 3.4346 0.0006 0.2829 1.0351
## factor(cratc)3346-5000 0.2421 0.1977 1.2246 0.2207 -0.1454 0.6297
## factor(cratc)>5000 -0.6026 0.4164 -1.4472 0.1478 -1.4186 0.2135
##
## 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.0497682705209593
C applications total (categorical)
# Create categories
count(nt, C_app)##
## C_app n
## 1 1.0 167
## 2 2.0 240
## 3 3.0 14
## 4 3.5 1
## 5 4.0 46
## 6 5.0 2
## 7 6.0 9
## 8 7.0 2
## 9 8.0 6
## 10 9.0 23
## 11 10.0 17
## 12 15.0 2
## 13 17.0 1
## 14 18.0 2
## 15 32.0 107
## 16 36.0 3
## 17 39.5 3
## 18 71.0 2
## 19 104.0 8nt$capc <- NA
for(i in 1:nrow(nt)) {
if(nt$C_app[i] == 1) {
nt$capc[i] <- "1"
} else if(nt$C_app[i] == 2) {
nt$capc[i] <- "2"
} else if(between(nt$C_app[i], 3, 6)) {
nt$capc[i] <- "3-6"
} else if(between(nt$C_app[i], 7, 10)) {
nt$capc[i] <- "7-10"
} else if(between(nt$C_app[i], 15, 22)) {
nt$capc[i] <- "15-22"
} else if(between(nt$C_app[i], 32, 40)) {
nt$capc[i] <- "32-40"
} else {
nt$capc[i] <- ">40"
}
}
nt$capc <- factor(nt$capc, levels = c("1", "2", "3-6", "7-10", "15-22", "32-40", ">40")))count(nt, capc)##
## capc n
## 1 1 167
## 2 2 240
## 3 3-6 72
## 4 7-10 48
## 5 15-22 5
## 6 32-40 113
## 7 >40 10any(is.na(nt$capc))## [1] FALSEnt.mv.capc <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt,
mods = ~ factor(capc) - 1)h <- c("C applications", "Estimate [95% CI]")
single.mod.cat(nt, "capc", nt.mv, nt.mv.capc, h)##
## Multivariate Meta-Analysis Model (k = 655; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0970 0.3115 53 no exp_ID
## sigma^2.2 0.0000 0.0001 655 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 648) = 875.5744, p-val < .0001
##
## Test of Moderators (coefficients 1:7):
## QM(df = 7) = 12.4477, p-val = 0.0868
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## factor(capc)1 0.1795 0.1040 1.7254 0.0844 -0.0244 0.3833 .
## factor(capc)2 -0.0644 0.0989 -0.6514 0.5148 -0.2582 0.1294
## factor(capc)3-6 -0.2460 0.1187 -2.0720 0.0383 -0.4786 -0.0133 *
## factor(capc)7-10 -0.0540 0.1496 -0.3608 0.7183 -0.3471 0.2392
## factor(capc)15-22 -0.3843 0.2710 -1.4183 0.1561 -0.9154 0.1468
## factor(capc)32-40 -0.2213 0.1444 -1.5325 0.1254 -0.5043 0.0617
## factor(capc)>40 -0.1259 0.2979 -0.4227 0.6725 -0.7097 0.4579
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0140172789583436
Months applying C
count(nt, C_app_tm)##
## C_app_tm n
## 1 0.00 167
## 2 1.00 6
## 3 2.00 1
## 4 3.50 1
## 5 4.00 176
## 6 5.00 1
## 7 6.00 5
## 8 7.00 11
## 9 9.00 39
## 10 10.00 2
## 11 12.00 5
## 12 13.00 3
## 13 14.00 1
## 14 15.00 4
## 15 15.25 1
## 16 15.50 4
## 17 16.00 47
## 18 18.00 6
## 19 20.00 2
## 20 26.00 4
## 21 27.00 12
## 22 29.00 5
## 23 37.00 2
## 24 38.25 107
## 25 41.00 23
## 26 76.00 2
## 27 80.00 3
## 28 84.00 4
## 29 148.00 8
## 30 205.00 3nt$capt <- NA
for(i in 1:nrow(nt)) {
if(nt$C_app_tm[i] == 0) {
nt$capt[i] <- "1 application total"
} else if(between(nt$C_app_tm[i], 0.001, 6)) {
nt$capt[i] <- "1-6"
} else if(between(nt$C_app_tm[i], 7, 12)) {
nt$capt[i] <- "7-12"
} else if(between(nt$C_app_tm[i], 13, 24)) {
nt$capt[i] <- "13-24"
} else if(between(nt$C_app_tm[i], 25, 36)) {
nt$capt[i] <- "25-36"
} else if(between(nt$C_app_tm[i], 37, 48)) {
nt$capt[i] <- "37-48"
} else {
nt$capt[i] <- ">48"
}
}
nt$capt <- factor(nt$capt, levels = c("1 application total", "1-6", "7-12", "13-24",
"25-36", "37-48", ">48"))count(nt, capt)##
## capt n
## 1 1 application total 167
## 2 1-6 190
## 3 7-12 57
## 4 13-24 68
## 5 25-36 21
## 6 37-48 132
## 7 >48 20any(is.na(nt$capt))## [1] FALSEnt.mv.capt <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt,
mods = ~ factor(capt) - 1)h <- c("Months applying C", "Estimate [95% CI]")
single.mod.cat(nt, "capt", nt.mv, nt.mv.capt, h)##
## Multivariate Meta-Analysis Model (k = 655; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0962 0.3101 53 no exp_ID
## sigma^2.2 0.0000 0.0001 655 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 648) = 866.5161, p-val < .0001
##
## Test of Moderators (coefficients 1:7):
## QM(df = 7) = 19.8604, p-val = 0.0059
##
## Model Results:
##
## estimate se zval pval ci.lb
## factor(capt)1 application total 0.2301 0.1056 2.1782 0.0294 0.0231
## factor(capt)1-6 -0.1934 0.0982 -1.9693 0.0489 -0.3859
## factor(capt)7-12 0.2619 0.1680 1.5593 0.1189 -0.0673
## factor(capt)13-24 -0.2551 0.1133 -2.2521 0.0243 -0.4772
## factor(capt)25-36 -0.5304 0.2255 -2.3523 0.0187 -0.9724
## factor(capt)37-48 -0.0723 0.1299 -0.5569 0.5776 -0.3268
## factor(capt)>48 -0.1842 0.2108 -0.8735 0.3824 -0.5974
## ci.ub
## factor(capt)1 application total 0.4372 *
## factor(capt)1-6 -0.0009 *
## factor(capt)7-12 0.5912
## factor(capt)13-24 -0.0331 *
## factor(capt)25-36 -0.0885 *
## factor(capt)37-48 0.1822
## factor(capt)>48 0.2291
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0224063183990693
Annual/perennial
any(is.na(nt$plant_anper))## [1] FALSEnt.mv.panp <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt,
mods = ~ factor(plant_anper) - 1)h <- c("Annual/perennial", "Estimate [95% CI]")
single.mod.cat(nt, "plant_anper", nt.mv, nt.mv.panp, h)##
## Multivariate Meta-Analysis Model (k = 655; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1034 0.3215 53 no exp_ID
## sigma^2.2 0.0004 0.0198 655 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 651) = 909.1067, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## QM(df = 4) = 10.6297, p-val = 0.0311
##
## Model Results:
##
## estimate se zval pval ci.lb
## factor(plant_anper)annual -0.3458 0.1155 -2.9940 0.0028 -0.5722
## factor(plant_anper)biennial -0.2311 0.3566 -0.6481 0.5169 -0.9301
## factor(plant_anper)perennial -0.0303 0.0608 -0.4991 0.6177 -0.1494
## factor(plant_anper)unknown -0.1615 0.1271 -1.2701 0.2041 -0.4106
## ci.ub
## factor(plant_anper)annual -0.1194 **
## factor(plant_anper)biennial 0.4678
## factor(plant_anper)perennial 0.0888
## factor(plant_anper)unknown 0.0877
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0115573747349454
Grass/forb/shrub
any(is.na(nt$plant_gfs))## [1] FALSEnt.mv.pgfs <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt,
mods = ~ factor(plant_gfs) - 1)h <- c("Grass/forb/shrub", "Estimate [95% CI]")
single.mod.cat(nt, "plant_gfs", nt.mv, nt.mv.pgfs, h)##
## Multivariate Meta-Analysis Model (k = 655; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1199 0.3463 53 no exp_ID
## sigma^2.2 0.0000 0.0001 655 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 650) = 918.8050, p-val < .0001
##
## Test of Moderators (coefficients 1:5):
## QM(df = 5) = 4.0902, p-val = 0.5365
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## factor(plant_gfs)forb -0.0256 0.0777 -0.3296 0.7417 -0.1779 0.1267
## factor(plant_gfs)graminoid -0.0683 0.0672 -1.0165 0.3094 -0.1999 0.0634
## factor(plant_gfs)shrub -0.0248 0.1214 -0.2045 0.8380 -0.2627 0.2130
## factor(plant_gfs)subshrub -0.2650 0.2280 -1.1627 0.2450 -0.7118 0.1817
## factor(plant_gfs)unknown -0.1700 0.1122 -1.5152 0.1297 -0.3900 0.0499
##
## factor(plant_gfs)forb
## factor(plant_gfs)graminoid
## factor(plant_gfs)shrub
## factor(plant_gfs)subshrub
## factor(plant_gfs)unknown
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.00443192828975222
Annual/perennial and grass/forb/shrub
any(is.na(nt$plant_apgfs))## [1] FALSEnt.mv.papgfs <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt,
mods = ~ factor(plant_apgfs) - 1)h <- c("Annual/perennial and grass/forb/shrub", "Estimate [95% CI]")
single.mod.cat(nt, "plant_apgfs", nt.mv, nt.mv.papgfs, h)##
## Multivariate Meta-Analysis Model (k = 655; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1175 0.3428 53 no exp_ID
## sigma^2.2 0.0000 0.0001 655 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 642) = 881.1929, p-val < .0001
##
## Test of Moderators (coefficients 1:13):
## QM(df = 13) = 37.2272, p-val = 0.0004
##
## Model Results:
##
## estimate se zval pval
## factor(plant_apgfs)annual forb -0.2320 0.1349 -1.7198 0.0855
## factor(plant_apgfs)annual graminoid -0.5032 0.3049 -1.6507 0.0988
## factor(plant_apgfs)annual unknown -1.0497 0.3415 -3.0735 0.0021
## factor(plant_apgfs)biennial forb -0.1921 0.3585 -0.5359 0.5920
## factor(plant_apgfs)perennial forb 0.0913 0.0853 1.0710 0.2842
## factor(plant_apgfs)perennial graminoid -0.0455 0.0696 -0.6534 0.5135
## factor(plant_apgfs)perennial shrub -0.1170 0.1259 -0.9297 0.3525
## factor(plant_apgfs)perennial subshrub -0.2922 0.2286 -1.2781 0.2012
## factor(plant_apgfs)perennial unknown -0.0703 0.1290 -0.5447 0.5859
## factor(plant_apgfs)unknown forb -0.4295 0.1811 -2.3718 0.0177
## factor(plant_apgfs)unknown graminoid -0.0435 0.2038 -0.2136 0.8308
## factor(plant_apgfs)unknown shrub 1.3462 0.4482 3.0037 0.0027
## factor(plant_apgfs)unknown unknown -0.0871 0.2680 -0.3249 0.7453
## ci.lb ci.ub
## factor(plant_apgfs)annual forb -0.4964 0.0324 .
## factor(plant_apgfs)annual graminoid -1.1008 0.0943 .
## factor(plant_apgfs)annual unknown -1.7191 -0.3803 **
## factor(plant_apgfs)biennial forb -0.8948 0.5105
## factor(plant_apgfs)perennial forb -0.0758 0.2584
## factor(plant_apgfs)perennial graminoid -0.1819 0.0910
## factor(plant_apgfs)perennial shrub -0.3637 0.1297
## factor(plant_apgfs)perennial subshrub -0.7402 0.1559
## factor(plant_apgfs)perennial unknown -0.3231 0.1825
## factor(plant_apgfs)unknown forb -0.7844 -0.0746 *
## factor(plant_apgfs)unknown graminoid -0.4431 0.3560
## factor(plant_apgfs)unknown shrub 0.4678 2.2247 **
## factor(plant_apgfs)unknown unknown -0.6122 0.4381
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.040533941179534
Plot size (categorical)
count(nt, plot)##
## plot n
## 1 7.853982e-03 2
## 2 2.500000e-01 15
## 3 1.000000e+00 89
## 4 1.300000e+00 2
## 5 1.500000e+00 4
## 6 2.000000e+00 2
## 7 2.250000e+00 7
## 8 3.750000e+00 223
## 9 4.000000e+00 88
## 10 6.000000e+00 12
## 11 9.000000e+00 10
## 12 1.600000e+01 20
## 13 2.500000e+01 48
## 14 5.000000e+01 6
## 15 7.500000e+01 6
## 16 1.000000e+02 107
## 17 1.600000e+02 8
## 18 NA 6nt$plotc <- NA
for(i in 1:nrow(nt)) {
if(is.na(nt$plot[i])) {
nt$plotc[i] <- "idk"
} else if(nt$plot[i] < 1) {
nt$plotc[i] <- "<1"
} else if(nt$plot[i] == 1) {
nt$plotc[i] <- "1"
} else if(between(nt$plot[i], 1.3, 3.75)) {
nt$plotc[i] <- "1.3-3.75"
} else if(nt$plot[i] == 4) {
nt$plotc[i] <- "4"
} else if(between(nt$plot[i], 6, 9)) {
nt$plotc[i] <- "6-9"
} else if(nt$plot[i] == 12) {
nt$plotc[i] <- "12"
} else if(between(nt$plot[i], 12.5, 16)) {
nt$plotc[i] <- "12.5-16"
} else if(nt$plot[i] == 25) {
nt$plotc[i] <- "25"
} else if(between(nt$plot[i], 28, 50)) {
nt$plotc[i] <- "28-50"
} else if(nt$plot[i] == 75) {
nt$plotc[i] <- "75"
} else {
nt$plotc[i] <- ">100"
}
}
nt$plotc <- factor(nt$plotc,
levels = c("<1", "1", "1.3-3.75", "4", "6-9", "12", "12.5-16",
"25", "28-50", "75", ">100"))count(nt, plotc)##
## plotc n
## 1 <1 17
## 2 1 89
## 3 1.3-3.75 238
## 4 4 88
## 5 6-9 22
## 6 12.5-16 20
## 7 25 48
## 8 28-50 6
## 9 75 6
## 10 >100 115
## 11 <NA> 6any(is.na(nt$plotc))## [1] TRUEnt.plotc.na <- nt %>%
filter(!is.na(plotc))
nt.mv.plotc <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt.plotc.na,
mods = ~ factor(plotc) - 1)
nt.mv.plotc.na <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt.plotc.na)h <- c("Plot size", "Estimate [95% CI]")
single.mod.cat(nt.plotc.na, "plotc", nt.mv.plotc.na, nt.mv.plotc, h)##
## Multivariate Meta-Analysis Model (k = 649; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0827 0.2875 50 no exp_ID
## sigma^2.2 0.0000 0.0001 649 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 639) = 824.0525, p-val < .0001
##
## Test of Moderators (coefficients 1:10):
## QM(df = 10) = 18.5659, p-val = 0.0461
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## factor(plotc)<1 0.5137 0.2551 2.0136 0.0441 0.0137 1.0137 *
## factor(plotc)1 -0.1142 0.1239 -0.9218 0.3566 -0.3570 0.1286
## factor(plotc)1.3-3.75 -0.2338 0.0918 -2.5473 0.0109 -0.4136 -0.0539 *
## factor(plotc)4 0.0708 0.1754 0.4039 0.6863 -0.2729 0.4146
## factor(plotc)6-9 -0.0528 0.2197 -0.2402 0.8101 -0.4833 0.3778
## factor(plotc)12.5-16 0.3573 0.2602 1.3731 0.1697 -0.1527 0.8672
## factor(plotc)25 0.0476 0.2125 0.2241 0.8227 -0.3688 0.4641
## factor(plotc)28-50 -0.5495 0.3687 -1.4904 0.1361 -1.2722 0.1731
## factor(plotc)75 0.0546 0.3103 0.1758 0.8604 -0.5535 0.6626
## factor(plotc)>100 -0.2317 0.1486 -1.5588 0.1190 -0.5229 0.0596
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0220336291807025
Seeding of native
count(nt, seeding_native)##
## seeding_native n
## 1 0.0 174
## 2 150.0 1
## 3 160.0 5
## 4 225.0 54
## 5 262.5 5
## 6 300.0 217
## 7 584.0 2
## 8 1300.0 6
## 9 3230.0 4
## 10 4000.0 4
## 11 8000.0 4
## 12 10000.0 15
## 13 NA 164nt <- nt %>%
mutate(seedn = rep(NA, nrow(nt)))
for(i in 1:nrow(nt)) {
if(is.na(nt$seeding_native[i])) {
nt$seedn[i] <- "native seeded"
} else if(nt$seeding_native[i] == 0) {
nt$seedn[i] <- "native not seeded"
} else {
nt$seedn[i] <- "native seeded"
}
}any(is.na(nt$seedn))## [1] FALSEnt.mv.seedn <- rma.mv(yi = yi,
V = vi,
random = ~ 1 | exp_ID / obs_ID,
data = nt,
mods = ~ factor(seedn) - 1)h <- c("Seeding of native", "Estimate [95% CI]")
single.mod.cat(nt, "seedn", nt.mv, nt.mv.seedn, h)##
## Multivariate Meta-Analysis Model (k = 655; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0972 0.3118 53 no exp_ID
## sigma^2.2 0.0000 0.0001 655 no exp_ID/obs_ID
##
## Test for Residual Heterogeneity:
## QE(df = 653) = 894.7390, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## QM(df = 2) = 10.2446, p-val = 0.0060
##
## Model Results:
##
## estimate se zval pval ci.lb
## factor(seedn)native not seeded -0.2840 0.0915 -3.1050 0.0019 -0.4632
## factor(seedn)native seeded 0.0335 0.0663 0.5059 0.6130 -0.0964
## ci.ub
## factor(seedn)native not seeded -0.1047 **
## factor(seedn)native seeded 0.1634
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [1] QM/QT = 0.0113202166450891