setwd("~/Google Drive/Agrosavia/colaboraciones/Laura")
datos4<-read.table("santamaria.csv", header=T, sep=',')
datos4$gen<-as.factor(datos4$gen)
datos4$forestal<-as.factor(datos4$forestal)
datos4$bloque<-as.factor(datos4$bloque)
attach(datos4)
library(ggplot2)
#Gráfica diámetro
ggplot(datos4, aes(semana, diam, group = gen, colour = gen)) +
facet_grid(~forestal) +
geom_smooth(method="lm", se=F) +
theme_classic() +
xlab ("Semana") +
ylab ("Diámetro") +
labs(colour = "Genotipo") +
theme_linedraw() +
theme(
plot.title = element_text(hjust = 0.5, size = 16),
strip.text = element_text(size = 16),
axis.title.y = element_text(size = 16),
axis.title.x = element_text(size = 16),
axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14)
)
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 34 rows containing non-finite values (stat_smooth).
# Gráfica altura
ggplot(datos4, aes(semana, alt, group = gen, colour = gen)) +
facet_grid(~forestal) +
geom_smooth(method="lm", se=F) +
theme_classic() +
xlab ("Semana") +
ylab ("Altura") +
labs(colour = "Genotipo") +
theme_linedraw() +
theme(
plot.title = element_text(hjust = 0.5, size = 16),
strip.text = element_text(size = 16),
axis.title.y = element_text(size = 16),
axis.title.x = element_text(size = 16),
axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14)
)
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 34 rows containing non-finite values (stat_smooth).
# Anova general
aov.diam<-aov(diam~semana*forestal*gen+bloque)
aov.alt<-aov(alt~semana*forestal*gen+bloque)
#Análisis para diámetro
library(nlme)
fit.compsym.diam <- gls(diam ~ semana*forestal*gen+bloque, data=datos4, corr=corCompSymm(, form= ~ 1 | gen),na.action=na.exclude)
fit.ar1.diam <- gls(diam ~ semana*forestal*gen+bloque, data=datos4, corr=corAR1(, form= ~ 1 | gen), na.action=na.exclude)
fit.ar1het.diam <- gls(diam ~ semana*forestal*gen+bloque, data=datos4, corr=corAR1(, form= ~ 1 | gen), weights=varIdent(form = ~ 1 | semana), na.action=na.exclude)
anova(fit.compsym.diam, fit.ar1.diam, fit.ar1het.diam) #compares the models
## Model df AIC BIC logLik Test L.Ratio p-value
## fit.compsym.diam 1 33 1929.753 2077.746 -931.8766
## fit.ar1.diam 2 33 1884.634 2032.627 -909.3171
## fit.ar1het.diam 3 36 1799.179 1960.626 -863.5895 2 vs 3 91.45515 <.0001
anova(fit.ar1.diam)
## Denom. DF: 655
## numDF F-value p-value
## (Intercept) 1 9452.403 <.0001
## semana 1 448.364 <.0001
## forestal 2 0.699 0.4977
## gen 4 2.206 0.0669
## bloque 1 1.840 0.1754
## semana:forestal 2 0.201 0.8176
## semana:gen 4 0.162 0.9574
## forestal:gen 8 6.682 <.0001
## semana:forestal:gen 8 1.318 0.2312
anova(fit.ar1het.diam)
## Denom. DF: 655
## numDF F-value p-value
## (Intercept) 1 10600.425 <.0001
## semana 1 858.726 <.0001
## forestal 2 0.590 0.5547
## gen 4 2.961 0.0193
## bloque 1 2.649 0.1041
## semana:forestal 2 0.361 0.6969
## semana:gen 4 0.312 0.8697
## forestal:gen 8 5.407 <.0001
## semana:forestal:gen 8 2.315 0.0188
# Análisis para altura
fit.compsym.alt <- gls(alt ~ semana*forestal*gen+bloque, data=datos4, corr=corCompSymm(, form= ~ 1 | gen),na.action=na.exclude)
fit.ar1.alt <- gls(alt ~ semana*forestal*gen+bloque, data=datos4, corr=corAR1(, form= ~ 1 | gen), na.action=na.exclude)
fit.ar1het.alt <- gls(alt ~ semana*forestal*gen+bloque, data=datos4, corr=corAR1(, form= ~ 1 | gen), weights=varIdent(form = ~ 1 | semana), na.action=na.exclude)
anova(fit.compsym.alt, fit.ar1.alt, fit.ar1het.alt) #compares the models
## Model df AIC BIC logLik Test L.Ratio p-value
## fit.compsym.alt 1 33 6643.458 6791.451 -3288.729
## fit.ar1.alt 2 33 6635.769 6783.762 -3284.885
## fit.ar1het.alt 3 36 6615.120 6776.567 -3271.560 2 vs 3 26.64967 <.0001
anova(fit.ar1.alt)
## Denom. DF: 655
## numDF F-value p-value
## (Intercept) 1 20018.297 <.0001
## semana 1 1240.223 <.0001
## forestal 2 3.123 0.0447
## gen 4 5.110 0.0005
## bloque 1 0.131 0.7177
## semana:forestal 2 0.132 0.8766
## semana:gen 4 0.689 0.5999
## forestal:gen 8 11.487 <.0001
## semana:forestal:gen 8 1.065 0.3858
anova(fit.ar1het.alt)
## Denom. DF: 655
## numDF F-value p-value
## (Intercept) 1 18288.746 <.0001
## semana 1 1544.340 <.0001
## forestal 2 3.503 0.0307
## gen 4 6.449 <.0001
## bloque 1 0.006 0.9377
## semana:forestal 2 0.294 0.7457
## semana:gen 4 0.730 0.5715
## forestal:gen 8 10.470 <.0001
## semana:forestal:gen 8 1.374 0.2043
#Tukey diámetro
library(multcompView)
interac.tuk.diam<-TukeyHSD(aov.diam, "forestal:gen", ordered = TRUE)
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana,
## forestal
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana,
## gen
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana,
## forestal, gen
fores.tuk.diam<-TukeyHSD(aov.diam, "forestal", ordered = TRUE)
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana,
## forestal
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana,
## gen
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana,
## forestal, gen
gen.tuk.diam<-TukeyHSD(aov.diam, "gen", ordered = TRUE)
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana,
## forestal
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana,
## gen
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana,
## forestal, gen
#Tukey altura
interac.tuk.alt<-TukeyHSD(aov.alt, "forestal:gen", ordered = TRUE)
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana,
## forestal
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana,
## gen
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana,
## forestal, gen
fores.tuk.alt<-TukeyHSD(aov.alt, "forestal", ordered = TRUE)
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana,
## forestal
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana,
## gen
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana,
## forestal, gen
gen.tuk.alt<-TukeyHSD(aov.alt, "gen", ordered = TRUE)
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana,
## forestal
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana,
## gen
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana,
## forestal, gen
#Etiquetas Tukey diámetro
#Genotipos
generate_label_df_gen_diam <- function(gen.tuk.diam, variable){
Tukey.levels <- gen.tuk.diam[[variable]][,4]
Tukey.labels <- data.frame(multcompLetters(Tukey.levels)['Letters'])
Tukey.labels$treatment=rownames(Tukey.labels)
Tukey.labels=Tukey.labels[order(Tukey.labels$treatment) , ]
return(Tukey.labels)
}
labels.gen.diam <- generate_label_df_gen_diam(gen.tuk.diam, "gen")
labels.gen.diam
## Letters treatment
## CCN51 a CCN51
## TCS01 ab TCS01
## TCS06 a TCS06
## TCS13 b TCS13
## TCS19 a TCS19
# Forestal
generate_label_df_forestal_diam <- function(fores.tuk.diam, variable){
Tukey.levels <- fores.tuk.diam[[variable]][,2]
Tukey.labels <- data.frame(multcompLetters(Tukey.levels)['Letters'])
Tukey.labels$treatment=rownames(Tukey.labels)
Tukey.labels=Tukey.labels[order(Tukey.labels$treatment) , ]
return(Tukey.labels)
}
labels.forestal.diam <- generate_label_df_forestal_diam(fores.tuk.diam, "forestal")
labels.forestal.diam
## Letters treatment
## Abarco b Abarco
## Roble c Roble
## Terminalia a Terminalia
# Interacción Forestal:Genotipo
generate_label_df_interac_diam <- function(interac.tuk.diam, variable){
Tukey.levels <- interac.tuk.diam[[variable]][,4]
Tukey.labels <- data.frame(multcompLetters(Tukey.levels)['Letters'])
Tukey.labels$treatment=rownames(Tukey.labels)
Tukey.labels=Tukey.labels[order(Tukey.labels$treatment) , ]
return(Tukey.labels)
}
labels.interac.diam <- generate_label_df_interac_diam(interac.tuk.diam, "forestal:gen")
labels.interac.diam
## Letters treatment
## Abarco:CCN51 cdef Abarco:CCN51
## Abarco:TCS01 abcd Abarco:TCS01
## Abarco:TCS06 abc Abarco:TCS06
## Abarco:TCS13 acdef Abarco:TCS13
## Abarco:TCS19 f Abarco:TCS19
## Roble:CCN51 abcd Roble:CCN51
## Roble:TCS01 abcde Roble:TCS01
## Roble:TCS06 acdef Roble:TCS06
## Roble:TCS13 abcd Roble:TCS13
## Roble:TCS19 acdef Roble:TCS19
## Terminalia:CCN51 def Terminalia:CCN51
## Terminalia:TCS01 acdef Terminalia:TCS01
## Terminalia:TCS06 ef Terminalia:TCS06
## Terminalia:TCS13 b Terminalia:TCS13
## Terminalia:TCS19 ab Terminalia:TCS19
#Etiquetas Tukey altura
#Genotipos
generate_label_df_gen_alt <- function(gen.tuk.alt, variable){
Tukey.levels <- gen.tuk.alt[[variable]][,4]
Tukey.labels <- data.frame(multcompLetters(Tukey.levels)['Letters'])
Tukey.labels$treatment=rownames(Tukey.labels)
Tukey.labels=Tukey.labels[order(Tukey.labels$treatment) , ]
return(Tukey.labels)
}
labels.gen.alt <- generate_label_df_gen_alt(gen.tuk.alt, "gen")
labels.gen.alt
## Letters treatment
## CCN51 ac CCN51
## TCS01 ab TCS01
## TCS06 c TCS06
## TCS13 ab TCS13
## TCS19 b TCS19
# Forestal
generate_label_df_forestal_alt <- function(fores.tuk.alt, variable){
Tukey.levels <- fores.tuk.alt[[variable]][,2]
Tukey.labels <- data.frame(multcompLetters(Tukey.levels)['Letters'])
Tukey.labels$treatment=rownames(Tukey.labels)
Tukey.labels=Tukey.labels[order(Tukey.labels$treatment) , ]
return(Tukey.labels)
}
labels.forestal.alt <- generate_label_df_forestal_alt(fores.tuk.alt, "forestal")
labels.forestal.alt
## Letters treatment
## Abarco b Abarco
## Roble b Roble
## Terminalia a Terminalia
# Interacción Forestal:Genotipo
generate_label_df_interac_alt <- function(interac.tuk.alt, variable){
Tukey.levels <- interac.tuk.alt[[variable]][,4]
Tukey.labels <- data.frame(multcompLetters(Tukey.levels)['Letters'])
Tukey.labels$treatment=rownames(Tukey.labels)
Tukey.labels=Tukey.labels[order(Tukey.labels$treatment) , ]
return(Tukey.labels)
}
labels.interac.alt <- generate_label_df_interac_alt(interac.tuk.alt, "forestal:gen")
labels.interac.alt
## Letters treatment
## Abarco:CCN51 f Abarco:CCN51
## Abarco:TCS01 abcde Abarco:TCS01
## Abarco:TCS06 abcd Abarco:TCS06
## Abarco:TCS13 cdef Abarco:TCS13
## Abarco:TCS19 acdef Abarco:TCS19
## Roble:CCN51 ab Roble:CCN51
## Roble:TCS01 abc Roble:TCS01
## Roble:TCS06 ef Roble:TCS06
## Roble:TCS13 cdef Roble:TCS13
## Roble:TCS19 abcde Roble:TCS19
## Terminalia:CCN51 def Terminalia:CCN51
## Terminalia:TCS01 def Terminalia:TCS01
## Terminalia:TCS06 f Terminalia:TCS06
## Terminalia:TCS13 b Terminalia:TCS13
## Terminalia:TCS19 ab Terminalia:TCS19
detach(datos4)