The effect of Tithonia diversifolia in a Brazilian silvopastoral system with Urochloa brizantha cv. Xaraés

---

Packages

library(dplyr)
library(ggplot2)
library(ggh4x)#grafico
library(geobr) #pacote do ibge para polignos do BR
library(elevatr) #Haster de relevo br
library(sp) #Polignos de Sao Paulo
library(sf) #vinculado ao ggplot2
library(cowplot)#combinar varios graficos
library(lme4)#modelo misto
library(lmerTest)#fornece p value da anova
library(emmeans)
library(multcompView)

citation("dplyr")

To cite package ‘dplyr’ in publications use:

  Wickham H, François R, Henry L, Müller K, Vaughan D (2023). _dplyr: A Grammar of Data Manipulation_.
  R package version 1.1.4, <https://CRAN.R-project.org/package=dplyr>.

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {dplyr: A Grammar of Data Manipulation},
    author = {Hadley Wickham and Romain François and Lionel Henry and Kirill Müller and Davis Vaughan},
    year = {2023},
    note = {R package version 1.1.4},
    url = {https://CRAN.R-project.org/package=dplyr},
  }

citation("ggplot2")

To cite ggplot2 in publications, please use

  H. Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.

A BibTeX entry for LaTeX users is

  @Book{,
    author = {Hadley Wickham},
    title = {ggplot2: Elegant Graphics for Data Analysis},
    publisher = {Springer-Verlag New York},
    year = {2016},
    isbn = {978-3-319-24277-4},
    url = {https://ggplot2.tidyverse.org},
  }

citation("ggh4x")

To cite package ‘ggh4x’ in publications use:

  van den Brand T (2024). _ggh4x: Hacks for 'ggplot2'_. R package version 0.2.8,
  <https://CRAN.R-project.org/package=ggh4x>.

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {ggh4x: Hacks for 'ggplot2'},
    author = {Teun {van den Brand}},
    year = {2024},
    note = {R package version 0.2.8},
    url = {https://CRAN.R-project.org/package=ggh4x},
  }

citation("geobr")

To cite package ‘geobr’ in publications use:

  Pereira R, Goncalves C (2024). _geobr: Download Official Spatial Data Sets of Brazil_. R package
  version 1.8.2, <https://CRAN.R-project.org/package=geobr>.

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {geobr: Download Official Spatial Data Sets of Brazil},
    author = {Rafael H. M. Pereira and Caio Nogueira Goncalves},
    year = {2024},
    note = {R package version 1.8.2},
    url = {https://CRAN.R-project.org/package=geobr},
  }

citation("elevatr")

To cite elevatr:

  Hollister, J.W. (2023). elevatr: Access Elevation Data from Various APIs. R package version 0.99.0.
  https://CRAN.R-project.org/package=elevatr/

A BibTeX entry for LaTeX users is

  @Manual{,
    author = {Jeffrey Hollister and Tarak Shah and Jakub Nowosad and Alec L. Robitaille and Marcus W. Beck and Mike Johnson},
    title = {elevatr: Access Elevation Data from Various APIs},
    year = {2023},
    note = {R package version 0.99.0},
    doi = {10.5281/zenodo.8335450},
    url = {https://github.com/jhollist/elevatr/},
  }

citation("sp")

To cite package sp in publications use:

  Pebesma E, Bivand R (2005). “Classes and methods for spatial data in R.” _R News_, *5*(2), 9-13.
  <https://CRAN.R-project.org/doc/Rnews/>.

  Bivand R, Pebesma E, Gomez-Rubio V (2013). _Applied spatial data analysis with R, Second edition_.
  Springer, NY. <https://asdar-book.org/>.

To see these entries in BibTeX format, use 'print(<citation>, bibtex=TRUE)', 'toBibtex(.)', or set
'options(citation.bibtex.max=999)'.

citation("sf")

To cite package sf in publications, please use:

  Pebesma, E., & Bivand, R. (2023). Spatial Data Science: With Applications in R. Chapman and Hall/CRC.
  https://doi.org/10.1201/9780429459016

  Pebesma, E., 2018. Simple Features for R: Standardized Support for Spatial Vector Data. The R Journal
  10 (1), 439-446, https://doi.org/10.32614/RJ-2018-009

To see these entries in BibTeX format, use 'print(<citation>, bibtex=TRUE)', 'toBibtex(.)', or set
'options(citation.bibtex.max=999)'.

citation("cowplot")

To cite package ‘cowplot’ in publications use:

  Wilke C (2024). _cowplot: Streamlined Plot Theme and Plot Annotations for 'ggplot2'_. R package
  version 1.1.3, <https://CRAN.R-project.org/package=cowplot>.

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {cowplot: Streamlined Plot Theme and Plot Annotations for 'ggplot2'},
    author = {Claus O. Wilke},
    year = {2024},
    note = {R package version 1.1.3},
    url = {https://CRAN.R-project.org/package=cowplot},
  }

citation("lme4")

To cite lme4 in publications use:

  Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015). Fitting Linear Mixed-Effects Models
  Using lme4. Journal of Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01.

A BibTeX entry for LaTeX users is

  @Article{,
    title = {Fitting Linear Mixed-Effects Models Using {lme4}},
    author = {Douglas Bates and Martin M{\"a}chler and Ben Bolker and Steve Walker},
    journal = {Journal of Statistical Software},
    year = {2015},
    volume = {67},
    number = {1},
    pages = {1--48},
    doi = {10.18637/jss.v067.i01},
  }

citation("lmerTest")

To cite lmerTest in publications use:

  Kuznetsova A, Brockhoff PB, Christensen RHB (2017). “lmerTest Package: Tests in Linear Mixed Effects
  Models.” _Journal of Statistical Software_, *82*(13), 1-26. doi:10.18637/jss.v082.i13
  <https://doi.org/10.18637/jss.v082.i13>.

A BibTeX entry for LaTeX users is

  @Article{,
    title = {{lmerTest} Package: Tests in Linear Mixed Effects Models},
    author = {Alexandra Kuznetsova and Per B. Brockhoff and Rune H. B. Christensen},
    journal = {Journal of Statistical Software},
    year = {2017},
    volume = {82},
    number = {13},
    pages = {1--26},
    doi = {10.18637/jss.v082.i13},
  }

citation("emmeans")

To cite package ‘emmeans’ in publications use:

  Lenth R (2024). _emmeans: Estimated Marginal Means, aka Least-Squares Means_. R package version
  1.10.0, <https://CRAN.R-project.org/package=emmeans>.

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {emmeans: Estimated Marginal Means, aka Least-Squares Means},
    author = {Russell V. Lenth},
    year = {2024},
    note = {R package version 1.10.0},
    url = {https://CRAN.R-project.org/package=emmeans},
  }

citation("multcompView")

To cite package ‘multcompView’ in publications use:

  Graves S, Piepho H, Dorai-Raj LSwhfS (2023). _multcompView: Visualizations of Paired Comparisons_. R
  package version 0.1-9, <https://CRAN.R-project.org/package=multcompView>.

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {multcompView: Visualizations of Paired Comparisons},
    author = {Spencer Graves and Hans-Peter Piepho and Luciano Selzer with help from Sundar Dorai-Raj},
    year = {2023},
    note = {R package version 0.1-9},
    url = {https://CRAN.R-project.org/package=multcompView},
  }

ATTENTION: This citation information has been auto-generated from the package DESCRIPTION file and may
need manual editing, see ‘help("citation")’.

Weather conditions

Data

#Weather
clima3.1=read.csv("D:\\Armazenamento\\DATA R\\Ensaio 3\\clima3.1.csv")
clima3.1$data = as.Date(as.character(clima3.1$data, format = "%Y%m$d"))
clima3.1 = arrange(clima3.1, data)
str(clima3.1)

'data.frame':   730 obs. of  3 variables:
 $ data : Date, format: "2022-01-01" "2022-01-01" "2022-01-02" "2022-01-02" ...
 $ group: chr  "Rainfall" "Radiation" "Rainfall" "Radiation" ...
 $ value: num  37.3 17 0.3 22.2 2.5 ...

print(clima3.1)

clima3.2=read.csv("D:\\Armazenamento\\DATA R\\Ensaio 3\\clima3.2.csv")
clima3.2$data = as.Date(as.character(clima3.2$data, format = "%Y%m$d"))
clima3.2 = arrange(clima3.2, data)
str(clima3.2)

'data.frame':   730 obs. of  3 variables:
 $ data : Date, format: "2022-01-01" "2022-01-01" "2022-01-02" "2022-01-02" ...
 $ group: chr  "RH" "Temperature" "RH" "Temperature" ...
 $ value: num  89.9 22.1 78.4 24.6 86.8 24.3 89.6 23.3 83.7 24.4 ...

print(clima3.2)

Plot

datebreaks = seq(as.Date("2022-01-01"), as.Date("2022-12-31"), by = "20 days")
p1=ggplot() +
  geom_bar(data = clima3.1, aes(data,value*1.3,fill=group),stat="identity",position = "dodge", width = 1)+
  geom_point(data = clima3.2, aes(data, value, shape=group), size=2)+
  geom_line(data = clima3.2, aes(data, value, group = group), color = "black")+
 scale_shape_manual(values = c(1, 16), name= NULL,labels = c("Relative humidity (%)","Temperature (ºC)"))+
  scale_fill_manual(values=c("gray28","darkblue"), name= NULL,labels=c("Radiation (MJ/m2 day)","Rainfall (mm)"))+
  scale_y_continuous(name= "Relative humidity (%); Temperature (ºC)",breaks = seq(0, 125, 5),   sec.axis=sec_axis(~./1.3,name="Radiation (MJ/m2 day); Rainfall (mm)",breaks=seq(0,98,4)))+coord_cartesian(ylim = c(0,120))+scale_x_date(name=NULL,breaks = datebreaks)
p1

p2=p1+theme(panel.background = element_rect(fill = "transparent", colour = "black"),
            legend.position = c(0.25, 0.88), 
            legend.direction="horizontal",
            legend.background = element_rect(fill = "transparent", size=0.5, linetype="solid",colour ="black"),
            legend.key.height = unit(0.4,"cm"),  
            legend.key.width = unit(0.4,"cm"),
            legend.text = element_text(size=14),
            legend.key=element_blank(),
            axis.line = element_line(colour = "black", size = 0.7, linetype = "solid"),
            axis.text.x = element_text(angle = 30, hjust = 1,size = 12),
            axis.title.y = element_text(size = 14),
            axis.text.y = element_text(size = 12),
            )
p2

Experimental site

Data

BRA2 = read_country()

Using year 2010

states = read_state()

Using year 2010
Warning: GDAL Error 1: database disk image is malformedWarning: GDAL Error 1: sqlite3_prepare_v2(SELECT COUNT(*) FROM sqlite_master WHERE name IN ('gpkg_metadata', 'gpkg_metadata_reference') AND type IN ('table', 'view')) failed: database disk image is malformedErro: Cannot open "C:\Users\Samsung\AppData\Local\Temp\RtmpItEjVq\31state_2010_simplified.gpkg"; The source could be corrupt or not supported. See `st_drivers()` for a list of supported formats.

País

Estado

p2.1 = ggplot()+
  geom_sf(data=SP, color="blacK", fill="wheat")+
  geom_raster(data = topo_df, aes(x = x, y = y, fill = elevation))+
  scale_fill_gradient(low = "yellow", high = "darkgreen", name = "Topography (m)")+
  geom_sf(data=piracicaba, color="blacK", fill="transparent")+
  theme_light()+
  labs(x="Longitude", y="Latitude")+
  ylim(-25.3,-19.3)+xlim(-52.8,-44.5)
p2.1

p2.2 = p2.1+
  theme(panel.background = element_rect(fill = "white"),
  axis.text = element_text(angle = 10, hjust = 0.6),
  axis.line = element_line(colour = "gray", size = 0.3, linetype = "solid"),
  panel.grid.major = element_line(color = "transparent"), 
  panel.grid.minor = element_line(color = "transparent"),
  legend.position = "NONE")+  
  guides(fill = guide_colorbar(title.position = "top", title.hjust = 0.5,title.vjust = 1.5))+
  annotate(geom="text", y=-19.4, x=-48.7,label="State of São Paulo (SP), Brazil", size=6, color="black")
p2.2
save_plot("P2.2.pdf", p2.2, ncol = 1, nrow = 1)

Cidade

Experimental site

Biomass production

Data

data1=read.csv("D:\\Armazenamento\\DATA R\\Ensaio 3\\assay3.1.csv")
data1$Ciclo=as.factor(data1$Ciclo)
data1$tratamentos=as.factor(data1$tratamentos)
data1$entrada=as.Date(data1$entrada, format="%d/%m/%Y")
print(data1)

str(data1)

'data.frame':   96 obs. of  80 variables:
 $ Ciclo                : Factor w/ 6 levels "Cycle 1","Cycle 2",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Piquete              : chr  "8B" "7B" "8A" "7A" ...
 $ tratamentos          : Factor w/ 2 levels "consorcio","exclusivo": 2 1 2 1 2 1 2 1 2 1 ...
 $ rep                  : int  1 1 2 2 3 3 4 4 5 5 ...
 $ entrada              : Date, format: "2022-01-21" "2022-01-21" "2022-01-24" "2022-01-24" ...
 $ saida                : chr  "23/01/2022" "23/01/2022" "26/01/2022" "26/01/2022" ...
 $ ocupacao             : int  3 3 3 3 3 3 2 3 3 3 ...
 $ descanso             : int  24 24 24 24 24 24 24 24 24 24 ...
 $ pregrazing.height    : num  30.6 30.9 33.5 32.2 31.8 ...
 $ postgrazing.height   : num  16.9 17.2 16.1 15.4 17.6 ...
 $ pregrazing.height.thi: num  NA NA NA NA NA NA NA NA NA NA ...
 $ area.thi.paddock     : num  0 8.02 0 6.58 0 7.74 0 7.45 0 3.84 ...
 $ area.thi.ha          : num  0 1114 0 914 0 ...
 $ n.plantas            : num  0 29 0 25 0 23 0 26 0 15 ...
 $ diametro.copa        : num  0 0.55 0 0.53 0 0.67 0 0.57 0 0.51 ...
 $ densidade.planta     : num  0 4028 0 3472 0 ...
 $ biomass.bra.paddock  : num  6.24 6.38 10.04 5.73 9.76 ...
 $ biomass.bra.ha       : num  866 886 1395 795 1356 ...
 $ biomass.thi.paddock  : num  0 1.68 0 0.43 0 2.87 0 3.35 0 2.86 ...
 $ biomass.thi.ha       : num  0 233.5 0 59.6 0 ...
 $ biomass.total.paddock: num  6.24 8.06 10.04 6.16 9.76 ...
 $ biomass.total.ha     : num  866 1120 1395 855 1356 ...
 $ resi.bra.paddock     : num  2.94 2.88 1.09 0 1.67 1.56 0 0 2.14 3.93 ...
 $ resi.bra.ha          : num  409 399 151 0 232 ...
 $ resi.thi.paddock     : num  0 0.34 0 0 0 1.93 0 0 0 0 ...
 $ resi.thi.ha          : num  0 47.9 0 0 0 ...
 $ acumulado.bra        : num  866 886 1395 795 1356 ...
 $ acumulado.thi        : num  0 233.5 0 59.6 0 ...
 $ acumulado.total      : num  866 1120 1395 855 1356 ...
 $ resi.total.paddock   : num  2.94 3.22 1.09 0 1.67 3.49 0 0 2.14 3.93 ...
 $ resi.total.ha        : num  409 447 151 0 232 ...
 $ bra.accumulation.rate: num  36.1 36.9 58.1 33.1 56.5 ...
 $ thi.accumulation.rate: num  NA 9.73 NA 2.48 NA ...
 $ bra.folhas           : num  641 571 661 600 584 ...
 $ bra.hastes           : num  148 237 165 185 233 ...
 $ bra.morto            : num  210 192 174 215 184 ...
 $ thi.folhas           : num  0 468 0 520 0 ...
 $ thi.hastes           : num  0 425 0 319 0 ...
 $ thi.morto            : num  0 107 0 162 0 ...
 $ bra.leaf.acu         : num  555 506 922 477 791 ...
 $ bra.stem.acu         : num  129 210 230 147 315 ...
 $ bra.dead.acu         : num  182 170 242 171 249 ...
 $ thi.leaf.acu         : num  0 109 0 31 0 ...
 $ thi.stem.acu         : num  0 99.3 0 19 0 ...
 $ thi.dead.acu         : num  0 24.95 0 9.64 0 ...
 $ total.leaf.acu       : num  555 615 922 508 791 ...
 $ total.stem.acu       : num  129 310 230 166 315 ...
 $ total.dead.acu       : num  182 195 242 180 249 ...
 $ n.animais            : num  5 5.33 5 5 4.67 5 2.33 5 3.67 4.33 ...
 $ pv.piquete           : num  143 149 143 140 135 ...
 $ pv.ha                : num  825 860 825 809 783 ...
 $ ua.ha                : num  1.83 1.91 1.83 1.8 1.74 1.8 0.91 1.8 1.4 1.61 ...
 $ DM.bra               : num  905 903 898 900 896 ...
 $ MM.bra               : num  127 124 128 129 123 ...
 $ MO.bra               : num  873 876 872 871 877 ...
 $ PB.bra               : num  118 107 169 137 105 ...
 $ aFDNom.bra           : num  547 608 546 562 582 ...
 $ aFDA.bra             : num  292 314 290 288 308 ...
 $ HEMC.bra             : num  254 294 255 274 275 273 272 270 278 277 ...
 $ DM.thi               : num  0 881 0 881 0 ...
 $ MM.thi               : num  0 164 0 176 0 ...
 $ MO.thi               : num  0 836 0 824 0 ...
 $ PB.thi               : num  0 124 0 144 0 ...
 $ aFDNom.thi           : num  0 356 0 352 0 ...
 $ aFDA.thi             : num  0 245 0 251 0 ...
 $ HEMC.thi             : num  0 110 0 101 0 89 0 108 0 94 ...
 $ total.DM.ac          : num  784 1005 1253 768 1214 ...
 $ total.MM.ac          : num  110 148 178 113 166 ...
 $ total.MO.ac          : num  756 972 1216 742 1189 ...
 $ total.PB.ac          : num  102 124 236 117 143 ...
 $ total.FDN.ac         : num  474 622 761 468 789 ...
 $ total.FDA.ac         : num  253 336 405 244 417 ...
 $ total.HEMC.ac        : num  220 286 356 224 373 ...
 $ d15N                 : num  NA NA NA NA NA NA NA NA NA NA ...
 $ X.N                  : num  NA NA NA NA NA NA NA NA NA NA ...
 $ d13C                 : num  NA NA NA NA NA NA NA NA NA NA ...
 $ X.C                  : num  NA NA NA NA NA NA NA NA NA NA ...
 $ C.N                  : num  NA NA NA NA NA NA NA NA NA NA ...
 $ C3.                  : num  NA NA NA NA NA NA NA NA NA NA ...
 $ C4.                  : num  NA NA NA NA NA NA NA NA NA NA ...

Rest time

psych::describeBy(data1$descanso)

Warning: no grouping variable requested

1- Pre-grazing canopy height (cm)

mod1 = lmer(pregrazing.height~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod1))

shapiro.test(resid(mod1))

    Shapiro-Wilk normality test

data:  resid(mod1)
W = 0.97071, p-value = 0.03009

mod1.1 = lmer(pregrazing.height^1.3~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod1.1))

shapiro.test(resid(mod1.1))

    Shapiro-Wilk normality test

data:  resid(mod1.1)
W = 0.97413, p-value = 0.05431

bartlett.test(resid(mod1.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod1.1) by tratamentos
Bartlett's K-squared = 0.0033781, df = 1, p-value = 0.9537

anova(mod1.1)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
tratamentos 32.729  32.729     1    89  0.3517 0.5546

medias1=emmeans(mod1,~tratamentos)
summary(medias1)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio     28.4 1.48 5.38     24.7     32.1
 exclusivo     28.1 1.48 5.38     24.3     31.8

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

2- Post-grazing canopy height (cm)

mod2 = lmer(postgrazing.height~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod2))

shapiro.test(resid(mod2))

    Shapiro-Wilk normality test

data:  resid(mod2)
W = 0.93578, p-value = 0.000149

bartlett.test(resid(mod2)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod2) by tratamentos
Bartlett's K-squared = 6.7228, df = 1, p-value = 0.009519

mod2.1 = lmer(postgrazing.height^0.05~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod2.1))

shapiro.test(resid(mod2.1))

    Shapiro-Wilk normality test

data:  resid(mod2.1)
W = 0.97381, p-value = 0.05141

bartlett.test(resid(mod2.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod2.1) by tratamentos
Bartlett's K-squared = 3.4331, df = 1, p-value = 0.0639

anova(mod2.1)

Type III Analysis of Variance Table with Satterthwaite's method
                Sum Sq    Mean Sq NumDF DenDF F value  Pr(>F)  
tratamentos 7.8291e-05 7.8291e-05     1    89  3.1252 0.08052 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias2=emmeans(mod2,~tratamentos)
summary(medias2)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio     13.9 1.24 5.12     10.7     17.1
 exclusivo     14.4 1.24 5.12     11.2     17.6

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

3- Herbage biomass (kg ha-1)

mod3 = lmer(biomass.bra.ha~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod3))

shapiro.test(resid(mod3))

    Shapiro-Wilk normality test

data:  resid(mod3)
W = 0.97777, p-value = 0.1023

bartlett.test(resid(mod3)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod3) by tratamentos
Bartlett's K-squared = 1.0228, df = 1, p-value = 0.3119

anova(mod3)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
tratamentos 616.66  616.66     1    89  0.0036 0.9522

medias3=emmeans(mod3,~tratamentos)
summary(medias3)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio     1251 81.1 9.32     1069     1434
 exclusivo     1246 81.1 9.32     1064     1429

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

4- Residual biomass (kg ha-1)

mod4 = lmer(resi.bra.ha~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod4))

shapiro.test(resid(mod4))

    Shapiro-Wilk normality test

data:  resid(mod4)
W = 0.85143, p-value = 2.208e-08

bartlett.test(resid(mod4)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod4) by tratamentos
Bartlett's K-squared = 0.10694, df = 1, p-value = 0.7437

mod4.1 = lmer(resi.bra.ha^0.3~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod4.1))

shapiro.test(resid(mod4.1)) 

    Shapiro-Wilk normality test

data:  resid(mod4.1)
W = 0.97707, p-value = 0.09052

bartlett.test(resid(mod4.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod4.1) by tratamentos
Bartlett's K-squared = 0.090194, df = 1, p-value = 0.7639

anova(mod4.1)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
tratamentos 6.9044  6.9044     1    89  1.3308 0.2518

medias4=emmeans(mod4,~tratamentos)
summary(medias4)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      152 48.7 6.86     36.0      267
 exclusivo      183 48.7 6.86     67.1      298

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

5- Herbage accumulation rate (kg ha-1 day-1)

mod5 = lmer(bra.accumulation.rate~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod5))

shapiro.test(resid(mod5))

    Shapiro-Wilk normality test

data:  resid(mod5)
W = 0.95934, p-value = 0.004587

bartlett.test(resid(mod5)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod5) by tratamentos
Bartlett's K-squared = 0.23164, df = 1, p-value = 0.6303

mod5.1 = lmer(bra.accumulation.rate^0.8~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod5.1))

shapiro.test(resid(mod5.1)) 

    Shapiro-Wilk normality test

data:  resid(mod5.1)
W = 0.98086, p-value = 0.1743

bartlett.test(resid(mod5.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod5.1) by tratamentos
Bartlett's K-squared = 2.4364e-07, df = 1, p-value = 0.9996

anova(mod5.1)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
tratamentos 1.7106  1.7106     1    89  0.1152 0.7351

medias5=emmeans(mod5,~tratamentos)
summary(medias5)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio     31.1 6.65 5.24     14.3     48.0
 exclusivo     30.4 6.65 5.24     13.5     47.2

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

6- Tithonia - Herbage biomass (kg ha-1)

mod6 = lmer(biomass.thi.ha~tratamentos+(1|Ciclo), data = data1)

medias6=emmeans(mod6,~tratamentos)
summary(medias6)

 tratamentos emmean  SE   df lower.CL upper.CL
 consorcio      933 157 6.59      558     1308
 exclusivo        0 157 6.59     -375      375

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

7- Tithonia - Residual biomass (kg ha-1)

mod7 = lmer(resi.thi.ha~tratamentos+(1|Ciclo), data = data1)

medias7=emmeans(mod7,~tratamentos)
summary(medias7)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      361 99.1 7.01      127      595
 exclusivo        0 99.1 7.01     -234      234

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

8- Tithonia - Plant density (N ha-1)

mod8 = lmer(densidade.planta~tratamentos+(1|Ciclo), data = data1)

medias8=emmeans(mod8,~tratamentos)
summary(medias8)

 tratamentos emmean SE df lower.CL upper.CL
 consorcio     2821 71 13     2668     2975
 exclusivo        0 71 13     -153      153

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

9 - Tithonia - Pre-grazing canopy height (cm)

mod9 = lmer(pregrazing.height.thi~tratamentos+(1|Ciclo), data = data1)

medias9=emmeans(mod9,~tratamentos)
summary(medias9)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      118 6.94 1.55     77.9    158.1
 exclusivo        0 6.94 1.55    -40.1     40.1

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

10- Total - Herbage biomass (kg ha-1)

mod10 = lmer(biomass.total.ha~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod10))

shapiro.test(resid(mod10))

    Shapiro-Wilk normality test

data:  resid(mod10)
W = 0.87096, p-value = 1.225e-07

bartlett.test(resid(mod10)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod10) by tratamentos
Bartlett's K-squared = 16.657, df = 1, p-value = 4.478e-05

mod10.1 = lmer(biomass.total.ha^0.3~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod10.1))

shapiro.test(resid(mod10.1)) 

    Shapiro-Wilk normality test

data:  resid(mod10.1)
W = 0.97673, p-value = 0.0854

bartlett.test(resid(mod10.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod10.1) by tratamentos
Bartlett's K-squared = 0.82181, df = 1, p-value = 0.3647

anova(mod10.1)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
tratamentos 50.947  50.947     1    89  45.957 1.274e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias10=emmeans(mod10,~tratamentos)
summary(medias10)

 tratamentos emmean  SE   df lower.CL upper.CL
 consorcio     2185 198 6.87     1714     2655
 exclusivo     1246 198 6.87      776     1717

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

11- Total - Residual biomass (kg ha-1)

mod11 = lmer(resi.total.ha~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod11))

shapiro.test(resid(mod11))

    Shapiro-Wilk normality test

data:  resid(mod11)
W = 0.73359, p-value = 6.617e-12

bartlett.test(resid(mod11)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod11) by tratamentos
Bartlett's K-squared = 28.357, df = 1, p-value = 1.009e-07

mod11.1 = lmer(resi.total.ha^0.4~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod11.1))

shapiro.test(resid(mod11.1)) 

    Shapiro-Wilk normality test

data:  resid(mod11.1)
W = 0.98258, p-value = 0.2329

bartlett.test(resid(mod11.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod11.1) by tratamentos
Bartlett's K-squared = 0.019999, df = 1, p-value = 0.8875

anova(mod11.1)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)    
tratamentos 434.28  434.28     1    89  21.475 1.22e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias11=emmeans(mod11,~tratamentos)
summary(medias11)

 tratamentos emmean  SE   df lower.CL upper.CL
 consorcio      512 106 6.92    261.3      763
 exclusivo      183 106 6.92    -68.3      434

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

Biomass plant-part composition

12- Leaf (g kg-1)

mod12 = lmer(bra.folhas~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod12))

shapiro.test(resid(mod12))

    Shapiro-Wilk normality test

data:  resid(mod12)
W = 0.98543, p-value = 0.3693

bartlett.test(resid(mod12)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod12) by tratamentos
Bartlett's K-squared = 0.49557, df = 1, p-value = 0.4815

anova(mod12)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)   
tratamentos  78769   78769     1    89  7.8581 0.00621 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias12=emmeans(mod12,~tratamentos)
summary(medias12)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      657 70.2 5.22      479      836
 exclusivo      715 70.2 5.22      536      893

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

13- Stem (g kg-1)

mod13 = lmer(bra.hastes~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod13))

shapiro.test(resid(mod13))

    Shapiro-Wilk normality test

data:  resid(mod13)
W = 0.91838, p-value = 1.679e-05

bartlett.test(resid(mod13)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod13) by tratamentos
Bartlett's K-squared = 4.7876, df = 1, p-value = 0.02867

mod13.1 = lmer(bra.hastes^0.3~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod13.1))

shapiro.test(resid(mod13.1)) 

    Shapiro-Wilk normality test

data:  resid(mod13.1)
W = 0.97675, p-value = 0.08566

bartlett.test(resid(mod13.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod13.1) by tratamentos
Bartlett's K-squared = 0.24263, df = 1, p-value = 0.6223

anova(mod13.1)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
tratamentos 3.2695  3.2695     1    89  9.1341 0.003276 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias13=emmeans(mod13,~tratamentos)
summary(medias13)

 tratamentos emmean SE   df lower.CL upper.CL
 consorcio      160 26 5.63     95.0      224
 exclusivo      123 26 5.63     58.1      187

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

14- Leaf:stem ratio

data1$leaf.stem.ratio=data1$bra.folhas/data1$bra.hastes

mod14 = lmer(leaf.stem.ratio~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod14))

shapiro.test(resid(mod14))

    Shapiro-Wilk normality test

data:  resid(mod14)
W = 0.90952, p-value = 6.048e-06

bartlett.test(resid(mod14)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod14) by tratamentos
Bartlett's K-squared = 14.138, df = 1, p-value = 0.0001699

mod14.1 = lmer(leaf.stem.ratio^0.3~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod14.1))

shapiro.test(resid(mod14.1)) 

    Shapiro-Wilk normality test

data:  resid(mod14.1)
W = 0.98914, p-value = 0.6241

bartlett.test(resid(mod14.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod14.1) by tratamentos
Bartlett's K-squared = 2.095, df = 1, p-value = 0.1478

anova(mod14.1)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
tratamentos 1.1427  1.1427     1    89  11.244 0.001175 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias14=emmeans(mod14,~tratamentos)
summary(medias14)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio     6.45 1.99 6.01     1.57     11.3
 exclusivo    10.08 1.99 6.01     5.20     14.9

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

15- Dead material (g kg-1)

mod15 = lmer(bra.morto~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod15))

shapiro.test(resid(mod15))

    Shapiro-Wilk normality test

data:  resid(mod15)
W = 0.92972, p-value = 6.769e-05

bartlett.test(resid(mod15)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod15) by tratamentos
Bartlett's K-squared = 3.6241, df = 1, p-value = 0.05695

mod15.1 = lmer(bra.morto^0.6~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod15.1))

shapiro.test(resid(mod15.1))

    Shapiro-Wilk normality test

data:  resid(mod15.1)
W = 0.97777, p-value = 0.1023

bartlett.test(resid(mod15.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod15.1) by tratamentos
Bartlett's K-squared = 1.7354, df = 1, p-value = 0.1877

anova(mod15.1)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
tratamentos 45.609  45.609     1    89  1.2948 0.2582

medias15=emmeans(mod15,~tratamentos)
summary(medias15)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      183 50.1 5.35     56.6      309
 exclusivo      163 50.1 5.35     36.3      289

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

16- Tithonia - Leaf (g kg-1)

mod16 = lmer(thi.folhas~tratamentos+(1|Ciclo), data = data1)

medias16=emmeans(mod16,~tratamentos)
summary(medias16)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      489 21.8 6.58    437.1    541.8
 exclusivo        0 21.8 6.58    -52.3     52.3

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

17- Tithonia - Stem (g kg-1)

mod17 = lmer(thi.hastes~tratamentos+(1|Ciclo), data = data1)

medias17=emmeans(mod17,~tratamentos)
summary(medias17)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      454 21.8 6.86    402.4    505.9
 exclusivo        0 21.8 6.86    -51.8     51.8

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

18- Tithonia - Leaf:stem ratio

data1$LSratio.thi=data1$thi.folhas/data1$thi.hastes
data1$LSratio.thi = replace(data1$LSratio.thi, is.na(data1$LSratio.thi), 0)

mod18 = lmer(LSratio.thi~tratamentos+(1|Ciclo), data = data1)

medias18=emmeans(mod18,~tratamentos)
summary(medias18)

 tratamentos emmean    SE   df lower.CL upper.CL
 consorcio     1.26 0.116 7.02    0.990    1.539
 exclusivo     0.00 0.116 7.02   -0.275    0.275

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

19- Tithonia - Dead material (g kg-1)

mod19 = lmer(thi.morto~tratamentos+(1|Ciclo), data = data1)

medias19=emmeans(mod19,~tratamentos)
summary(medias19)

 tratamentos emmean  SE   df lower.CL upper.CL
 consorcio     56.4 4.2 13.6    47.39    65.46
 exclusivo      0.0 4.2 13.6    -9.04     9.04

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

20- Total - Leaf accumulation (kg ha-1)

mod20 = lmer(total.leaf.acu~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod20))

shapiro.test(resid(mod20))

    Shapiro-Wilk normality test

data:  resid(mod20)
W = 0.96454, p-value = 0.01065

bartlett.test(resid(mod20)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod20) by tratamentos
Bartlett's K-squared = 0.58922, df = 1, p-value = 0.4427

mod20.1 = lmer(total.leaf.acu^0.8~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod20.1))

shapiro.test(resid(mod20.1))

    Shapiro-Wilk normality test

data:  resid(mod20.1)
W = 0.97559, p-value = 0.06998

bartlett.test(resid(mod20.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod20.1) by tratamentos
Bartlett's K-squared = 0.13441, df = 1, p-value = 0.7139

anova(mod20.1)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
tratamentos  92407   92407     1    89  19.817 2.462e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias20=emmeans(mod20,~tratamentos)
summary(medias20)

 tratamentos emmean  SE   df lower.CL upper.CL
 consorcio     1207 157 5.56      816     1598
 exclusivo      889 157 5.56      498     1280

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

21- total - Stem accumulation (kg ha-1)

mod21 = lmer(total.stem.acu~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod21))

shapiro.test(resid(mod21))

    Shapiro-Wilk normality test

data:  resid(mod21)
W = 0.75476, p-value = 2.309e-11

bartlett.test(resid(mod21)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod21) by tratamentos
Bartlett's K-squared = 47.572, df = 1, p-value = 5.301e-12

mod21.1 = lmer(total.stem.acu^0.3~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod21.1))

shapiro.test(resid(mod21.1))

    Shapiro-Wilk normality test

data:  resid(mod21.1)
W = 0.97661, p-value = 0.08353

bartlett.test(resid(mod21.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod21.1) by tratamentos
Bartlett's K-squared = 3.3034, df = 1, p-value = 0.06914

anova(mod21.1)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
tratamentos 142.71  142.71     1    89  81.843 3.006e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias21=emmeans(mod21,~tratamentos)
summary(medias21)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      716 99.5 7.91    486.5      946
 exclusivo      162 99.5 7.91    -68.1      392

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

22- total - Dead material accumulation (kg ha-1)

mod22 = lmer(total.dead.acu~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod22))

shapiro.test(resid(mod22))

    Shapiro-Wilk normality test

data:  resid(mod22)
W = 0.8398, p-value = 8.489e-09

bartlett.test(resid(mod22)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod22) by tratamentos
Bartlett's K-squared = 2.7119, df = 1, p-value = 0.0996

mod22.1 = lmer(total.dead.acu^0.3~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod22.1))

shapiro.test(resid(mod22.1))

    Shapiro-Wilk normality test

data:  resid(mod22.1)
W = 0.97945, p-value = 0.1368

bartlett.test(resid(mod22.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod22.1) by tratamentos
Bartlett's K-squared = 1.9703, df = 1, p-value = 0.1604

anova(mod22.1)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
tratamentos 7.7089  7.7089     1    89  12.672 0.0005981 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias22=emmeans(mod22,~tratamentos)
summary(medias22)

 tratamentos emmean   SE  df lower.CL upper.CL
 consorcio      261 57.5 5.7    118.7      404
 exclusivo      195 57.5 5.7     52.9      338

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

Nutrive value

23- DM (g kg-1)

mod23 = lmer(DM.bra~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod23))

shapiro.test(resid(mod23))

    Shapiro-Wilk normality test

data:  resid(mod23)
W = 0.99278, p-value = 0.889

bartlett.test(resid(mod23)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod23) by tratamentos
Bartlett's K-squared = 0.19701, df = 1, p-value = 0.6571

anova(mod23)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
tratamentos 1.4925  1.4925     1    89  0.1023 0.7498

medias23=emmeans(mod23,~tratamentos)
summary(medias23)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      923 8.17 5.02      902      944
 exclusivo      923 8.17 5.02      902      944

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

24- Ash (g kg-1)

mod24 = lmer(MM.bra~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod24))

shapiro.test(resid(mod24))

    Shapiro-Wilk normality test

data:  resid(mod24)
W = 0.9799, p-value = 0.1478

bartlett.test(resid(mod24)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod24) by tratamentos
Bartlett's K-squared = 0.17572, df = 1, p-value = 0.6751

anova(mod24)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
tratamentos 320.36  320.36     1    89  3.2179 0.07623 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias24=emmeans(mod24,~tratamentos)
summary(medias24)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      108 2.88 6.52      101      114
 exclusivo      111 2.88 6.52      104      118

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

25- CP (g kg-1)

mod25 = lmer(PB.bra~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod25))

shapiro.test(resid(mod25))

    Shapiro-Wilk normality test

data:  resid(mod25)
W = 0.96829, p-value = 0.01994

bartlett.test(resid(mod25)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod25) by tratamentos
Bartlett's K-squared = 4.871, df = 1, p-value = 0.02731

#Transformation
mod25.1 = lmer(PB.bra^0.7~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod25.1))

shapiro.test(resid(mod25.1))

    Shapiro-Wilk normality test

data:  resid(mod25.1)
W = 0.97815, p-value = 0.1092

bartlett.test(resid(mod25.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod25.1) by tratamentos
Bartlett's K-squared = 3.6096, df = 1, p-value = 0.05745

anova(mod25.1)

Type III Analysis of Variance Table with Satterthwaite's method
             Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
tratamentos 0.23864 0.23864     1    89  0.0244 0.8762

medias25=emmeans(mod25,~tratamentos)
summary(medias25)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      112 8.33 5.54     90.9      133
 exclusivo      113 8.33 5.54     91.8      133

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

26- aNDFom (g kg-1)

mod26 = lmer(aFDNom.bra~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod26))

shapiro.test(resid(mod26))

    Shapiro-Wilk normality test

data:  resid(mod26)
W = 0.97954, p-value = 0.1388

bartlett.test(resid(mod26)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod26) by tratamentos
Bartlett's K-squared = 3.1871, df = 1, p-value = 0.07422

anova(mod26)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
tratamentos 7866.3  7866.3     1    89  9.1269 0.003288 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias26=emmeans(mod26,~tratamentos)
summary(medias26)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      596 8.06 6.72      577      616
 exclusivo      578 8.06 6.72      559      597

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

27- ADFom (g kg-1)

mod27 = lmer(aFDA.bra~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod27))

shapiro.test(resid(mod27))

    Shapiro-Wilk normality test

data:  resid(mod27)
W = 0.96655, p-value = 0.01488

bartlett.test(resid(mod27)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod27) by tratamentos
Bartlett's K-squared = 3.6496, df = 1, p-value = 0.05608

mod27.1 = lmer(aFDA.bra^0.5~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod27.1))

shapiro.test(resid(mod27.1))

    Shapiro-Wilk normality test

data:  resid(mod27.1)
W = 0.97569, p-value = 0.07116

bartlett.test(resid(mod27.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod27.1) by tratamentos
Bartlett's K-squared = 2.7552, df = 1, p-value = 0.09694

anova(mod27.1)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
tratamentos 2.4852  2.4852     1    89  6.4005 0.01317 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias27=emmeans(mod27,~tratamentos)
summary(medias27)

 tratamentos emmean   SE  df lower.CL upper.CL
 consorcio      302 8.38 5.8      282      323
 exclusivo      291 8.38 5.8      270      311

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

28- Hemicelulose (g kg-1)

mod28 = lmer(HEMC.bra~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod28))

shapiro.test(resid(mod28))

    Shapiro-Wilk normality test

data:  resid(mod28)
W = 0.98437, p-value = 0.3126

bartlett.test(resid(mod28)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod28) by tratamentos
Bartlett's K-squared = 2.0862, df = 1, p-value = 0.1486

anova(mod28)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
tratamentos 1053.4  1053.4     1    89  4.3565 0.03973 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias28=emmeans(mod28,~tratamentos)
summary(medias28)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      294 6.37 5.68      278      310
 exclusivo      288 6.37 5.68      272      303

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

29- Tithonia - DM (g kg-1)

mod29 = lmer(DM.thi~tratamentos+(1|Ciclo), data = data1)

medias29=emmeans(mod29,~tratamentos)
summary(medias29)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      912 5.77 6.11    898.1    926.3
 exclusivo        0 5.77 6.11    -14.1     14.1

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

30- Tithonia - Ash (g kg-1)

mod30 = lmer(MM.thi~tratamentos+(1|Ciclo), data = data1)

medias30=emmeans(mod30,~tratamentos)
summary(medias30)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      136 5.68 6.71    122.4    149.5
 exclusivo        0 5.68 6.71    -13.5     13.5

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

31- Tithonia - CP (g kg-1)

mod31 = lmer(PB.thi~tratamentos+(1|Ciclo), data = data1)

medias31=emmeans(mod31,~tratamentos)
summary(medias31)

 tratamentos emmean   SE  df lower.CL upper.CL
 consorcio      169 10.3 6.2      144      194
 exclusivo        0 10.3 6.2      -25       25

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

32- Tithonia - aNDFom (g kg-1)

mod32 = lmer(aFDNom.thi~tratamentos+(1|Ciclo), data = data1)

medias32=emmeans(mod32,~tratamentos)
summary(medias32)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      379 7.68 7.87    361.2    396.7
 exclusivo        0 7.68 7.87    -17.8     17.8

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

33- Tithonia - ADFom (g kg-1)

mod33 = lmer(aFDA.thi~tratamentos+(1|Ciclo), data = data1)

medias33=emmeans(mod33,~tratamentos)
summary(medias33)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      271 6.99 7.67    254.5    287.0
 exclusivo        0 6.99 7.67    -16.2     16.2

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

34- Tithonia - Hemicellulose (g kg-1)

mod34 = lmer(HEMC.thi~tratamentos+(1|Ciclo), data = data1)

medias34=emmeans(mod34,~tratamentos)
summary(medias34)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      108 2.93 10.6   101.75   114.71
 exclusivo        0 2.93 10.6    -6.48     6.48

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

35- total - DM (kg ha-1)

mod35 = lmer(total.DM.ac~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod35))

shapiro.test(resid(mod35))

    Shapiro-Wilk normality test

data:  resid(mod35)
W = 0.86943, p-value = 1.066e-07

bartlett.test(resid(mod35)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod35) by tratamentos
Bartlett's K-squared = 16.128, df = 1, p-value = 5.919e-05

mod35.1 = lmer(total.DM.ac^0.3~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod35.1))

shapiro.test(resid(mod35.1))

    Shapiro-Wilk normality test

data:  resid(mod35.1)
W = 0.97679, p-value = 0.08616

bartlett.test(resid(mod35.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod35.1) by tratamentos
Bartlett's K-squared = 0.73388, df = 1, p-value = 0.3916

anova(mod35.1)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
tratamentos 48.223  48.223     1    89  45.392 1.542e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias35=emmeans(mod35,~tratamentos)
summary(medias35)

 tratamentos emmean  SE   df lower.CL upper.CL
 consorcio     2020 193 6.67     1559     2482
 exclusivo     1152 193 6.67      690     1614

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

36- total - Ash (kg ha-1)

mod36 = lmer(total.MM.ac~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod36))

shapiro.test(resid(mod36))

    Shapiro-Wilk normality test

data:  resid(mod36)
W = 0.85716, p-value = 3.594e-08

bartlett.test(resid(mod36)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod36) by tratamentos
Bartlett's K-squared = 17.279, df = 1, p-value = 3.228e-05

mod36.1 = lmer(total.MM.ac^0.3~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod36.1))

shapiro.test(resid(mod36.1))

    Shapiro-Wilk normality test

data:  resid(mod36.1)
W = 0.98393, p-value = 0.291

bartlett.test(resid(mod36.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod36.1) by tratamentos
Bartlett's K-squared = 0.46817, df = 1, p-value = 0.4938

anova(mod36.1)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
tratamentos 15.076  15.076     1    89   54.56 7.736e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias36=emmeans(mod36,~tratamentos)
summary(medias36)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      247 22.1 6.72    194.6      300
 exclusivo      139 22.1 6.72     85.9      191

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

37- total - CP (kg ha-1)

mod37 = lmer(total.PB.ac~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod37))

shapiro.test(resid(mod37))

    Shapiro-Wilk normality test

data:  resid(mod37)
W = 0.89023, p-value = 7.795e-07

bartlett.test(resid(mod37)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod37) by tratamentos
Bartlett's K-squared = 12.533, df = 1, p-value = 0.0003998

mod37.1 = lmer(total.PB.ac^0.3~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod37.1))

shapiro.test(resid(mod37.1))

    Shapiro-Wilk normality test

data:  resid(mod37.1)
W = 0.98358, p-value = 0.2746

bartlett.test(resid(mod37.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod37.1) by tratamentos
Bartlett's K-squared = 0.30755, df = 1, p-value = 0.5792

anova(mod37.1)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
tratamentos 20.825  20.825     1    89  58.156 2.525e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias37=emmeans(mod37,~tratamentos)
summary(medias37)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      269 20.7 7.71    221.2      317
 exclusivo      140 20.7 7.71     92.5      189

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

38- total - aNDFom (kg ha-1)

mod38 = lmer(total.FDN.ac~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod38))

shapiro.test(resid(mod38))

    Shapiro-Wilk normality test

data:  resid(mod38)
W = 0.91583, p-value = 1.245e-05

bartlett.test(resid(mod38)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod38) by tratamentos
Bartlett's K-squared = 9.9732, df = 1, p-value = 0.001588

mod38.1 = lmer(total.FDN.ac^0.3~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod38.1))

shapiro.test(resid(mod38.1))

    Shapiro-Wilk normality test

data:  resid(mod38.1)
W = 0.99205, p-value = 0.8428

bartlett.test(resid(mod38.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod38.1) by tratamentos
Bartlett's K-squared = 0.16163, df = 1, p-value = 0.6877

anova(mod38.1)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
tratamentos 22.581  22.581     1    89   30.65 3.087e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias38=emmeans(mod38,~tratamentos)
summary(medias38)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio     1121 93.7 7.14      900     1342
 exclusivo      724 93.7 7.14      503      945

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

39- total - ADFom (kg ha-1)

mod39 = lmer(total.FDA.ac~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod39))

shapiro.test(resid(mod39))

    Shapiro-Wilk normality test

data:  resid(mod39)
W = 0.90909, p-value = 5.766e-06

bartlett.test(resid(mod39)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod39) by tratamentos
Bartlett's K-squared = 17.591, df = 1, p-value = 2.738e-05

mod39.1 = lmer(total.FDA.ac^0.3~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod39.1))

shapiro.test(resid(mod39.1))

    Shapiro-Wilk normality test

data:  resid(mod39.1)
W = 0.98804, p-value = 0.5414

bartlett.test(resid(mod39.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod39.1) by tratamentos
Bartlett's K-squared = 1.2356, df = 1, p-value = 0.2663

anova(mod39.1)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
tratamentos 26.376  26.376     1    89  43.463 2.974e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias39=emmeans(mod39,~tratamentos)
summary(medias39)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      650 55.9 7.29      519      782
 exclusivo      365 55.9 7.29      233      496

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

40- total - Hemicellulose (kg ha-1)

mod40 = lmer(total.HEMC.ac~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod40))

shapiro.test(resid(mod40))

    Shapiro-Wilk normality test

data:  resid(mod40)
W = 0.91058, p-value = 6.816e-06

bartlett.test(resid(mod40)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod40) by tratamentos
Bartlett's K-squared = 2.8813, df = 1, p-value = 0.08961

mod40.1 = lmer(total.HEMC.ac^0.3~tratamentos+(1|Ciclo), data = data1)
hist(resid(mod40.1))

shapiro.test(resid(mod40.1))

    Shapiro-Wilk normality test

data:  resid(mod40.1)
W = 0.98857, p-value = 0.5804

bartlett.test(resid(mod40.1)~tratamentos, data=data1)

    Bartlett test of homogeneity of variances

data:  resid(mod40.1) by tratamentos
Bartlett's K-squared = 0.17538, df = 1, p-value = 0.6754

anova(mod40.1)

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
tratamentos 5.7717  5.7717     1    89  14.477 0.0002598 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

medias40=emmeans(mod40,~tratamentos)
summary(medias40)

 tratamentos emmean   SE   df lower.CL upper.CL
 consorcio      471 38.8 7.15      380      562
 exclusivo      359 38.8 7.15      268      451

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

Animal parameters

41- Stocking rate (UA ha-1)

mod41 = aov(ua.ha~tratamentos*Ciclo, data = data1)
anova(mod41)

Analysis of Variance Table

Response: ua.ha
                  Df Sum Sq Mean Sq F value  Pr(>F)    
tratamentos        1  0.524  0.5236  3.4455 0.06693 .  
Ciclo              5 57.582 11.5165 75.7800 < 2e-16 ***
tratamentos:Ciclo  5  0.383  0.0765  0.5037 0.77267    
Residuals         84 12.766  0.1520                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

tukey = TukeyHSD(mod41)
print(tukey)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = ua.ha ~ tratamentos * Ciclo, data = data1)

$tratamentos
                          diff       lwr        upr     p adj
exclusivo-consorcio -0.1477083 -0.305952 0.01053535 0.0669292

$Ciclo
                     diff         lwr         upr     p adj
Cycle 2-Cycle 1  0.104375 -0.29760651  0.50635651 0.9737791
Cycle 3-Cycle 1  0.964375  0.56239349  1.36635651 0.0000000
Cycle 4-Cycle 1  0.487500  0.08551849  0.88948151 0.0084025
Cycle 5-Cycle 1  1.819375  1.41739349  2.22135651 0.0000000
Cycle 6-Cycle 1  1.966250  1.56426849  2.36823151 0.0000000
Cycle 3-Cycle 2  0.860000  0.45801849  1.26198151 0.0000003
Cycle 4-Cycle 2  0.383125 -0.01885651  0.78510651 0.0707726
Cycle 5-Cycle 2  1.715000  1.31301849  2.11698151 0.0000000
Cycle 6-Cycle 2  1.861875  1.45989349  2.26385651 0.0000000
Cycle 4-Cycle 3 -0.476875 -0.87885651 -0.07489349 0.0106668
Cycle 5-Cycle 3  0.855000  0.45301849  1.25698151 0.0000003
Cycle 6-Cycle 3  1.001875  0.59989349  1.40385651 0.0000000
Cycle 5-Cycle 4  1.331875  0.92989349  1.73385651 0.0000000
Cycle 6-Cycle 4  1.478750  1.07676849  1.88073151 0.0000000
Cycle 6-Cycle 5  0.146875 -0.25510651  0.54885651 0.8935307

$`tratamentos:Ciclo`
                                        diff         lwr         upr     p adj
exclusivo:Cycle 1-consorcio:Cycle 1 -0.22875 -0.88406024  0.42656024 0.9896282
consorcio:Cycle 2-consorcio:Cycle 1 -0.04875 -0.70406024  0.60656024 1.0000000
exclusivo:Cycle 2-consorcio:Cycle 1  0.02875 -0.62656024  0.68406024 1.0000000
consorcio:Cycle 3-consorcio:Cycle 1  0.89750  0.24218976  1.55281024 0.0008457
exclusivo:Cycle 3-consorcio:Cycle 1  0.80250  0.14718976  1.45781024 0.0047978
consorcio:Cycle 4-consorcio:Cycle 1  0.54125 -0.11406024  1.19656024 0.2091416
exclusivo:Cycle 4-consorcio:Cycle 1  0.20500 -0.45031024  0.86031024 0.9958452
consorcio:Cycle 5-consorcio:Cycle 1  1.77875  1.12343976  2.43406024 0.0000000
exclusivo:Cycle 5-consorcio:Cycle 1  1.63125  0.97593976  2.28656024 0.0000000
consorcio:Cycle 6-consorcio:Cycle 1  1.93000  1.27468976  2.58531024 0.0000000
exclusivo:Cycle 6-consorcio:Cycle 1  1.77375  1.11843976  2.42906024 0.0000000
consorcio:Cycle 2-exclusivo:Cycle 1  0.18000 -0.47531024  0.83531024 0.9986869
exclusivo:Cycle 2-exclusivo:Cycle 1  0.25750 -0.39781024  0.91281024 0.9740006
consorcio:Cycle 3-exclusivo:Cycle 1  1.12625  0.47093976  1.78156024 0.0000079
exclusivo:Cycle 3-exclusivo:Cycle 1  1.03125  0.37593976  1.68656024 0.0000590
consorcio:Cycle 4-exclusivo:Cycle 1  0.77000  0.11468976  1.42531024 0.0083773
exclusivo:Cycle 4-exclusivo:Cycle 1  0.43375 -0.22156024  1.08906024 0.5366498
consorcio:Cycle 5-exclusivo:Cycle 1  2.00750  1.35218976  2.66281024 0.0000000
exclusivo:Cycle 5-exclusivo:Cycle 1  1.86000  1.20468976  2.51531024 0.0000000
consorcio:Cycle 6-exclusivo:Cycle 1  2.15875  1.50343976  2.81406024 0.0000000
exclusivo:Cycle 6-exclusivo:Cycle 1  2.00250  1.34718976  2.65781024 0.0000000
exclusivo:Cycle 2-consorcio:Cycle 2  0.07750 -0.57781024  0.73281024 0.9999997
consorcio:Cycle 3-consorcio:Cycle 2  0.94625  0.29093976  1.60156024 0.0003289
exclusivo:Cycle 3-consorcio:Cycle 2  0.85125  0.19593976  1.50656024 0.0020057
consorcio:Cycle 4-consorcio:Cycle 2  0.59000 -0.06531024  1.24531024 0.1190902
exclusivo:Cycle 4-consorcio:Cycle 2  0.25375 -0.40156024  0.90906024 0.9766989
consorcio:Cycle 5-consorcio:Cycle 2  1.82750  1.17218976  2.48281024 0.0000000
exclusivo:Cycle 5-consorcio:Cycle 2  1.68000  1.02468976  2.33531024 0.0000000
consorcio:Cycle 6-consorcio:Cycle 2  1.97875  1.32343976  2.63406024 0.0000000
exclusivo:Cycle 6-consorcio:Cycle 2  1.82250  1.16718976  2.47781024 0.0000000
consorcio:Cycle 3-exclusivo:Cycle 2  0.86875  0.21343976  1.52406024 0.0014523
exclusivo:Cycle 3-exclusivo:Cycle 2  0.77375  0.11843976  1.42906024 0.0078637
consorcio:Cycle 4-exclusivo:Cycle 2  0.51250 -0.14281024  1.16781024 0.2807227
exclusivo:Cycle 4-exclusivo:Cycle 2  0.17625 -0.47906024  0.83156024 0.9989167
consorcio:Cycle 5-exclusivo:Cycle 2  1.75000  1.09468976  2.40531024 0.0000000
exclusivo:Cycle 5-exclusivo:Cycle 2  1.60250  0.94718976  2.25781024 0.0000000
consorcio:Cycle 6-exclusivo:Cycle 2  1.90125  1.24593976  2.55656024 0.0000000
exclusivo:Cycle 6-exclusivo:Cycle 2  1.74500  1.08968976  2.40031024 0.0000000
exclusivo:Cycle 3-consorcio:Cycle 3 -0.09500 -0.75031024  0.56031024 0.9999976
consorcio:Cycle 4-consorcio:Cycle 3 -0.35625 -1.01156024  0.29906024 0.7981386
exclusivo:Cycle 4-consorcio:Cycle 3 -0.69250 -1.34781024 -0.03718976 0.0289731
consorcio:Cycle 5-consorcio:Cycle 3  0.88125  0.22593976  1.53656024 0.0011498
exclusivo:Cycle 5-consorcio:Cycle 3  0.73375  0.07843976  1.38906024 0.0152150
consorcio:Cycle 6-consorcio:Cycle 3  1.03250  0.37718976  1.68781024 0.0000575
exclusivo:Cycle 6-consorcio:Cycle 3  0.87625  0.22093976  1.53156024 0.0012628
consorcio:Cycle 4-exclusivo:Cycle 3 -0.26125 -0.91656024  0.39406024 0.9710724
exclusivo:Cycle 4-exclusivo:Cycle 3 -0.59750 -1.25281024  0.05781024 0.1084764
consorcio:Cycle 5-exclusivo:Cycle 3  0.97625  0.32093976  1.63156024 0.0001811
exclusivo:Cycle 5-exclusivo:Cycle 3  0.82875  0.17343976  1.48406024 0.0030152
consorcio:Cycle 6-exclusivo:Cycle 3  1.12750  0.47218976  1.78281024 0.0000077
exclusivo:Cycle 6-exclusivo:Cycle 3  0.97125  0.31593976  1.62656024 0.0002002
exclusivo:Cycle 4-consorcio:Cycle 4 -0.33625 -0.99156024  0.31906024 0.8514723
consorcio:Cycle 5-consorcio:Cycle 4  1.23750  0.58218976  1.89281024 0.0000007
exclusivo:Cycle 5-consorcio:Cycle 4  1.09000  0.43468976  1.74531024 0.0000172
consorcio:Cycle 6-consorcio:Cycle 4  1.38875  0.73343976  2.04406024 0.0000000
exclusivo:Cycle 6-consorcio:Cycle 4  1.23250  0.57718976  1.88781024 0.0000008
consorcio:Cycle 5-exclusivo:Cycle 4  1.57375  0.91843976  2.22906024 0.0000000
exclusivo:Cycle 5-exclusivo:Cycle 4  1.42625  0.77093976  2.08156024 0.0000000
consorcio:Cycle 6-exclusivo:Cycle 4  1.72500  1.06968976  2.38031024 0.0000000
exclusivo:Cycle 6-exclusivo:Cycle 4  1.56875  0.91343976  2.22406024 0.0000000
exclusivo:Cycle 5-consorcio:Cycle 5 -0.14750 -0.80281024  0.50781024 0.9997985
consorcio:Cycle 6-consorcio:Cycle 5  0.15125 -0.50406024  0.80656024 0.9997432
exclusivo:Cycle 6-consorcio:Cycle 5 -0.00500 -0.66031024  0.65031024 1.0000000
consorcio:Cycle 6-exclusivo:Cycle 5  0.29875 -0.35656024  0.95406024 0.9269758
exclusivo:Cycle 6-exclusivo:Cycle 5  0.14250 -0.51281024  0.79781024 0.9998559
exclusivo:Cycle 6-consorcio:Cycle 6 -0.15625 -0.81156024  0.49906024 0.9996493

tukey.cld = multcompLetters4(mod41, tukey)
print(tukey.cld)

$tratamentos
$tratamentos$Letters
consorcio exclusivo 
      "a"       "a" 

$tratamentos$LetterMatrix
             a
consorcio TRUE
exclusivo TRUE


$Ciclo
Cycle 6 Cycle 5 Cycle 3 Cycle 4 Cycle 2 Cycle 1 
    "a"     "a"     "b"     "c"    "cd"     "d" 

$`tratamentos:Ciclo`
consorcio:Cycle 6 consorcio:Cycle 5 exclusivo:Cycle 6 exclusivo:Cycle 5 consorcio:Cycle 3 exclusivo:Cycle 3 
              "a"               "a"               "a"               "a"               "b"              "bc" 
consorcio:Cycle 4 exclusivo:Cycle 4 exclusivo:Cycle 2 consorcio:Cycle 1 consorcio:Cycle 2 exclusivo:Cycle 1 
            "bcd"             "cde"              "de"              "de"              "de"               "e" 

data_summary = group_by(data1, tratamentos, Ciclo) %>%
  summarise(data=mean(as.Date(entrada)),
            mean=mean(ua.ha), 
            herbage=mean(biomass.bra.ha), 
            sd=sd(ua.ha)) %>%
  arrange(desc(mean))

`summarise()` has grouped output by 'tratamentos'. You can override using the `.groups` argument.

print(data_summary)

scale = max(data_summary$herbage)/max(data_summary$mean)
scale

[1] 416.2487

plot1=ggplot(data_summary, aes(x =interaction(Ciclo, data), y=mean))+ 
  scale_x_discrete(NULL,guide = "axis_nested")+
  geom_bar(aes(y=mean, fill=tratamentos), stat = "identity", position = "dodge", alpha = 0.5)+
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd, group= tratamentos), size=0.3, position = position_dodge(0.9), width = 0.25)+
  geom_line(aes(y=herbage/(scale-150),group= tratamentos), size=0.3, linetype =1)+
  coord_cartesian(ylim = c(0,8))+
   geom_point(aes(y=herbage/(scale-150), shape = tratamentos), size = 3)+
  scale_y_continuous(name= "Stocking rate (UA/ha)",breaks = seq(0, 8, 1),   sec.axis=sec_axis(~.*(scale-150), name="Grass HB (kg/ha)",    breaks=seq(0,2200,150)))+
scale_shape_manual(values = c(1, 15), name= "Grass HB:",labels = c("SPS", "EP"))+
scale_fill_manual(name = "Stocking rate:",values=c("gray", "gray28"),labels=c("SPS","EP"))
plot1

plot2=plot1+
  theme(
  panel.background = element_rect(fill = "transparent",colour = "black",size = 0.3),
  legend.key=element_blank(),
  legend.background = element_rect(fill = "transparent", size=0.3, linetype="solid",colour ="black"),
  legend.position = c(0.25, 0.85), legend.direction="horizontal",
  legend.key.size = unit(0.5,"cm"), 
  legend.key.height = unit(0.5,"cm"),  
  legend.key.width = unit(0.5,"cm"), 
  legend.title = element_text(size=10), 
  legend.text = element_text(size=9),
  axis.title.y = element_text(size = 12),
  axis.title.x = element_text(size = 12),
  axis.text.x = element_text(size = 10),
  axis.text.y = element_text(size = 10),
  axis.line = element_line(colour = "black", size = 0.3,  linetype = "solid"))+
annotate("text", size=4, x=5.5, y=8,label="Stocking rate P values")+
annotate("text",size=3.5, x=5.5, y=7.5,label= "    Treat = 0.066")+
annotate("text",size=3.5, x=5.5, y=7,  label= "    Cycles = <.001")+
annotate("text",size=3.5, x=5.5, y=6.5,label= "      T*C = 0.773")+
annotate("text", size=5, x=1, y=2,label="d")+
annotate("text",size=5, x=2, y=2,label= "cd")+
annotate("text",size=5, x=3, y=2.9,  label= "b")+
annotate("text",size=5, x=4, y=2.5,label= "c")+
annotate("text", size=5, x=5, y=3.8,label="a")+
annotate("text",size=5, x=6, y=3.8,label= "a")
plot2

save_plot("Stocking rate.pdf", plot2, ncol = 1, nrow = 1) 

42 - Herbage biomass by cycles

mod42 = aov(biomass.bra.ha~tratamentos*Ciclo, data = data1)
anova(mod42)

Analysis of Variance Table

Response: biomass.bra.ha
                  Df   Sum Sq Mean Sq F value Pr(>F)  
tratamentos        1      617     617  0.0034 0.9534  
Ciclo              5  2301019  460204  2.5661 0.0328 *
tratamentos:Ciclo  5   114143   22829  0.1273 0.9858  
Residuals         84 15064509  179339                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

tukey1 = TukeyHSD(mod42)
print(tukey1)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = biomass.bra.ha ~ tratamentos * Ciclo, data = data1)

$tratamentos
                         diff       lwr      upr     p adj
exclusivo-consorcio -5.068958 -176.9713 166.8334 0.9533789

$Ciclo
                      diff        lwr      upr     p adj
Cycle 2-Cycle 1 -329.17250 -765.85065 107.5056 0.2494602
Cycle 3-Cycle 1 -169.46688 -606.14502 267.2113 0.8667283
Cycle 4-Cycle 1 -188.31125 -624.98940 248.3669 0.8067805
Cycle 5-Cycle 1   82.54812 -354.13002 519.2263 0.9937527
Cycle 6-Cycle 1   91.32312 -345.35502 528.0013 0.9900419
Cycle 3-Cycle 2  159.70562 -276.97252 596.3838 0.8931425
Cycle 4-Cycle 2  140.86125 -295.81690 577.5394 0.9346594
Cycle 5-Cycle 2  411.72062  -24.95752 848.3988 0.0761673
Cycle 6-Cycle 2  420.49562  -16.18252 857.1738 0.0658913
Cycle 4-Cycle 3  -18.84438 -455.52252 417.8338 0.9999954
Cycle 5-Cycle 3  252.01500 -184.66315 688.6931 0.5466169
Cycle 6-Cycle 3  260.79000 -175.88815 697.4681 0.5085891
Cycle 5-Cycle 4  270.85938 -165.81877 707.5375 0.4656232
Cycle 6-Cycle 4  279.63438 -157.04377 716.3125 0.4290969
Cycle 6-Cycle 5    8.77500 -427.90315 445.4531 0.9999999

$`tratamentos:Ciclo`
                                          diff        lwr       upr     p adj
exclusivo:Cycle 1-consorcio:Cycle 1  -97.56000  -809.4327  614.3127 0.9999987
consorcio:Cycle 2-consorcio:Cycle 1 -370.24750 -1082.1202  341.6252 0.8400204
exclusivo:Cycle 2-consorcio:Cycle 1 -385.65750 -1097.5302  326.2152 0.8016733
consorcio:Cycle 3-consorcio:Cycle 1 -201.48125  -913.3539  510.3914 0.9982783
exclusivo:Cycle 3-consorcio:Cycle 1 -235.01250  -946.8852  476.8602 0.9934425
consorcio:Cycle 4-consorcio:Cycle 1 -274.82000  -986.6927  437.0527 0.9772229
exclusivo:Cycle 4-consorcio:Cycle 1 -199.36250  -911.2352  512.5102 0.9984347
consorcio:Cycle 5-consorcio:Cycle 1  -13.95750  -725.8302  697.9152 1.0000000
exclusivo:Cycle 5-consorcio:Cycle 1   81.49375  -630.3789  793.3664 0.9999998
consorcio:Cycle 6-consorcio:Cycle 1   69.95375  -641.9189  781.8264 1.0000000
exclusivo:Cycle 6-consorcio:Cycle 1   15.13250  -696.7402  727.0052 1.0000000
consorcio:Cycle 2-exclusivo:Cycle 1 -272.68750  -984.5602  439.1852 0.9785243
exclusivo:Cycle 2-exclusivo:Cycle 1 -288.09750  -999.9702  423.7752 0.9677217
consorcio:Cycle 3-exclusivo:Cycle 1 -103.92125  -815.7939  607.9514 0.9999975
exclusivo:Cycle 3-exclusivo:Cycle 1 -137.45250  -849.3252  574.4202 0.9999556
consorcio:Cycle 4-exclusivo:Cycle 1 -177.26000  -889.1327  534.6127 0.9994709
exclusivo:Cycle 4-exclusivo:Cycle 1 -101.80250  -813.6752  610.0702 0.9999980
consorcio:Cycle 5-exclusivo:Cycle 1   83.60250  -628.2702  795.4752 0.9999997
exclusivo:Cycle 5-exclusivo:Cycle 1  179.05375  -532.8189  890.9264 0.9994184
consorcio:Cycle 6-exclusivo:Cycle 1  167.51375  -544.3589  879.3864 0.9996908
exclusivo:Cycle 6-exclusivo:Cycle 1  112.69250  -599.1802  824.5652 0.9999941
exclusivo:Cycle 2-consorcio:Cycle 2  -15.41000  -727.2827  696.4627 1.0000000
consorcio:Cycle 3-consorcio:Cycle 2  168.76625  -543.1064  880.6389 0.9996680
exclusivo:Cycle 3-consorcio:Cycle 2  135.23500  -576.6377  847.1077 0.9999623
consorcio:Cycle 4-consorcio:Cycle 2   95.42750  -616.4452  807.3002 0.9999990
exclusivo:Cycle 4-consorcio:Cycle 2  170.88500  -540.9877  882.7577 0.9996260
consorcio:Cycle 5-consorcio:Cycle 2  356.29000  -355.5827 1068.1627 0.8709482
exclusivo:Cycle 5-consorcio:Cycle 2  451.74125  -260.1314 1163.6139 0.6009067
consorcio:Cycle 6-consorcio:Cycle 2  440.20125  -271.6714 1152.0739 0.6386875
exclusivo:Cycle 6-consorcio:Cycle 2  385.38000  -326.4927 1097.2527 0.8024005
consorcio:Cycle 3-exclusivo:Cycle 2  184.17625  -527.6964  896.0489 0.9992430
exclusivo:Cycle 3-exclusivo:Cycle 2  150.64500  -561.2277  862.5177 0.9998897
consorcio:Cycle 4-exclusivo:Cycle 2  110.83750  -601.0352  822.7102 0.9999951
exclusivo:Cycle 4-exclusivo:Cycle 2  186.29500  -525.5777  898.1677 0.9991582
consorcio:Cycle 5-exclusivo:Cycle 2  371.70000  -340.1727 1083.5727 0.8365888
exclusivo:Cycle 5-exclusivo:Cycle 2  467.15125  -244.7214 1179.0239 0.5499850
consorcio:Cycle 6-exclusivo:Cycle 2  455.61125  -256.2614 1167.4839 0.5881418
exclusivo:Cycle 6-exclusivo:Cycle 2  400.79000  -311.0827 1112.6627 0.7601566
exclusivo:Cycle 3-consorcio:Cycle 3  -33.53125  -745.4039  678.3414 1.0000000
consorcio:Cycle 4-consorcio:Cycle 3  -73.33875  -785.2114  638.5339 0.9999999
exclusivo:Cycle 4-consorcio:Cycle 3    2.11875  -709.7539  713.9914 1.0000000
consorcio:Cycle 5-consorcio:Cycle 3  187.52375  -524.3489  899.3964 0.9991055
exclusivo:Cycle 5-consorcio:Cycle 3  282.97500  -428.8977  994.8477 0.9716842
consorcio:Cycle 6-consorcio:Cycle 3  271.43500  -440.4377  983.3077 0.9792617
exclusivo:Cycle 6-consorcio:Cycle 3  216.61375  -495.2589  928.4864 0.9967304
consorcio:Cycle 4-exclusivo:Cycle 3  -39.80750  -751.6802  672.0652 1.0000000
exclusivo:Cycle 4-exclusivo:Cycle 3   35.65000  -676.2227  747.5227 1.0000000
consorcio:Cycle 5-exclusivo:Cycle 3  221.05500  -490.8177  932.9277 0.9961013
exclusivo:Cycle 5-exclusivo:Cycle 3  316.50625  -395.3664 1028.3789 0.9380047
consorcio:Cycle 6-exclusivo:Cycle 3  304.96625  -406.9064 1016.8389 0.9517543
exclusivo:Cycle 6-exclusivo:Cycle 3  250.14500  -461.7277  962.0177 0.9890594
exclusivo:Cycle 4-consorcio:Cycle 4   75.45750  -636.4152  787.3302 0.9999999
consorcio:Cycle 5-consorcio:Cycle 4  260.86250  -451.0102  972.7352 0.9847363
exclusivo:Cycle 5-consorcio:Cycle 4  356.31375  -355.5589 1068.1864 0.8708988
consorcio:Cycle 6-consorcio:Cycle 4  344.77375  -367.0989 1056.6464 0.8935897
exclusivo:Cycle 6-consorcio:Cycle 4  289.95250  -421.9202 1001.8252 0.9661897
consorcio:Cycle 5-exclusivo:Cycle 4  185.40500  -526.4677  897.2777 0.9991948
exclusivo:Cycle 5-exclusivo:Cycle 4  280.85625  -431.0164  992.7289 0.9732116
consorcio:Cycle 6-exclusivo:Cycle 4  269.31625  -442.5564  981.1889 0.9804646
exclusivo:Cycle 6-exclusivo:Cycle 4  214.49500  -497.3777  926.3677 0.9969993
exclusivo:Cycle 5-consorcio:Cycle 5   95.45125  -616.4214  807.3239 0.9999990
consorcio:Cycle 6-consorcio:Cycle 5   83.91125  -627.9614  795.7839 0.9999997
exclusivo:Cycle 6-consorcio:Cycle 5   29.09000  -682.7827  740.9627 1.0000000
consorcio:Cycle 6-exclusivo:Cycle 5  -11.54000  -723.4127  700.3327 1.0000000
exclusivo:Cycle 6-exclusivo:Cycle 5  -66.36125  -778.2339  645.5114 1.0000000
exclusivo:Cycle 6-consorcio:Cycle 6  -54.82125  -766.6939  657.0514 1.0000000

tukey.cld1 = multcompLetters4(mod42, tukey1)
print(tukey.cld1)

$tratamentos
$tratamentos$Letters
consorcio exclusivo 
      "a"       "a" 

$tratamentos$LetterMatrix
             a
consorcio TRUE
exclusivo TRUE


$Ciclo
$Ciclo$Letters
Cycle 6 Cycle 5 Cycle 1 Cycle 3 Cycle 4 Cycle 2 
    "a"     "a"     "a"     "a"     "a"     "a" 

$Ciclo$LetterMatrix
           a
Cycle 6 TRUE
Cycle 5 TRUE
Cycle 1 TRUE
Cycle 3 TRUE
Cycle 4 TRUE
Cycle 2 TRUE


$`tratamentos:Ciclo`
$`tratamentos:Ciclo`$Letters
exclusivo:Cycle 5 consorcio:Cycle 6 exclusivo:Cycle 6 consorcio:Cycle 1 consorcio:Cycle 5 exclusivo:Cycle 1 
              "a"               "a"               "a"               "a"               "a"               "a" 
exclusivo:Cycle 4 consorcio:Cycle 3 exclusivo:Cycle 3 consorcio:Cycle 4 consorcio:Cycle 2 exclusivo:Cycle 2 
              "a"               "a"               "a"               "a"               "a"               "a" 

$`tratamentos:Ciclo`$LetterMatrix
                     a
exclusivo:Cycle 5 TRUE
consorcio:Cycle 6 TRUE
exclusivo:Cycle 6 TRUE
consorcio:Cycle 1 TRUE
consorcio:Cycle 5 TRUE
exclusivo:Cycle 1 TRUE
exclusivo:Cycle 4 TRUE
consorcio:Cycle 3 TRUE
exclusivo:Cycle 3 TRUE
consorcio:Cycle 4 TRUE
consorcio:Cycle 2 TRUE
exclusivo:Cycle 2 TRUE

data.season = data1 %>%
  group_by(Ciclo, tratamentos) %>%
  summarise(
    data=mean(as.Date(entrada)),
    bra = mean(biomass.bra.ha),
    thi = mean(biomass.thi.ha),
    sd=sd(biomass.bra.ha))

`summarise()` has grouped output by 'Ciclo'. You can override using the `.groups` argument.

print(data.season)

str(data.season)

gropd_df [12 × 6] (S3: grouped_df/tbl_df/tbl/data.frame)
 $ Ciclo      : Factor w/ 6 levels "Cycle 1","Cycle 2",..: 1 1 2 2 3 3 4 4 5 5 ...
 $ tratamentos: Factor w/ 2 levels "consorcio","exclusivo": 1 2 1 2 1 2 1 2 1 2 ...
 $ data       : Date[1:12], format: "2022-01-31" "2022-01-31" "2022-02-28" "2022-02-28" ...
 $ bra        : num [1:12] 1383 1286 1013 998 1182 ...
 $ thi        : num [1:12] 455 0 149 0 580 ...
 $ sd         : num [1:12] 494 289 309 296 151 ...
 - attr(*, "groups")= tibble [6 × 2] (S3: tbl_df/tbl/data.frame)
  ..$ Ciclo: Factor w/ 6 levels "Cycle 1","Cycle 2",..: 1 2 3 4 5 6
  ..$ .rows: list<int> [1:6] 
  .. ..$ : int [1:2] 1 2
  .. ..$ : int [1:2] 3 4
  .. ..$ : int [1:2] 5 6
  .. ..$ : int [1:2] 7 8
  .. ..$ : int [1:2] 9 10
  .. ..$ : int [1:2] 11 12
  .. ..@ ptype: int(0) 
  ..- attr(*, ".drop")= logi TRUE

data.season = data.season %>% mutate_if(is.numeric, ~na_if(., 0))

`mutate_if()` ignored the following grouping variables:

plot1.1=ggplot(data.season, aes(x =interaction(Ciclo, data), y=bra))+ 
  scale_x_discrete(NULL,guide = "axis_nested")+
  geom_bar(aes(y=bra, fill=tratamentos), stat = "identity", position = "dodge", alpha = 0.5)+
 # geom_errorbar(aes(ymin=bra-sd, ymax=bra+sd, group= tratamentos), size=0.3, position = position_dodge(0.9), width = 0.25)+
  geom_line(aes(y=thi, group=tratamentos), size=0.3, linetype =1)+
  geom_point(aes(y=thi, shape = tratamentos), size = 3)+
  scale_y_continuous(name= "HB production (kg/ha)",breaks = seq(0, 2600, 200))+
scale_shape_manual(values = c(1,NA), name= "Tithonia HB:",labels = c("SPS",""))+
scale_fill_manual(name = "Grass HB:",values=c("gray", "gray28"),labels=c("SPS","EP"))+
coord_cartesian(ylim = c(0,2600))
plot1.1

plot2.1=plot1.1+
  theme(
  panel.background = element_rect(fill = "transparent",colour = "black",size = 0.3),
  legend.key=element_blank(),
  legend.background = element_rect(fill = "transparent", size=0.3, linetype="solid",colour ="black"),
  legend.position = c(0.25, 0.85), legend.direction="horizontal",
  legend.key.size = unit(0.5,"cm"), 
  legend.key.height = unit(0.5,"cm"),  
  legend.key.width = unit(0.5,"cm"), 
  legend.title = element_text(size=10), 
  legend.text = element_text(size=9),
  axis.title.y = element_text(size = 12),
  axis.title.x = element_text(size = 12),
  axis.text.x = element_text(size = 10),
  axis.text.y = element_text(size = 10),
  axis.line = element_line(colour = "black", size = 0.3,  linetype = "solid"))+
annotate("text", size=4, x=5.5, y=2500,label="Stocking rate P values")+
annotate("text",size=3.5, x=5.5, y=2400,label= "    Treat = 0.953")+
annotate("text",size=3.5, x=5.5, y=2300,  label= "    Cycles = 0.033")+
annotate("text",size=3.5, x=5.5, y=2200,label= "      T*C = 0.986")
plot2.1

43 - Plot isotopic

d13c=data1[c(36:48,81:96),]
print(d13c)

media = aggregate(C3. ~ tratamentos, data = d13c, mean)

plot3.1 = ggplot(data=d13c, aes(x=tratamentos, y=C3., fill=tratamentos))+
  geom_boxplot(outlier.shape = NA,size=0.3)+
  coord_flip()+
  geom_text(data = media, aes(x = tratamentos, y = C3., label = sprintf("%.2f", C3.)),
            vjust = -4, color = "black", size = 5, position = position_dodge(width = 0.75))+
  scale_x_discrete(name= "Treatments",labels=c("consorcio" = "SPS","exclusivo" = "EP"))+
  scale_y_continuous(name= expression(paste("",C[3]," in feces (%)")),breaks = seq(0, 20, 1))+
  scale_fill_manual(name = "Treatments: ",values=c("gray", "gray40"),labels=c("SPS","EP"))
plot3.1

plot3.2=plot3.1+
  theme(panel.background = element_rect(fill = "transparent",colour = "black",size = 0.3),
  axis.title.y = element_text(size = 12),
  axis.title.x = element_text(size = 12),
  axis.text.x = element_text(size = 10),
  axis.text.y = element_text(size = 10),
  axis.line = element_line(colour = "black", size = 0.3,  linetype = "solid"),
  legend.key=element_blank(),
  legend.background = element_rect(fill = "transparent", size=0.3, linetype="solid",colour ="black"),
  legend.position = c(0.75, 0.9), legend.direction="horizontal",
  legend.title = element_text(size=12), 
  legend.text = element_text(size=10))
plot3.2

save_plot("isotopic.pdf", plot3.2, ncol = 1, nrow = 1)