setwd("~/Library/CloudStorage/GoogleDrive-icarounam@gmail.com/Mi unidad/Agrosavia/Colaboraciones/Fabricio/Revisión/data")
library(ggplot2)
library(emmeans)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
datos<-read.table("cafor_dbh.csv", header=T, sep=',')
datos$Shade<-as.factor(datos$Shade)
datos$Bloque<-as.factor(datos$Bloque)
datos$Shade<-as.factor(datos$Shade)
attach(datos)
#General anova
aov.diam.shoot<-aov(Mg.ha.carbon.shoot~Shade+Bloque)
aov.diam.root<-aov(Mg.ha.carbon.root~Shade+Bloque)
aov.diam.total<-aov(Mg.ha.carbon.total~Shade+Bloque)
summary(aov.diam.shoot)
## Df Sum Sq Mean Sq F value Pr(>F)
## Shade 2 3151.6 1575.8 57.644 1.34e-13 ***
## Bloque 2 54.7 27.4 1.001 0.375
## Residuals 49 1339.5 27.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov.diam.root)
## Df Sum Sq Mean Sq F value Pr(>F)
## Shade 2 0.20973 0.10487 59.041 8.87e-14 ***
## Bloque 2 0.00145 0.00073 0.409 0.666
## Residuals 49 0.08703 0.00178
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov.diam.total)
## Df Sum Sq Mean Sq F value Pr(>F)
## Shade 2 3203 1601.5 57.708 1.32e-13 ***
## Bloque 2 55 27.6 0.993 0.378
## Residuals 49 1360 27.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Mg.ha.carbon.shoot
#Shade
contrast <- emmeans(aov.diam.shoot, ~Shade)
ex <- expression(Estimated~marginal~mean(Mg~C~ha^-1))
MgC.shoot<-plot(contrast, comparisons = TRUE, xlab =ex, ylab = "Shade tree")
MgC.shoot
ggsave("MgCshoot.tiff", plot= MgC.shoot, width = 16, height = 12, units= c("cm"), dpi = 1000)
medias.Shade.shoot <- emmeans(aov.diam.shoot, pairwise ~ Shade)
medias.Shade.shoot
## $emmeans
## Shade emmean SE df lower.CL upper.CL
## C. pyriformis 6.82 1.23 49 4.35 9.3
## T. rosea 5.12 1.23 49 2.65 7.6
## T. superba 22.11 1.23 49 19.64 24.6
##
## Results are averaged over the levels of: Bloque
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## C. pyriformis - T. rosea 1.7 1.74 49 0.975 0.5959
## C. pyriformis - T. superba -15.3 1.74 49 -8.773 <.0001
## T. rosea - T. superba -17.0 1.74 49 -9.748 <.0001
##
## Results are averaged over the levels of: Bloque
## P value adjustment: tukey method for comparing a family of 3 estimates
cld_gen.shoot <-multcomp::cld(contrast, alpha = 0.05, Letters = LETTERS, reversed=T)
cld_gen.shoot
## Shade emmean SE df lower.CL upper.CL .group
## T. superba 22.11 1.23 49 19.64 24.6 A
## C. pyriformis 6.82 1.23 49 4.35 9.3 B
## T. rosea 5.12 1.23 49 2.65 7.6 B
##
## Results are averaged over the levels of: Bloque
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.
#Mg.ha.carbon.root
#Shade
contrast <- emmeans(aov.diam.root, ~Shade)
MgC.root<-plot(contrast, comparisons = TRUE, xlab = ex, ylab = "Shade tree")
MgC.root
ggsave("MgCroot.tiff", plot= MgC.root, width = 16, height = 12, units= c("cm"), dpi = 1000)
medias.Shade.root <- emmeans(aov.diam.root, pairwise ~ Shade)
medias.Shade.root
## $emmeans
## Shade emmean SE df lower.CL upper.CL
## C. pyriformis 0.200 0.00993 49 0.180 0.220
## T. rosea 0.168 0.00993 49 0.148 0.188
## T. superba 0.313 0.00993 49 0.293 0.333
##
## Results are averaged over the levels of: Bloque
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## C. pyriformis - T. rosea 0.0317 0.014 49 2.256 0.0719
## C. pyriformis - T. superba -0.1135 0.014 49 -8.078 <.0001
## T. rosea - T. superba -0.1452 0.014 49 -10.334 <.0001
##
## Results are averaged over the levels of: Bloque
## P value adjustment: tukey method for comparing a family of 3 estimates
cld_gen.root <-multcomp::cld(contrast, alpha = 0.05, Letters = LETTERS, reversed=T)
cld_gen.root
## Shade emmean SE df lower.CL upper.CL .group
## T. superba 0.313 0.00993 49 0.293 0.333 A
## C. pyriformis 0.200 0.00993 49 0.180 0.220 B
## T. rosea 0.168 0.00993 49 0.148 0.188 B
##
## Results are averaged over the levels of: Bloque
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.
#Mg.ha.carbon.total
#Shade
contrast <- emmeans(aov.diam.total, ~Shade)
MgC.total<-plot(contrast, comparisons = TRUE, xlab = ex, ylab = "Shade tree")
MgC.total
ggsave("MgCtotal.tiff", plot= MgC.total, width = 16, height = 12, units= c("cm"), dpi = 1000)
medias.Shade.total <- emmeans(aov.diam.total, pairwise ~ Shade)
medias.Shade.total
## $emmeans
## Shade emmean SE df lower.CL upper.CL
## C. pyriformis 7.02 1.24 49 4.53 9.52
## T. rosea 5.29 1.24 49 2.80 7.79
## T. superba 22.43 1.24 49 19.93 24.92
##
## Results are averaged over the levels of: Bloque
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## C. pyriformis - T. rosea 1.73 1.76 49 0.986 0.5891
## C. pyriformis - T. superba -15.40 1.76 49 -8.772 <.0001
## T. rosea - T. superba -17.13 1.76 49 -9.758 <.0001
##
## Results are averaged over the levels of: Bloque
## P value adjustment: tukey method for comparing a family of 3 estimates
cld_gen.total <-multcomp::cld(contrast, alpha = 0.05, Letters = LETTERS, reversed=T)
cld_gen.total
## Shade emmean SE df lower.CL upper.CL .group
## T. superba 22.43 1.24 49 19.93 24.92 A
## C. pyriformis 7.02 1.24 49 4.53 9.52 B
## T. rosea 5.29 1.24 49 2.80 7.79 B
##
## Results are averaged over the levels of: Bloque
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.