setwd("~/Library/CloudStorage/GoogleDrive-icarounam@gmail.com/Mi unidad/Agrosavia/Colaboraciones/Laura/Intercrop")
dat<-read.table("maiz2M.csv", header=T, sep=",")
dat$trat<-as.factor(dat$trat)
dat$bloq<-as.factor(dat$bloq)
dat$forestal<-as.factor(dat$forestal)
dat$mainplot<-as.factor(dat$mainplot)
library(lmerTest)
## Loading required package: lme4
## Loading required package: Matrix
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.1.2
##
## 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
library(agricolae)
library(emmeans)
library(lattice)
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
library(lsmeans)
## The 'lsmeans' package is now basically a front end for 'emmeans'.
## Users are encouraged to switch the rest of the way.
## See help('transition') for more information, including how to
## convert old 'lsmeans' objects and scripts to work with 'emmeans'.
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.1.2
## Metodo 1 sin bloques (usando modelos mixtos, los errores no se ven en la tabla de ANOVA pero el modelo de Split plot, ya esta incluido)
mod <- lmer(prod ~ forestal * trat + (1|mainplot),
data=dat)
mod %>%
VarCorr() %>%
as.data.frame() %>%
select(grp, vcov)
## grp vcov
## 1 mainplot 1492056
## 2 Residual 1594391
summary(mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prod ~ forestal * trat + (1 | mainplot)
## Data: dat
##
## REML criterion at convergence: 1858.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5373 -0.6589 -0.0974 0.3818 3.5275
##
## Random effects:
## Groups Name Variance Std.Dev.
## mainplot (Intercept) 1492056 1221
## Residual 1594391 1263
## Number of obs: 113, groups: mainplot, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4520.808 908.707 3.287 4.975 0.01252 *
## forestalRoble -2766.115 1293.904 3.377 -2.138 0.11198
## forestalTerminalia -2647.644 1285.105 3.287 -2.060 0.12347
## tratT3 -3025.166 410.841 104.036 -7.363 4.38e-11 ***
## forestalRoble:tratT3 1892.330 595.043 104.081 3.180 0.00194 **
## forestalTerminalia:tratT3 2601.158 572.913 104.019 4.540 1.52e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) frstlR frstlT tratT3 frR:T3
## forestalRbl -0.702
## forstlTrmnl -0.707 0.497
## tratT3 -0.214 0.150 0.151
## frstlRbl:T3 0.147 -0.236 -0.104 -0.690
## frstlTrm:T3 0.153 -0.108 -0.217 -0.717 0.495
mod %>% anova(ddf="Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## forestal 3600873 1800437 2 3.00 1.1292 0.4309
## trat 65162363 65162363 1 104.06 40.8698 4.694e-09 ***
## forestal:trat 34695914 17347957 2 104.05 10.8806 5.118e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Metodo 1 con bloques
mod2 <- lmer(prod ~ bloq + forestal * trat + (1|bloq:forestal),
data=dat)
mod2 %>%
VarCorr() %>%
as.data.frame() %>%
select(grp, vcov)
## grp vcov
## 1 bloq:forestal 1847318
## 2 Residual 1594339
summary(mod2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prod ~ bloq + forestal * trat + (1 | bloq:forestal)
## Data: dat
##
## REML criterion at convergence: 1841.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5154 -0.6635 -0.1174 0.3623 3.5347
##
## Random effects:
## Groups Name Variance Std.Dev.
## bloq:forestal (Intercept) 1847318 1359
## Residual 1594339 1263
## Number of obs: 113, groups: bloq:forestal, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4902.485 1151.345 2.124 4.258 0.04580 *
## bloqB2 -763.354 1135.288 2.008 -0.672 0.57037
## forestalRoble -2772.780 1424.613 2.213 -1.946 0.17863
## forestalTerminalia -2647.644 1416.599 2.164 -1.869 0.19288
## tratT3 -3028.248 410.849 104.027 -7.371 4.22e-11 ***
## forestalRoble:tratT3 1900.674 595.110 104.034 3.194 0.00186 **
## forestalTerminalia:tratT3 2604.240 572.915 104.018 4.546 1.48e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) bloqB2 frstlR frstlT tratT3 frR:T3
## bloqB2 -0.493
## forestalRbl -0.614 0.005
## forstlTrmnl -0.615 0.000 0.497
## tratT3 -0.172 0.007 0.136 0.137
## frstlRbl:T3 0.123 -0.014 -0.214 -0.095 -0.690
## frstlTrm:T3 0.123 -0.005 -0.098 -0.196 -0.717 0.495
mod2 %>% anova(ddf="Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## bloq 720796 720796 1 2.000 0.4521 0.5706
## forestal 2944298 1472149 2 1.981 0.0058 0.9942
## trat 65093046 65093046 1 104.042 40.8276 4.769e-09 ***
## forestal:trat 34818896 17409448 2 104.023 10.9195 4.957e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Medias para forestal
emm_forestal<-emmeans(mod2, ~ forestal, lmer.df = "kenward-roger")
## NOTE: Results may be misleading due to involvement in interactions
multcomp::cld(emm_forestal, adjust="none", Letters=LETTERS, reversed=T)
## forestal emmean SE df lower.CL upper.CL .group
## Abarco 3007 983 2.00 -1229 7242 A
## Terminalia 1661 982 1.99 -2589 5911 A
## Roble 1184 985 2.02 -3022 5391 A
##
## Results are averaged over the levels of: bloq, trat
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## 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.
## Medias para trat
emm_IS<-emmeans(mod2, ~ trat, lmer.df = "kenward-roger")
## NOTE: Results may be misleading due to involvement in interactions
multcomp::cld(emm_IS, adjust="none", Letters=LETTERS, reversed=T)
## trat emmean SE df lower.CL upper.CL .group
## T1 2714 581 2.19 410 5018 A
## T3 1187 580 2.17 -1124 3499 B
##
## Results are averaged over the levels of: bloq, forestal
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## 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.
## Medias interacción forestal*trat
emm_int<-emmeans(mod2, ~ trat*forestal,lmer.df = "kenward-roger")
mult<-multcomp::cld(emm_int, adjust="none", Letters=LETTERS, reversed=T)
mult
## trat forestal emmean SE df lower.CL upper.CL .group
## T1 Abarco 4521 1002 2.15 494 8547 A C
## T1 Terminalia 1873 1002 2.15 -2153 5900 ABCD
## T1 Roble 1748 1013 2.25 -2174 5670 AB
## T3 Abarco 1493 1006 2.20 -2489 5474 B D
## T3 Terminalia 1449 1002 2.15 -2577 5476 ABCD
## T3 Roble 620 1004 2.17 -3385 4626 CD
##
## Results are averaged over the levels of: bloq
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## 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.
## Medias comparando tratamientos en cada forestal
emm_trat_for<-emmeans(mod2, pairwise~trat|forestal,lmer.df = "kenward-roger")
multcomp::cld(emm_trat_for, adjust="none", Letters=LETTERS, reversed=T)
## forestal = Abarco:
## trat emmean SE df lower.CL upper.CL .group
## T1 4521 1002 2.15 494 8547 A
## T3 1493 1006 2.20 -2489 5474 B
##
## forestal = Roble:
## trat emmean SE df lower.CL upper.CL .group
## T1 1748 1013 2.25 -2174 5670 A
## T3 620 1004 2.17 -3385 4626 B
##
## forestal = Terminalia:
## trat emmean SE df lower.CL upper.CL .group
## T1 1873 1002 2.15 -2153 5900 A
## T3 1449 1002 2.15 -2577 5476 A
##
## Results are averaged over the levels of: bloq
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## 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.
##Guardando base de datos comparación interacción
write.csv(mult, "~/Library/CloudStorage/GoogleDrive-icarounam@gmail.com/Mi unidad/Agrosavia/Colaboraciones/Laura/Intercrop/mulprod.csv" )
datgraf<-read.table("mulprod.csv", header=T, sep=",")
## Gráfica
graf<-ggplot(datgraf, aes(x=trat, y=emmean, group=forestal, color=forestal)) +
geom_line()+
geom_pointrange(aes(ymin=emmean-SE, ymax=emmean+SE))
graf<-graf + labs(x = "Intercropping system", y = "yield",
tag = "A")
graf
