Loading packages

library(plyr)
library(gridExtra)
library(tidyr)
library(dplyr)
library(data.table)
library(tidyverse)
library(ggpubr)
library(ggsci)
library(knitr)
library(kableExtra)
library(htmltools)

Data Manipulation

Convert values to CWM trait values

I’m using the datasheet and species proportional basal area (transformed from ~80% to correspond to 100%)

setwd("C:/Marina/Artigos/Nutrient demand")
data<-read.table("DADOS.txt", h=T)
NPP<-read.table("NPP_sami.txt", h=T)
data1 <- data %>% unite ("site.sp", Sp, site, sep= ".")
data1.cwm <- group_by(data1, site.sp) %>%   
  mutate_if(          
    is.numeric, funs(weighted.mean(.,AB_trans, na.rm=T))
  )

## `mutate_if()` ignored the following grouping variables:
## Column `site.sp`

data<- separate(data1.cwm, site.sp, c("Sp", "site"))

Nutrient resorption

Calculating resorption based on the formula

\[RE = (1- \frac{Nut_{old}}{Nut_{mat}} MLCF)*100\] , where MLCF is the green to senesced leaf Ca ratio.

I measured resorption with cwm values (among site) and for untransformed values (among species)

resorption.cwm<-transform(data,
                      resN = ((1-((N_old/N_mat)*(Ca_mat/Ca_old)))*100)) %>% transform(resP = ((1-((P_old/P_mat)*(Ca_mat/Ca_old)))*100)) %>%transform(resK = ((1-((K_old/K_mat)*(Ca_mat/Ca_old)))*100)) %>%transform(resMg = ((1-((Mg_old/Mg_mat)*(Ca_mat/Ca_old)))*100))

resorption<-transform(data1,
                      N = ((1-((N_old/N_mat)*(Ca_mat/Ca_old)))*100)) %>% transform(P = ((1-((P_old/P_mat)*(Ca_mat/Ca_old)))*100)) %>% transform(K = ((1-((K_old/K_mat)*(Ca_mat/Ca_old)))*100)) %>% transform(Mg = ((1-((Mg_old/Mg_mat)*(Ca_mat/Ca_old)))*100)) %>% separate(site.sp, c("Sp", "site"))

Figure 2 - Nutrient resorption

The first two panels (species values - untransformed data)

RES<-resorption[c(1:2,41:44)]

RES<- RES %>% gather("Nutrient","mean", 3:6) %>% mutate(site = revalue(site, c("CT"= "Savanna", "CD"="Transition Forest"))) 
+704

## [1] 704

g1<-ggplot(data=subset(RES, site=="Transition Forest"), aes(x=Sp, y=mean, fill=Nutrient))+
  labs(y= "Resorption (%)", x="Transition Forest species") + ylim(-50,100)+
  scale_fill_npg()+
  geom_boxplot() +
  theme_classic()+labs(tag = "A")+
  guides(fill=FALSE)

g2<-ggplot(data=subset(RES, site=="Savanna"), aes(x=Sp, y=mean, fill=Nutrient))+
  labs(y= "Resorption (%)", x="Savanna species") +ylim(-50,100)+
  geom_boxplot()+
  scale_fill_npg()+ labs(tag = "B")+
  theme_classic()+
  guides(fill=FALSE)

The third panel (site mean) with CWM

RES.cwm<-resorption.cwm[c(1:2,41:44)]
RES.cwm<- RES.cwm %>% gather("Nutrient","mean", 3:6) %>% mutate(site = revalue(site, c("CT"= "Savanna", "CD"="Transition Forest"))) %>% mutate(Nutrient = revalue(Nutrient, c("resN"= "N", "resP"="P","resK"="K", "resMg"="Mg")))

RES.cwm$Nutrient<-factor(RES.cwm$Nutrient, levels= c("N","P","K","Mg"))

g3<-ggplot(data=RES.cwm, aes(x=factor(site), y=mean, fill=Nutrient)) +
  labs(y= "Resorption (%)", x="Site") + ylim(-20,100)+
  geom_boxplot()+
  scale_fill_npg()+ labs(tag = "C")+
  theme_classic()

Figure 2

Figure 3

Calculating nutrient demand

Data provided by Sami (for first panel (NPP) in Figure 3)

NPP$component<-factor(NPP$component, levels= c("Canopy","Wood","Fine_roots"))
NPP<- NPP %>% mutate(component = revalue(component, c("Fine_roots"= "Fine roots")))
NPP<- NPP %>% mutate(site = revalue(site, c("CT"= "Savanna", "CD"="Transition Forest")))

library(RColorBrewer)
my_palette <- brewer.pal(name="Greens",n=4)[2:5]

g.NPP<-ggplot(data=NPP, aes(x=site, y=value, fill=component))+
  labs(y= expression(paste("Mg C  ", ha^-1, year^-1)), x=" ") +
  geom_bar(stat="identity", width=0.5, position=position_dodge())+
  geom_errorbar(aes(ymin=value-sem, ymax=value+sem), width=.2,
  position=position_dodge(.5),stat="identity") +
  scale_fill_manual(values=my_palette)+
  theme_classic()+
  guides(fill=FALSE)

Organizing dataset

For calculating nutrient demand and error propagation (combining CWM values and NPP values)

final<- resorption.cwm[c(1:2, 5:9,23:27, 32:36, 41:44)]

error<-ddply(final, c("site"),
         numcolwise(sd), na.rm=T) %>% mutate(site = revalue(site, c("CT"= "Savanna", "CD"="Transition Forest")))
medias<- ddply(final, c("site"),
               numcolwise(mean), na.rm=T) %>% mutate(site = revalue(site, c("CT"= "Savanna", "CD"="Transition Forest")))

NPP_medias<- spread(NPP[1:3], component, value)
NPP_error<- spread(NPP[c(1:2,4)], component, sem)
medias1<- merge(medias, NPP_medias)
error1<- merge(error, NPP_error)

Nutrient demand for different components

Demand

canopy_res<-transform(medias1,
                      N = (Canopy/(39/N_mat)*resN)) %>% transform(P = (Canopy/(39/P_mat)*resP))%>% transform(K = (Canopy/(39/K_mat)*resK)) %>% transform(Ca = (Canopy/(39/Ca_mat)*0))%>% transform(Mg = (Canopy/(39/Mg_mat)*resMg)) #assuming C 39%

canopy_resoprtion<- canopy_res[c(1, 24:28)]

canopy<-transform(canopy_res,
                      N = (Canopy/(50/(N_mat*100))-N)) %>% transform(P = (Canopy/(50/(P_mat*100))-P))%>% transform(K = (Canopy/(50/(K_mat*100))-K)) %>% transform(Ca = (Canopy/(50/(Ca_mat*100))-Ca))%>% transform(Mg = (Canopy/(50/(Mg_mat*100))-Mg))#assuming C 50%

canopy<- canopy[c(1, 24:28)]

wood<-transform(medias1,
                N = (Wood/(48/(N_SW*100)))) %>% transform(P = (Wood/(48/(P_SW*100))))%>% transform(K = (Wood/(48/(K_SW*100)))) %>% transform(Ca = (Wood/(48/(Ca_SW*100))))%>% transform(Mg = (Wood/(48/(Mg_SW*100))))#assuming C 48%

wood<- wood[c(1, 24:28)]

root<-transform(medias1,
                N = (`Fine roots`/(50/(N_root*100)))) %>% transform(P = (Fine.roots/(50/(P_root*100))))%>% transform(K = (Fine.roots/(50/(K_root*100)))) %>% transform(Ca = (Fine.roots/(50/(Ca_root*100))))%>% transform(Mg = (Fine.roots/(50/(Mg_root*100))))#assuming C 50%

root<- root[c(1, 24:28)]

bio<-rbind(canopy_resoprtion, canopy, wood, root) %>% transform(component = c("canopy_res", "canopy_res", "canopy", "canopy", "wood", "wood", "root", "root"))
kable(bio, digits=2) %>%  kable_styling(full_width = F, bootstrap_options = c("striped", "hover"))

site	N	P	K	Ca	Mg	component
Savanna	42.35	2.63	17.38	0.00	4.61	canopy_res
Transition Forest	83.09	5.41	33.65	0.00	5.55	canopy_res
Savanna	28.94	1.12	10.04	28.07	12.04	canopy
Transition Forest	83.43	2.97	18.64	44.42	20.83	canopy
Savanna	23.51	4.11	15.31	8.86	4.28	wood
Transition Forest	15.93	0.73	6.55	6.81	3.80	wood
Savanna	70.99	4.54	27.27	19.55	12.05	root
Transition Forest	77.83	2.87	20.78	23.64	9.41	root

Error Propagation

Propagating the error for nutrient demand calculations using quadrature approach

canopy_res.e<-transform(error1,
                      N.e = ((sqrt((N_mat/medias1$N_mat)^2+(Canopy/medias1$Canopy)^2+(resN/medias1$resN)^2))))%>% transform(P.e = (sqrt((P_mat/medias1$P_mat)^2+(Canopy/medias1$Canopy)^2+(resP/medias1$resP)^2)))%>% transform(K.e = (sqrt((K_mat/medias1$K_mat)^2+(Canopy/medias1$Canopy)^2+(resK/medias1$resK)^2))) %>% transform(Ca.e = 0)%>% transform(Mg.e = (sqrt((Mg_mat/medias1$Mg_mat)^2+(Canopy/medias1$Canopy)^2+(resMg/medias1$resMg)^2)))
canopy_res.sem<- canopy_res.e[c(1, 24:28)]

canopy.e<-transform(error1,
                        N.e = (sqrt((N_mat/medias1$N_mat)^2+(Canopy/medias1$Canopy)^2)))%>% transform(P.e = (sqrt((P_mat/medias1$P_mat)^2+(Canopy/medias1$Canopy)^2)))%>% transform(K.e = (sqrt((K_mat/medias1$K_mat)^2+(Canopy/medias1$Canopy)^2))) %>% transform(Ca.e = (sqrt((Ca_mat/medias1$Ca_mat)^2+(Canopy/medias1$Canopy)^2)))%>% transform(Mg.e = (sqrt((Mg_mat/medias1$Mg_mat)^2+(Canopy/medias1$Canopy)^2)))
canopy.sem<- canopy.e[c(1, 24:28)]

wood.e<-transform(error1,
                    N.e = (sqrt((N_SW/medias1$N_SW)^2+(Wood/medias1$Wood)^2)))%>% transform(P.e = (sqrt((P_SW/medias1$P_SW)^2+(Wood/medias1$Wood)^2)))%>% transform(K.e = (sqrt((K_SW/medias1$K_SW)^2+(Wood/medias1$Wood)^2))) %>% transform(Ca.e = (sqrt((Ca_SW/medias1$Ca_SW)^2+(Wood/medias1$Wood)^2)))%>% transform(Mg.e = (sqrt((Mg_SW/medias1$Mg_SW)^2+(Wood/medias1$Wood)^2)))
wood.sem<- wood.e[c(1, 24:28)]

root.e<-transform(error1,
                  N.e = (sqrt((N_root/medias1$N_root)^2+(`Fine roots`/medias1$`Fine roots`)^2)))%>% transform(P.e = (sqrt((P_root/medias1$P_root)^2+(Fine.roots/medias1$`Fine roots`)^2)))%>% transform(K.e = (sqrt((K_root/medias1$K_root)^2+(Fine.roots/medias1$`Fine roots`)^2))) %>% transform(Ca.e = (sqrt((Ca_root/medias1$Ca_root)^2+(Fine.roots/medias1$`Fine roots`)^2)))%>% transform(Mg.e = (sqrt((Mg_root/medias1$Mg_root)^2+Fine.roots/medias1$`Fine roots`)^2))
root.sem<- root.e[c(1, 24:28)]

bio.e<-rbind(canopy_res.sem, canopy.sem, wood.sem, root.sem) 
kable(bio.e, digits=2) %>%  kable_styling(full_width = F, bootstrap_options = c("striped", "hover"))

site	N.e	P.e	K.e	Ca.e	Mg.e
Savanna	0.48	0.51	0.63	0.00	1.33
Transition Forest	0.40	0.38	0.52	0.00	0.91
Savanna	0.24	0.35	0.45	0.56	0.42
Transition Forest	0.22	0.27	0.37	0.80	0.35
Savanna	0.34	0.62	0.43	0.99	0.71
Transition Forest	0.44	0.46	0.40	0.99	0.80
Savanna	0.25	0.30	0.25	0.26	0.25
Transition Forest	0.27	0.26	0.29	0.33	0.27

Other panels for Figure 3

bio.t<- cbind(bio, bio.e)
bio.t<-bio.t[2:13]
bio.t$component<-factor(bio$component, levels= c("canopy_res","canopy","wood","root")) 
bio<- mutate(bio.t, component = revalue(component, c("canopy_res"= "Canopy Resorption", "canopy"="Canopy", "wood"="Wood", "root"="Fine roots")))

g.N<-ggplot(data=bio, aes(x=site, y=N, fill=component))+
  labs(y= expression(paste("kg N  ", ha^-1, year^-1)), x=" ") +
  geom_bar(stat="identity", width=0.5, position=position_dodge())+
  geom_errorbar(aes(ymin=N-(N.e), ymax=N+(N.e)), width=.2,
                position=position_dodge(.5),stat="identity") +
  scale_fill_brewer(palette = "Greens")+
  theme_classic()+ ylim(0,200)+
  theme(legend.position=c(0.5,1), legend.direction= "horizontal", legend.title=element_blank())

g.P<-ggplot(data=bio, aes(x=site, y=P, fill=component))+
  labs(y= expression(paste("kg P  ", ha^-1, year^-1)), x=" ") +
  geom_bar(stat="identity", width=0.5, position=position_dodge())+
  geom_errorbar(aes(ymin=P-(P.e), ymax=P+(P.e)), width=.2,
                position=position_dodge(.5),stat="identity") +
  scale_fill_brewer(palette = "Greens")+
  theme_classic()+guides(fill=FALSE)

g.K<-ggplot(data=bio, aes(x=site, y=K, fill=component))+
  labs(y= expression(paste("kg K  ", ha^-1, year^-1)), x=" ") +
  geom_bar(stat="identity", width=0.5, position=position_dodge())+
  geom_errorbar(aes(ymin=K-(K.e), ymax=K+(K.e)), width=.2,
                position=position_dodge(.5),stat="identity") +
  scale_fill_brewer(palette = "Greens")+
  theme_classic()+guides(fill=FALSE)

g.Ca<-ggplot(data=bio, aes(x=site, y=Ca, fill=component))+
  labs(y= expression(paste("kg Ca  ", ha^-1, year^-1)), x=" ") +
  geom_bar(stat="identity", width=0.5, position=position_dodge())+
  geom_errorbar(aes(ymin=Ca-(Ca.e), ymax=Ca+(Ca.e)), width=.2,
                position=position_dodge(.5),stat="identity") +
  scale_fill_brewer(palette = "Greens")+
  theme_classic()+guides(fill=FALSE)

g.Mg<-ggplot(data=bio, aes(x=site, y=Mg, fill=component))+
  labs(y= expression(paste("kg Mg  ", ha^-1, year^-1)), x=" ") +
  geom_bar(stat="identity", width=0.5, position=position_dodge())+
  geom_errorbar(aes(ymin=Mg-(Mg.e), ymax=Mg+(Mg.e)), width=.2,
                position=position_dodge(.5),stat="identity") +
  scale_fill_brewer(palette = "Greens")+
  theme_classic()+guides(fill=FALSE)

Figure 3

Figure 4

Efficiencies calculation

N<- spread(bio.t[c(1, 6:7)], component, N)
P<- spread(bio.t[c(2, 6:7)], component, P)
K<- spread(bio.t[c(3, 6:7)], component, K)
Ca<- spread(bio.t[c(4, 6:7)], component, Ca)
Mg<- spread(bio.t[c(5, 6:7)], component, Mg)

uptake1<-rbind(N,P,K,Ca,Mg) %>% transform(Nutrient = c("N", "N", "P", "P", "K", "K", "Ca", "Ca", "Mg", "Mg"))

Uptake<- uptake1 %>% group_by(site, Nutrient) %>% transform(uptake = canopy+wood+root,
                                                            demand = canopy_res+canopy+wood+root)
Uptake$NPP<- ifelse(Uptake$site == "Transition Forest", (10.51323), (8.797829))


Uptake_eff<- Uptake %>% group_by(site, Nutrient) %>% transform(uptake_eff = ((NPP/uptake)*1000),
                                                            use_eff = ((NPP/demand)*1000))#from Mg to Kg

Error calculation

N.e<- spread(bio.t[6:8], component, N.e)
P.e<- spread(bio.t[c(6:7, 9)], component, P.e)
K.e<- spread(bio.t[c(6:7, 10)], component, K.e)
Ca.e<- spread(bio.t[c(6:7, 11)], component, Ca.e)
Mg.e<- spread(bio.t[c(6:7, 12)], component, Mg.e)

uptake1.e<-rbind(N.e,P.e,K.e,Ca.e,Mg.e) %>% transform(Nutrient = c("N.e", "N.e", "P.e", "P.e", "K.e", "K.e", "Ca.e", "Ca.e", "Mg.e", "Mg.e"))

Uptake.e<- uptake1.e %>% group_by(site, Nutrient) %>% transform(uptake.e = (sqrt(canopy^2+wood^2+root^2)*1000),
                                                                demand.e = (sqrt(canopy_res^2+canopy^2+wood^2+root^2))*1000) 

Uptake.e$NPP.e<-ifelse(Uptake.e$site == "Transition Forest", 0.96397, 0.930912)

Uptake.t<-cbind(Uptake, Uptake.e)[c(6:10, 16:18)]

Uptake_eff.e<- Uptake.t %>% group_by(site, Nutrient) %>% transform(uptake_eff.e = sqrt(((NPP.e/NPP)^2)+((uptake.e/uptake)^2)), use_eff.e = sqrt(((NPP.e/NPP)^2)+((demand.e/demand)^2)))

NUE<-cbind(Uptake_eff, Uptake_eff.e)[c(1, 10:12, 20,21)]

NUE$Nutrient<-factor(NUE$Nutrient, levels= c("N","P","K", "Ca", "Mg"))

Figure 4

g.uptake<-ggplot(data=NUE, aes(x=site, y=uptake_eff, fill=Nutrient))+
  labs(y= expression(paste("Uptake efficiency (kg C per kg [Nut])")), x="Site") +
  geom_bar(stat="identity", position=position_dodge())+
  geom_errorbar(aes(ymin=uptake_eff-uptake_eff.e, ymax=uptake_eff+uptake_eff.e), width=.2,
                position=position_dodge(.5),stat="identity") +
  facet_wrap(~Nutrient, ncol= 5, scales = "free_y")+
  scale_fill_npg()+
  theme_classic()
g.uptake

g.use<-ggplot(data=NUE, aes(x=site, y=use_eff, fill=Nutrient))+
  labs(y= expression(paste("Use efficiency (kg C per kg [Nut])")), x="Site") +
  geom_bar(stat="identity", position=position_dodge())+
  geom_errorbar(aes(ymin=use_eff-use_eff.e, ymax=use_eff+use_eff.e), width=.2,
                position=position_dodge(.5),stat="identity") +
  facet_wrap(~Nutrient, ncol= 5, scales = "free_y")+
  scale_fill_npg()+
  theme_classic()
g.use

fig4<-ggarrange(g.uptake, g.use, ncol=1, nrow=2, common.legend = TRUE, legend="right")
fig4

Table 1

Proportional contribution (Table 1)

Calculating values

prop<- uptake1 %>% group_by(site, Nutrient) %>%  transform (sum.t = (canopy+canopy_res+wood+root)) %>% 
                                                 transform  (total_canopy.p = ((canopy+canopy_res)/sum.t)*100,
                                                            resorption.p = (canopy_res/(canopy+canopy_res))*100,
                                                            wood.p = (wood/sum.t)*100,
                                                            root.p = (root/ sum.t)*100) %>%  
                                                  transform (sum.p = (total_canopy.p+wood.p+root.p))

Error propagation

uptake.t<-cbind(prop, uptake1.e)[c(1:11, 14:17)]
e.t<- uptake.t %>% group_by(site, Nutrient) %>% transform (total.e = sqrt((canopy.1^2+canopy_res.1^2+wood.1^2+root.1^2)),
                                                           canopy.e.t = sqrt((canopy.1^2+canopy_res.1^2)),
                                                           sum.canopy = ((canopy+canopy_res)))

prop.e<- e.t %>% group_by(site, Nutrient) %>%  transform(total_canopy.e = sqrt((canopy.e.t/canopy)^2+(total.e/sum.t)^2)*100,
                                                           resorption.e = sqrt((canopy_res.1/canopy_res)^2+(canopy.e.t/sum.canopy)^2)*100,
                                                           wood.e = sqrt((wood.1/wood)^2+(total.e/sum.t)^2)*100,
                                                           root.e = sqrt((root.1/root)^2+ (total.e/sum.t)^2)*100)

table1<-prop.e[c(1,6,8:11, 19:22)]
kable(table1, digits=2)%>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed", "responsive"))

site	Nutrient	total_canopy.p	resorption.p	wood.p	root.p	total_canopy.e	resorption.e	wood.e	root.e
Savanna	N	43.00	59.41	14.18	42.82	1.91	1.37	1.49	0.54
Transition Forest	N	63.98	49.90	6.12	29.90	0.61	0.55	2.78	0.43
Savanna	P	30.26	70.02	33.14	36.60	55.62	25.56	16.86	10.04
Transition Forest	P	69.97	64.55	6.11	23.92	16.83	8.96	63.40	10.97
Savanna	K	39.17	63.40	21.87	38.96	7.78	4.57	3.08	1.60
Transition Forest	K	65.68	64.35	8.23	26.10	3.57	1.97	6.21	1.74
Savanna	Ca	49.71	0.00	15.68	34.61	2.87	NaN	11.33	2.45
Transition Forest	Ca	59.33	0.00	9.09	31.58	2.51	NaN	14.66	2.24
Savanna	Mg	50.49	27.68	12.98	36.54	12.57	30.15	17.24	5.23
Transition Forest	Mg	66.63	21.05	9.60	23.77	5.70	16.80	21.40	4.38

Error Propagation

Marina Scalon

May 30, 2020

Loading packages

Data Manipulation

Convert values to CWM trait values

Nutrient resorption

Figure 2 - Nutrient resorption

The first two panels (species values - untransformed data)

The third panel (site mean) with CWM

Figure 2

Figure 3

Calculating nutrient demand

Organizing dataset

Nutrient demand for different components

Demand

Error Propagation

Other panels for Figure 3

Figure 3

Figure 4

Efficiencies calculation

Error calculation

Figure 4

Table 1

Proportional contribution (Table 1)

Error propagation