setwd("~/Google Drive/Agrosavia/colaboraciones/Laura")
datos<-read.table("campohermoso.csv", header=T, sep=',')
datos$gen<-as.factor(datos$gen)
datos$forestal<-as.factor(datos$forestal)
datos$bloque<-as.factor(datos$bloque)
attach(datos)
library(ggplot2)
#Gráfica diámetro
ggplot(datos, 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 333 rows containing non-finite values (stat_smooth).
# Gráfica altura
ggplot(datos, 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 332 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=datos, corr=corCompSymm(, form= ~ 1 | gen),na.action=na.exclude)
fit.ar1.diam <- gls(diam ~ semana*forestal*gen+bloque, data=datos, corr=corAR1(, form= ~ 1 | gen), na.action=na.exclude)
fit.ar1het.diam <- gls(diam ~ semana*forestal*gen+bloque, data=datos, 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 34 1870.932 2026.389 -901.4659
## fit.ar1.diam 2 34 1858.949 2014.407 -895.4745
## fit.ar1het.diam 3 37 1857.064 2026.238 -891.5319 2 vs 3 7.885065 0.0484
anova(fit.ar1.diam)
## Denom. DF: 715
## numDF F-value p-value
## (Intercept) 1 10145.383 <.0001
## semana 1 110.023 <.0001
## forestal 2 47.989 <.0001
## gen 4 1.885 0.1113
## bloque 2 41.483 <.0001
## semana:forestal 2 0.833 0.4351
## semana:gen 4 0.476 0.7530
## forestal:gen 8 5.778 <.0001
## semana:forestal:gen 8 0.646 0.7389
anova(fit.ar1het.diam)
## Denom. DF: 715
## numDF F-value p-value
## (Intercept) 1 10189.583 <.0001
## semana 1 143.762 <.0001
## forestal 2 48.642 <.0001
## gen 4 1.778 0.1313
## bloque 2 39.717 <.0001
## semana:forestal 2 0.890 0.4110
## semana:gen 4 0.623 0.6461
## forestal:gen 8 6.062 <.0001
## semana:forestal:gen 8 0.808 0.5955
# Análisis para altura
fit.compsym.alt <- gls(alt ~ semana*forestal*gen+bloque, data=datos, corr=corCompSymm(, form= ~ 1 | gen),na.action=na.exclude)
fit.ar1.alt <- gls(alt ~ semana*forestal*gen+bloque, data=datos, corr=corAR1(, form= ~ 1 | gen), na.action=na.exclude)
fit.ar1het.alt <- gls(alt ~ semana*forestal*gen+bloque, data=datos, 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 34 7321.417 7476.922 -3626.709
## fit.ar1.alt 2 34 7316.113 7471.618 -3624.056
## fit.ar1het.alt 3 37 7301.865 7471.091 -3613.932 2 vs 3 20.24765 2e-04
anova(fit.ar1het.alt)
## Denom. DF: 716
## numDF F-value p-value
## (Intercept) 1 6620.156 <.0001
## semana 1 191.354 <.0001
## forestal 2 38.975 <.0001
## gen 4 7.450 <.0001
## bloque 2 33.728 <.0001
## semana:forestal 2 0.831 0.4359
## semana:gen 4 1.823 0.1225
## forestal:gen 8 8.072 <.0001
## semana:forestal:gen 8 0.598 0.7796
#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 ab TCS06
## TCS13 ab TCS13
## TCS19 b 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 a Abarco
## Roble a 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 d Abarco:CCN51
## Abarco:TCS01 d Abarco:TCS01
## Abarco:TCS06 abcd Abarco:TCS06
## Abarco:TCS13 d Abarco:TCS13
## Abarco:TCS19 d Abarco:TCS19
## Roble:CCN51 acd Roble:CCN51
## Roble:TCS01 acd Roble:TCS01
## Roble:TCS06 d Roble:TCS06
## Roble:TCS13 d Roble:TCS13
## Roble:TCS19 abc Roble:TCS19
## Terminalia:CCN51 cd Terminalia:CCN51
## Terminalia:TCS01 abc Terminalia:TCS01
## Terminalia:TCS06 ab 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 b CCN51
## TCS01 a TCS01
## TCS06 a TCS06
## TCS13 a TCS13
## TCS19 a 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 a Abarco
## Roble a 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 cdef Abarco:TCS01
## Abarco:TCS06 abcd Abarco:TCS06
## Abarco:TCS13 acde Abarco:TCS13
## Abarco:TCS19 def Abarco:TCS19
## Roble:CCN51 abcd Roble:CCN51
## Roble:TCS01 cde Roble:TCS01
## Roble:TCS06 ef Roble:TCS06
## Roble:TCS13 cde Roble:TCS13
## Roble:TCS19 acd Roble:TCS19
## Terminalia:CCN51 cde Terminalia:CCN51
## Terminalia:TCS01 abc Terminalia:TCS01
## Terminalia:TCS06 ab Terminalia:TCS06
## Terminalia:TCS13 ab Terminalia:TCS13
## Terminalia:TCS19 b Terminalia:TCS19
detach(datos)