Water_stelization_diversity
YQF
2025/2/24
1 Library
2 Function
Ftable <- function(model, reorder = F, order, percent = F){
table <- summary(model)$fixed
table2 <-summary(model, prob = 0.9)$fixed %>% as_tibble()
names <- rownames(table)
table <- table %>%
as_tibble() %>%
dplyr::rename('L95%CL' = `l-95% CI`, 'U95%CL' = `u-95% CI`) %>%
dplyr::select(-Rhat, -Bulk_ESS, -Tail_ESS)
table$`L90%CL` <- table2$`l-90% CI`
table$`U90%CL` <- table2$`u-90% CI`
if (percent == T){
table <- table %>%
mutate(percent = exp(as.numeric(Estimate)) - 1)
}
table <-table %>%
mutate(Estimate = ifelse(`L95%CL` * `U95%CL` > 0,
paste(sprintf('%.3f', round(Estimate, 4)), '*', sep = ''),
ifelse(`L90%CL` * `U90%CL` > 0,
paste(sprintf('%.3f', round(Estimate, 4)), '..', sep = ''), round(Estimate, 4))),
rowname = names) %>%
mutate(
rowname = str_replace(rowname, 'Intercept', 'Intercept (Grand Mean)'),
rowname = str_replace(rowname, 'Diversity1', 'Diversity'),
rowname = str_replace(rowname, 'Water1', 'Water'),
rowname = str_replace(rowname, 'Soil1', 'Soil'),
rowname = gsub('(^)([[:alpha:]])', '\\1\\U\\2',rowname, perl = T) # capitalize)
) %>%
dplyr::rename('SE ' = Est.Error) %>%
column_to_rownames(var = 'rowname')
if (reorder == T){
table <- table[order,]
}
table <- table %>%
kable(digits = 4, escape = F, table.attr = "style = \"color: black;\"") %>%
kable_styling(position = "left")%>%
footnote(general = " `L95%CL` * `U95%CL` > 0,significant,`L90%CL` * `U90%CL` > 0,marginally significant ")
return(table)
}plot.resid<-function(model,breaks=30){
plot(resid(model)~fitted(model), xlab = "Fitted values", ylab="Residuals")
hist(resid(model),breaks=breaks)
qqnorm(resid(model))
qqline(resid(model))
}
cbrt <- function(x) { sign(x)*abs(x)^(1/3)}
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
conf.interval=.95, .drop=TRUE) {
library(plyr)
# New version of length which can handle NA's: if na.rm==T, don't count them
length2 <- function (x, na.rm=FALSE) {
if (na.rm) sum(!is.na(x))
else length(x)
}
# This does the summary. For each group's data frame, return a vector with
# N, mean, and sd
datac <- ddply(data, groupvars, .drop=.drop,
.fun = function(xx, col) {
c(N = length2(xx[[col]], na.rm=na.rm),
mean = mean (xx[[col]], na.rm=na.rm),
sd = sd (xx[[col]], na.rm=na.rm)
)
},
measurevar
)
# Rename the "mean" column
#datac <- plyr::rename(datac, c("mean" = measurevar))
names(datac)[names(datac) == 'mean'] <- measurevar
datac$se <- datac$sd / sqrt(datac$N) # Calculate standard error of the mean
# Confidence interval multiplier for standard error
# Calculate t-statistic for confidence interval:
# e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
ciMult <- suppressWarnings(qt(conf.interval/2 + .5, datac$N-1))
datac$ci <- datac$se * ciMult
return(datac)
}theme <- theme(text = element_text(size=22,colour = 'black',family='serif'),
axis.text.x = element_text(size=22,colour = 'black'),
axis.text.y = element_text(size=22,colour = 'black'),
legend.text = element_text(size=22),
strip.text.x = element_text(size = 22, colour = "black"),
axis.line.x=element_line(linetype=1,color="black",size=0.5),
axis.line.y=element_line(linetype=1,color="black",size=0.5),
axis.ticks = element_line(colour = "black",lineend = "square",size = 0.5),
axis.ticks.length.x = unit(+.2,"cm"),
axis.ticks.length.y=unit(+.2,"cm"),
panel.background = element_rect(fill = "transparent"),
panel.grid.minor = element_blank(),panel.grid.major = element_blank(),
legend.background = element_rect(fill="transparent"),
legend.title = element_blank(),
panel.border = element_rect(color = "black", size = 0.5, fill = NA),
plot.background = element_rect(),
legend.key = element_rect(fill = "transparent"),
plot.title = element_text(hjust = 0.5)
)3 Data preparation
data_germination <-read.csv("germination.csv",stringsAsFactors=TRUE,header = T,check.name=FALSE,fileEncoding="latin1")
data_germination <- data_germination %>%
filter(!is.na(data_germination$HarvestNum)) %>%
droplevels()
data_biomass <-read.csv("biomass.csv",stringsAsFactors=TRUE,header = T,check.name=FALSE,fileEncoding="latin1")
data_biomass <- data_biomass %>%
filter(!is.na(data_biomass $Target_biomass),!is.na(data_biomass $Native_Biomass)) %>%
droplevels()4 Target_biomass
4.1 model selection
4.1.1 data transform selection
u<-gl(1,1,length(data_biomass$Target_biomass))
dat1_1 <- lme(Target_biomass ~ Diversity*Water*Soil,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
# weights=varIdent(form= ~ 1|Target_Species),
data= data_biomass, method="REML",
control=list(maxIter=100,msMaxIter=100,niterEM=100))
dat1_2 <- lme(sqrt(Target_biomass) ~ Diversity*Water*Soil,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
# weights=varIdent(form= ~ 1|Target_Species),
data= data_biomass, method="REML",
control=list(maxIter=100,msMaxIter=100,niterEM=100))
dat1_3 <- lme(log(Target_biomass) ~ Diversity*Water*Soil,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
data=data_biomass, method="REML",
control=list(maxIter=100,msMaxIter=100,niterEM=100))
dat1_4 <- lme((Target_biomass)^(1/3) ~ Diversity*Water*Soil,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
# weights=varIdent(form= ~ 1|Target_Species),
data= data_biomass, method="REML",
control=list(maxIter=100,msMaxIter=100,niterEM=100))
dat1_5 <- lme(cbrt(Target_biomass) ~ Diversity*Water*Soil,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
# weights=varIdent(form= ~ 1|Target_Species),
data= data_biomass, method="REML",
control=list(maxIter=100,msMaxIter=100,niterEM=100))opm1a <- par(mfrow = c(3,3), mar = c(4,4,2,2))
plot.resid(dat1_1)
plot.resid(dat1_2)
plot.resid(dat1_3)
## df AIC
## dat1_1 12 2717.0090
## dat1_2 12 1038.6427
## dat1_3 12 1534.0391
## dat1_4 12 392.4718
## dat1_5 12 392.4718
4.1.2 variance structure selection
dat1_3_01 <- lme(log(Target_biomass) ~ Diversity*Water*Soil,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
weights=varIdent(form=~1|Native_Communities),
data=data_biomass, method="REML",
control=list(maxIter=100,msMaxIter=100,niterEM=100))
dat1_3_02 <- lme(log(Target_biomass) ~ Diversity*Water*Soil,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
weights=varIdent(form=~1|Target_Species),
data=data_biomass, method="REML",
control=list(maxIter=100,msMaxIter=100,niterEM=100))
dat1_3_03 <- lme(log(Target_biomass) ~ Diversity*Water*Soil,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
weights=varComb(varIdent(form=~1|Target_Species),
varIdent(form=~1|Native_Communities)),
data=data_biomass, method="REML",
control=list(maxIter=200,msMaxIter=200,niterEM=200))
AIC(dat1_3_01,dat1_3_02,dat1_3_03)## df AIC
## dat1_3_01 27 1542.677
## dat1_3_02 20 1492.435
## dat1_3_03 35 1507.446opm1a <- par(mfrow = c(4,3), mar = c(4,4,2,2))
plot.resid(dat1_3)
plot.resid(dat1_3_01)
plot.resid(dat1_3_02)
plot.resid(dat1_3_03)
## df AIC
## dat1_3 12 1534.039
## dat1_3_01 27 1542.677
## dat1_3_02 20 1492.435
## dat1_3_03 35 1507.446## Warning in model.matrix.default(formula(terms(mf)), mf, contrasts.arg =
## object$contrasts): variable 'Target_Family' is absent, its contrast will be
## ignored## Warning in model.matrix.default(formula(terms(mf)), mf, contrasts.arg =
## object$contrasts): variable 'Target_Species' is absent, its contrast will be
## ignored## Warning in model.matrix.default(formula(terms(mf)), mf, contrasts.arg =
## object$contrasts): variable 'Native_Communities' is absent, its contrast will
## be ignored## R2m R2c
## [1,] 0.02801812 0.5653814## Linear mixed-effects model fit by REML
## Data: data_biomass
## AIC BIC logLik
## 1492.435 1579.136 -726.2177
##
## Random effects:
## Composite Structure: Blocked
##
## Block 1: Target_FamilyAmaranthaceae, Target_FamilyAsteraceae, Target_FamilyFabaceae, Target_FamilyOxalidaceae, Target_FamilyPlantaginaceae, Target_FamilyPoaceae, Target_FamilySolanaceae
## Formula: ~Target_Family - 1 | u
## Structure: Multiple of an Identity
## Target_FamilyAmaranthaceae Target_FamilyAsteraceae
## StdDev: 0.0009649158 0.0009649158
## Target_FamilyFabaceae Target_FamilyOxalidaceae
## StdDev: 0.0009649158 0.0009649158
## Target_FamilyPlantaginaceae Target_FamilyPoaceae
## StdDev: 0.0009649158 0.0009649158
## Target_FamilySolanaceae
## StdDev: 0.0009649158
##
## Block 2: Target_SpeciesAmaranthus_retroflexus, Target_SpeciesBidens_frondosa, Target_SpeciesBidens_pilosa, Target_SpeciesDatura_stramonium, Target_SpeciesLolium_perenne, Target_SpeciesMedicago_sativa, Target_SpeciesOxalis_corymbosa, Target_SpeciesPlantago_virginica, Target_SpeciesSesbania_cannabina
## Formula: ~Target_Species - 1 | u
## Structure: Multiple of an Identity
## Target_SpeciesAmaranthus_retroflexus Target_SpeciesBidens_frondosa
## StdDev: 1.234345 1.234345
## Target_SpeciesBidens_pilosa Target_SpeciesDatura_stramonium
## StdDev: 1.234345 1.234345
## Target_SpeciesLolium_perenne Target_SpeciesMedicago_sativa
## StdDev: 1.234345 1.234345
## Target_SpeciesOxalis_corymbosa Target_SpeciesPlantago_virginica
## StdDev: 1.234345 1.234345
## Target_SpeciesSesbania_cannabina
## StdDev: 1.234345
##
## Block 3: Native_CommunitiesH_NC10, Native_CommunitiesH_NC11, Native_CommunitiesH_NC12, Native_CommunitiesH_NC13, Native_CommunitiesH_NC14, Native_CommunitiesH_NC15, Native_CommunitiesH_NC16, Native_CommunitiesH_NC9, Native_CommunitiesL_NC1, Native_CommunitiesL_NC2, Native_CommunitiesL_NC3, Native_CommunitiesL_NC4, Native_CommunitiesL_NC5, Native_CommunitiesL_NC6, Native_CommunitiesL_NC7, Native_CommunitiesL_NC8
## Formula: ~Native_Communities - 1 | u
## Structure: Multiple of an Identity
## Native_CommunitiesH_NC10 Native_CommunitiesH_NC11
## StdDev: 0.349276 0.349276
## Native_CommunitiesH_NC12 Native_CommunitiesH_NC13
## StdDev: 0.349276 0.349276
## Native_CommunitiesH_NC14 Native_CommunitiesH_NC15
## StdDev: 0.349276 0.349276
## Native_CommunitiesH_NC16 Native_CommunitiesH_NC9
## StdDev: 0.349276 0.349276
## Native_CommunitiesL_NC1 Native_CommunitiesL_NC2 Native_CommunitiesL_NC3
## StdDev: 0.349276 0.349276 0.349276
## Native_CommunitiesL_NC4 Native_CommunitiesL_NC5 Native_CommunitiesL_NC6
## StdDev: 0.349276 0.349276 0.349276
## Native_CommunitiesL_NC7 Native_CommunitiesL_NC8 Residual
## StdDev: 0.349276 0.349276 1.153673
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Target_Species
## Parameter estimates:
## Sesbania_cannabina Bidens_frondosa Amaranthus_retroflexus
## 1.0000000 0.6236618 0.9002894
## Bidens_pilosa Oxalis_corymbosa Medicago_sativa
## 0.6051503 0.5248881 0.8491004
## Datura_stramonium Lolium_perenne Plantago_virginica
## 0.7170674 0.7306140 0.4581722
## Fixed effects: log(Target_biomass) ~ Diversity * Water * Soil
## Value Std.Error DF t-value
## (Intercept) -0.4691664 0.4389378 564 -1.068868
## DiversityLow 0.5110880 0.2150412 564 2.376698
## WaterNormal 0.0045424 0.1259590 564 0.036062
## SoilSterilized -0.3611878 0.1274647 564 -2.833631
## DiversityLow:WaterNormal -0.1047982 0.1777919 564 -0.589443
## DiversityLow:SoilSterilized -0.2049368 0.1788618 564 -1.145783
## WaterNormal:SoilSterilized 0.7270786 0.1791858 564 4.057680
## DiversityLow:WaterNormal:SoilSterilized -0.1094487 0.2521827 564 -0.434006
## p-value
## (Intercept) 0.2856
## DiversityLow 0.0178
## WaterNormal 0.9712
## SoilSterilized 0.0048
## DiversityLow:WaterNormal 0.5558
## DiversityLow:SoilSterilized 0.2524
## WaterNormal:SoilSterilized 0.0001
## DiversityLow:WaterNormal:SoilSterilized 0.6645
## Correlation:
## (Intr) DvrstL WtrNrm SlStrl DvL:WN
## DiversityLow -0.245
## WaterNormal -0.142 0.291
## SoilSterilized -0.141 0.287 0.490
## DiversityLow:WaterNormal 0.101 -0.412 -0.708 -0.347
## DiversityLow:SoilSterilized 0.101 -0.409 -0.349 -0.713 0.495
## WaterNormal:SoilSterilized 0.100 -0.204 -0.703 -0.711 0.498
## DiversityLow:WaterNormal:SoilSterilized -0.071 0.290 0.499 0.505 -0.705
## DvL:SS WtN:SS
## DiversityLow
## WaterNormal
## SoilSterilized
## DiversityLow:WaterNormal
## DiversityLow:SoilSterilized
## WaterNormal:SoilSterilized 0.507
## DiversityLow:WaterNormal:SoilSterilized -0.709 -0.711
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.0973936 -0.6221299 0.1050973 0.6977027 2.6060249
##
## Number of Observations: 572
## Number of Groups: 14.2 likelihood Ratio test
u<-gl(1,1,length(data_biomass$Target_biomass))
Three_Ref_1_1 <- lme(log(Target_biomass) ~ Diversity*Water*Soil,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
weights=varIdent(form=~1|Target_Species),
data=data_biomass, method="ML",
control=list(maxIter=200,msMaxIter=200,niterEM=200))
Three_1_1 <- drop1(Three_Ref_1_1,test="Chisq")[-1,]
# Two-way interaction
Two_Ref_1_1 <- update(Three_Ref_1_1,.~.- Diversity:Water:Soil)
Two_1_1 <- drop1(Two_Ref_1_1,test="Chisq")[-1,]
#Main effect
Main_Ref_1_1 <- update(Two_Ref_1_1,.~.- (Diversity:Water+
Diversity:Soil+
Water:Soil))
Main_1_1 <- drop1(Main_Ref_1_1,test="Chisq")[-1,]4.3 Table
table<-rbind(Main_1_1,Two_1_1,Three_1_1)
table$" "=ifelse(table$`Pr(>Chi)`<0.05,"*",
ifelse(table$`Pr(>Chi)`<0.1,"..",""))#add significant levels
kable(table,digits=4)| Df | AIC | LRT | Pr(>Chi) | ||
|---|---|---|---|---|---|
| Diversity | 1 | 1500.389 | 3.1391 | 0.0764 | .. |
| Water | 1 | 1515.460 | 18.2103 | 0.0000 | * |
| Soil | 1 | 1501.261 | 4.0109 | 0.0452 | * |
| Diversity:Water | 1 | 1473.460 | 1.5960 | 0.2065 | |
| Diversity:Soil | 1 | 1475.896 | 4.0324 | 0.0446 | * |
| Water:Soil | 1 | 1498.158 | 26.2943 | 0.0000 | * |
| Diversity:Water:Soil | 1 | 1473.864 | 0.1901 | 0.6628 |
4.4 Fig_A_pred
data_biomass$Target_pred <- exp(predict(dat1_3_02))
Fig_A <- summarySE(data_biomass, measurevar="Target_pred", groupvars=c("Target_Species","Diversity","Water","Soil"))
Fig_A_01 <- summarySE(Fig_A, measurevar="Target_pred", groupvars=c("Diversity","Water","Soil"))4.4.1 SoilxWaterxDiversity(ns)
pdf('./fig/Target SoilxWaterxDiversity.pdf', height = 15/2.54, width = 15/2.54)
pd <- position_dodge(.3)
Fig_A_01$Soil <- ordered(Fig_A_01$Soil, levels = c('Living', 'Sterilized'))
Fig_A_01$Diversity <- ordered(Fig_A_01$Diversity, levels = c('Low', 'High'))
Fig_A_01$Water <- ordered(Fig_A_01$Water, levels = c('Normal', 'Drought'))
p <- ggplot(Fig_A_01, aes(x=Water, y=Target_pred, color=Soil)) +
ylim(1,6) +
geom_errorbar(aes(ymin=Target_pred-se, ymax=Target_pred+se),
size=2, width=0,position = pd) +
geom_line(colour="darkgrey",linetype="dashed", size=1,
aes(group=Soil), position = pd)+
geom_point(size = 8,position = pd) +
facet_grid( ~Diversity ) +
labs(title="",x="", y = " Alien target species
aboveground biomass")
#remove the background
p <- p + theme_bw()
p <- p + scale_color_manual(breaks = c("Living", "Sterilized"),
values=c("#bcd4e8ff", "#64879aff"))
#change the direction of x-aix labels
ptswd <- p + theme+
theme(#legend.position=c(0.96, 0.5),
legend.background=element_rect(fill=NA),
legend.title=element_text(size=17),
legend.text=element_text(size=17),
legend.spacing.y = unit(0.2, "cm")
)+ theme(legend.title=element_blank())
ptswd <- ptswd + theme(legend.position = c( .85 , .85 ))+coord_cartesian(clip = 'off',ylim = c(0,3))+
theme(axis.line.x=element_line(linetype=1,color="black",size=0),
axis.line.y=element_line(linetype=1,color="black",size=0),
axis.ticks = element_line(colour = "black",lineend = "square",size = 0.5))+
scale_y_continuous(breaks=seq(0,3,1))## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.## png
## 2
## Diversity Water Soil N Target_pred sd se ci
## 1 High Drought Living 9 1.0928244 1.0362062 0.3454021 0.7964986
## 2 High Drought Sterilized 9 0.7588686 0.7245332 0.2415111 0.5569255
## 3 High Normal Living 9 1.1009584 1.0429808 0.3476603 0.8017060
## 4 High Normal Sterilized 9 1.5828053 1.5008016 0.5002672 1.1536182
## 5 Low Drought Living 9 1.8746300 1.7775072 0.5925024 1.3663130
## 6 Low Drought Sterilized 9 1.0642675 1.0091288 0.3363763 0.7756850
## 7 Low Normal Living 9 1.6958016 1.6079437 0.5359812 1.2359749
## 8 Low Normal Sterilized 9 1.7854328 1.6929312 0.5643104 1.30130214.4.2 SoilxWater(*)
pdf('./fig/Target SoilxWater.pdf', height = 15/2.54, width = 15/2.54)
pd <- position_dodge(.3)
Fig_A_01 <- summarySE(Fig_A, measurevar="Target_pred", groupvars=c("Soil","Water"))
#Fig_A_01
Fig_A_01$Soil <- ordered(Fig_A_01$Soil, levels = c('Living', 'Sterilized'))
Fig_A_01$Water <- ordered(Fig_A_01$Water, levels = c('Normal', 'Drought'))
p <- ggplot(Fig_A_01, aes(x=Water, y=Target_pred, color=Soil)) +
ylim(0,2) +
#scale_y_continuous(limits = c(0,4), breaks = c(0,2,4)) +
geom_errorbar(aes(ymin=Target_pred-se, ymax=Target_pred+se),
size=2, width=0,position = pd) +
geom_line(colour="darkgrey",linetype="dashed", size=1,
aes(group=Soil), position = pd)+
geom_point(size = 8,position = pd) +
#facet_grid( ~ Competition) +
labs(title="",x="", y = " Alien target species
aboveground biomass")
#remove the background
p <- p + theme_bw()
p <- p + scale_color_manual(breaks = c("Living", "Sterilized"),
values=c("#bcd4e8ff", "#64879aff"))
#change the direction of x-aix labels
ptsw <- p + theme+
theme(#legend.position=c(0.96, 0.5),
legend.background=element_rect(fill=NA),
legend.title=element_text(size=17),
legend.text=element_text(size=17),
legend.spacing.y = unit(0.2, "cm")
)+ theme(legend.title=element_blank())
ptsw <- ptsw + theme(legend.position = c( .85 , .85 ))+coord_cartesian(clip = 'off',ylim = c(0.5,2.5))+
theme(axis.line.x=element_line(linetype=1,color="black",size=0),
axis.line.y=element_line(linetype=1,color="black",size=0),
axis.ticks = element_line(colour = "black",lineend = "square",size = 0.5))
ptsw
dev.off()## png
## 2
## Soil Water N Target_pred sd se ci
## 1 Living Drought 18 1.483727 1.4676217 0.3459218 0.7298311
## 2 Living Normal 18 1.398380 1.3499152 0.3181781 0.6712970
## 3 Sterilized Drought 18 0.911568 0.8665694 0.2042524 0.4309348
## 4 Sterilized Normal 18 1.684119 1.5554862 0.3666316 0.77352514.4.3 SoilxDiversity(*)
pdf('./fig/Target SoilxDiversity.pdf', height = 15/2.54, width = 15/2.54)
pd <- position_dodge(.4)
Fig_A_01 <- summarySE(Fig_A, measurevar="Target_pred", groupvars=c("Soil","Diversity"))
#Fig_A_01
Fig_A_01$Soil <- ordered(Fig_A_01$Soil, levels = c('Living', 'Sterilized'))
Fig_A_01$Diversity <- ordered(Fig_A_01$Diversity, levels = c('Low', 'High'))
p <- ggplot(Fig_A_01, aes(x=Diversity, y=Target_pred, color=Soil)) +
geom_errorbar(aes(ymin=Target_pred-se, ymax=Target_pred+se),
size=2, width=0,position = pd) +
geom_line(colour="darkgrey",linetype="dashed", size=1,
aes(group=Soil), position = pd)+
geom_point(size = 8,position = pd) +
#facet_grid( ~ Competition) +
labs(title="",x="", y = " Alien target species
aboveground biomass")
#remove the background
p <- p + theme_bw()
p <- p + scale_color_manual(breaks = c("Living", "Sterilized"),
values=c("#bcd4e8ff", "#64879aff"))
#change the direction of x-aix labels
ptsd <- p + theme+
theme(#legend.position=c(0.96, 0.5),
legend.background=element_rect(fill=NA),
legend.title=element_text(size=17),
legend.text=element_text(size=17),
legend.spacing.y = unit(0.2, "cm")
)+ theme(legend.title=element_blank())
ptsd <- ptsd + theme(legend.position = c( .85 , .85 ))+coord_cartesian(clip = 'off',ylim = c(0.5,2.5))+
theme(axis.line.x=element_line(linetype=1,color="black",size=0),
axis.line.y=element_line(linetype=1,color="black",size=0),
axis.ticks = element_line(colour = "black",lineend = "square",size = 0.5))
ptsd
dev.off()## png
## 2
## Soil Diversity N Target_pred sd se ci
## 1 Living High 18 1.096891 1.008568 0.2377217 0.5015490
## 2 Living Low 18 1.785216 1.646816 0.3881582 0.8189421
## 3 Sterilized High 18 1.170837 1.219300 0.2873917 0.6063434
## 4 Sterilized Low 18 1.424850 1.401999 0.3304543 0.69719764.4.4 Soil(*)
## Soil N Target_pred sd se ci
## 1 Living 36 1.441054 1.390382 0.2317304 0.4704377
## 2 Sterilized 36 1.297844 1.301314 0.2168857 0.44030144.4.5 Water(*)
## Water N Target_pred sd se ci
## 1 Drought 36 1.197648 1.222747 0.2037911 0.4137180
## 2 Normal 36 1.541250 1.442672 0.2404453 0.48812995 Native Community Biomass
5.1 model selection
5.1.1 data transform selection
u<-gl(1,1,length(data_biomass$Native_Biomass))
dat1_1 <- lme(Native_Biomass ~ Diversity*Water*Soil,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
data= data_biomass, method="REML",
control=list(maxIter=100,msMaxIter=100,niterEM=100))
dat1_2 <- lme(sqrt(Native_Biomass) ~ Diversity*Water*Soil,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
data= data_biomass, method="REML",
control=list(maxIter=100,msMaxIter=100,niterEM=100))
dat1_3 <- lme(log(Native_Biomass) ~ Diversity*Water*Soil,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
data=data_biomass, method="REML",
control=list(maxIter=100,msMaxIter=100,niterEM=100))
dat1_4 <- lme((Native_Biomass)^(1/4) ~ Diversity*Water*Soil,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
data= data_biomass, method="REML",
control=list(maxIter=100,msMaxIter=100,niterEM=100))
dat1_5 <- lme(cbrt(Native_Biomass) ~ Diversity*Water*Soil,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
data= data_biomass, method="REML",
control=list(maxIter=100,msMaxIter=100,niterEM=100))opm1a <- par(mfrow = c(3,3), mar = c(4,4,2,2))
plot.resid(dat1_1)
plot.resid(dat1_2)
plot.resid(dat1_3)
## df AIC
## dat1_1 12 3776.8278
## dat1_2 12 1257.7658
## dat1_3 12 608.5731
## dat1_4 12 -289.3374
## dat1_5 12 280.3179
5.1.2 variance structure selection
dat1_2_01 <- lme(sqrt(Native_Biomass) ~ Diversity*Water*Soil,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
weights=varIdent(form=~1|Native_Communities),
data=data_biomass, method="REML",
control=list(maxIter=100,msMaxIter=100,niterEM=100))
dat1_2_02 <- lme(sqrt(Native_Biomass) ~ Diversity*Water*Soil,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
weights=varIdent(form=~1|Target_Species),
data=data_biomass, method="REML",
control=list(maxIter=100,msMaxIter=100,niterEM=100))
dat1_2_03 <- lme(sqrt(Native_Biomass) ~ Diversity*Water*Soil,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
weights=varComb(varIdent(form=~1|Target_Species),
varIdent(form=~1|Native_Communities)),
data=data_biomass, method="REML",
control=list(maxIter=200,msMaxIter=200,niterEM=200))opm1a <- par(mfrow = c(4,3), mar = c(4,4,2,2))
plot.resid(dat1_2)
plot.resid(dat1_2_01)
plot.resid(dat1_2_02)
plot.resid(dat1_2_03)
## df AIC
## dat1_2 12 1257.766
## dat1_2_01 27 1141.175
## dat1_2_02 20 1271.806
## dat1_2_03 35 1153.640## Warning in model.matrix.default(formula(terms(mf)), mf, contrasts.arg =
## object$contrasts): variable 'Target_Family' is absent, its contrast will be
## ignored## Warning in model.matrix.default(formula(terms(mf)), mf, contrasts.arg =
## object$contrasts): variable 'Target_Species' is absent, its contrast will be
## ignored## Warning in model.matrix.default(formula(terms(mf)), mf, contrasts.arg =
## object$contrasts): variable 'Native_Communities' is absent, its contrast will
## be ignored## R2m R2c
## [1,] 0.4072467 0.8995693## Linear mixed-effects model fit by REML
## Data: data_biomass
## AIC BIC logLik
## 1141.175 1258.222 -543.5876
##
## Random effects:
## Composite Structure: Blocked
##
## Block 1: Target_FamilyAmaranthaceae, Target_FamilyAsteraceae, Target_FamilyFabaceae, Target_FamilyOxalidaceae, Target_FamilyPlantaginaceae, Target_FamilyPoaceae, Target_FamilySolanaceae
## Formula: ~Target_Family - 1 | u
## Structure: Multiple of an Identity
## Target_FamilyAmaranthaceae Target_FamilyAsteraceae
## StdDev: 4.448578e-05 4.448578e-05
## Target_FamilyFabaceae Target_FamilyOxalidaceae
## StdDev: 4.448578e-05 4.448578e-05
## Target_FamilyPlantaginaceae Target_FamilyPoaceae
## StdDev: 4.448578e-05 4.448578e-05
## Target_FamilySolanaceae
## StdDev: 4.448578e-05
##
## Block 2: Target_SpeciesAmaranthus_retroflexus, Target_SpeciesBidens_frondosa, Target_SpeciesBidens_pilosa, Target_SpeciesDatura_stramonium, Target_SpeciesLolium_perenne, Target_SpeciesMedicago_sativa, Target_SpeciesOxalis_corymbosa, Target_SpeciesPlantago_virginica, Target_SpeciesSesbania_cannabina
## Formula: ~Target_Species - 1 | u
## Structure: Multiple of an Identity
## Target_SpeciesAmaranthus_retroflexus Target_SpeciesBidens_frondosa
## StdDev: 0.1071631 0.1071631
## Target_SpeciesBidens_pilosa Target_SpeciesDatura_stramonium
## StdDev: 0.1071631 0.1071631
## Target_SpeciesLolium_perenne Target_SpeciesMedicago_sativa
## StdDev: 0.1071631 0.1071631
## Target_SpeciesOxalis_corymbosa Target_SpeciesPlantago_virginica
## StdDev: 0.1071631 0.1071631
## Target_SpeciesSesbania_cannabina
## StdDev: 0.1071631
##
## Block 3: Native_CommunitiesH_NC10, Native_CommunitiesH_NC11, Native_CommunitiesH_NC12, Native_CommunitiesH_NC13, Native_CommunitiesH_NC14, Native_CommunitiesH_NC15, Native_CommunitiesH_NC16, Native_CommunitiesH_NC9, Native_CommunitiesL_NC1, Native_CommunitiesL_NC2, Native_CommunitiesL_NC3, Native_CommunitiesL_NC4, Native_CommunitiesL_NC5, Native_CommunitiesL_NC6, Native_CommunitiesL_NC7, Native_CommunitiesL_NC8
## Formula: ~Native_Communities - 1 | u
## Structure: Multiple of an Identity
## Native_CommunitiesH_NC10 Native_CommunitiesH_NC11
## StdDev: 1.00826 1.00826
## Native_CommunitiesH_NC12 Native_CommunitiesH_NC13
## StdDev: 1.00826 1.00826
## Native_CommunitiesH_NC14 Native_CommunitiesH_NC15
## StdDev: 1.00826 1.00826
## Native_CommunitiesH_NC16 Native_CommunitiesH_NC9
## StdDev: 1.00826 1.00826
## Native_CommunitiesL_NC1 Native_CommunitiesL_NC2 Native_CommunitiesL_NC3
## StdDev: 1.00826 1.00826 1.00826
## Native_CommunitiesL_NC4 Native_CommunitiesL_NC5 Native_CommunitiesL_NC6
## StdDev: 1.00826 1.00826 1.00826
## Native_CommunitiesL_NC7 Native_CommunitiesL_NC8 Residual
## StdDev: 1.00826 1.00826 0.4579523
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Native_Communities
## Parameter estimates:
## H_NC10 H_NC13 H_NC11 L_NC6 L_NC7 H_NC14 L_NC8 H_NC16
## 1.0000000 0.7805098 1.6956536 0.9460205 1.2312976 1.4192632 1.7469138 1.9776915
## L_NC3 H_NC12 H_NC9 L_NC2 L_NC1 L_NC5 L_NC4 H_NC15
## 1.4297715 0.8637591 0.6025466 0.7322874 2.3902410 2.6287027 1.7431994 0.8126567
## Fixed effects: sqrt(Native_Biomass) ~ Diversity * Water * Soil
## Value Std.Error DF t-value
## (Intercept) 3.799702 0.3623941 564 10.484999
## DiversityLow -0.616659 0.5125311 564 -1.203163
## WaterNormal 1.010217 0.0712004 564 14.188369
## SoilSterilized 0.404447 0.0707640 564 5.715439
## DiversityLow:WaterNormal -0.257252 0.1191396 564 -2.159246
## DiversityLow:SoilSterilized 0.162719 0.1188794 564 1.368773
## WaterNormal:SoilSterilized 0.874227 0.1003815 564 8.709039
## DiversityLow:WaterNormal:SoilSterilized 0.124408 0.1683032 564 0.739190
## p-value
## (Intercept) 0.0000
## DiversityLow 0.2294
## WaterNormal 0.0000
## SoilSterilized 0.0000
## DiversityLow:WaterNormal 0.0313
## DiversityLow:SoilSterilized 0.1716
## WaterNormal:SoilSterilized 0.0000
## DiversityLow:WaterNormal:SoilSterilized 0.4601
## Correlation:
## (Intr) DvrstL WtrNrm SlStrl DvL:WN
## DiversityLow -0.700
## WaterNormal -0.097 0.068
## SoilSterilized -0.097 0.068 0.495
## DiversityLow:WaterNormal 0.058 -0.116 -0.598 -0.296
## DiversityLow:SoilSterilized 0.058 -0.116 -0.294 -0.595 0.498
## WaterNormal:SoilSterilized 0.068 -0.048 -0.709 -0.705 0.424
## DiversityLow:WaterNormal:SoilSterilized -0.041 0.082 0.423 0.420 -0.708
## DvL:SS WtN:SS
## DiversityLow
## WaterNormal
## SoilSterilized
## DiversityLow:WaterNormal
## DiversityLow:SoilSterilized
## WaterNormal:SoilSterilized 0.420
## DiversityLow:WaterNormal:SoilSterilized -0.706 -0.596
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.536593953 -0.739997053 0.001810396 0.791518172 2.607282388
##
## Number of Observations: 572
## Number of Groups: 15.1.3 likelihood Ratio test
u<-gl(1,1,length(data_biomass$Native_Biomass))
Three_Ref_1_1 <- lme(log(Native_Biomass) ~ Diversity*Water*Soil,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
weights=varIdent(form=~1|Native_Communities),
data=data_biomass, method="ML",
control=list(maxIter=200,msMaxIter=200,niterEM=200))
Three_1_1 <- drop1(Three_Ref_1_1,test="Chisq")[-1,]
# Two-way interaction
Two_Ref_1_1 <- update(Three_Ref_1_1,.~.- Diversity:Water:Soil)
Two_1_1 <- drop1(Two_Ref_1_1,test="Chisq")[-1,]
#Main effect
Main_Ref_1_1 <- update(Two_Ref_1_1,.~.- (Diversity:Water+
Diversity:Soil+
Water:Soil))
Main_1_1 <- drop1(Main_Ref_1_1,test="Chisq")[-1,]5.2 Table
table<-rbind(Main_1_1,Two_1_1,Three_1_1)
table$" "=ifelse(table$`Pr(>Chi)`<0.05,"*",
ifelse(table$`Pr(>Chi)`<0.1,"..",""))#add significant levels
kable(table,digits=4)| Df | AIC | LRT | Pr(>Chi) | ||
|---|---|---|---|---|---|
| Diversity | 1 | 281.9394 | 2.2740 | 0.1316 | |
| Water | 1 | 656.6364 | 376.9710 | 0.0000 | * |
| Soil | 1 | 562.7526 | 283.0872 | 0.0000 | * |
| Diversity:Water | 1 | 224.8159 | 0.0347 | 0.8522 | |
| Diversity:Soil | 1 | 232.5204 | 7.7392 | 0.0054 | * |
| Water:Soil | 1 | 277.3800 | 52.5988 | 0.0000 | * |
| Diversity:Water:Soil | 1 | 226.7812 | 0.0932 | 0.7602 |
5.3 Fig_B_pred
data_biomass$native_pred <- (predict(dat1_2))^2
Fig_B <- summarySE(data_biomass, measurevar="native_pred", groupvars=c("Target_Species","Diversity","Water","Soil"))5.3.1 SoilxWaterxDiversity(ns)
pdf('./fig/NC SoilxWaterxDiversity.pdf', height = 15/2.54, width = 15/2.54)
pd <- position_dodge(.3)
Fig_B_01 <- summarySE(Fig_B, measurevar="native_pred", groupvars=c("Diversity","Water","Soil"))
Fig_B_01$Soil <- ordered(Fig_B_01$Soil, levels = c('Living', 'Sterilized'))
Fig_B_01$Diversity <- ordered(Fig_B_01$Diversity, levels = c('Low', 'High'))
Fig_B_01$Water <- ordered(Fig_B_01$Water, levels = c('Normal', 'Drought'))
p <- ggplot(Fig_B_01, aes(x=Water, y=native_pred, color=Soil)) +
ylim(1,6) +
geom_errorbar(aes(ymin=native_pred-se, ymax=native_pred+se),
size=2, width=0,position = pd) +
geom_line(colour="darkgrey",linetype="dashed", size=1,
aes(group=Soil), position = pd)+
geom_point(size = 8,position = pd) +
facet_grid( ~Diversity ) +
labs(title="",x="", y = " Native community
aboveground biomass")
#remove the background
p <- p + theme_bw()
p <- p + scale_color_manual(breaks = c("Living", "Sterilized"),
values=c("#bcd4e8ff", "#64879aff"))
#change the direction of x-aix labels
pncswd <- p + theme+
theme(#legend.position=c(0.96, 0.5),
legend.background=element_rect(fill=NA),
legend.title=element_text(size=17),
legend.text=element_text(size=17),
legend.spacing.y = unit(0.2, "cm")
)+ theme(legend.title=element_blank())
pncswd <- pncswd + theme(legend.position = c( .85 , .85 ))+
theme(axis.line.x=element_line(linetype=1,color="black",size=0),
axis.line.y=element_line(linetype=1,color="black",size=0),
axis.ticks = element_line(colour = "black",lineend = "square",size = 0.5))+coord_cartesian(clip = 'off',ylim = c(10,40))+
scale_y_continuous(breaks=seq(10,40,10))## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.## png
## 2
## Diversity Water Soil N native_pred sd se ci
## 1 High Drought Living 9 15.47692 0.8852927 0.2950976 0.6804962
## 2 High Drought Sterilized 9 20.03762 1.3724776 0.4574925 1.0549797
## 3 High Normal Living 9 20.91376 1.0931625 0.3643875 0.8402791
## 4 High Normal Sterilized 9 36.51148 1.3669065 0.4556355 1.0506974
## 5 Low Drought Living 9 12.79916 0.7621779 0.2540593 0.5858618
## 6 Low Drought Sterilized 9 17.36461 0.9035038 0.3011679 0.6944945
## 7 Low Normal Living 9 15.59275 0.8514438 0.2838146 0.6544777
## 8 Low Normal Sterilized 9 28.98129 1.1896826 0.3965609 0.91447105.3.2 SoilxWater(*)
pdf('./fig/NC SoilxWater.pdf', height = 15/2.54, width = 15/2.54)
pd <- position_dodge(.3)
Fig_B_01 <- summarySE(Fig_B, measurevar="native_pred", groupvars=c("Soil","Water"))
#Fig_B_01
Fig_B_01$Soil <- ordered(Fig_B_01$Soil, levels = c('Living', 'Sterilized'))
Fig_B_01$Water <- ordered(Fig_B_01$Water, levels = c('Normal', 'Drought'))
p <- ggplot(Fig_B_01, aes(x=Water, y=native_pred, color=Soil)) +
#scale_y_continuous(limits = c(0,4), breaks = c(0,2,4)) +
geom_errorbar(aes(ymin=native_pred-se, ymax=native_pred+se),
size=2, width=0,position = pd) +
geom_line(colour="darkgrey",linetype="dashed", size=1,
aes(group=Soil), position = pd)+
geom_point(size = 8,position = pd) +
#facet_grid( ~ Competition) +
labs(title="",x="", y = " Native community
aboveground biomass")
#remove the background
p <- p + theme_bw()
p <- p + scale_color_manual(breaks = c("Living", "Sterilized"),
values=c("#bcd4e8ff", "#64879aff"))
#change the direction of x-aix labels
pncsw <- p + theme+
theme(#legend.position=c(0.96, 0.5),
legend.background=element_rect(fill=NA),
legend.title=element_text(size=17),
legend.text=element_text(size=17),
legend.spacing.y = unit(0.2, "cm")
)+ theme(legend.title=element_blank())
pncsw <- pncsw + theme(legend.position = c( .85 , .85 ))+coord_cartesian(clip = 'off',ylim = c(10,40))+
theme(axis.line.x=element_line(linetype=1,color="black",size=0),
axis.line.y=element_line(linetype=1,color="black",size=0),
axis.ticks = element_line(colour = "black",lineend = "square",size = 0.5))
pncsw
dev.off()## png
## 2
## Soil Water N native_pred sd se ci
## 1 Living Drought 18 14.13804 1.593812 0.3756652 0.7925842
## 2 Living Normal 18 18.25325 2.897960 0.6830556 1.4411214
## 3 Sterilized Drought 18 18.70111 1.778178 0.4191206 0.8842671
## 4 Sterilized Normal 18 32.74639 4.068801 0.9590255 2.02336705.3.3 SoilxDiversity(*)
pdf('./fig/NC SoilxDiversity.pdf', height = 15/2.54, width = 15/2.54)
pd <- position_dodge(.4)
Fig_B_01 <- summarySE(Fig_B, measurevar="native_pred", groupvars=c("Soil","Diversity"))
#Fig_B_01
Fig_B_01$Soil <- ordered(Fig_B_01$Soil, levels = c('Living', 'Sterilized'))
Fig_B_01$Diversity <- ordered(Fig_B_01$Diversity, levels = c('Low', 'High'))
p <- ggplot(Fig_B_01, aes(x=Diversity, y=native_pred, color=Soil)) +
geom_errorbar(aes(ymin=native_pred-se, ymax=native_pred+se),
size=2, width=0,position = pd) +
geom_line(colour="darkgrey",linetype="dashed", size=1,
aes(group=Soil), position = pd)+
geom_point(size = 8,position = pd) +
#facet_grid( ~ Competition) +
labs(title="",x="", y = " Native community
aboveground biomass")
#remove the background
p <- p + theme_bw()
p <- p + scale_color_manual(breaks = c("Living", "Sterilized"),
values=c("#bcd4e8ff", "#64879aff"))
#change the direction of x-aix labels
pncsd <- p + theme+
theme(#legend.position=c(0.96, 0.5),
legend.background=element_rect(fill=NA),
legend.title=element_text(size=17),
legend.text=element_text(size=17),
legend.spacing.y = unit(0.2, "cm")
)+ theme(legend.title=element_blank())
pncsd <- pncsd + theme(legend.position = c( .85 , .85 ))+coord_cartesian(clip = 'off',ylim = c(10,35))+
theme(axis.line.x=element_line(linetype=1,color="black",size=0),
axis.line.y=element_line(linetype=1,color="black",size=0),
axis.ticks = element_line(colour = "black",lineend = "square",size = 0.5))+
scale_y_continuous(breaks=seq(10,35,5))
pncsd
dev.off()## png
## 2
## Soil Diversity N native_pred sd se ci
## 1 Living High 18 18.19534 2.958998 0.6974426 1.4714753
## 2 Living Low 18 14.19595 1.637170 0.3858847 0.8141455
## 3 Sterilized High 18 28.27455 8.579263 2.0221517 4.2663673
## 4 Sterilized Low 18 23.17295 6.063955 1.4292878 3.01553375.3.4 Soil(*)
## Soil N native_pred sd se ci
## 1 Living 36 16.19565 3.109284 0.518214 1.052030
## 2 Sterilized 36 25.72375 7.765523 1.294254 2.6274755.3.5 Water(*)
## Water N native_pred sd se ci
## 1 Drought 36 16.41958 2.850220 0.4750366 0.9643756
## 2 Normal 36 25.49982 8.132237 1.3553728 2.75155316 Aboveground biomass proportion of alien target species
data_biomass <- data_biomass %>%
mutate(target.pro=(Target_biomass)/(Target_biomass+Native_Biomass)) %>%
droplevels()6.1 model selection
6.1.1 data transform selection
u<-gl(1,1,length(data_biomass$target.pro))
# no transform
m_Above_1<- lme(target.pro~ Soil*Water*Diversity,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
control=list(maxIter=200,msMaxIter=200,niterEM=200),
data=data_biomass,method="REML")
# log
m_Above_2<- lme(log(target.pro)~ Soil*Water*Diversity,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
control=list(maxIter=200,msMaxIter=200,niterEM=200),
data=data_biomass,method="REML")
# sqrt
m_Above_3<- lme(sqrt(target.pro)~ Soil*Water*Diversity,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
control=list(maxIter=200,msMaxIter=200,niterEM=200),
data=data_biomass,method="REML")
m_Above_4<- lme(cbrt(target.pro)~ Soil*Water*Diversity,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
control=list(maxIter=200,msMaxIter=200,niterEM=200),
data=data_biomass,method="REML")
# ^(1/3)
m_Above_5<- lme((target.pro)^(1/3)~ Soil*Water*Diversity,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
control=list(maxIter=200,msMaxIter=200,niterEM=200),
data=data_biomass,method="REML")
# logit
m_Above_6<- lme(logit(target.pro)~ Soil*Water*Diversity,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
control=list(maxIter=200,msMaxIter=200,niterEM=200),
data=data_biomass,method="REML")
opm1a <- par(mfrow = c(3,3), mar = c(4,4,2,2))
plot.resid(m_Above_1)
plot.resid(m_Above_2)# log
plot.resid(m_Above_3)# ^(sqrt)
## df AIC
## m_Above_1 12 -658.1533
## m_Above_2 12 1521.0183
## m_Above_3 12 -662.0432
## m_Above_4 12 -762.2615
## m_Above_5 12 -762.2615
## m_Above_6 12 1687.5812logit
opm1a <- par(mfrow = c(2,3), mar = c(4,4,2,2))
plot.resid(m_Above_2)# log
plot.resid(m_Above_6)# logit
## df AIC
## m_Above_2 12 1521.018
## m_Above_6 12 1687.581log transform is better(Residuals better)
6.1.2 variance structure selection
# Target_species
m_Above_2vf <- lme(log(target.pro)~ Soil*Water*Diversity,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
weights=varComb(varIdent(form= ~1|Target_Species)), control=list(maxIter=200,msMaxIter=200,niterEM=200),
data=data_biomass,method="REML")
opm1a <- par(mfrow = c(2,3), mar = c(4,4,2,2))
plot.resid(m_Above_2)
plot.resid(m_Above_2vf)
## df AIC
## m_Above_2 12 1521.018
## m_Above_2vf 20 1493.813m_Above_2vf is better(residuals is similiar,AIC is smaller)
## Warning in model.matrix.default(formula(terms(mf)), mf, contrasts.arg =
## object$contrasts): variable 'Target_Family' is absent, its contrast will be
## ignored## Warning in model.matrix.default(formula(terms(mf)), mf, contrasts.arg =
## object$contrasts): variable 'Target_Species' is absent, its contrast will be
## ignored## Warning in model.matrix.default(formula(terms(mf)), mf, contrasts.arg =
## object$contrasts): variable 'Native_Communities' is absent, its contrast will
## be ignored## R2m R2c
## [1,] 0.05838564 0.6556214## Linear mixed-effects model fit by REML
## Data: data_biomass
## AIC BIC logLik
## 1493.813 1580.514 -726.9066
##
## Random effects:
## Composite Structure: Blocked
##
## Block 1: Target_FamilyAmaranthaceae, Target_FamilyAsteraceae, Target_FamilyFabaceae, Target_FamilyOxalidaceae, Target_FamilyPlantaginaceae, Target_FamilyPoaceae, Target_FamilySolanaceae
## Formula: ~Target_Family - 1 | u
## Structure: Multiple of an Identity
## Target_FamilyAmaranthaceae Target_FamilyAsteraceae
## StdDev: 0.0009518097 0.0009518097
## Target_FamilyFabaceae Target_FamilyOxalidaceae
## StdDev: 0.0009518097 0.0009518097
## Target_FamilyPlantaginaceae Target_FamilyPoaceae
## StdDev: 0.0009518097 0.0009518097
## Target_FamilySolanaceae
## StdDev: 0.0009518097
##
## Block 2: Target_SpeciesAmaranthus_retroflexus, Target_SpeciesBidens_frondosa, Target_SpeciesBidens_pilosa, Target_SpeciesDatura_stramonium, Target_SpeciesLolium_perenne, Target_SpeciesMedicago_sativa, Target_SpeciesOxalis_corymbosa, Target_SpeciesPlantago_virginica, Target_SpeciesSesbania_cannabina
## Formula: ~Target_Species - 1 | u
## Structure: Multiple of an Identity
## Target_SpeciesAmaranthus_retroflexus Target_SpeciesBidens_frondosa
## StdDev: 1.195501 1.195501
## Target_SpeciesBidens_pilosa Target_SpeciesDatura_stramonium
## StdDev: 1.195501 1.195501
## Target_SpeciesLolium_perenne Target_SpeciesMedicago_sativa
## StdDev: 1.195501 1.195501
## Target_SpeciesOxalis_corymbosa Target_SpeciesPlantago_virginica
## StdDev: 1.195501 1.195501
## Target_SpeciesSesbania_cannabina
## StdDev: 1.195501
##
## Block 3: Native_CommunitiesH_NC10, Native_CommunitiesH_NC11, Native_CommunitiesH_NC12, Native_CommunitiesH_NC13, Native_CommunitiesH_NC14, Native_CommunitiesH_NC15, Native_CommunitiesH_NC16, Native_CommunitiesH_NC9, Native_CommunitiesL_NC1, Native_CommunitiesL_NC2, Native_CommunitiesL_NC3, Native_CommunitiesL_NC4, Native_CommunitiesL_NC5, Native_CommunitiesL_NC6, Native_CommunitiesL_NC7, Native_CommunitiesL_NC8
## Formula: ~Native_Communities - 1 | u
## Structure: Multiple of an Identity
## Native_CommunitiesH_NC10 Native_CommunitiesH_NC11
## StdDev: 0.7397163 0.7397163
## Native_CommunitiesH_NC12 Native_CommunitiesH_NC13
## StdDev: 0.7397163 0.7397163
## Native_CommunitiesH_NC14 Native_CommunitiesH_NC15
## StdDev: 0.7397163 0.7397163
## Native_CommunitiesH_NC16 Native_CommunitiesH_NC9
## StdDev: 0.7397163 0.7397163
## Native_CommunitiesL_NC1 Native_CommunitiesL_NC2 Native_CommunitiesL_NC3
## StdDev: 0.7397163 0.7397163 0.7397163
## Native_CommunitiesL_NC4 Native_CommunitiesL_NC5 Native_CommunitiesL_NC6
## StdDev: 0.7397163 0.7397163 0.7397163
## Native_CommunitiesL_NC7 Native_CommunitiesL_NC8 Residual
## StdDev: 0.7397163 0.7397163 1.067537
##
## Combination of variance functions:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Target_Species
## Parameter estimates:
## Sesbania_cannabina Bidens_frondosa Amaranthus_retroflexus
## 1.0000000 0.5554415 0.8449116
## Bidens_pilosa Oxalis_corymbosa Medicago_sativa
## 0.6158151 0.6746881 0.9071252
## Datura_stramonium Lolium_perenne Plantago_virginica
## 0.7655588 0.8580667 0.5450924
## Fixed effects: log(target.pro) ~ Soil * Water * Diversity
## Value Std.Error DF t-value
## (Intercept) -3.271339 0.4850196 564 -6.744756
## SoilSterilized -0.622242 0.1270064 564 -4.899298
## WaterNormal -0.260743 0.1260335 564 -2.068838
## DiversityLow 0.713711 0.3905349 564 1.827523
## SoilSterilized:WaterNormal 0.376081 0.1789103 564 2.102063
## SoilSterilized:DiversityLow -0.170235 0.1784743 564 -0.953835
## WaterNormal:DiversityLow 0.113602 0.1777833 564 0.638993
## SoilSterilized:WaterNormal:DiversityLow -0.183333 0.2519004 564 -0.727798
## p-value
## (Intercept) 0.0000
## SoilSterilized 0.0000
## WaterNormal 0.0390
## DiversityLow 0.0681
## SoilSterilized:WaterNormal 0.0360
## SoilSterilized:DiversityLow 0.3406
## WaterNormal:DiversityLow 0.5231
## SoilSterilized:WaterNormal:DiversityLow 0.4670
## Correlation:
## (Intr) SlStrl WtrNrm DvrstL SlS:WN
## SoilSterilized -0.128
## WaterNormal -0.128 0.491
## DiversityLow -0.403 0.158 0.160
## SoilSterilized:WaterNormal 0.091 -0.710 -0.704 -0.113
## SoilSterilized:DiversityLow 0.091 -0.712 -0.349 -0.226 0.505
## WaterNormal:DiversityLow 0.091 -0.348 -0.709 -0.226 0.499
## SoilSterilized:WaterNormal:DiversityLow -0.064 0.504 0.500 0.160 -0.710
## SlS:DL WtN:DL
## SoilSterilized
## WaterNormal
## DiversityLow
## SoilSterilized:WaterNormal
## SoilSterilized:DiversityLow
## WaterNormal:DiversityLow 0.495
## SoilSterilized:WaterNormal:DiversityLow -0.708 -0.706
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.0846888 -0.5982175 0.1152389 0.6619631 2.4445957
##
## Number of Observations: 572
## Number of Groups: 16.2 likelihood ratio test
Three_Ref <- lme(logit(target.pro)~ Soil*Water*Diversity,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
weights=varComb(varIdent(form= ~1|Target_Species)), control=list(maxIter=200,msMaxIter=200,niterEM=200),
data=data_biomass,method="ML")
Three <- drop1(Three_Ref ,test="Chisq")[-1,]
# Two-way interaction
Two_Ref <- lme(logit(target.pro)~ Soil+Water+Diversity+Soil:Water + Soil:Diversity + Water:Diversity,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
weights=varComb(varIdent(form= ~1|Target_Species)), control=list(maxIter=600,msMaxIter=600,niterEM=600),
data=data_biomass,method="ML")
Two <- drop1(Two_Ref ,test="Chisq")[-1,]
#Main effect
Main_Ref <- lme(logit(target.pro)~ Soil+Water+Diversity,
random=list(u=pdBlocked(list(pdIdent(form=~Target_Family-1),
pdIdent(form=~Target_Species-1),
pdIdent(form=~Native_Communities-1)
))),
weights=varComb(varIdent(form= ~1|Target_Species)), control=list(maxIter=600,msMaxIter=600,niterEM=600),
data=data_biomass,method="ML")
Main <- drop1(Main_Ref,test="Chisq")[-1,]6.3 LRT table
table<-rbind(Main,Two,Three)
table$" "=ifelse(table$`Pr(>Chi)`<0.05,"*",
ifelse(table$`Pr(>Chi)`<0.1,"..",""))#add significant levels
kable(table,digits=4)| Df | AIC | LRT | Pr(>Chi) | ||
|---|---|---|---|---|---|
| Soil | 1 | 1702.278 | 57.3473 | 0.0000 | * |
| Water | 1 | 1645.849 | 0.9180 | 0.3380 | |
| Diversity | 1 | 1647.844 | 2.9133 | 0.0879 | .. |
| Soil:Water | 1 | 1644.226 | 3.3919 | 0.0655 | .. |
| Soil:Diversity | 1 | 1647.529 | 6.6948 | 0.0097 | * |
| Water:Diversity | 1 | 1640.858 | 0.0234 | 0.8783 | |
| Soil:Water:Diversity | 1 | 1642.834 | 0.2420 | 0.6228 |
6.4 Fig_C_pred
data_biomass$target.pro_pred <- (invlogit(predict(m_Above_2vf)))
Fig_C <- summarySE(data_biomass, measurevar="target.pro_pred", groupvars=c("Target_Species","Diversity","Water","Soil"))6.4.1 SoilxWaterxDiversity(ns)
Fig_C_01 <- summarySE(Fig_C, measurevar="target.pro_pred", groupvars=c("Soil",
"Water","Diversity"))
Fig_C_01$Soil <- ordered(Fig_C_01$Soil, levels = c('Living', 'Sterilized'))
Fig_C_01$Diversity <- ordered(Fig_C_01$Diversity, levels = c('Low', 'High'))
Fig_C_01$Water <- ordered(Fig_C_01$Water, levels = c('Normal', 'Drought'))
pd <- position_dodge(.3)
p <- ggplot(Fig_C_01, aes(x=Water, y=target.pro_pred, color=Soil)) +
ylim(0,.2) +
geom_errorbar(aes(ymin=target.pro_pred-se, ymax=target.pro_pred+se),
size=2, width=0,position = pd) +
geom_line(colour="darkgrey",linetype="dashed", size=1,
aes(group=Soil), position = pd)+
geom_point(size = 8,position = pd) +
facet_grid( ~Diversity ) +
labs(title="",x="", y = " Alien target species
aboveground biomass proportion")
#remove the background
p <- p + theme_bw()
p <- p + scale_color_manual(breaks = c("Living", "Sterilized"),
values=c("#bcd4e8ff", "#64879aff"))
#change the direction of x-aix labels
ptpo1 <- p + theme+
theme(#legend.position=c(0.96, 0.5),
legend.background=element_rect(fill=NA),
legend.title=element_text(size=17),
legend.text=element_text(size=17),
legend.spacing.y = unit(0.2, "cm")
)+ theme(legend.title=element_blank())
ptpo1 <- ptpo1 + theme(legend.position = c( .85 , .85 ))+coord_cartesian(clip = 'off',ylim = c(0,.2))
ptpo1
## Soil Water Diversity N target.pro_pred sd se
## 1 Living Drought High 9 0.06024069 0.04876416 0.01625472
## 2 Living Drought Low 9 0.12644049 0.09106166 0.03035389
## 3 Living Normal High 9 0.04806338 0.03951444 0.01317148
## 4 Living Normal Low 9 0.11325972 0.08301178 0.02767059
## 5 Sterilized Drought High 9 0.03388085 0.02878152 0.00959384
## 6 Sterilized Drought Low 9 0.06750379 0.05278942 0.01759647
## 7 Sterilized Normal High 9 0.03799939 0.03160735 0.01053578
## 8 Sterilized Normal Low 9 0.07013766 0.05463438 0.01821146
## ci
## 1 0.03748345
## 2 0.06999619
## 3 0.03037349
## 4 0.06380851
## 5 0.02212343
## 6 0.04057754
## 7 0.02429556
## 8 0.041995706.4.2 SoilxDiversity(*)
pdf('./fig/TP SoilxDiversity.pdf', height = 15/2.54, width = 25/2.54)
pd <- position_dodge(width = .4)
Fig_C_01 <- summarySE(Fig_C, measurevar="target.pro_pred", groupvars=c("Soil","Diversity"))
Fig_C_01$Soil <- ordered(Fig_C_01$Soil, levels = c('Living', 'Sterilized'))
Fig_C_01$Diversity <- ordered(Fig_C_01$Diversity, levels = c('Low', 'High'))
p <- ggplot(Fig_C_01, aes(x=Diversity, y=target.pro_pred, color=Soil)) +
#scale_y_continuous(limits = c(0,4), breaks = c(0,2,4)) +
geom_errorbar(aes(ymin=target.pro_pred-se, ymax=target.pro_pred+se),
size=2, width=0,position = pd) +
geom_line(colour="darkgrey",linetype="dashed", size=1,
aes(group=Soil), position = pd)+
geom_point(size = 8,position = pd) +
#facet_grid( ~ Competition) +
labs(title="",x="", y = " Alien target species
aboveground biomass proportion")
#remove the background
p <- p + theme_bw()
p <- p + scale_color_manual(breaks = c("Living", "Sterilized"),
values=c("#bcd4e8ff", "#64879aff"))
#change the direction of x-aix labels
ptposd <- p + theme+
theme(#legend.position=c(0.96, 0.5),
legend.background=element_rect(fill=NA),
legend.title=element_text(size=17),
legend.text=element_text(size=17),
legend.spacing.y = unit(0.2, "cm")
)+ theme(legend.title=element_blank()) + theme(legend.position = c( .85 , .9 ))+coord_cartesian(clip = 'off',ylim = c(0.02,0.16))+
# annotate("text", x = 0.5, y = 35, label = expression(S ~ '*' ),colour = 'black',family='serif', hjust = 0, size = 7.8)+
# annotate("text", x = 1.5, y = 35, label = expression(D ~ 'ns'),colour = 'black',family='serif', hjust = 0, size = 7.8)+
# annotate("text", x = 0.5, y = 32, label = expression(S:D ~ '..'), colour = 'black',family='serif',hjust = 0, size = 7.8)+
scale_y_continuous(breaks=seq(0.02,0.16,0.07))
ptposd
dev.off()## png
## 2
## Soil Diversity N target.pro_pred sd se ci
## 1 Living High 18 0.05415204 0.04350927 0.01025523 0.02163665
## 2 Living Low 18 0.11985011 0.08479986 0.01998752 0.04216998
## 3 Sterilized High 18 0.03594012 0.02940141 0.00692998 0.01462098
## 4 Sterilized Low 18 0.06882072 0.05213349 0.01228798 0.025925386.4.3 SoilxWater(marginal significant)
pdf('./fig/TP SoilxWater.pdf', height = 15/2.54, width = 15/2.54)
pd <- position_dodge(.3)
Fig_C_01 <- summarySE(Fig_C, measurevar="target.pro_pred", groupvars=c("Soil","Water"))
#Fig_C_01
Fig_C_01$Soil <- ordered(Fig_C_01$Soil, levels = c('Living', 'Sterilized'))
Fig_C_01$Water <- ordered(Fig_C_01$Water, levels = c('Normal', 'Drought'))
p <- ggplot(Fig_C_01, aes(x=Water, y=target.pro_pred, color=Soil)) +
#scale_y_continuous(limits = c(0,4), breaks = c(0,2,4)) +
geom_errorbar(aes(ymin=target.pro_pred-se, ymax=target.pro_pred+se),
size=2, width=0,position = pd) +
geom_line(colour="darkgrey",linetype="dashed", size=1,
aes(group=Soil), position = pd)+
geom_point(size = 8,position = pd) +
#facet_grid( ~ Competition) +
labs(title="",x="", y = " Alien target species
aboveground biomass proportion")
#remove the background
p <- p + theme_bw()
p <- p + scale_color_manual(breaks = c("Living", "Sterilized"),
values=c("#bcd4e8ff", "#64879aff"))
#change the direction of x-aix labels
p <- p + theme+
theme(#legend.position=c(0.96, 0.5),
legend.background=element_rect(fill=NA),
legend.title=element_text(size=17),
legend.text=element_text(size=17),
legend.spacing.y = unit(0.2, "cm")
)+ theme(legend.title=element_blank())
ptposw <- p + theme(legend.position = c( .85 , .85 ))+coord_cartesian(clip = 'off',ylim = c(0.03,.13))
ptposw
dev.off()## png
## 2
## Soil Water N target.pro_pred sd se ci
## 1 Living Drought 18 0.09334059 0.07862129 0.01853122 0.03909745
## 2 Living Normal 18 0.08066155 0.07143335 0.01683700 0.03552297
## 3 Sterilized Drought 18 0.05069232 0.04472665 0.01054217 0.02224204
## 4 Sterilized Normal 18 0.05406853 0.04634869 0.01092449 0.023048666.4.4 Soil(*)
## Soil N target.pro_pred sd se ci
## 1 Living 36 0.08700107 0.07431119 0.012385198 0.02514329
## 2 Sterilized 36 0.05238042 0.04492217 0.007487029 0.015199486.4.5 Water(ns)
## Water N target.pro_pred sd se ci
## 1 Drought 36 0.07201646 0.06664622 0.01110770 0.02254984
## 2 Normal 36 0.06736504 0.06085825 0.01014304 0.020591477 germination
## Factor levels
data_germination$Diversity<-factor(data_germination$Diversity, levels=c("High","Low"))
data_germination$Water<-factor(data_germination$Water, c("Normal","Drought"))
data_germination$Soil<-factor(data_germination$Soil, levels=c("Living","Sterilized"))
## Apply sum contrust
contrasts(data_germination$Diversity) <- contr.sum(2)
contrasts(data_germination$Water) <- contr.sum(2)
contrasts(data_germination$Soil) <- contr.sum(2)7.1 Brms model
# prior = c(set_prior ("normal(0,1)",class="b"),
# set_prior ("normal(0,10)",class="Intercept"))
# modelgermination<-brm(HarvestNum~Diversity*Water*Soil+
# (1|TargetFamily/TargetSpecies)+(1|NativeCommunities),
# prior = prior,
# control = list(adapt_delta = 0.9999, max_treedepth = 20),
# data=data_germination,family ="negbinomial2",seed = 123,
# iter = 3000, cores = 4)
# saveRDS(modelgermination, file = "modelgermination.rds")
modelgermination<-readRDS(file = "modelgermination.rds")## Family: negbinomial2
## Links: mu = log; sigma = identity
## Formula: HarvestNum ~ Diversity * Water * Soil + (1 | TargetFamily/TargetSpecies) + (1 | NativeCommunities)
## Data: data_germination (Number of observations: 564)
## Draws: 4 chains, each with iter = 3000; warmup = 1500; thin = 1;
## total post-warmup draws = 6000
##
## Multilevel Hyperparameters:
## ~NativeCommunities (Number of levels: 16)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.10 0.06 0.01 0.23 1.00 1999 2411
##
## ~TargetFamily (Number of levels: 7)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.85 0.62 0.04 2.30 1.00 1405 2523
##
## ~TargetFamily:TargetSpecies (Number of levels: 9)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.07 0.44 0.45 2.12 1.00 2097 3258
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 1.53 0.57 0.41 2.68 1.00 2557
## Diversity1 -0.07 0.05 -0.16 0.03 1.00 5294
## Water1 0.36 0.04 0.28 0.43 1.00 9543
## Soil1 0.16 0.04 0.08 0.23 1.00 9734
## Diversity1:Water1 0.05 0.04 -0.03 0.12 1.00 8426
## Diversity1:Soil1 0.03 0.04 -0.04 0.11 1.00 9829
## Water1:Soil1 -0.21 0.04 -0.29 -0.14 1.00 9043
## Diversity1:Water1:Soil1 -0.08 0.04 -0.15 -0.00 1.00 9506
## Tail_ESS
## Intercept 3030
## Diversity1 4213
## Water1 4362
## Soil1 3914
## Diversity1:Water1 3899
## Diversity1:Soil1 4893
## Water1:Soil1 4748
## Diversity1:Water1:Soil1 4398
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.53 0.05 0.43 0.63 1.00 7735 4496
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
7.2 Table
| Estimate | SE | L95%CL | U95%CL | L90%CL | U90%CL | |
|---|---|---|---|---|---|---|
| Intercept (Grand Mean) | 1.533* | 0.5686 | 0.4127 | 2.6838 | 0.6366 | 2.4512 |
| Diversity | -0.0672 | 0.0482 | -0.1647 | 0.0285 | -0.1466 | 0.0128 |
| Water | 0.359* | 0.0383 | 0.2836 | 0.4346 | 0.2970 | 0.4232 |
| Soil | 0.157* | 0.0379 | 0.0823 | 0.2316 | 0.0951 | 0.2182 |
| Diversity:Water | 0.0487 | 0.0390 | -0.0298 | 0.1240 | -0.0162 | 0.1121 |
| Diversity:Soil | 0.0309 | 0.0384 | -0.0439 | 0.1050 | -0.0325 | 0.0951 |
| Water:Soil | -0.213* | 0.0382 | -0.2886 | -0.1380 | -0.2762 | -0.1496 |
| Diversity:Water:Soil | -0.076* | 0.0374 | -0.1496 | -0.0029 | -0.1377 | -0.0145 |
| Note: | ||||||
L95%CL * U95%CL >
0,significant,L90%CL * U90%CL >
0,marginally significant
|
7.3 Figures
7.4 SoilxWaterxDiversity(*)
pdf('./fig/seeldings SoilxWaterxDiversity.pdf', height = 15/2.54, width = 15/2.54)
pd <- position_dodge(.3)
d_AboveT_Fig1_mean <- summarySE(d_AboveT_Fig1, measurevar="HarvestNum", groupvars=c("Soil",
"Water","Diversity"))
d_AboveT_Fig1_mean$Soil <- ordered(d_AboveT_Fig1_mean$Soil, levels = c('Living', 'Sterilized'))
d_AboveT_Fig1_mean$Diversity <- ordered(d_AboveT_Fig1_mean$Diversity, levels = c('Low', 'High'))
d_AboveT_Fig1_mean$Water <- ordered(d_AboveT_Fig1_mean$Water, levels = c('Normal', 'Drought'))
p <- ggplot(d_AboveT_Fig1_mean, aes(x=Water, y=HarvestNum, color=Soil)) +
geom_errorbar(aes(ymin=HarvestNum-se, ymax=HarvestNum+se),
size=2, width=0,position = pd) +
geom_line(colour="darkgrey",linetype="dashed", size=1,
aes(group=Soil), position = pd)+
geom_point(size = 8,position = pd) +
facet_grid( ~Diversity ) +
labs(title="",x="", y = " survival number of
alien target seedlings")
#remove the background
p <- p + theme_bw()
p <- p + scale_color_manual(breaks = c("Living", "Sterilized"),
values=c("#bcd4e8ff", "#64879aff"))
#change the direction of x-aix labels
psuswd <- p + theme+
theme(#legend.position=c(0.96, 0.5),
legend.background=element_rect(fill=NA),
legend.title=element_text(size=17),
legend.text=element_text(size=17),
legend.spacing.y = unit(0.2, "cm")
)+ theme(legend.title=element_blank())
psuswd <- psuswd + theme(legend.position = c( .85 , .85 ))+
theme(axis.line.x=element_line(linetype=1,color="black",size=0),
axis.line.y=element_line(linetype=1,color="black",size=0),
axis.ticks = element_line(colour = "black",lineend = "square",size = 0.5))+coord_cartesian(clip = 'off',ylim = c(0,15))+
scale_y_continuous(breaks=seq(0,15,5))
psuswd
dev.off()## png
## 2
## Soil Water Diversity N HarvestNum sd se ci
## 1 Living Normal High 9 8.736111 7.472461 2.4908203 5.743842
## 2 Living Normal Low 9 8.865079 6.094475 2.0314917 4.684628
## 3 Living Drought High 9 6.180556 4.005422 1.3351406 3.078840
## 4 Living Drought Low 9 6.746032 4.980925 1.6603083 3.828678
## 5 Sterilized Normal High 9 10.857143 8.335131 2.7783770 6.406949
## 6 Sterilized Normal Low 9 10.210317 7.487206 2.4957354 5.755176
## 7 Sterilized Drought High 9 2.571429 2.604435 0.8681451 2.001946
## 8 Sterilized Drought Low 9 4.472222 4.465126 1.4883753 3.4322007.5 SoilxWater(*)
pdf('./fig/seeldings SoilxWater.pdf', height = 15/2.54, width = 15/2.54)
pd <- position_dodge(.3)
d_AboveT_Fig1_mean <- summarySE(d_AboveT_Fig1, measurevar="HarvestNum", groupvars=c("Soil",
"Water"))
d_AboveT_Fig1_mean$Soil <- ordered(d_AboveT_Fig1_mean$Soil, levels = c('Living', 'Sterilized'))
d_AboveT_Fig1_mean$Water <- ordered(d_AboveT_Fig1_mean$Water, levels = c('Normal', 'Drought'))
p <- ggplot(d_AboveT_Fig1_mean, aes(x=Water, y=HarvestNum, color=Soil)) +
geom_errorbar(aes(ymin=HarvestNum-se, ymax=HarvestNum+se),
size=2, width=0,position = pd) +
geom_line(colour="darkgrey",linetype="dashed", size=1,
aes(group=Soil), position = pd)+
geom_point(size = 8,position = pd) +
#facet_grid( ~Diversity ) +
labs(title="",x="", y = " survival number of
alien target seedlings")
#remove the background
p <- p + theme_bw()
p <- p + scale_color_manual(breaks = c("Living", "Sterilized"),
values=c("#bcd4e8ff", "#64879aff"))
#change the direction of x-aix labels
psuswd <- p + theme+
theme(#legend.position=c(0.96, 0.5),
legend.background=element_rect(fill=NA),
legend.title=element_text(size=17),
legend.text=element_text(size=17),
legend.spacing.y = unit(0.2, "cm")
)+ theme(legend.title=element_blank())
psuswd <- psuswd + theme(legend.position = c( .85 , .85 ))+
theme(axis.line.x=element_line(linetype=1,color="black",size=0),
axis.line.y=element_line(linetype=1,color="black",size=0),
axis.ticks = element_line(colour = "black",lineend = "square",size = 0.5))+coord_cartesian(clip = 'off',ylim = c(0,15))+
scale_y_continuous(breaks=seq(0,15,5))
psuswd
dev.off()## png
## 2
## Soil Water N HarvestNum sd se ci
## 1 Living Normal 18 8.800595 6.615122 1.5591992 3.289623
## 2 Living Drought 18 6.463294 4.394268 1.0357389 2.185218
## 3 Sterilized Normal 18 10.533730 7.693175 1.8132988 3.825726
## 4 Sterilized Drought 18 3.521825 3.678411 0.8670098 1.8292317.6 Soil(*)
d_AboveT_Fig1_mean <- summarySE(d_AboveT_Fig1, measurevar="HarvestNum", groupvars=c("Soil"))
d_AboveT_Fig1_mean## Soil N HarvestNum sd se ci
## 1 Living 36 7.631944 5.660255 0.9433759 1.915155
## 2 Sterilized 36 7.027778 6.925455 1.1542425 2.3432378 Picture in article
pdf('./fig/Figure.main2factor.pdf', height = 35/2.54, width = 45/2.54)
p1<-ptsw
p2<-pncsw
p3<-ptposw
p4<-ptsd
p5<-pncsd
p6<-ptposd
ggarrange(p1, p2,p3,p4,p5,p6,ncol = 3, nrow = 2,
labels = c("A","B","C","D","E","F"),
font.label = list(size = 22, face = "bold"))
dev.off()## png
## 2#pdf('./fig/main.pdf', height = 15/2.54, width = 65/2.54)
#pdf('./fig/main1.pdf', height = 15/2.54, width = 55/2.54)
pdf('./fig/Figure.all.pdf', height = 15/2.54, width = 55/2.54)
p1<-ptswd
p2<-pncswd
p3<-ptpo1
ggarrange(p1, p2, p3, ncol = 3, nrow = 1,
widths = c(0.6, 0.6, 0.6),
labels = c("A","B","C"),
font.label = list(size = 22, face = "bold"))
dev.off()## png
## 29 Session information
## R version 4.4.2 (2024-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=Chinese (Simplified)_China.utf8
## [2] LC_CTYPE=Chinese (Simplified)_China.utf8
## [3] LC_MONETARY=Chinese (Simplified)_China.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=Chinese (Simplified)_China.utf8
##
## time zone: Asia/Shanghai
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] LaplacesDemon_16.1.6 kableExtra_1.4.0 brms_2.21.0
## [4] Rcpp_1.0.12 gcookbook_2.0 magrittr_2.0.3
## [7] reshape2_1.4.4 car_3.1-2 carData_3.0-5
## [10] performance_0.12.1 funrar_1.5.0 phytools_2.3-0
## [13] maps_3.4.2 ape_5.8 htmltools_0.5.8.1
## [16] ggpubr_0.6.0 emmeans_1.10.3 multcomp_1.4-25
## [19] TH.data_1.1-2 MASS_7.3-61 survival_3.7-0
## [22] mvtnorm_1.2-5 Rmisc_1.5.1 plyr_1.8.9
## [25] lattice_0.22-6 MuMIn_1.48.4 knitr_1.48
## [28] lme4_1.1-35.3 Matrix_1.7-1 nlme_3.1-166
## [31] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1
## [34] dplyr_1.1.4 purrr_1.0.2 readr_2.1.5
## [37] tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.1
## [40] tidyverse_2.0.0 xlsx_0.6.5 rstudioapi_0.16.0
##
## loaded via a namespace (and not attached):
## [1] tensorA_0.36.2.1 jsonlite_1.8.8 estimability_1.5.1
## [4] farver_2.1.2 nloptr_2.0.3 rmarkdown_2.27
## [7] vctrs_0.6.5 minqa_1.2.7 rstatix_0.7.2
## [10] distributional_0.4.0 curl_5.2.1 DEoptim_2.2-8
## [13] broom_1.0.6 sass_0.4.9 StanHeaders_2.32.10
## [16] bslib_0.7.0 sandwich_3.1-0 zoo_1.8-12
## [19] cachem_1.1.0 igraph_2.0.3 lifecycle_1.0.4
## [22] iterators_1.0.14 pkgconfig_2.0.3 R6_2.5.1
## [25] fastmap_1.2.0 digest_0.6.37 numDeriv_2016.8-1.1
## [28] colorspace_2.1-0 labeling_0.4.3 clusterGeneration_1.3.8
## [31] fansi_1.0.6 timechange_0.3.0 abind_1.4-5
## [34] compiler_4.4.2 withr_3.0.0 doParallel_1.0.17
## [37] backports_1.5.0 inline_0.3.19 optimParallel_1.0-2
## [40] highr_0.11 QuickJSR_1.3.1 pkgbuild_1.4.4
## [43] ggsignif_0.6.4 scatterplot3d_0.3-44 loo_2.8.0
## [46] tools_4.4.2 glue_1.7.0 quadprog_1.5-8
## [49] grid_4.4.2 checkmate_2.3.1 generics_0.1.3
## [52] gtable_0.3.5 tzdb_0.4.0 hms_1.1.3
## [55] xml2_1.3.6 utf8_1.2.4 foreach_1.5.2
## [58] pillar_1.9.0 posterior_1.6.0 rJava_1.0-11
## [61] splines_4.4.2 tidyselect_1.2.1 gridExtra_2.3
## [64] V8_4.4.2 svglite_2.1.3 stats4_4.4.2
## [67] xfun_0.44 expm_0.999-9 bridgesampling_1.1-2
## [70] matrixStats_1.3.0 rstan_2.32.6 stringi_1.8.4
## [73] yaml_2.3.8 boot_1.3-31 evaluate_0.24.0
## [76] codetools_0.2-20 xlsxjars_0.6.1 cli_3.6.3
## [79] RcppParallel_5.1.8 systemfonts_1.1.0 xtable_1.8-4
## [82] munsell_0.5.1 jquerylib_0.1.4 coda_0.19-4.1
## [85] parallel_4.4.2 rstantools_2.4.0 bayesplot_1.11.1
## [88] Brobdingnag_1.2-9 phangorn_2.11.1 viridisLite_0.4.2
## [91] scales_1.3.0 insight_0.20.2 combinat_0.0-8
## [94] rlang_1.1.4 cowplot_1.1.3 fastmatch_1.1-4
## [97] mnormt_2.1.1