setwd("~/Downloads")
dat<-read.table("bdfrijol.csv", header=T, sep=",")
dat$trat<-as.factor(dat$trat)
dat$bloq<-as.factor(dat$bloq)
dat$forestal<-as.factor(dat$forestal)
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'.
## Metodo 1 (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 + trat:forestal +
bloq + (1|bloq:forestal),
data=dat)
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 1128339 564170 2 2 1.7462 0.36414
## trat 1465168 1465168 1 111 4.5349 0.03542 *
## bloq 704615 704615 1 2 2.1809 0.27776
## forestal:trat 957692 478846 2 111 1.4821 0.23162
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Medias para trat
mean_comparisons <- mod %>%
emmeans(specs = ~ trat,
lmer.df = "kenward-roger")
## NOTE: Results may be misleading due to involvement in interactions
mean_comparisons # adjusted means
## trat emmean SE df lower.CL upper.CL
## T2 1635 111 3.26 1296 1973
## T3 1414 111 3.26 1075 1752
##
## Results are averaged over the levels of: forestal, bloq
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Método 2 (método clásico usando paquete Agricolae)
mod2<-with(dat,sp.plot(bloq,forestal,trat,prod))
##
## ANALYSIS SPLIT PLOT: prod
## Class level information
##
## forestal : Abarco Roble Terminalia
## trat : T3 T2
## bloq : B1 B2
##
## Number of observations: 120
##
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## bloq 1 2532957 2532957 7.8759 0.005944 **
## forestal 2 4056162 2028081 1.7462 0.364140
## Ea 2 2322856 1161428 1.0985 0.337072
## trat 1 1465168 1465168 3.8942 0.142976
## forestal:trat 2 957692 478846 1.2727 0.397907
## Eb 3 1128729 376243 1.0985 0.337072
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## cv(a) = 70.7 %, cv(b) = 40.2 %, Mean = 1524.337
## No hay significancia para las fuentes de variación principales (solo el bloque) pero se corre el ejemplo para las medias de todas las fuentes de variación
#Medias para forestal
gla<-mod2$gl.a; glb<-mod2$gl.b; glc<-mod2$gl.c
Ea<-mod2$Ea; Eb<-mod2$Eb; Ec<-mod2$Ec
par(mfrow=c(1,3),cex=0.6)
out1<-with(dat,LSD.test(prod,forestal,gla,Ea,console=TRUE))
##
## Study: prod ~ forestal
##
## LSD t Test for prod
##
## Mean Square Error: 1161428
##
## forestal, means and individual ( 95 %) CI
##
## prod std r LCL UCL Min Max
## Abarco 1585.751 691.9981 40 852.5850 2318.917 670.00 3496.53
## Roble 1712.430 613.4921 40 979.2638 2445.596 779.66 3492.00
## Terminalia 1274.831 500.9478 40 541.6645 2007.997 460.39 2478.96
##
## Alpha: 0.05 ; DF Error: 2
## Critical Value of t: 4.302653
##
## least Significant Difference: 1036.854
##
## Treatments with the same letter are not significantly different.
##
## prod groups
## Roble 1712.430 a
## Abarco 1585.751 a
## Terminalia 1274.831 a
#Medias para tratamientos
out2<-with(dat,LSD.test(prod,trat,glb,Eb,console=TRUE))
##
## Study: prod ~ trat
##
## LSD t Test for prod
##
## Mean Square Error: 376242.9
##
## trat, means and individual ( 95 %) CI
##
## prod std r LCL UCL Min Max
## T2 1634.835 724.5785 60 1382.824 1886.846 460.39 3496.53
## T3 1413.840 500.1066 60 1161.829 1665.851 519.55 2788.77
##
## Alpha: 0.05 ; DF Error: 3
## Critical Value of t: 3.182446
##
## least Significant Difference: 356.3975
##
## Treatments with the same letter are not significantly different.
##
## prod groups
## T2 1634.835 a
## T3 1413.840 a
#Medias para la interacción forestal*tratamientos
out3<-with(dat, LSD.test(prod, forestal:trat, glb, Eb, console = TRUE))
##
## Study: prod ~ forestal:trat
##
## LSD t Test for prod
##
## Mean Square Error: 376242.9
##
## forestal:trat, means and individual ( 95 %) CI
##
## prod std r LCL UCL Min Max
## Abarco:T2 1619.435 833.0952 20 1182.9385 2055.930 670.00 3496.53
## Abarco:T3 1552.068 535.2473 20 1115.5720 1988.564 820.67 2680.46
## Roble:T2 1948.201 648.3788 20 1511.7055 2384.697 1149.06 3492.00
## Roble:T3 1476.659 484.9033 20 1040.1625 1913.154 779.66 2788.77
## Terminalia:T2 1336.869 564.6155 20 900.3730 1773.365 460.39 2478.96
## Terminalia:T3 1212.793 433.8340 20 776.2965 1649.288 519.55 2200.37
##
## Alpha: 0.05 ; DF Error: 3
## Critical Value of t: 3.182446
##
## least Significant Difference: 617.2985
##
## Treatments with the same letter are not significantly different.
##
## prod groups
## Roble:T2 1948.201 a
## Abarco:T2 1619.435 ab
## Abarco:T3 1552.068 ab
## Roble:T3 1476.659 ab
## Terminalia:T2 1336.869 ab
## Terminalia:T3 1212.793 b